Assemble matrix of a subdomain

[TaikoOsO]

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

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

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.

[dokken]

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.

[TaikoOsO]

Thank you for your reply,

a example is shown here:

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.

[dokken]

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

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.

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.

[TaikoOsO]

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:

# 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()

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.

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?

[dokken]

keep_diagonal should not fix this issue by it-self, as it only keeps zero-entries in the sparsity pattern (ref: Bitbucket)

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

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 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()

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…