I am trying to implement fluid-rigid body interaction with a gear-like structure moving in fluids and here I am trying to move the mesh around. I modify the code from https://fenicsproject.org/qa/13470/ale-move-class-and-meshes/
from dolfin import * from mshr import * import matplotlib as plt def solve_linear_elasticity(mesh, boundaries, d): c = Constant(d) V = VectorFunctionSpace(mesh, "Lagrange", 1) u = TrialFunction(V) v = TestFunction(V) E, nu = 10.0, 5.0 mu = E/(2.0*(1.0 + nu)) lmbda = E*nu/((1.0 + nu)*(1.0 -2.0*nu)) sigma = 2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(2) F = inner(sigma, grad(v))*dx a, L = lhs(F), rhs(F) bcs = [DirichletBC(V, Constant((0.0, 0.0)), boundaries, 1), DirichletBC(V.sub(0), c, boundaries, 2)] displacement = Function(V) solve(a==L, displacement, bcs) return displacement def update_mesh(mesh, displacement, boundaries): new_mesh = Mesh(mesh) new_boundaries = MeshFunction("size_t", new_mesh, 1) new_boundaries.set_values(boundaries.array()) ALE.move(new_mesh, displacement) return new_mesh, new_boundaries # Original mesh # Generate mesh domain_vertices = [Point(20, 23.125), Point(20, 21.25), Point(17.2925, 21.5625), Point(18.9175, 20.625), Point(17.2925, 18.4375), Point(18.9175, 19.375), Point(20, 16.875), Point(20, 18.75), Point(22.705, 18.4375), Point(21.0825, 19.375), Point(22.705, 21.5625), Point(20.995, 20.8275)] inclusion = Polygon(domain_vertices) mesh = generate_mesh(Rectangle(Point(0, 0), Point(40, 40))-inclusion, 10) class on_inclusion(SubDomain): def inside(self, x, on_boundary): return (on_boundary and not(near(x, 0) or near(x, 40) or near(x, 0) or near(x, 40))) # Initialize mesh function for boundary domains boundaries = MeshFunction('size_t', mesh, 1, mesh.domains()) boundaries.set_all(1) inc = on_inclusion() inc.mark(boundaries, 2) displacement = solve_linear_elasticity(mesh, boundaries, 5) mesh, boundaries = update_mesh(mesh, displacement, boundaries) plt.figure plot(mesh)
Before the moving it looks like this (the gear structure is in the middle)
After the moving the meshes becomes this
It seems the cells are entangled. Is there anyway to prevent this?
Besides, I would like to include the rotation of the object. Then the boundary condition would vary at each point on the gear, and I wonder how should I specify it?