# Extract dof index on different cell-subdomains **URL:** \ **Category:** dolfinx\ **Created:** [March 7, 2022, 3:25pm UTC](https://fenicsproject.discourse.group/t/extract-dof-index-on-different-cell-subdomains/7795 "2022-03-07T15:25:18Z")\ **Posts on this page:** 4\ **Page:** 1
**Author:** ![YYWayne](https://avatars.discourse-cdn.com/v4/letter/y/ad7895/32.png) [@YYWayne](https://fenicsproject.discourse.group/u/YYWayne)\ **Post date:** [March 7, 2022, 3:25pm UTC](https://fenicsproject.discourse.group/t/extract-dof-index-on-different-cell-subdomains/7795/1 "2022-03-07T15:25:18Z")
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) ![1646664708(1)](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/e/ee93d1de2dc2371b161dc25dad43938e1df43dc6.png) ![1646664633(1)](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/4d1f2828309c12493a4441280041a2273d3aff91.png) ![1646664662(1)](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/542c60823f376c782b3b1b353ba5c8a577be1b2f.png) 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. ```auto 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) ``` ![1646666366(1)](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/2/2ec42d0c79389db18ca2ad2544ba362bce557f82.png) 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. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [March 7, 2022, 5:53pm UTC](https://fenicsproject.discourse.group/t/extract-dof-index-on-different-cell-subdomains/7795/2 "2022-03-07T17:53:42Z")
You had some typos in your code extracting matrix data; > [@YYWayne](#): > > ```auto > 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]))) > > ``` should be ```auto Matrix1 = np.real(A1.getValues(range(A1.getSize()[0]), range(A1.getSize()[1]))) Matrix2 = np.real(A2.getValues(range(A2.getSize()[0]), range(A2.getSize()[1]))) Matrix = np.real(A.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]))) ``` which gives the picture: ![sp](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/0/00dd2d09d19f4e2fdbd50c2daf6cad1fe62de6f5.png) ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [March 7, 2022, 5:55pm UTC](https://fenicsproject.discourse.group/t/extract-dof-index-on-different-cell-subdomains/7795/3 "2022-03-07T17:55:45Z")
And it is a lot easier to spot the structure if you do skewed split, say at 0.3 ```python def marker1(x): return x[0] <= 0.3+1e-7 indices1 = mesh.locate_entities(msh, 2, marker1) def marker2(x): return x[0] >= 0.3-1e-7 ``` giving ![sp](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/c/c795603ffa5c7c14d8c066cfb2a9d0cb9a7549ee.png) ---
**Author:** ![YYWayne](https://avatars.discourse-cdn.com/v4/letter/y/ad7895/32.png) [@YYWayne](https://fenicsproject.discourse.group/u/YYWayne)\ **Post date:** [March 8, 2022, 1:33am UTC](https://fenicsproject.discourse.group/t/extract-dof-index-on-different-cell-subdomains/7795/4 "2022-03-08T01:33:18Z")
Thanks Dokken. Skewed split does give a easier spot.