Integrating an expression over an internal boundary

[sbhasan]

Hi,

I am using dolfinx 0.4.1 installed through conda. My geometry is a square box containing two concentric circles inside. I want to integrate an algebraic expression (which will involve the solution as well) over the circumference of each circle i.e. internal boundaries. The example I am trying to follow is from the following thread:

How to define right hand side L as integral over smaller circle?

There are a few issues in your code. There is a difference between interior facets (a facet connected to two cells), and an exterior facet, (a facet connected to one cell). An interior facet integral is done by using the "dS"-measure, while an exterior facet is using the "ds"-measure. Your mesh is quite coarse, and this the circle is not exactly at the radius you are trying to use for marking. You can get the correct tolerance by reducing the tolerance. In the following code I’ve fixed these …

Therefore, my code looks like this

dS = ufl.Measure(“dS”, model.mesh, subdomain_data=boundary_index)
S = dolfinx.fem.assemble(dS) # test by evaluating the circumference

However, I could not find the counterpart of MeshFunction() used in the referred thread in the new dolfinx I am using. Moreover, I also could not understand how to distinguish between the circumferences of the two circles separately. I would appreciate if someone could show me how it is done.

Thanks a lot

[dokken]

You should use MeshTags. See for instance:
https://jorgensd.github.io/dolfinx-tutorial/chapter3/subdomains.html?highlight=meshtags#convert-msh-files-to-xdmf-using-meshio
and
https://jorgensd.github.io/dolfinx-tutorial/chapter3/robin_neumann_dirichlet.html

[sbhasan]

Thanks for your quick reply!

In my version (0.4.1) dolfinx.MeshTags does not appear to be a valid attribute. Further, I am sadly still very dumb about the finer details. Say I have the ct and ft objects that contain mesh and edge information. How would I feed either of them into ufl.Measure or fem.assemble so as to perform the integration?

[nate]

There was an interface change back in 0.3.0 or so:

In [3]: dolfinx.mesh.meshtags?                                                  
Signature:
dolfinx.mesh.meshtags(
    mesh: 'Mesh',
    dim: 'int',
    indices: 'np.ndarray',
    values: 'np.ndarray',
) -> 'MeshTagsMetaClass'
Docstring:
Create a MeshTags object that associates data with a subset of mesh entities.

Args:
    mesh: The mesh
    dim: Topological dimension of the mesh entity
    indices: Entity indices (local to process)
    values: The corresponding value for each entity

Returns:
    A MeshTags object

Note:
    The type of the returned MeshTags is inferred from the type of
    ``values``.

[sbhasan]

Thanks a lot Nate for the heads up!

So now this is how I am trying to do it:

# we already have ft instance computed as described in the tutorial shared above
bnd_tag = 6 # physical tag assigned to the boundary/edge we want to compute the integral on
arc_indices = ft.indices[ft.values == bnd_tag]
arc_tags = dolfinx.mesh.meshtags(model.mesh, 1, arc_indices, 6)
dS = ufl.Measure("dS", domain=model.mesh, subdomain_data=arc_tags)
print(dolfinx.fem.assemble.assemble_scalar(dS))

In the print statement, I used assemble_scalar as I could not find assemble() function itself. However, this too gives me the following error:

Traceback (most recent call last):
File “/install/miniconda3/envs/fenics/lib/python3.9/site-packages/IPython/core/interactiveshell.py”, line 3398, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File “”, line 1, in <cell line: 1>
dolfinx.cpp.fem.assemble_scalar(dS)
TypeError: assemble_scalar(): incompatible function arguments. The following argument types are supported:
1. (M: dolfinx::fem::Form, constants: numpy.ndarray[numpy.float64], coefficients: Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float64]]) → float
2. (M: dolfinx::fem::Form, constants: numpy.ndarray[numpy.float32], coefficients: Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float32]]) → float
3. (M: dolfinx::fem::Form<std::complex >, constants: numpy.ndarray[numpy.complex128], coefficients: Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.complex128]]) → complex
Invoked with: Measure(‘interior_facet’, subdomain_id=‘everywhere’, domain=Mesh(VectorElement(FiniteElement(‘Lagrange’, triangle, 1, variant=‘equispaced’), dim=2, variant=‘equispaced’), 0), subdomain_data=<dolfinx.mesh.MeshTagsMetaClass object at 0x7f132a6e71d0>)
Did you forget to #include <pybind11/stl.h>? Or <pybind11/complex.h>,
<pybind11/functional.h>, <pybind11/chrono.h>, etc. Some automatic
conversions are optional and require extra headers to be included
when compiling your pybind11 module.

Obviously assemble_scalar was a guess, trying assemble_vector also leads to some other error. I would appreciate if anyone could pull me out of this

[dokken]

I would suggest you read through the tutorial: The FEniCSx tutorial — FEniCSx tutorial
as it explains how to use assemble_scalar, assemble_vector etc, as you need to call dolfinx.fem.form(1*dS) prior to inputting into any of the aforementioned functions.

[sbhasan]

That was a super comment! Thank you so much for your tutorials and being ever so ready to direct to them!

I am able to calculate integrals in individual subdomains as well as on facets now!