How to optimise re-meshing, or an alternative method for simulating a moving boundary through a pore network?

Hi everyone,

I’ll get straight to the point with my problem and example code, then add a bit more information at the bottom about why I’m doing this for those interested in the motivation behind it.

I’m solving a moving boundary problem iteratively that I’m re-meshing for on each iteration using mshr’s Polygon mesh function and the new position of my moving boundary. The re-meshing every iteration is proving to be very time consuming, especially due to the high resolution I need on my moving boundary. Is there a way to optimise this?

Ideally I’d like my code to run in parallel, however it takes even longer re-meshing in parallel. I’ve also tried using gmsh, but that is also slower (though I may not be doing it the best way). Any ideas?

My example code:

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
import time

rank = MPI.comm_world.rank
mesh_file = File('MeshFolder/moving_mesh.pvd')

#Set inital position of moving boundary with high resolution
points_y = np.linspace(0.1, -0.1, num=int(5001), endpoint=True).reshape(-2,1)
points_x = np.zeros(len(points_y)).reshape(-2,1)
points = np.concatenate((points_x,points_y),axis=1)

#First and last points have alpha=0 so they hold the corner positions
alpha = np.ones(len(points))
alpha[[0,len(points)-1]] = 0

#Create initial mesh
domain_vertices = [Point(0,-5),Point(5,-5),Point(5,5),Point(0,5)]
for idx in range(len(points)):
polygon_domain = Polygon(domain_vertices)
mesh = generate_mesh(polygon_domain,100)
mesh_file << (mesh,0)

for count in range(5):
    t0 = time.time()
    #Calculate new boundary position
    points[:,0] -= 0.1*alpha
    #Create new mesh
    domain_vertices = [Point(0,-5),Point(5,-5),Point(5,5),Point(0,5)]
    for idx in range(len(points)):
    polygon_domain = Polygon(domain_vertices)
    mesh = generate_mesh(polygon_domain,100)
    mesh_file << (mesh,count+1)
    #Display time taken
    if rank == 0:
        print(count, time.time()-t0)

I’m developing a front-tracking numerical method for simulating a fluid flowing through a pore network. I start with a random pore network (Fig1 - red is solid, green is empty space) with a liquid above it. I have some equations I solve in the liquid which affect the movement of the liquid. I start by meshing the liquid phase then iteratively solve the equations in the liquid, then move the parts of the boundary touching the green regions. As this results in my mesh distorting and I need the moving boundary to be able to split and move in different directions and re-join as it moves through the pore network, I then create a new mesh and interpolate the solution onto it. In order to accurately model the movement through the pore network I need a very high resolution in this region (Fig2).

It may make sense to use an arbitrary Lagrangian–Eulerian (ALE) formulation (or even a fully-Lagrangian formulation, for low Reynolds numbers) to deform the mesh until some criterion for excessive distortion or topology change triggers a re-mesh. This would cut down on the number of steps requiring re-meshing. Given a desired displacement field u (which must be a Function in VectorFunctionSpace(mesh,"CG",1)), you can use the function ALE.move(mesh,u) to deform mesh by that displacement.


Hi kamensky, thanks for your reply. Unfortuately this is something I’ve already tried unsuccesfully as I can’t come up with a function that gets the nodes on the pore walls to slide along them, they always end up moving into the solid or away from the side into the liquid.

An alternative I’ve thought of is to create a routine that slides all the pore wall boundary nodes based on the movement of the moving boundary, then use ALE.move(mesh,boundary), then when they get over a certain distance apart add new nodes in and remesh to maintain high resolution along the boundary. However, I fear sliding all the boundary nodes (due to the number of them and needing to be really accurate in matching the pore wall) will end up taking even longer than remeshing, but maybe this is my best option


I don’t know if it will be useful, but there is a class of methods called phase field models where you can control the different materials in your model by introducing a new variable that takes different values depending on which material it is. The interface is diffusive and you only need to mesh the domain once. Typically, they are used to simulate the evolution of microstructures or solidification, but nowadays many disciplines are starting to use them.
Here are some of the papers that deal with these models:
An introduction to phase-field modeling of microstructure evolution
Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model

I hope you find it useful,

1 Like