Construct a physical domain with the help of transformation of a reference domain


If I have defined a reference geometry via Gmsh file in the code, then how can I define this transformation to get a physical domain?

\boldsymbol{Reference~ domain ~and~ Transformation~ is ~defined~ as}

The reference domain is defined as \hat{\Omega} given by the rectangle [0,1] \times[-1,1] . A fracture in the reference domain is then positioned in the middle of the rectangle, defined by

\hat{x}=200(0.05-\hat{y})(0.05+\hat{y}), \quad \hat{y} \in[-0.05,0.05] .

The physical domain \Omega, with more realistic geometry, is defined as a transformation of the reference domain \hat{\Omega} by the mapping

\left[\begin{array}{l} x \\ y \end{array}\right]=\left[\begin{array}{c} \hat{x} \\ 5 \cos \left(\frac{\hat{x}+\hat{y}}{100}\right) \cos \left(\frac{\pi \hat{x}+\hat{y}}{100}\right)^2+\hat{y} / 2-\hat{x} / 10 \end{array}\right] .

Please provide some example similar to this. Thank you in advance.

@niravshah241 has a library for mesh motion at GitHub - niravshah241/MDFEniCSx: Mesh Deformation using FEniCSx which may be of interest. If I were to summarize it in one line of code, it would be MDFEniCSx/demo/0_fundamental_deformation/ at main · niravshah241/MDFEniCSx · GitHub

Sir, I didn’t understand.
one line of code is

mesh.geometry.x[:, :mesh.geometry.dim] += \
    uh.x.array.reshape(reference_coordinates.shape[0], gdim)

How can I define this tranformation map ?

You need to create a perturbation function (uh) in the example above, that is the difference between the deformed domain and the reference domain (if you want to use +=), or you can simply use

mesh.geometry.x[:, :mesh.geometry.dim] = your_transformation_function(mesh.geometry.x)
where the my_transformation_function(x) takes in an x that is (3, num_points), i.e. each column will be the x,y,z coordinate of a node in the mesh.

You could also see: Changing the shape in the surface normal directions - #6 by dokken for more code explaining different deformation techniques.

Could you try this?

# "mesh" is Mesh on the reference domain
# Store reference mesh for visualiation
with, "reference_mesh.xdmf", "w") as mesh_file_xdmf:

# Store the reference mesh coordinates
reference_coordinates = mesh.geometry.x.copy()

# Define a VectorFunctionSpace with same degree as the mesh degree
V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", mesh.geometry.cmaps[0].degree))

# Define and interpolate deformation function
uh = dolfinx.fem.Function(V)
def deformation_func(x):
    return (x[0], 5. * np.cos((x[0] + x[1])/100.) * (np.cos((np.pi * x[0] + x[1])/100.))**2 + x[1]/2 - x[0]/10)


# Overwrite reference mesh coordinates
mesh.geometry.x[:, :mesh.geometry.dim] = uh.x.array.reshape(reference_coordinates.shape[0], mesh.geometry.dim)

# Store deformed mesh for visualiation
with, "mesh.xdmf", "w") as mesh_file_xdmf:
1 Like

Thank you, sir. I am not using DolfinX. I am trying to write in Dolfin. Is the way I am writing fine? The syntax works but when I plot the deformed mesh, it seems a weired plot.

from dolfin import *
import matplotlib.pyplot as plt

mesh = Mesh('Stokes.xml')
subdomains = MeshFunction("size_t", mesh, "Stokes_physical_region.xml")
bdry = MeshFunction("size_t", mesh, "Stokes_facet_region.xml")

reference_coordinates = mesh.coordinates()

# Define a VectorFunctionSpace with the same degree as the mesh degree
V = VectorFunctionSpace(mesh, "Lagrange", mesh.geometry().dim())
uh = Function(V)
# Define and interpolate deformation function
class deformation_func(UserExpression):
      def eval(self, values, x):
             values[0] = x[0]
             values[1] = 5*cos((x[0] + x[1])/100) * (cos((pi * x[0] + x[1])/100))**2 + x[1]/2 - x[0]/10
      def value_shape(self):
          return (2,)
u_exact = deformation_func(degree=2)


# Overwrite reference mesh coordinates
mesh.coordinates()[:, :mesh.geometry().dim()] = uh.compute_vertex_values().reshape(reference_coordinates.shape[0], mesh.geometry().dim())


Can you try below if you are using dolfin?

# Overwrite reference mesh coordinates
mesh.coordinates()[:, :mesh.geometry().dim()] = np.vstack((uh.compute_vertex_values()[:mesh.coordinates().shape[0]], uh.compute_vertex_values()[mesh.coordinates().shape[0]:])).T
1 Like