Mixing QUGaR and FEniCSx assemblies: when to call Mat.assemble()?

dolfinx
1

Hello everyone,

I am currently working with a library called QUGaR (GitHub - pantolin/qugar: Quadratures for Unfitted GeometRies), which allows the generation of quadrature rules for immersed problems.

In my formulation, I have a volumetric term that relies on QUGaR, but I also have a penalty term defined on internal facets that I would like to implement in the classical way using ufl.dS.

The issue is that I cannot mix the QUGaR form (a) and the standard FEniCSx form (a2) in a single assembly, so I have to assemble them separately.

Here is how I did it. I was wondering whether I need to call self._A.assemble() after each dolfinx.petsc.assemble, or only once after all the dolfinx.petsc.assemble calls.

_init_

_a = form_custom(
        a,
        dtype=PETSc.ScalarType,  # type: ignore
        form_compiler_options=form_compiler_options,
        jit_options=jit_options,
)
self._a = typing.cast(CustomForm | list[CustomForm], _a)

self._A = create_matrix(self._a)  # type: ignore
_L = form_custom(
        L,
        dtype=PETSc.ScalarType,  # type: ignore
        form_compiler_options=form_compiler_options,
        jit_options=jit_options,
 )
 self._L = typing.cast(CustomForm | list[CustomForm], _L)

    # Assemble additional FEniCSx bilinear form if provided
  self._a2 = None
  if a2 is not None:
        self._a2 = dolfinx.fem.form(
            a2,
            dtype=PETSc.ScalarType,  # type: ignore
            form_compiler_options=form_compiler_options,
            jit_options=jit_options,
  )

solve

    self._A.zeroEntries()

    self.A_coeffs = _pack_coefficients(self._a)
    dolfinx.fem.petsc.assemble_matrix(self._A, self._a, bcs=self.bcs, coeffs=self.A_coeffs)
    
    if self._a2 is not None:
        # Assemble a2 directly into self._A (adds to existing values)
        dolfinx.fem.petsc.assemble_matrix(self._A, self._a2, bcs=self.bcs)
    
    # Finalize matrix assembly (only once, after all contributions)
    self._A.assemble()

form_custom and _pack_coefficient are fonctions provided by qugar who allow to compute the term for cutted elements.

2

I’ve not used QUGaR, but A.assemble() is typically the last thing you call to complete the matrix set-up. A small MWE with analytical expression would give you your definitive answer though…

3

Sometimes one needs to “flush” assembly to go from add mode to insert mode and vice versa (MatAssemblyBegin — PETSc v3.24.3-356-g763f2bca3976 documentation)
This is for instance done in:
dolfinx_mpc/python/src/dolfinx_mpc/numba/assemble_matrix.py at 54fece302ee3c46f0f6af078f8fb7b494be73492 · jorgensd/dolfinx_mpc · GitHub