Integration Measure over small dimension

Hello community,

I have two questions, which are the same but in two different levels.
Let suppose a cubic domain \Omega, and a subdomain of it : \Gamma, corresponding to 3 joined surfaces of the cube.
Let (p, u) two trial functions defined respectively on \Omega and \Gamma. And (v, w) the test functions, also defined respectively on \Omega and \Gamma.
My first question concerns the integration by part and its implementation of such a term:
\int_\Gamma \partial_\tau^2u\cdot w = \int_{\partial\Gamma}\partial_\tau u\cdot w-\int_\Gamma \partial_\tau u\cdot \partial_\tau w
As u and w are defined in the same function space, of the same domain. The integration over the subdomain \Gamma is performed by

ns = CellNormal(submesh)
dx1 = Measure("dx", domain=submesh)
d_tau_u = tangential_proj(grad(u), ns)
d_tau_w = tangential_proj(grad(w), ns)
a1 = inner(d_tau_u, d_tau_w) * dx1

And my question is, can the border term (I’m not sure about the english translation of the first term) \int_{\partial\Gamma}\partial_\tau u\cdot w be computed with a ds1 = Measure("ds", domain=submesh) ? I’m pretty sure it’s a good way to do it but I would really enjoy a confirmation from you… it would give something like :

ns = CellNormal(submesh)
ds1 = Measure("ds", domain=submesh)
d_tau_u = tangential_proj(grad(u), ns)
a2 = inner(d_tau_u, w) * ds1

Now, in a more complexe way. For the same integration by part as previously, but this time with a mix of the functions :
\int_\Gamma \partial_\tau^2p\cdot w = \int_{\partial\Gamma}\partial_\tau p\cdot w-\int_\Gamma \partial_\tau p\cdot \partial_\tau w
The fact that the terms involve p, the integration operant must be defined on the mesh. So once again, the term \int_\Gamma \partial_\tau p\cdot \partial_\tau w is computed with:

n = FacetNormal(mesh)
ds = Measure("ds", domain=mesh, subdomain_data=mesh_bc_tags)
d_tau_p = tangential_proj(grad(p), n)
d_tau_w = tangential_proj(grad(w), n)
a1 = inner(d_tau_p, d_tau_w) * ds(3) # the tag 3 corresponds to the joined surfaces

and the question is, how do I obtain a Measure over the boundary of the boudary… I thought about using dP = Measure("dP", domain=mesh, subdomain_data=mesh_bc_tags), but how can I be sure I’m integration only on the verties of the tagged surface and not over the whole domain ? Should firstly I tag the verties I’ll integrate over ?
For the term \int_{\partial\Gamma}\partial_\tau p\cdot w, would something like this work :

n = FacetNormal(mesh)
dP = Measure("dP", domain=mesh, subdomain_data=mesh_bc_tags)
d_tau_p = tangential_proj(grad(p), n)
a2 = inner(d_tau_p, w) * dP(3) 

The question is not asked in the proper way for this forum. Do not hesitate to ask question for additional code, I’m not sure if a MWE would be relevant, but I can produce one if it’s needed.

Thank you very much.
Pierre.

May be relevant Integrate a field quantity along a line in 3D - #8 by Hoseong_Jeong .

I would be interested if you find a way to do line integrations. I’m not sure “dP“ can be used as it represents the vertices, not the edges…

First of all, you need to specify:

  1. Which version of DOLFIN you are using, (legacy, DOLFINx, and what version of any of those).
  2. How is submesh created, as there are various ways to do that, depending on the installation.
  3. By not providing a small, reproducible example, for instance with a cube and a surface of the cube, one has to make alot of assumptions on your definitions, such as tangetial_proj.

It seems like your question is boiling down to:

  • Can I integrate over the edges the “exterior” edges of a 3D mesh, which would be viewed as the exterior boundary of a 2D submesh?

whose short answer is, no; not at the moment.

There is progress towards this, in for instance:

and just recently @schnellerhase has extended DOLFINx to exterior vertex integrals, ref: Vertex integrals by schnellerhase · Pull Request #3726 · FEniCS/dolfinx · GitHub

Hi Dokken,

Thank you for the reply,
I work with the last version DOLFINx, v0.10.0.0. I create my geometry and the mesh with the python API of gmsh, I tag 3 joined face (e.g. for a cube, Top - Right - Back) and I extract from that a submesh with submesh, entity_map = msh.create_submesh(final_mesh, fdim, facet_tags.find(3))[0:2]. The tangetial_proj function is the tangential projection of a vector according to the normal vector of the surface, as proposed in a FEniCSx post that I don’t find anymore. Basically :

def tangential_proj(u, n):
    
    return (ufl.Identity(n.ufl_shape[0]) - ufl.outer(n, n)) * u

If you don’t mind, I keep this post open untilI a solution to integrate over the edge of a 3D mesh is found, and an exemple can be provided. I’ll update this topic with that.

Does it still faisable to integrate over a contour of a surface with Measure("ds", domain=submesh) ?

I’m working on a MWE to provide here.

Many thanks,
Pierre.

This should be possible, if you only have functons that live on the submesh, it can’t use quantities fromt he parent mesh, as of now.

I just merged the ridge integrals into FFCx. Will see if I can get a version supported in DOLFINx next week.

Very good to know it solves part of my problem. For the ridge integrals, good luck on that !

Once again, many thanks.

Hi again,
Here is a MWE that could work if it was possible to integrate over the edge of a surface of a volume domain. Actually the geometry is the one I use, and might not be relevant for everyone.. And I hope this is not too much case specific. The error provided by this code comes after.

from mpi4py import MPI
import gmsh

from dolfinx.io import gmshio
import ufl
from ufl import (
    TestFunction,
    TrialFunction,
    Measure,
    inner,
    grad,
)
from basix.ufl import element
from dolfinx.fem import functionspace, form, petsc

import gmsh
from mpi4py import MPI
import dolfinx.io.gmshio as gmshio

import gmsh
from mpi4py import MPI
import dolfinx.io.gmshio as gmshio
import dolfinx.mesh as dmesh
import numpy as np


def cubic_domain(side_box = 0.11, radius = 0.1, lc = 8e-3, model_name = "Cubic"):

    gmsh.initialize()
    comm = MPI.COMM_WORLD
    model_rank = 0
    model = gmsh.model()
    gmsh.model.add(model_name)

    # Definition of the points
    p1 = gmsh.model.geo.addPoint(0, 0, 0, lc)
    p2 = gmsh.model.geo.addPoint(side_box, 0, 0, lc)
    p3 = gmsh.model.geo.addPoint(side_box, side_box, 0, lc)
    p4 = gmsh.model.geo.addPoint(0, side_box, 0, lc)

    p5 = gmsh.model.geo.addPoint(0, 0, side_box, lc)
    p6 = gmsh.model.geo.addPoint(side_box, 0, side_box, lc)
    p7 = gmsh.model.geo.addPoint(side_box, side_box, side_box, lc)
    p8 = gmsh.model.geo.addPoint(0, side_box, side_box, lc)

    p9 = gmsh.model.geo.addPoint(side_box, radius, 0, lc)
    p10 = gmsh.model.geo.addPoint(side_box, 0, radius, lc)

    # Definition of the lines
    l1 = gmsh.model.geo.addLine(p1, p4)
    l2 = gmsh.model.geo.addLine(p4, p3)
    l3 = gmsh.model.geo.addLine(p3, p9)
    l4 = gmsh.model.geo.addLine(p9, p2)
    l5 = gmsh.model.geo.addLine(p2, p1)

    l6 = gmsh.model.geo.addLine(p5, p6)
    l7 = gmsh.model.geo.addLine(p6, p7)
    l8 = gmsh.model.geo.addLine(p7, p8)
    l9 = gmsh.model.geo.addLine(p8, p5)

    l10 = gmsh.model.geo.addLine(p2, p10)
    l11 = gmsh.model.geo.addLine(p10, p6)
    l12 = gmsh.model.geo.addLine(p3, p7)
    l13 = gmsh.model.geo.addLine(p4, p8)
    l14 = gmsh.model.geo.addLine(p1, p5)

    # definition of the quarter circle
    c15 = gmsh.model.geo.addCircleArc(p9, p2, p10)

    # Curve loops
    cl1 = [l1, l2, l3, l4, l5]
    cl2 = [l6, l7, l8, l9]
    cl3 = [-l3, l12, -l7, -l11, -c15]
    cl4 = [c15, -l10, -l4]
    cl5 = [-l2, l13, -l8, -l12]
    cl6 = [-l5, l10, l11, -l6, -l14]
    cl7 = [-l1, l14, -l9, -l13]

    # Surfaces
    s1 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl1)])
    s2 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl2)])
    s3 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl3)])
    s4 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl4)])
    s5 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl5)])
    s6 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl6)])
    s7 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl7)])

    # Surface loop
    sl1 = gmsh.model.geo.addSurfaceLoop([s1, s2, s5, s3, s4, s6, s7])

    # Definition of the volume
    v1 = gmsh.model.geo.addVolume([sl1])

    gmsh.model.geo.synchronize()

    # Definition of the physical identities
    gmsh.model.addPhysicalGroup(1, [l12, -l7, -l6, -l14, l1, l2], tag=111)
    gmsh.model.addPhysicalGroup(2, [s4], tag=1)  # Neumann
    gmsh.model.addPhysicalGroup(2, [s7, s5, s2], tag=3)
    gmsh.model.addPhysicalGroup(3, [v1], tag=1)  # Volume physique

    # Generation of the 3D mesh
    gmsh.model.mesh.generate(3)

    # Save the mesh
    gmsh.write(model_name + ".msh")

    # Affect the mesh, volumic tags and surfacic tag to variables
    mesh_data = gmshio.model_to_mesh(model, comm, model_rank)
    #print(f'mesh_data : {mesh_data}')
    final_mesh = mesh_data.mesh
    cell_tags = mesh_data.cell_tags
    facet_tags = mesh_data.facet_tags
    ridge_tags = mesh_data.ridge_tags   
    

    # Close gmsh
    gmsh.finalize()
    
    tdim = final_mesh.topology.dim
    fdim = tdim - 1

    xref = [side_box, 0, 0]

    

    submesh, entity_map = dmesh.create_submesh(final_mesh, fdim, facet_tags.find(3))[0:2]
    entity_maps_mesh = [entity_map]

    #mesh_info = [final_mesh, cell_tags, facet_tags, xref]
    mesh_info = [final_mesh, cell_tags, facet_tags, xref, ridge_tags] # this line is new
    submesh_info = [submesh, entity_maps_mesh]

    return mesh_info, submesh_info


# Example usage:
if __name__ == "__main__":
    # Generate a cube mesh and retrieve mesh data
    mesh_info, submesh_info = cubic_domain()
    
    final_mesh   = mesh_info[0]
    cell_tags    = mesh_info[1]
    facet_tags   = mesh_info[2]
    xref         = mesh_info[3]
    ridge_tags   = mesh_info[4]
    
    submesh          = submesh_info[0]
    entity_maps_mesh = submesh_info[1]

    # Define volume (dx) and surface (ds) measures
    dx = Measure("dx", domain=final_mesh)
    ds = Measure("ds", domain=final_mesh, subdomain_data=facet_tags)
    dr = Measure("dr", domain=final_mesh, subdomain_data=ridge_tags)
    dx1 = Measure("dx", domain=submesh)
    ds1 = Measure("ds", domain=submesh)

    # Define a Lagrange element of degree 3
    P1 = element("Lagrange", final_mesh.basix_cell(), 2)
    P = functionspace(final_mesh, P1)
    Q1 = element("Lagrange", submesh.basix_cell(), 2)
    Q = functionspace(submesh, Q1)

    # Define trial and test functions
    p, v = TrialFunction(P), TestFunction(P)
    u, w = TrialFunction(Q), TestFunction(Q)

    # Define the weak form
    a_00_0 = inner(p, v) * dx
    a_01_0 = inner(u, v) * ds(3)
    a_10_0 = inner(grad(p)[0] +grad(p)[1] + grad(p)[2], w) * dr(111)
    a_10_1 = inner(grad(p), grad(w)) * ds(3)
    a_11_0 = inner(grad(u)[0] +grad(u)[1] + grad(u)[2], w) * ds1
    a_11_1 = inner(grad(u), grad(w)) * dx1

    # Assemble the total form
    A_00 = form(a_00_0)
    A_01 = form(a_01_0, entity_maps=entity_maps_mesh)
    A_10 = form(a_10_0 + a_10_1, entity_maps=entity_maps_mesh)
    A_11 = form(a_11_0 + a_11_0)
    print("end")

The error :

  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/root/ABC_FEniCSx_classical/mwe_FEniCSx.py", line 179, in <module>
    A_10 = form(a_10_0 + a_10_1, entity_maps=entity_maps_mesh)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-dolfinx-main-iqvj2zu4tzntqaicsg2qyaehn6q5jb4r/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 438, in form
    return _create_form(form)
           ^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-dolfinx-main-iqvj2zu4tzntqaicsg2qyaehn6q5jb4r/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 430, in _create_form
    return _form(form)
           ^^^^^^^^^^^
  File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-dolfinx-main-iqvj2zu4tzntqaicsg2qyaehn6q5jb4r/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 339, in _form
    ufcx_form, module, code = jit.ffcx_jit(
                              ^^^^^^^^^^^^^
  File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-dolfinx-main-iqvj2zu4tzntqaicsg2qyaehn6q5jb4r/lib/python3.12/site-packages/dolfinx/jit.py", line 61, in mpi_jit
    return local_jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-dolfinx-main-iqvj2zu4tzntqaicsg2qyaehn6q5jb4r/lib/python3.12/site-packages/dolfinx/jit.py", line 216, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 226, in compile_forms
    raise e
  File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 206, in compile_forms
    impl = _compile_objects(
           ^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 331, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/compiler.py", line 113, in compile_ufl_objects
    ir = compute_ir(analysis, _object_names, _prefix, options, visualise)
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/ir/representation.py", line 189, in compute_ir
    _compute_integral_ir(
  File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/ir/representation.py", line 262, in _compute_integral_ir
    entity_type = _entity_types[itg_data.integral_type]
                  ~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^
KeyError: 'ridge'

The merge has been done, I then use the version 0.10.0 of dolfinx, so I post an update here to have some help on the use of the new Measure. My MWE hasn’t changes, but here is the error I have when I run the MWE code :

Traceback (most recent call last):
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/root/ABC_FEniCSx_classical/mwe_FEniCSx.py", line 179, in <module>
    A_10 = form(a_10_0 + a_10_1, entity_maps=entity_maps_mesh)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 439, in form
    return _create_form(form)
           ^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 431, in _create_form
    return _form(form)
           ^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 354, in _form
    ufcx_form, module, code = jit.ffcx_jit(
                              ^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/dolfinx/jit.py", line 61, in mpi_jit
    return local_jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/dolfinx/jit.py", line 216, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 226, in compile_forms
    raise e
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 206, in compile_forms
    impl = _compile_objects(
           ^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 331, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/ffcx/compiler.py", line 113, in compile_ufl_objects
    ir = compute_ir(analysis, _object_names, _prefix, options, visualise)
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/ffcx/ir/representation.py", line 170, in compute_ir
    _compute_integral_ir(
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/ffcx/ir/representation.py", line 408, in _compute_integral_ir
    integral_ir = compute_integral_ir(
                  ^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/ffcx/ir/integral.py", line 141, in compute_integral_ir
    mt_table_reference = build_optimized_tables(
                         ^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/ffcx/ir/elementtables.py", line 500, in build_optimized_tables
    get_ffcx_table_values(
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/ffcx/ir/elementtables.py", line 154, in get_ffcx_table_values
    tbl = component_element.tabulate(deriv_order, entity_points)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/basix/ufl.py", line 604, in tabulate
    tables = self._element.tabulate(nderivs, points)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/basix/ufl.py", line 452, in tabulate
    tab = self._element.tabulate(nderivs, points)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/fenicsx-env/.spack-env/view/lib/python3.12/site-packages/basix/finite_element.py", line 90, in tabulate
    return self._e.tabulate(n, x)
           ^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Point dim (1) does not match element dim (2).

Ridge integrals are not fully exposed in DOLFINx yet, as it is introduced in this PR (Unification of one-sided sub-entity integrals by jorgensd · Pull Request #3900 · FEniCS/dolfinx · GitHub).
However, the issue your are observing seems to be unrelated to that. I’ll try to boil it down to a minimal reproducible example.