Hii
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.
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 dolfinx.io.XDMFFile(mesh.comm, "reference_mesh.xdmf", "w") as mesh_file_xdmf:
mesh_file_xdmf.write_mesh(mesh)
# 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)
uh.interpolate(deformation_func)
# 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 dolfinx.io.XDMFFile(mesh.comm, "mesh.xdmf", "w") as mesh_file_xdmf:
mesh_file_xdmf.write_mesh(mesh)
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)
uh.interpolate(u_exact)
# Overwrite reference mesh coordinates
mesh.coordinates()[:, :mesh.geometry().dim()] = uh.compute_vertex_values().reshape(reference_coordinates.shape[0], mesh.geometry().dim())
plot(mesh)
plt.show()
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