Hi, I am looking for extracting the dofs index on each subdomain (divided according to geometry coordinates, or partitioned by METIS in mpirun), so I can probably assemble sub-block matrix by looking into entries of original big matrix. Finally additive Schwarz preconditioner can be applied.
(R1 and R2 after permutation according to dof sequence are what im trying to get, figure1 is a example for 1d case when dofs are arranged sequentially)
I created a UnitSquare with celltag 1 when x<=0.5 and 2 when x>=0.5. I tried assemble the matrix in cases that integral with dx(1), dx(2) or both, but matrix patterns are same.
import dolfinx
import numpy as np
import matplotlib.pyplot as plt
import ufl
from mpi4py import MPI
from dolfinx import mesh
msh = mesh.create_unit_square(MPI.COMM_WORLD, 4, 4, cell_type = mesh.CellType.quadrilateral)
tdim = msh.topology.dim
# cell tag
num_cells = msh.topology.index_map(tdim).size_local + \
msh.topology.index_map(tdim).num_ghosts
marker1 = lambda x: x[0]<=0.5+1e-7
indices1 = mesh.locate_entities(msh, 2, marker1)
marker2 = lambda x: x[0]>=0.5-1e-7
indices2 = mesh.locate_entities(msh, 2, marker2)
cell_marker = np.zeros(num_cells, dtype=np.int32)
cell_marker[indices1] = 1
cell_marker[indices2] = 2
cell_tag = dolfinx.mesh.MeshTags(msh, dim=2, indices=np.arange(num_cells),
values=cell_marker)
dx = ufl.Measure('dx', subdomain_data=cell_tag, domain=msh)
V = dolfinx.fem.FunctionSpace(msh,("CG",1))
uD = dolfinx.fem.Function(V)
uD.interpolate(lambda x: x[0]-x[0]+1)
u,v = ufl.TrialFunction(V), ufl.TestFunction(V)
msh.topology.create_connectivity(tdim-1, tdim)
boundary_facets = np.flatnonzero(dolfinx.mesh.compute_boundary_facets(msh.topology)) # substitutive: facet_marker + np.where
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, tdim-1, boundary_facets)
bc = dolfinx.fem.dirichletbc(uD, boundary_dofs)
bcs = [bc]
a1 = ufl.inner(u,v) * dx(1)
a2 = ufl.inner(u,v) * dx(2)
a = ufl.inner(u,v) * dx(1) + ufl.inner(u,v) * dx(2)
a_block = dolfinx.fem.form([[a1,None], [None,a2]])
A1 = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a1), bcs=bcs)
A1.assemble()
A2 = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a2), bcs=bcs)
A2.assemble()
A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a), bcs=bcs)
A.assemble()
A_block = dolfinx.fem.petsc.assemble_matrix_block(a_block, bcs=bcs)
A_block.assemble()
print('norm:',A1.norm(),A2.norm(),A.norm(),A_block.norm())
Matrix1 = np.real(A1.getValues(range(A1.getSize()[0]), range(A1.getSize()[1])))
Matrix2 = np.real(A1.getValues(range(A2.getSize()[0]), range(A2.getSize()[1])))
Matrix = np.real(A1.getValues(range(A.getSize()[0]), range(A.getSize()[1])))
Matrix_block = np.real(A_block.getValues(range(A_block.getSize()[0]), range(A_block.getSize()[1])))
fig, axs = plt.subplots(2,2)
fig.set_size_inches(15,15)
plt.subplot(2,2,1)
plt.spy(Matrix1)
plt.subplot(2,2,2)
plt.spy(Matrix2)
plt.subplot(2,2,3)
plt.spy(Matrix)
plt.subplot(2,2,4)
plt.spy(Matrix_block)
Why matrix patterns are same after assembling? I must have misunderstood the function of assemble_matrix(). And how to achieve my goal (extrct dof index(different to mesh.dofmap when fem order>1) and assemble subdomains into a block diagnol form)? Great thanks.
