How to solve a coupled multiphysics problem using one mesh with subdomains?
4 months ago by
I am a new user of FEniCS trying to solve a multiphysics problem with its help. I have generated a simple mesh from two subdomains with a conforming interface between them, and created a FunctionSpace on the whole mesh using a MixedElement (composed of a Vector- and FiniteElement). Now I want to solve for two dependent variables (scalar and vector) in the first subdomain and a single scalar dependent variable in the second subdomain (the scalar variable is common for the two subdomains=the same variable), with a coupling condition on the interface. How can I do this in FEniCS? If I add any (non-physical, unwanted) component with the vector variable (even multiplied by a very small number) to the bilinear form in the second subdomain, I am able to get some (wrong) results. Otherwise, the solution is nan everywhere except where the boundary conditions were applied. Eliminating superfluous DOFs by setting
does not work correctly either. What is the proper way for solving such problems? What is worth mentioning here is that I have found a quite recent presentation (Link, page 5, 7th slide) tackling similar issue to some extent, but the key function FunctionSpaceProduct that in my opinion will serve the discussed purpose is not available in my distribution of dolfin (2017.2.0, ubuntu). Am I doing something wrong?
DirichletBC(W.sub(0), Constant((0.,)*ndim), SecondDomainIdentifier, method='pointwise')
A minimal working example is attached below.
Thank you very much for your time and help! Best regards.
# run using dolfin v2017.2.0 from dolfin import * from mshr import * def build_mesh(l1, l2, h, nl1, nl2, nh): domain_1 = Rectangle(Point(.0, .0), Point(l1, h)) # https://fenicsproject.org/qa/11220/how-to-merge-concatenate-add-meshed-in-one-mesh/ domain_2 = Rectangle(Point(l1, h), Point(l1+l2, .0)) domain_tot = domain_1 + domain_2 domain_tot.set_subdomain(1, domain_1) # https://fenicsproject.org/qa/9512/consisting-geometries-boundary-between-geometries-contained/ domain_tot.set_subdomain(2, domain_2) m = generate_mesh(domain_tot, 25) # FIXME: should use nl1, nl2, nh to control the mesh resolution return m #===== BORDERS ===== # Geometric parameters H = 4.E-3 # height L2 = 24.E-3 # second thickness L1 = L2 # first thickness #===== MESH ===== # Create mesh mesh = build_mesh(L1, L2, H, 10, 10, 8) #n = FacetNormal(mesh) # FIXME: this normal vector is not consistent on interior facets as it is probably defined cell-wise V = VectorFunctionSpace(mesh, "CG", mesh.topology().dim(), dim=2) n = project(Constant((-1., 0.)), V) # normal vectors pointing in the negative x-direction (outside the second domain at the interface) ndim = mesh.geometry().dim() #===== SUBDOMAINS ===== # Define subdomain markers subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), mesh.domains()) # Define integration measure for domain integrals dx = Measure('dx', domain=mesh, subdomain_data=subdomains) #===== FUNCTION SPACES ===== Uh = VectorElement('P', mesh.ufl_cell(), 2) Ph = FiniteElement('P', mesh.ufl_cell(), 1) W = FunctionSpace(mesh, MixedElement([Uh, Ph])) (u, p) = TrialFunctions(W) (du, dp) = TestFunctions(W) #===== BOUNDARY CONDITIONS ===== u_bc = Constant((0.,)*ndim) p_bc = Constant(1.) # Define boundaries tol = 1.E-14 def Left(x, on_boundary): return on_boundary and x < tol def Right(x, on_boundary): return on_boundary and near(x, L1 + L2, tol) def Symmetry(x, on_boundary): return on_boundary and near(x, 0., tol) or near(x, H, tol) def Interface(x): return near(x, L1, tol) class SecondDomain(SubDomain): # https://fenicsproject.org/qa/9049/how-to-create-a-dirichletbc-on-a-single-vertex-or-dof/ def inside(self, x, on_boundary): return x >= L1 - tol SecondDomainIdentifier = SecondDomain() bcs = [DirichletBC(W.sub(0), u_bc, Left), DirichletBC(W.sub(0).sub(1), Constant(0.), Symmetry), DirichletBC(W.sub(1), p_bc, Right)\ # , DirichletBC(W.sub(0), Constant((0.,)*ndim), SecondDomainIdentifier, method='pointwise')\ ] # Create mesh function over cell facets boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim()-1) boundary_subdomains.set_all(0) # 0 for interior cell boundaries interface_boundary = AutoSubDomain(Interface) interface_boundary.mark(boundary_subdomains, 1) # 1 for interface boundary # Define integration measure for boundary integrals dS = Measure('dS', domain=mesh, subdomain_data=boundary_subdomains) # measure for internal boundaries #===== PROBLEM DEFINITION ===== a = p*dp*dx(2) #a += inner(u, du)*dx(2) # TODO: non-physical component a += ( p*dp + inner(u, du) )*dx(1) a_interface = ( avg(dp)*inner(u('+'),n) + avg(p)*inner(du('+'),n) )*dS(1) L = Constant(0.)*dp*dx(1) A = assemble(a) # component matrix A_interface = assemble(a_interface) # interface coupling matrix A += A_interface # governing matrix b = assemble(L) # right-hand-side vector for bc in bcs: bc.apply(A, b) #===== SOLUTION ===== up = Function(W) solve(A, up.vector(), b, 'lu')
Community: FEniCS Project
Please login to add an answer/comment or follow this question.