6.2. Laplace smoothed mesh
Since the Lagrangian coordinates are initialized with the undeformed initial Eulerian coordinates, we can define the displacement from the initial configuration as \(\vec{d}=\vec{x}-\vec{\xi}\). One can smooth this displacement by solving a Laplace equation for \(\vec{d}\), i.e. \(\nabla_\xi^2\vec{d}=0\), where \(\nabla_\xi\) denotes the derivatives with respect to the Lagrangian coordinates. Thereby, any deformation that is imposed e.g. on the boundaries, will be smooth out along the mesh.
The weak formulation with test function \(\vec{\chi}\) reads
Let use hence define this equation class:
from pyoomph import *
from pyoomph.expressions import *
class LaplaceSmoothedMesh(Equations):
def define_fields(self):
# let the mesh coordinates become a variables, approximated with second order Lagrange basis functions
self.activate_coordinates_as_dofs(coordinate_space="C2")
def define_residuals(self):
x,xtest=var_and_test("mesh") # Eulerian mesh coordinates
xi=var("lagrangian") # fixed Lagrangian coordinates
d=x-xi # displacement
# Weak formulation: gradients and integrals are carried out with respect to the Lagrangian coordinates
self.add_residual(weak(grad(d,lagrangian=True), grad(xtest, lagrangian=True),lagrangian=True) )
We do not define fields, but we activate the mesh coordinates as dependent variables with the call activate_coordinates_as_dofs(). You can pass an argument coordinate_space to select the space. If further fields are added, the coordinate space must at least comprise the highest space of all defined fields, i.e. we cannot have a coordinate_space of "C1" and defining other fields on the space "C2". If the argument is omitted, the coordinate space will be automatically determined by the highest space of all added fields. The rest is trivial, except the usage of the variables "mesh" and "lagrangian" and the keyword arguments lagrangian=True to the grad() and weak() calls.
As an example problem, let us deform a rectangular mesh by prescribing Dirichlet boundary conditions at the interfaces and let the internal mesh relax based on the Laplace smoothing:
class LaplaceSmoothProblem(Problem):
def define_problem(self):
self.initial_adaption_steps=0
self.add_mesh(RectangularQuadMesh(N=6))
eqs=LaplaceSmoothedMesh()
eqs+=MeshFileOutput()
eqs+=DirichletBC(mesh_x=0,mesh_y=True)@"left" # fix the mesh at x=0 and keep y in place
eqs+=DirichletBC(mesh_x=True,mesh_y=0)@"bottom" # fix the mesh at y=0 and keep x in place
xi=var("lagrangian") # Lagrangian coordinate
eqs+=DirichletBC(mesh_x=1+0.5*xi[1])@"right" # linear slope at the left
eqs+=DirichletBC(mesh_y=1+0.25*xi[0]*(1-xi[0]))@"top" # quadratic deformation at the top
eqs+=SpatialErrorEstimator(mesh=1) # Adapt where large deformations are present
self.add_equations(eqs@"domain")
if __name__=="__main__":
with LaplaceSmoothProblem() as problem:
problem.output()
problem.solve(spatial_adapt=4)
problem.output_at_increased_time()
A few new things occur here. First, we set the property initial_adaption_steps of the problem class to 0. This controls the initial adaption, i.e. the adaption steps taken after the first solve. We deactivate this to get the middle mesh in Fig. 6.1. If this is not set, but a SpatialErrorEstimator is present, pyoomph will already adapt with respect to the initial condition. Then, the DirichletBC terms have values that are set to True instead to some value. This will fix the value of the variable at the interface, but it will not influence its value. Thereby, we can e.g. fix the \(y\)-coordinates of the "left" interface. Finally, note that we use the Lagrangian coordinate to prescribe the deformation in the DirichletBC term. We cannot use the Eulerian coordinate (i.e. var("mesh") or var("coordinate")) here, since these are now unknowns. Dirichlet boundary conditions may only depend on independent variables.
Finally, the SpatialErrorEstimator will refine the mesh where the deformation is rather discontinuous (cf. right panel in Fig. 6.1).
Fig. 6.1 Laplace smoothing: (left) undeformed mesh. (center) mesh after applying the Dirichlet boundary conditions that deform the mesh at the interfaces. (right) Relaxed mesh.