7.4. Mesh reconstruction upon large deformations
It happens frequently that a moving mesh deforms strongly so that the elements are not well suited for solving equations on it. In that case, it might be necessary to regenerate the mesh by fitting the interfaces with splines, recreate a new mesh and interpolate all defined fields, including their history values for time stepping, from the previous mesh to the new mesh. pyoomph invokes this with the RemeshWhen class, which must be added to bulk domains.
How the new mesh is created is controlled by the property remesher of the MeshTemplate. By default, a GmshTemplate just calls its own define_geometry() method again, i.e. you describe the new geometry yourself and thereby keep full control over the mesh. This approach, which allows full control over the local mesh size, is discussed in Section 7.6.
Alternatively, the geometry of the new mesh can be reconstructed automatically from the deformed one by a Remesher2d instance from the module pyoomph.meshes.remesher, which is what we do here. It fits the current boundaries by splines and estimates suitable mesh sizes from the current element sizes, but it only works for two-dimensional meshes. To use it, we have to import the module and set the remesher attribute to a Remesher2d, which requires the mesh template as first argument for the constructor:
from laplace_smoothed_mesh import *
from pyoomph.meshes.remesher import *
class RemeshingProblem(Problem):
def __init__(self):
super(RemeshingProblem, self).__init__()
self.remeshing=True # shall we remesh or not
self.remesh_options=RemeshingOptions(max_expansion=2,min_expansion=0.3,min_quality_decrease=0.2) # when to remesh
def define_problem(self):
# Create a mesh and add a remesher
mesh=RectangularQuadMesh(N=6)
mesh.remesher=Remesher2d(mesh)
The RemeshingOptions control when a mesh reconstruction shall be invoked. The parameters max_expansion and min_expansion give a threshold for which remeshing is invoked whenever an element has grown or shrunken above max_expansion or below min_expansion with respect to its initial size. This is determined based on the Cartesian area of each element. Furthermore, a min_quality_decrease can be set to the RemeshingOptions, since the area of an element can still remain quite close to its initial area, but the element becomes strongly anisotropic, e.g. by expanding in one direction while collapsing in the perpendicular direction.
Fig. 7.3 Without remeshing (top), the mesh can strongly deform. With remeshing (bottom), one can prevent this (at the cost of some computational time and small interpolation errors). One can furthermore control the local mesh size, e.g. to ensure a fine mesh at the sharp corner.
After that, the equations are assembled as usual. Here, we fix all boundaries except the "right" one, which will move with time:
# Add the mesh and use the Lagrange smoothed mesh
self.add_mesh(mesh)
eqs=LaplaceSmoothedMesh()
eqs+=MeshFileOutput()
# Fix some interfaces
eqs+=DirichletBC(mesh_x=0,mesh_y=True)@"left"
eqs+=DirichletBC(mesh_x=True,mesh_y=0)@"bottom"
eqs+=DirichletBC(mesh_y=1)@"top"
# Moving boundary
xi=var("lagrangian")
eqs+=DirichletBC(mesh_x=1+0.5*xi[1]*var("time"))@"right" # move the right interface with time
To automatically invoke remeshing, we just must add a RemeshWhen object to the domain, which takes RemeshingOptions as argument. When there are no further instructions added, pyoomph tries to estimate the local mesh resolution based on the previous resolution. This does not always work well and sometimes it is beneficial to instruct pyoomph to use specific mesh sizes instead. To that end, one can add RemeshMeshSize objects to interfaces and corners, by what the local mesh size (nondimensional, in terms of the spatial scale) can be set:
# Remeshing based on the options
eqs+=RemeshWhen(self.remesh_options)
# optional: setting particular sizes at interfaces or corners
eqs+=RemeshMeshSize(size=0.2)@"left" # size of 0.2 at the left interface
eqs += RemeshMeshSize(size=0.02) @ "right/top" # size of 0.02 at the top right corner
self.add_equations(eqs@"domain")
To run the simulation, there is nothing specifically to be done:
if __name__=="__main__":
with RemeshingProblem() as problem:
# problem.remesh_options.active=False we can deactivate remeshing by the remesh options as well
problem.run(10,outstep=True,startstep=0.5,maxstep=0.5)
As usual, the RemeshingOptions can be modified before the run() statement. In particular, one can deactivate it by setting active=False. A comparison without and with remeshing is shown in Fig. 7.3.
One additional option is to add an instance of the class SetLagrangianToEulerianAfterSolve to the equations of the moving mesh domain. This will set the Lagrangian coordinates to the Eulerian coordinates after each successful time step. Thereby, the mesh displacement in (7.1) will be zero at the beginning of the next time step.
Note
We really recommend to use the advanced remeshing methods described in Section 7.6 for real cases.