Imposing a discontinuity in the solution, on an internal interface

General
1

In an attempt to solve Ohm’s law (as well as the common heat equation) in a domain with 2 subdomains with different material properties, I want to simulate an interface resistance (electrical and thermal). Physically, the interfacial electrical resistance is directly proportional to the drop in voltage across the interface, or equivalently in the thermal problem, the interfacial thermal resistance is directly proportional to the temperature drop at the interface.

I am wondering how I could impose a discontinuity in the solution, at a given (internal) surface (for 3D problems)/curve(2D problems). Neumann’s boundary condition won’t make it.

2

There are a few options:

  1. Use discontinuous elements in the whole mesh, and use a DG formulation.
  2. Use the Meshview and mixed assembly framework by Cecile Daversin: GitHub - cdaversin/mixed-dimensional-examples: Code to reproduce numerical examples presented in "Abstractions and automated algorithms for mixed-dimensional finite element methods" (2019)
  3. create two separate meshes and use DOLFINx_MPC or some similar framework to enforce a jump at these boundaries,
1 Like
3

Thank you very much @dokken. After some inspection, I think I’ll attempt to go for 1).

I have found some old lecture notes regarding an implementation of DG: https://fenicsproject.org/pub/course/lectures/2013-11-logg-fcc/lecture_10_discontinuous_galerkin.pdf. I tried to adapt it to my code (I had to change CellSize by CellDiameter), unfortunately the solution I get is totally unchanged by the extra term I add in the weak form.
From a working code, I added:

from ufl import avg, FacetNormal, CellDiameter, jump
n = FacetNormal ( mesh )
h = CellDiameter ( mesh )
dS = Measure("dS", domain=mesh, subdomain_data=ft)
discontinuous_term = 10.0 / avg(h) * dot ( jump (u , n ) , jump (volt , n ) ) * dS
F = dot(grad(u), sigma_tensor * grad(volt))*dx + u * J * ds(inj_current_curve)# + discontinuous_term

Where “volt” is the unknown function I am solving for, and “u” is the test function. Am I doing something wrong?

My full code, up to that point, is:

from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh("meshes/rectangles_smaller_contacts_with_internal_interface.msh", MPI.COMM_WORLD, gdim=2)
proc = MPI.COMM_WORLD.rank 
import meshio
def create_mesh(mesh, cell_type, prune_z=True):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

if proc == 0:
    # Read in mesh
    msh = meshio.read("meshes/rectangles_smaller_contacts_with_internal_interface.msh")
   
    # Create and save one file for the mesh, and one file for the facets 
    triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
    line_mesh = create_mesh(msh, "line", prune_z=True)
    meshio.write("meshes/mesh_smaller_contacts_with_internal_interface.xdmf", triangle_mesh)
    meshio.write("meshes/mt_smaller_contacts_with_internal_interface.xdmf", line_mesh)
# Now read the mesh files produced above.
from dolfinx.io import XDMFFile
with XDMFFile(MPI.COMM_WORLD, "meshes/mesh_smaller_contacts_with_internal_interface.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, "meshes/mt_smaller_contacts_with_internal_interface.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological)
from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
                 dx, grad, inner, as_tensor, inv)
import numpy as np
from ufl import TensorElement
from dolfinx.fem import TensorFunctionSpace, FunctionSpace
# Define a tensor element (also possible to define a tensor space)
# Dokken's solution below
tensor_el = TensorElement("DG", mesh.ufl_cell(), 0, shape=(2, 2))

T = FunctionSpace(mesh, tensor_el)
rho = Function(T)
sigma_tensor = Function(T)

def rho_A(x):
    tensor = np.array([[1, 0], [0, 1.]])
    values = np.repeat(tensor, x.shape[1])
    return values.reshape(tensor.shape[0]*tensor.shape[1], x.shape[1])


def rho_B(x):
    tensor = np.array([[10, 0], [0, 10.]])
    values = np.repeat(tensor, x.shape[1])
    return values.reshape(tensor.shape[0]*tensor.shape[1], x.shape[1])

def sigma_A(x):
    tensor = np.linalg.inv(np.array([[1, 0], [0, 1.]]))
    values = np.repeat(tensor, x.shape[1])
    return values.reshape(tensor.shape[0]*tensor.shape[1], x.shape[1])

def sigma_B(x):
    tensor = np.linalg.inv(np.array([[1, 0], [0, 1.]]))
    values = np.repeat(tensor, x.shape[1])
    return values.reshape(tensor.shape[0]*tensor.shape[1], x.shape[1])

rho.interpolate(rho_A, cells=ct.find(11))
rho.interpolate(rho_B,cells=ct.find(12))

sigma_tensor.interpolate(sigma_A, cells=ct.find(11))
sigma_tensor.interpolate(sigma_B,cells=ct.find(12))

V = FunctionSpace(mesh, ("Lagrange", 1))
volt = Function(V)
u = TestFunction(V)
from ufl import avg, FacetNormal, CellDiameter, jump
#from dolfinx.functions.specialfunctions import CellSize
n = FacetNormal ( mesh )
h = CellDiameter ( mesh )
# Reminder of the mesh:
# Physical Surface("material A", 11) = {2};
# Physical Surface("material B", 12) = {1};
# Physical Curve("left_mat_A", 13) = {8};
# Physical Curve("right_mat_B", 14) = {2};
# Physical Curve("Internal interface", 34) = {25};
inj_current_curve = 13
vanish_voltage_curve = 14 # This corresponds to the curve the current leaves the material.
V_current_contact_out = 0.0 # Voltage value of the curve where the current leaves the material.
reading_voltage_surface_0 = 13
reading_voltage_surface_1 = 14
surface_Joule_effect_internal_interface = 34
# We define the boundary conditions
from petsc4py.PETSc import ScalarType
u_bc = Function(V)
#left_facets = ft.find(inj_current_curve)
left_facets = ft.find(vanish_voltage_curve)
left_dofs = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
bcs = [dirichletbc(ScalarType(V_current_contact_out), left_dofs, V)]
from dolfinx.fem import assemble_scalar, form
from ufl import Measure
dx = Measure("dx", domain=mesh,subdomain_data=ct)
ds = Measure("ds", domain=mesh, subdomain_data=ft)
dS = Measure("dS", domain=mesh, subdomain_data=ft)
the_current = 2.0 # Current, in amperes.
J = the_current / assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))
print(J)
print('Length of curve where current is injected', assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh))))
print('Length of curve where current leaves the wire', assemble_scalar(form(1 * ds(vanish_voltage_curve, domain=mesh))))
print('Length of curve of the interface', assemble_scalar(form(1 * dS(surface_Joule_effect_internal_interface, domain=mesh))))
from ufl import dot
# Weak form.
discontinuous_term = 10.0 / avg(h) * dot ( jump (u , n ) , jump (volt , n ) ) * dS(surface_Joule_effect_internal_interface)
F = dot(grad(u), sigma_tensor * grad(volt))*dx + u * J * ds(inj_current_curve) + discontinuous_term

Where my file.msh is

$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
5
1 13 "left_mat_A"
1 14 "right_mat_B"
1 34 "Internal interface"
2 11 "material A"
2 12 "material B"
$EndPhysicalNames
$Entities
10 11 2 0
15 0 0 0 0 
18 0 2 0 0 
19 -1 0 0 0 
20 -1 2 0 0 
21 -1 1.4 0 0 
22 -1 1.1 0 0 
23 1 0 0 0 
24 1 0.1 0 0 
25 1 0.2 0 0 
26 1 2 0 0 
20 -1.0000001 1.0999999 -1e-07 -0.9999999000000001 1.4000001 1e-07 1 13 2 21 -22 
21 0.9999999000000001 0.0999999 -1e-07 1.0000001 0.2000001 1e-07 1 14 2 25 -24 
25 -1e-07 -9.999999994736442e-08 -1e-07 1e-07 2.0000001 1e-07 1 34 2 18 -15 
26 -1.0000001 -1e-07 -1e-07 9.999999994736442e-08 1e-07 1e-07 0 2 19 -15 
27 -1.0000001 1.9999999 -1e-07 9.999999994736442e-08 2.0000001 1e-07 0 2 18 -20 
28 -1.0000001 1.3999999 -1e-07 -0.9999999000000001 2.0000001 1e-07 0 2 20 -21 
29 -1.0000001 -9.999999994736442e-08 -1e-07 -0.9999999000000001 1.1000001 1e-07 0 2 22 -19 
30 -9.999999994736442e-08 -1e-07 -1e-07 1.0000001 1e-07 1e-07 0 2 15 -23 
31 0.9999999000000001 -1.000000000028756e-07 -1e-07 1.0000001 0.1000001 1e-07 0 2 23 -24 
32 0.9999999000000001 0.1999999 -1e-07 1.0000001 2.0000001 1e-07 0 2 25 -26 
33 -9.999999994736442e-08 1.9999999 -1e-07 1.0000001 2.0000001 1e-07 0 2 26 -18 
1 -9.999999994736442e-08 -9.999999994736442e-08 -1e-07 1.0000001 2.0000001 1e-07 1 12 6 30 31 -21 32 33 25 
2 -1.0000001 -9.999999994736442e-08 -1e-07 9.999999994736442e-08 2.0000001 1e-07 1 11 6 26 -25 27 28 20 29 
$EndEntities
$Nodes
20 68 1 68
0 15 0 1
1
0 0 0
0 18 0 1
2
0 2 0
0 19 0 1
3
-1 0 0
0 20 0 1
4
-1 2 0
0 21 0 1
5
-1 1.4 0
0 22 0 1
6
-1 1.1 0
0 23 0 1
7
1 0 0
0 24 0 1
8
1 0.1 0
0 25 0 1
9
1 0.2 0
0 26 0 1
10
1 2 0
1 20 0 0
1 21 0 0
1 25 0 7
11
12
13
14
15
16
17
0 1.75 0
0 1.5 0
0 1.25 0
0 0.9999999999999998 0
0 0.75 0
0 0.5 0
0 0.25 0
1 26 0 3
18
19
20
-0.75 0 0
-0.4999999999999999 0 0
-0.25 0 0
1 27 0 3
21
22
23
-0.2500000000000001 2 0
-0.5000000000000001 2 0
-0.75 2 0
1 30 0 3
24
25
26
0.2500000000000001 0 0
0.5000000000000001 0 0
0.75 0 0
1 32 0 1
27
1 1.35290261948941 0
1 33 0 3
28
29
30
0.75 2 0
0.4999999999999999 2 0
0.25 2 0
2 1 0 16
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
0.2165063509461096 0.625 0
0.2129253918967202 1.627990015724681 0
0.6847799766951608 1.782467503425071 0
0.2245557494800761 1.115829697898795 0
0.3761500219765536 0.2122076945141456 0
0.6252300043953107 0.1824415389028292 0
0.4882046775611123 0.4333071134732195 0
0.530637335813078 0.8279338736844883 0
0.2104334777704233 0.3670424373543549 0
0.1943398872478528 0.8637527143166566 0
0.3855208408169349 1.795130069965427 0
0.5114034636855107 1.521164648922431 0
0.5172113347772932 1.22568945855073 0
0.2443493233066001 1.373445636849439 0
0.1817398161387643 0.1817398161387643 0
0.1817398161387641 1.818260183861236 0
2 2 0 22
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
-0.2023532308534801 1.370696913165413 0
-0.6249999999999999 0.2165063509461096 0
-0.208019165177519 0.629230046878511 0
-0.6250000000000001 1.78349364905389 0
-0.6491443534004766 1.231401242067281 0
-0.6058347963465431 0.8513089777016829 0
-0.3978471061286779 1.013711996672085 0
-0.1999165963785483 1.623132796884349 0
-0.2035784395509943 0.3768550297896387 0
-0.4547027214775108 0.506582648977719 0
-0.4501452254496273 1.507808568681515 0
-0.36780027299449 1.788782533080165 0
-0.825 1.745454371537976 0
-0.7098579157700208 1.533631566268133 0
-0.3691701628612115 0.2136139743087052 0
-0.7371075035648108 0.5348795955251022 0
-0.1899667471104287 0.8786838711583045 0
-0.1840244701869089 1.12444789757689 0
-0.1817398161387642 1.818260183861236 0
-0.181739816138764 0.181739816138764 0
-0.3712741072481359 0.7759035082776604 0
-0.3767028772038342 1.249613323632637 0
$EndNodes
$Elements
5 121 1 121
1 20 1 1
1 5 6 
1 21 1 1
2 9 8 
1 25 1 8
3 2 11 
4 11 12 
5 12 13 
6 13 14 
7 14 15 
8 15 16 
9 16 17 
10 17 1 
2 1 2 50
11 33 42 27 
12 42 43 27 
13 42 44 43 
14 43 44 34 
15 25 35 24 
16 8 9 26 
17 12 32 11 
18 37 39 35 
19 38 43 34 
20 24 45 1 
21 1 45 17 
22 11 46 2 
23 2 46 30 
24 17 39 16 
25 34 40 38 
26 30 41 29 
27 32 44 42 
28 13 44 12 
29 31 39 37 
30 16 31 15 
31 29 33 28 
32 32 42 41 
33 15 40 14 
34 14 34 13 
35 12 44 32 
36 37 38 31 
37 41 42 33 
38 16 39 31 
39 29 41 33 
40 38 40 31 
41 35 45 24 
42 34 44 13 
43 32 46 11 
44 26 36 25 
45 31 40 15 
46 36 37 35 
47 14 40 34 
48 17 45 39 
49 30 46 41 
50 25 36 35 
51 9 36 26 
52 39 45 35 
53 41 46 32 
54 7 8 26 
55 27 38 9 
56 9 37 36 
57 9 38 37 
58 10 33 27 
59 28 33 10 
60 27 43 38 
2 2 2 61
61 5 60 59 
62 18 48 3 
63 57 58 50 
64 48 62 3 
65 50 60 57 
66 52 53 51 
67 51 60 5 
68 48 61 56 
69 54 58 57 
70 11 54 12 
71 12 47 13 
72 23 50 22 
73 12 54 47 
74 56 61 55 
75 16 55 17 
76 6 52 51 
77 15 49 16 
78 22 58 21 
79 50 58 22 
80 2 65 11 
81 21 65 2 
82 17 66 1 
83 1 66 20 
84 20 61 19 
85 23 59 50 
86 49 55 16 
87 57 60 51 
88 19 48 18 
89 59 60 50 
90 54 57 47 
91 49 56 55 
92 53 68 51 
93 47 64 13 
94 4 59 23 
95 6 51 5 
96 19 61 48 
97 56 62 48 
98 56 67 52 
99 15 63 49 
100 14 63 15 
101 13 64 14 
102 11 65 54 
103 52 67 53 
104 57 68 47 
105 49 67 56 
106 55 66 17 
107 52 62 56 
108 58 65 21 
109 51 68 57 
110 54 65 58 
111 14 64 63 
112 20 66 61 
113 61 66 55 
114 47 68 64 
115 63 67 49 
116 63 64 53 
117 53 67 63 
118 64 68 53 
119 3 62 6 
120 6 62 52 
121 5 59 4 
$EndElements
4

I would start with something way simpler than this. Consider a Poisson Equation, where the solution should have a jump at an internal boundary. This problem can be posed on a built in mesh, and the variational forms should be quite straightforward.

I simply do not have time to parse longer pieces of code. Hopefully someone else comes along and has the time for it:)

1 Like
5

I see… I’ll see if I can focus on a simpler example then.
I suspect my mistake resides in the expression for “discontinuous_term”, i.e. in the extra term in the weak form due to the discontinuity in the solution.

6

@dokken Would you happen to have an idea on whether I should use “DG” elements rather than “Lagrange” when I define the spaces for the solution (volt)?

In other words, changing the line in

V = FunctionSpace(mesh, ("Lagrange", 1))
volt = Function(V)
u = TestFunction(V)

to

V = FunctionSpace(mesh, ("DG", 1))
volt = Function(V)
u = TestFunction(V)

I tried, but if I do so, Newton’s method fails to converge…

7

If you switch to DG, you need to use a DG variational formulation, see
https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/undocumented/dg-poisson/demo_dg-poisson.py
or a slightly outdated version for DOLFINx that describes some of the math

2 Likes
8

I see… then I guess this explains why I didn’t notice any impact on the solution when I modified the weak form, when I wasn’t using DG elements. And I think I know why Newton’s method fails to converge: my weak form is very likely wrong. I will study your examples and come back to you if I am stuck at deriving the weak form with the terms involving the jumps.

9

Hello dokken, after skimming through the references, I have a few questions/doubts of understanding regarding the derivation of the weak form.

Do I need to use Nitsche terms? It seems like so, but then I found examples on the Internet where they didn’t seem to use those terms: https://people.sc.fsu.edu/~jburkardt/fenics_src/dpg_laplace/dpg_laplace.py

Those Nitsche terms only apply to the external boundaries, right? They do not apply on the internal boundary where there is a discontinuity in the solution?

What are exactly, alpha and gamma? They seem to be related to the jump of the discontinuity in the solution, but I do not understand why there would be more than 1 constant involved. Why 2?

The equation I am trying to solve (Ohm’s law) is actually very simple, even compared to the Laplace equation. From \vec J = \sigma \vec E (Ohm’s law), taking the divergence and integrating by parts, one reaches the eq. -\hat n \cdot \vec J v |_{\partial \Omega} - \int_{\Omega} \sigma \nabla V \cdot \nabla v dx= 0 where the unknown is V (that’s the voltage), the test function is v, \sigma is the conductivity which is either a scalar or a tensor (essentially similar to \kappa the thermal conductivity in the examples in Fenicsx tutorial).

The regular continuous element’s weak form is just that equation. However, I am still unsure about how to derive it correctly when there are DG elements involved… I am still going through your examples, in particular the one of your website.

10

Nitsche terms are used to enforce Dirichlet boundary conditions weakly (it can be used for CG or DG methods). The integrals are only over the exterior facets (ds).

These are parameters to enforce coersivity. I would suggest reading some introductory literature to DG or papers such as
https://onlinelibrary.wiley.com/doi/epdf/10.1002/nme.2863
http://www.math.chalmers.se/Math/Research/Preprints/2010/25.pdf
or chapter 30.5.3 of the fenics book