# Assemble matrix of a subdomain

**URL:** <https://fenicsproject.discourse.group/t/assemble-matrix-of-a-subdomain/11464>\
**Category:** dolfinx\
**Created:** [June 11, 2023, 9:51am UTC](https://fenicsproject.discourse.group/t/assemble-matrix-of-a-subdomain/11464 "2023-06-11T09:51:15Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![TaikoOsO](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/taikooso/32/4559_2.png) [@TaikoOsO](https://fenicsproject.discourse.group/u/TaikoOsO)\
**Post date:** [June 11, 2023, 9:51am UTC](https://fenicsproject.discourse.group/t/assemble-matrix-of-a-subdomain/11464/1 "2023-06-11T09:51:15Z")

</div>

Dear community,

I am trying to conver my dolfin code to dolfinx, however, I dont know how to assemble matrix of a subdomain in dolfinx.I know how to [assemble matrix over subdomain in dolfin](https://fenicsproject.discourse.group/t/efficient-matrix-assembly-over-subdomain/7000/2)

When I tried to write similar code to define the ufl.measure(‘dx’) in dolfinx:

```auto
subdomain_indices = np.where(Xnp>1)[0]
subdomain_values = [999] * len(subdomain_indices)
subdomain = mesh.meshtags(self.domain, self.tdim, subdomain_indices,subdomain_values)
dx_solid = ufl.Measure('dx', domain=self.domain, subdomain_data=subdomain,subdomain_id=999)

a = fem.form(ufl.inner(stress(_v), strain(_dv)) * dx_solid)

```

It shows an erroro happen when using fem.form that

> incompatible constructor arguments.

In this example, `Xnp` is a numpy array get from a `Xe=fem.Function(self.Ve)` using `Xe.x.array`. It is used to define subdomain.

Would you mind helping me to solve the problem?

I’d appreciate your kind help.

---

<div class="post-metadata">

**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:** [June 11, 2023, 9:53am UTC](https://fenicsproject.discourse.group/t/assemble-matrix-of-a-subdomain/11464/2 "2023-06-11T09:53:51Z")

</div>

Please make a minimal reproducible example, as a snippet of your code without all definitions required to reproduce the error, makes it near impossible helping you.

---

<div class="post-metadata">

**Author:** ![TaikoOsO](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/taikooso/32/4559_2.png) [@TaikoOsO](https://fenicsproject.discourse.group/u/TaikoOsO)\
**Post date:** [June 11, 2023, 10:42am UTC](https://fenicsproject.discourse.group/t/assemble-matrix-of-a-subdomain/11464/3 "2023-06-11T10:42:18Z")

</div>

Thank you for your reply,

a example is shown here:

```auto
from dolfinx import mesh, fem, la, io
import numpy as np
from petsc4py import PETSc
import ufl
from mpi4py import MPI

domain = mesh.create_rectangle( MPI.COMM_WORLD, np.array([[0.,0.],[80,50]]),[160, 100],cell_type=mesh.CellType.quadrilateral)
Vu = fem.VectorFunctionSpace(domain, ("CG", 1))
_v, _dv = ufl.TrialFunction(Vu), ufl.TestFunction(Vu)

Ve = fem.FunctionSpace(domain, ("DG", 0))
Xe = fem.Function(Ve)
Xe.vector.set(1.)

# boundary condition
leftEdge = mesh.locate_entities_boundary(domain, 1, lambda ijk: np.isclose(ijk[0], 0))
boundary_dofs = fem.locate_dofs_topological(Vu, 1, leftEdge)
uD = np.array((0,) * 2, dtype=PETSc.ScalarType)
bc = [fem.dirichletbc(uD, boundary_dofs, Vu)]

# point force
point = mesh.locate_entities_boundary(domain, 0, lambda ijk: np.isclose(ijk.T, [80,25,0]).all(axis=1))
pointForce_dofs = fem.locate_dofs_topological(Vu.sub(1), entity_dim=0,entities=point)
point_value = -1

lambda_ = 1.0 * 0.3 / ((1.0 + 0.3) * (1.0 - 2.0 * 0.3))
mu = 1.0 / (2.0 * (1.0 + 0.3))

def epsilon(u):
    return ufl.sym(ufl.grad(u))
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)

# define sub domain 
Xnp = Xe.x.array.copy()
subdomain_indices = np.where(Xnp>0.5)[0]
subdomain_values = [999] * len(subdomain_indices)
subdomain = mesh.meshtags(domain, 2, subdomain_indices, subdomain_values)
dx_999 = ufl.Measure('dx', domain=domain, subdomain_data=subdomain,subdomain_id=999)
        
k = ufl.inner(sigma(_v), epsilon(_dv)) * dx_999
f = fem.Constant(domain, PETSc.ScalarType((0,) * 2))
L = ufl.dot(f, _v) * ufl.dx

# assembly
a, L = fem.form(k), fem.form(L)
# assemble matrix and vector
A = fem.petsc.assemble_matrix(a, bcs=bc)
A.assemble()
b = fem.petsc.assemble_vector(L)
# set boundary conditions
fem.petsc.apply_lifting(b, [a], bcs=[bc])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(b, bc)
# apply point force
b.array[pointForce_dofs] = -1

# solver
opts = PETSc.Options()
opts['ksp_type'] = 'preonly'
opts['pc_type'] = 'cholesky'
solver = PETSc.KSP().create(domain.comm)
solver.setFromOptions()
    
# solve
uh = fem.Function(Vu)
solver.setOperators(A)
solver.solve(b, uh.vector)
uh.x.scatter_forward()

obj = 0.5 * uh.vector.dot(b)
print(obj)

```

if change `dx_999` to `ufl.dx` for defining `k`, the code works.

---

<div class="post-metadata">

**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:** [June 11, 2023, 1:15pm UTC](https://fenicsproject.discourse.group/t/assemble-matrix-of-a-subdomain/11464/4 "2023-06-11T13:15:28Z")

</div>

Next time, please provide the full error message as well,

```bash
Traceback (most recent call last):
  File "/root/shared/mwe234.py", line 46, in <module>
    a, L = fem.form(k), fem.form(L)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 176, in form
    return _create_form(form)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 171, in _create_form
    return _form(form)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 165, in _form
    return formcls(ufcx_form, V, coeffs, constants, subdomains, mesh, module.ffi, code)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 51, in __init__
    super(). __init__ (ffi.cast("uintptr_t", ffi.addressof(self._ufcx_form)),
TypeError: __init__ (): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.Form_float64(spaces: List[dolfinx::fem::FunctionSpace], integrals: Dict[dolfinx::fem::IntegralType, Tuple[List[Tuple[int, object]], dolfinx.cpp.mesh.MeshTags_int32]], coefficients: List[dolfinx.cpp.fem.Function_float64], constants: List[dolfinx.cpp.fem.Constant_float64], need_permutation_data: bool, mesh: dolfinx.cpp.mesh.Mesh = None)
    2. dolfinx.cpp.fem.Form_float64(form: int, spaces: List[dolfinx::fem::FunctionSpace], coefficients: List[dolfinx.cpp.fem.Function_float64], constants: List[dolfinx.cpp.fem.Constant_float64], subdomains: Dict[dolfinx::fem::IntegralType, dolfinx.cpp.mesh.MeshTags_int32], mesh: dolfinx.cpp.mesh.Mesh)

Invoked with: <cdata 'uintptr_t' 140175546442848>, [<dolfinx.cpp.fem.FunctionSpace object at 0x7f7d375b53b0>, <dolfinx.cpp.fem.FunctionSpace object at 0x7f7d375b53b0>], [], [], {<IntegralType.cell: 0>: <dolfinx.mesh.MeshTagsMetaClass object at 0x7f7d373722a0>, <IntegralType.exterior_facet: 1>: None, <IntegralType.interior_facet: 2>: None, <IntegralType.vertex: 3>: None}, <dolfinx.mesh.Mesh object at 0x7f7d39d81800>

```

This says that the input to the form constructor is wrong. This is because of the way you have defined your meshtag/subdomain marker.

```python
subdomain_indices = np.where(Xnp>0.5)[0].astype(np.int32)
subdomain_values = np.full_like(subdomain_indices, 999, dtype=np.int32)
subdomain = mesh.meshtags(domain, 2, subdomain_indices, subdomain_values)
dx_999 = ufl.Measure('dx', domain=domain, subdomain_data=subdomain,subdomain_id=999)
        
k = ufl.inner(sigma(_v), epsilon(_dv)) * dx_999

```

Runs, as we require that the indices and values in a mesh-tag object, that is passed to C++ uses `int32`.

However, I would like to note that by assembling the matrix on only a part of the domain, the matrix contains rows with only zeros.

---

<div class="post-metadata">

**Author:** ![TaikoOsO](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/taikooso/32/4559_2.png) [@TaikoOsO](https://fenicsproject.discourse.group/u/TaikoOsO)\
**Post date:** [June 12, 2023, 9:12am UTC](https://fenicsproject.discourse.group/t/assemble-matrix-of-a-subdomain/11464/5 "2023-06-12T09:12:17Z")

</div>

Thank you very much. The code works when Xnp==1 for all cells.However, when I tried to create a subdoamin of the whole model, for example, `Xnp[0:500] = 0`:

```auto
# define sub domain 
Xnp = Xe.x.array.copy()
Xnp[0:500] = 0
subdomain_indices = np.where(Xnp>0.5)[0].astype(np.int32)
subdomain_values = np.full_like(subdomain_indices, 999, dtype=np.int32)
subdomain = mesh.meshtags(domain, 2, subdomain_indices, subdomain_values)
dx_999 = ufl.Measure('dx', domain=domain, subdomain_data=subdomain,subdomain_id=999)
        
k = ufl.inner(sigma(_v), epsilon(_dv)) * dx_999

```

Then an error happend in **petsc4py.PETSc.KSP.solve()**

```auto
Traceback (most recent call last)
---> 69 solver.solve(b, uh.vector)

File PETSc/KSP.pyx:403, in petsc4py.PETSc.KSP.solve()

Error: error code 73
[0] KSPSolve() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:1078
[0] KSPSolve_Private() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:843
[0] KSPSetUp() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:407
[0] PCSetUp() at /usr/local/petsc/src/ksp/pc/interface/precon.c:990
[0] PCSetUp_Cholesky() at /usr/local/petsc/src/ksp/pc/impls/factor/cholesky/cholesky.c:84
[0] MatCholeskyFactorSymbolic() at /usr/local/petsc/src/mat/interface/matrix.c:3139
[0] MatCholeskyFactorSymbolic_SeqAIJ() at /usr/local/petsc/src/mat/impls/aij/seq/aijfact.c:2765
[0] Object is in wrong state
[0] Matrix is missing diagonal entry 0

```

When I used FEniCS, I use the codes shown below to assemble the matrix and they works well.

```auto
subdomain = MeshFunction("size_t", domain, domain.topology().dim(), 0)
subdomain.array()[:] = np.where(Xnp > 0.5, 1, 0)
dx_sub = Measure('dx', domain=domain, subdomain_data=subdomain)
_a = inner(sigma(_v), epsilon(_dv)) * dx_sub(1)

A = assemble(_a, keep_diagonal=True)

```

Would you mind helping me to solve the problem in dolfinx? I was wondering if there is function like `assemble` with `keep_diagonal` in dolfinx?

---

<div class="post-metadata">

**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:** [June 12, 2023, 9:25am UTC](https://fenicsproject.discourse.group/t/assemble-matrix-of-a-subdomain/11464/6 "2023-06-12T09:25:56Z")

</div>

`keep_diagonal` should not fix this issue by it-self, as it only keeps zero-entries in the sparsity pattern (ref: [Bitbucket](https://bitbucket.org/fenics-project/dolfin/src/8060a84e8ada0df018772497534092727b06bd56/python/dolfin/fem/assembling.py?at=master#lines-113:116))

You would have to call `ident_zeros` in legacy DOLFIN to resolve that issue (which would add identity rows in these places).  
For this to be done with dolfinx, see for instance:

> [@Dirichlet Boundary Conditions on Subdomains in FEniCSx](https://fenicsproject.discourse.group/t/dirichlet-boundary-conditions-on-subdomains-in-fenicsx/9716/2):
>
> I guess you are referring to A.ident\_zeros() in legacy dolfin? As this simply calls petsc, you can still do this in DOLFINx, as illustrated in the following code: (The following code is a refined snippet of @CastriMik code at [thomas-fermi-demo/demo\_tf\_mwe.py at mpi-issue · mikics/thomas-fermi-demo · GitHub](https://github.com/mikics/thomas-fermi-demo/blob/mpi-issue/demo_tf_mwe.py#L191) which avoid global communciation) problem.\_A.zeroEntries() fem.petsc.assemble\_matrix(problem.\_A, problem.\_a, bcs=problem.bcs) problem.\_A.assemble() # Get diagonal of assembled A matrix di…

or

> [@Modify matrix diagonal -- dolfinx version for A.ident\_zeros()](https://fenicsproject.discourse.group/t/modify-matrix-diagonal-dolfinx-version-for-a-ident-zeros/9918/3):
>
> I would suggest the following: import time import dolfinx import numpy as np import ufl from dolfinx import cpp, fem from mpi4py import MPI from petsc4py import PETSc N = 75 degree = 1 mesh = dolfinx.mesh.create\_unit\_cube(MPI.COMM\_WORLD, N, N, N) V = dolfinx.fem.FunctionSpace(mesh, ("CG", degree)) u = ufl.TrialFunction(V) v = ufl.TestFunction(V) print( f"N={N}, Num dofs global {V.dofmap.index\_map.size\_global\*V.dofmap.index\_map\_bs}") # form is defined on the boundary =\> matrix wil have man…