2D and 1D variational terms in 3D simulation

variational formulation
, ,
1

I am trying to solve a following problem:

\int_\Omega \nabla u \nabla v + \int_{\Gamma_i} 100 \nabla_t u \nabla_t v = \int_{\partial \Gamma_i} f(u) + \int_{\partial \Omega} v q n

First integral is the standard integral over the domain. Second integral is over a surface within the domain \Omega, mesh elements are aligned with the surface and \nabla_t is a 2D gradient within the surface. Third integral is a boundary condition on the internal surface’s boundary and fourth is a standard Neumann BC on the domain’s boundary.

When I include the second term in my simulation, there is no effect on the results. What could be the possible problem?

As far as I know, currently there is no easy option to include 1D integral as ufl. Measure only supports dim-1 integrals. Can it be done with the “custom” Measure type?

Following is the code I am using:

from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, 
                         locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities, create_submesh

from ufl import (TestFunction, TrialFunction, FacetNormal,Measure, Dx,
                 dx, grad, inner)
from mpi4py import MPI


with XDMFFile(MPI.COMM_WORLD, "mesh/t10_tetra.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
with XDMFFile(MPI.COMM_WORLD, "mesh/t10_surf.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim - 2)
with XDMFFile(MPI.COMM_WORLD, "mesh/t10_line.xdmf", "r") as xdmf:
    lt = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

left_marker=20006
right_marker=20004
V = functionspace(mesh, ("Lagrange", 1))
right_facets = ft.find(right_marker)
right_dofs = locate_dofs_topological(V, mesh.topology.dim - 1, right_facets)
bcs = [dirichletbc(default_scalar_type(0), right_dofs, V)]

ds = Measure("ds", domain=mesh, subdomain_data=ft)
u, v = TrialFunction(V), TestFunction(V)
n = FacetNormal(mesh)
a = inner(grad(u), grad(v)) * dx + default_scalar_type(100)*(inner(Dx(u,0),Dx(v,0)))*ds(40001) ###second term has no effect
L = Constant(mesh, default_scalar_type(0)) * v * dx + inner(v,default_scalar_type(100))*ds(left_marker)#  

ds1 = Measure("ds", domain=mesh, subdomain_data=lt) 
L2 = 0.01*inner(v,default_scalar_type(100))*ds1(10001) ### ds only works for dim-1

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

with XDMFFile(mesh.comm, "out.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(uh)

and the mesh can be found here:
https://mega.nz/folder/vvxGXKQI#L8VEvztkVvdi1IF5UZrlMQ
Marker 40001 is the surface inside the boundary.

2

You are then using the wrong integration-measure

means the exterior facets, i.e. those only connected to a single cell.
What you want to use is Measure("dS", domain=mesh, subdomain_data=ft) to integrate over interior facets.

3

Thank you for the fast reply. I should have added that I tried this solution and I get an error:

ValueError: Discontinuous type Argument must be restricted.

from within ufl/algorithms/apply_restrictions.py, when setting up the solver.

The fragment looks like this:

ds = Measure("dS", domain=mesh, subdomain_data=ft)
u, v = TrialFunction(V), TestFunction(V)
n = FacetNormal(mesh)
a = inner(grad(u), grad(v)) * dx + default_scalar_type(100)*(inner(Dx(u,0),Dx(v,0)))*ds(40001) ###second term has no effect
L = Constant(mesh, default_scalar_type(0)) * v * dx
4

That does make sense, as Dx(u, 0) is not continuous across the surface. See for instance: Add one-sided integration example as example of custom integration · Issue #158 · jorgensd/dolfinx-tutorial · GitHub
on how to do a one-sided integration.
If you do not want a one-sided integration, you would need to use jump or avg operators.

5

Thank you, I will check the link and try to make sense of it, as currently I don’t quite understand why would the continuity be important :slight_smile:

Can I also ask you to comment on the second issue, with the 1D integral?

6

The derivative of u in the normal direction is not continuous at an inter-element interface.
Therefore you do not get a single result by calling Dx(u, 0) at the interface between 2 boundaries.

For 1D integrals (ridge integrals) within a 3D geometry; This is a work in progress at Codim 2 coupling by jorgensd · Pull Request #731 · FEniCS/ffcx · GitHub

7
  1. I think I understand, you mean that on the two sides of the integration domain gradients are different. But I should clarify that here we do not have this problem, as the gradients are calculated with respect to tangential components of the integration domain (and mesh elements). In this simplified example, all elements on the internal surface lie in the plane XY, meaning that derivatives of u with respect to x and y are independent of the normal component. For continuous FEM they are the same on both sides of the surface. In the general case mesh elements are also always aligned with the internal surface and the normal component of the gradient is then easily removed with the integral defined as:

\int_{\Gamma_i} (\nabla u - n (dot(n,\nabla u)))*(\nabla v - n (dot(n,\nabla v)) ds

  1. This is great news, I will try to deal with the other features of my problem while the functionality is implemented.
9

Note that UFL is a symbolic language, at the time of creating the expression there is no notion of what is the tangential direction at the surface. If you know that the expression is continuous, you can choose either side, by adding a (“+”) or (“-”) restriction to the relevant expressions

10

This is exactly what I was looking for, thank you very much for all the help!

The final script:

from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, 
                         locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities, create_submesh

from ufl import (TestFunction, TrialFunction, FacetNormal,Measure, Dx,
                 dx, grad, inner)
from mpi4py import MPI

with XDMFFile(MPI.COMM_WORLD, "mesh/t10_tetra.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")
tdim = mesh.topology.dim
fdim = mesh.topology.dim - 1
ldim = mesh.topology.dim - 2
mesh.topology.create_connectivity(tdim, fdim)
with XDMFFile(MPI.COMM_WORLD, "mesh/t10_surf.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(fdim, ldim)
with XDMFFile(MPI.COMM_WORLD, "mesh/t10_line.xdmf", "r") as xdmf:
    lt = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(fdim, tdim)

left_marker=20006
right_marker=20004
mid_marker=40001
V = functionspace(mesh, ("Lagrange", 1))
right_facets = ft.find(right_marker)
right_dofs = locate_dofs_topological(V, fdim, right_facets)
bcs = [dirichletbc(default_scalar_type(0), right_dofs, V)]

# dS for internal surface integration
dS = Measure("dS", domain=mesh, subdomain_data=ft)
# ds for boundary integration
ds = Measure("ds", domain=mesh, subdomain_data=ft)
u, v = TrialFunction(V), TestFunction(V)
n = FacetNormal(mesh)
#
a = inner(grad(u), grad(v)) * dx + default_scalar_type(100)*(inner(grad(u)('+')-n('+')*inner(grad(u)('+'),n('+')),grad(v)('+')-n('+')*inner(grad(v)('+'),n('+'))))*dS(mid_marker)
L = Constant(mesh, default_scalar_type(0)) * v * dx - inner(v,default_scalar_type(100))*ds(left_marker)# 

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

with XDMFFile(mesh.comm, "out.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(uh)