For a fluid dynamics system with the ficticious domain method, I need to assemble fluid and solid systems on two overlapped meshes seperately and then sum them up via interpolation:
# 1. Assemble seperately A x = b # fluid part As x = bs # solid part # 2. Sum them up via interpolation A += I.T As I b += I.T bs # 3. Solve solve(A, x, b)
I is the interpolation matrix created by
The function space is created by the mixed element, where vector element for velocity and finite element for pressure are needed. But only “first order finite element + serial run” scenario works, other scenarios such as “higher order finite element” and “parallel run” or their combination fail with Segmentation fault (signal 11) error. Here is a simple demo:
# question.py from dolfin import * import petsc4py petsc4py.init() # Order: 1 or 2 V_order = 1 # Change here to 2 then fail Q_order = 1 # Fluid mesh w = 1 mesh = RectangleMesh(Point(0, 0), Point(w, w), 64, 64, diagonal='crossed') # Cell/Solid mesh immersed in fluid r = w / 4. cell = RectangleMesh(Point(w/2-r/2, w/2-r/2), Point(w/2+r/2, w/2+r/2), 25, 25, diagonal='crossed') ## Fluid FunctionSpace V_VE = VectorElement("Lagrange", mesh.ufl_cell(), V_order) V_FE = FiniteElement("Lagrange", mesh.ufl_cell(), Q_order) # ===> switch following V_ME = MixedElement([V_VE, V_FE]) # vector + scalar: u + p, or with the equivalent following one: # V_ME = MixedElement([V_FE, V_FE, V_FE]) # scalar + scalar + scalar: ux + uy + p V_FS = FunctionSpace(mesh, V_ME) # Cell FunctionSpace + interpolation matrix VC_VE = VectorElement("Lagrange", cell.ufl_cell(), V_order) VC_FE = FiniteElement("Lagrange", cell.ufl_cell(), Q_order) # ===> switch following VC_ME = MixedElement([VC_VE, VC_FE]) # vector + scalar: u + p, or with the equivalent following one: # VC_ME = MixedElement([V_FE, V_FE, V_FE]) # scalar + scalar + scalar: ux + uy + p VC_FS = FunctionSpace(cell, VC_ME) # Interpolation matrix Ih = PETScDMCollection.create_transfer_matrix(V_FS, VC_FS) # <class 'dolfin.cpp.la.PETScMatrix'>, ns * n
It only works with
python question.py or
mpirun -n 1 question.py. Changing finite element order or running with more processes cause error. The only way I can do is to decompose the vector element to several finite elements, but this loses the generalization and increases the workload. How can I workaround this?