# Mesh moving / ALE in DOLFINx — example or official API

**URL:** <https://fenicsproject.discourse.group/t/mesh-moving-ale-in-dolfinx-example-or-official-api/18323>\
**Category:** General\
**Tags:** mesh, dolfinx\
**Created:** [October 30, 2025, 2:29am UTC](https://fenicsproject.discourse.group/t/mesh-moving-ale-in-dolfinx-example-or-official-api/18323 "2025-10-30T02:29:20Z")\
**Posts on this page:** 1\
**Showing post:** 9

<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:** [November 7, 2025, 2:32pm UTC](https://fenicsproject.discourse.group/t/mesh-moving-ale-in-dolfinx-example-or-official-api/18323/9 "2025-11-07T14:32:13Z")

</div>

Using a “+” restriction is not correct on an interior facet (as it is arbitrary).  
To do a proper, one-sided, oriented integral, see for instance:

> <https://github.com/jorgensd/dolfinx-tutorial/issues/158#issue-2015212119>
>
> \`\`\`python
> \# Computing a onesided integral over an interior facet with consisten…t orientation
> \# Enabled by: https://github.com/FEniCS/dolfinx/pull/2269
> \# Copyright 2023 Jørgen S. Dokken
> \# SPDX: MIT
> 
> import dolfinx
> import numpy as np
> import ufl
> from mpi4py import MPI
> 
> mesh = dolfinx.mesh.create\_unit\_square(
> MPI.COMM\_WORLD, 10, 10, ghost\_mode=dolfinx.mesh.GhostMode.shared\_facet)
> 
> \# Create connectivties required for defining integration entities
> tdim = mesh.topology.dim
> fdim = tdim - 1
> mesh.topology.create\_entities(fdim)
> mesh.topology.create\_connectivity(fdim, tdim)
> mesh.topology.create\_connectivity(tdim, fdim)
> 
> 
> \# Get number of cells on process
> cell\_map = mesh.topology.index\_map(tdim)
> num\_cells = cell\_map.size\_local + cell\_map.num\_ghosts
> 
> \# Create markers for each size of the interface
> cell\_values = np.ones(num\_cells, dtype=np.int32)
> cell\_values\[dolfinx.mesh.locate\_entities(
> mesh, tdim, lambda x: x\[0\] \<= 0.5+1e-13)\] = 2
> ct = dolfinx.mesh.meshtags(mesh, tdim, np.arange(
> num\_cells, dtype=np.int32), cell\_values)
> 
> facet\_map = mesh.topology.index\_map(fdim)
> num\_facets = facet\_map.size\_local + facet\_map.num\_ghosts
> 
> \# Create facet markers
> facet\_values = np.ones(num\_facets, dtype=np.int32)
> facet\_values\[dolfinx.mesh.locate\_entities(
> mesh, fdim, lambda x: np.isclose(x\[0\], 0.5))\] = 2
> ft = dolfinx.mesh.meshtags(mesh, fdim, np.arange(
> num\_facets, dtype=np.int32), facet\_values)
> 
> \# Give a set of facets marked with a value (in this case 2), get a consistent orientation for an interior integral
> facets\_to\_integrate = ft.find(2)
> 
> f\_to\_c = mesh.topology.connectivity(fdim, tdim)
> c\_to\_f = mesh.topology.connectivity(tdim, fdim)
> \# Compute integration entities for a single facet of a cell.
> \# Each facet is represented as a tuple (cell\_index, local\_facet\_index), where cell\_index is local to process
> \# local\_facet\_index is the local indexing of a facet for a given cell
> integration\_entities = \[\]
> for i, facet in enumerate(facets\_to\_integrate):
> # Only loop over facets owned by the process to avoid duplicate integration
> if facet \>= facet\_map.size\_local:
> continue
> # Find cells connected to facet
> cells = f\_to\_c.links(facet)
> # Get value of cells
> marked\_cells = ct.values\[cells\]
> # Get the cell marked with 2
> correct\_cell = np.flatnonzero(marked\_cells == 2)
> 
> assert len(correct\_cell) == 1
> # Get local index of facet
> local\_facets = c\_to\_f.links(cells\[correct\_cell\[0\]\])
> local\_index = np.flatnonzero(local\_facets == facet)
> assert len(local\_index) == 1
> 
> # Append integration entities
> integration\_entities.append(cells\[correct\_cell\[0\]\])
> integration\_entities.append(local\_index\[0\])
> 
> \# Create custom integration measure for one-sided integrals
> breakpoint()
> ds = ufl.Measure("ds", domain=mesh, subdomain\_data=\[
> (8, np.asarray(integration\_entities, dtype=np.int32))\])
> n = ufl.FacetNormal(mesh)
> x = ufl.SpatialCoordinate(mesh)
> \# Exact integral is \[y/2\*\*2\]\_0^1= 1/2
> L = ufl.dot(ufl.as\_vector((x\[1\], 0)), n)\*ds(8)
> L\_compiled = dolfinx.fem.form(L)
> 
> print(
> f"Correct integral: {mesh.comm.allreduce(dolfinx.fem.assemble\_scalar(L\_compiled), op=MPI.SUM)}")
> 
> 
> \# Create reference implementation where we use a restricted two-sided integral with no notion of orientation
> dS = ufl.Measure("dS", domain=mesh, subdomain\_data=ft, subdomain\_id=2)
> n = ufl.FacetNormal(mesh)
> x = ufl.SpatialCoordinate(mesh)
> L2 = ufl.dot(ufl.as\_vector((x\[1\], 0)), n("+"))\*dS
> L2\_compiled = dolfinx.fem.form(L2)
> print(
> f"Wrong integral: {mesh.comm.allreduce(dolfinx.fem.assemble\_scalar(L2\_compiled), op=MPI.SUM)}")
> \`\`\`
> Source: https://fenicsproject.discourse.group/t/wrong-facetnormal-vector-on-internal-boundaries/12887/2?u=dokken

and

> [@Correct way to compute internal interface fluxes](https://fenicsproject.discourse.group/t/correct-way-to-compute-internal-interface-fluxes/18343/2):
>
> The key take-away is that: In most finite element methods, the jump of a quantity across a surface does not require a specific orientation. For instance, in DG methods, you get jump integrals across all internal facets of the grids, but mathematically the formulation doesn’t care which cell is positive and which one is negative. DOLFINx follows the same notion, that for an interior facet integral (dS) which cell comes first (+) and which one is last (-) doesn’t matter if put in the correct va…

---

_[View the full topic](https://fenicsproject.discourse.group/t/mesh-moving-ale-in-dolfinx-example-or-official-api/18323)._