# Hyperelastic deformation with applied pressue on internal surface

**URL:** <https://fenicsproject.discourse.group/t/hyperelastic-deformation-with-applied-pressue-on-internal-surface/14589>\
**Category:** General\
**Tags:** mesh, boundary-conditions, variational-form\
**Created:** [May 9, 2024, 4:45am UTC](https://fenicsproject.discourse.group/t/hyperelastic-deformation-with-applied-pressue-on-internal-surface/14589 "2024-05-09T04:45:08Z")\
**Posts on this page:** 1\
**Showing post:** 2

<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:** [May 9, 2024, 4:50am UTC](https://fenicsproject.discourse.group/t/hyperelastic-deformation-with-applied-pressue-on-internal-surface/14589/2 "2024-05-09T04:50:25Z")

</div>

> [@Chris\_S](#):
>
> ```auto
> T_intr = P_intr * n # Internal traction due to pressure
> F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T_intr) * ds(2)
> 
> ```

`ds` is an integration measure over exterior facets, thus those that are connected to only one cell. By adding a new domain in the circle this becomes and interior facet.

I would follow the notes in:

> <https://github.com/jorgensd/dolfinx-tutorial/issues/158>
>
> \`\`\`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

on how to do a one-sided interior facet integral.

---

_[View the full topic](https://fenicsproject.discourse.group/t/hyperelastic-deformation-with-applied-pressue-on-internal-surface/14589)._