I am attempting to solve Ohm’s law in a rectangular shape, with 2 different materials next to each other. I have created the file.msh using gmsh, here is the file:
$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
4
1 13 "left_mat_A"
1 14 "right_mat_B"
2 11 "material A"
2 12 "material B"
$EndPhysicalNames
$Entities
8 10 2 0
1 0 0 0 0
2 1 0 0 0
3 1 2 0 0
4 0 2 0 0
5 -1 0 0 0
6 0 0 0 0
7 0 2 0 0
8 -1 2 0 0
1 -9.999999994736442e-08 -1e-07 -1e-07 1.0000001 1e-07 1e-07 0 2 1 -2
2 0.9999999000000001 -9.999999994736442e-08 -1e-07 1.0000001 2.0000001 1e-07 1 14 2 2 -3
3 -9.999999994736442e-08 1.9999999 -1e-07 1.0000001 2.0000001 1e-07 0 2 3 -4
4 -1e-07 -9.999999994736442e-08 -1e-07 1e-07 2.0000001 1e-07 0 2 4 -1
5 -1.0000001 -1e-07 -1e-07 9.999999994736442e-08 1e-07 1e-07 0 2 5 -6
6 -1e-07 -9.999999994736442e-08 -1e-07 1e-07 2.0000001 1e-07 0 2 6 -7
7 -1.0000001 1.9999999 -1e-07 9.999999994736442e-08 2.0000001 1e-07 0 2 7 -8
8 -1.0000001 -9.999999994736442e-08 -1e-07 -0.9999999000000001 2.0000001 1e-07 1 13 2 8 -5
9 -1.0000001 -9.999999994736442e-08 -1e-07 -0.9999999000000001 2.0000001 1e-07 0 2 8 -5
10 0.9999999000000001 -9.999999994736442e-08 -1e-07 1.0000001 2.0000001 1e-07 0 2 3 -2
1 -9.999999994736442e-08 -9.999999994736442e-08 -1e-07 1.0000001 2.0000001 1e-07 1 12 4 1 2 3 4
2 -1.0000001 -9.999999994736442e-08 -1e-07 9.999999994736442e-08 2.0000001 1e-07 1 11 4 5 6 7 8
$EndEntities
$Nodes
18 90 1 90
0 1 0 1
1
0 0 0
0 2 0 1
2
1 0 0
0 3 0 1
3
1 2 0
0 4 0 1
4
0 2 0
0 5 0 1
5
-1 0 0
0 6 0 1
6
0 0 0
0 7 0 1
7
0 2 0
0 8 0 1
8
-1 2 0
1 1 0 3
9
10
11
0.2500000000000001 0 0
0.5000000000000001 0 0
0.75 0 0
1 2 0 7
12
13
14
15
16
17
18
1 0.2500000000000002 0
1 0.5000000000000001 0
1 0.7500000000000001 0
1 1 0
1 1.25 0
1 1.5 0
1 1.75 0
1 3 0 3
19
20
21
0.75 2 0
0.4999999999999999 2 0
0.25 2 0
1 4 0 7
22
23
24
25
26
27
28
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 5 0 3
29
30
31
-0.75 0 0
-0.4999999999999999 0 0
-0.25 0 0
1 6 0 7
32
33
34
35
36
37
38
0 0.2500000000000002 0
0 0.5000000000000001 0
0 0.7500000000000001 0
0 1 0
0 1.25 0
0 1.5 0
0 1.75 0
1 7 0 3
39
40
41
-0.2500000000000001 2 0
-0.5000000000000001 2 0
-0.75 2 0
1 8 0 7
42
43
44
45
46
47
48
-1 1.75 0
-1 1.5 0
-1 1.25 0
-1 0.9999999999999998 0
-1 0.75 0
-1 0.5 0
-1 0.25 0
2 1 0 21
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
0.2500000000000001 0.25 0
0.2500000000000001 0.5 0
0.25 0.75 0
0.25 0.9999999999999999 0
0.25 1.25 0
0.25 1.5 0
0.25 1.75 0
0.5000000000000001 0.25 0
0.5 0.5 0
0.5 0.75 0
0.5 1 0
0.5 1.25 0
0.4999999999999999 1.5 0
0.4999999999999999 1.75 0
0.75 0.2500000000000002 0
0.75 0.5 0
0.75 0.7500000000000001 0
0.75 1 0
0.75 1.25 0
0.75 1.5 0
0.75 1.75 0
2 2 0 21
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
-0.75 0.25 0
-0.75 0.5 0
-0.75 0.75 0
-0.75 0.9999999999999999 0
-0.75 1.25 0
-0.75 1.5 0
-0.75 1.75 0
-0.4999999999999999 0.25 0
-0.4999999999999999 0.5 0
-0.5 0.75 0
-0.5 1 0
-0.5 1.25 0
-0.5 1.5 0
-0.5000000000000001 1.75 0
-0.25 0.2500000000000002 0
-0.25 0.5 0
-0.25 0.7500000000000001 0
-0.25 1 0
-0.25 1.25 0
-0.2500000000000001 1.5 0
-0.2500000000000001 1.75 0
$EndNodes
$Elements
4 144 1 144
1 2 1 8
1 2 12
2 12 13
3 13 14
4 14 15
5 15 16
6 16 17
7 17 18
8 18 3
1 8 1 8
9 8 42
10 42 43
11 43 44
12 44 45
13 45 46
14 46 47
15 47 48
16 48 5
2 1 2 64
17 1 9 28
18 28 9 49
19 28 49 27
20 27 49 50
21 27 50 26
22 26 50 51
23 26 51 25
24 25 51 52
25 25 52 24
26 24 52 53
27 24 53 23
28 23 53 54
29 23 54 22
30 22 54 55
31 22 55 4
32 4 55 21
33 9 10 49
34 49 10 56
35 49 56 50
36 50 56 57
37 50 57 51
38 51 57 58
39 51 58 52
40 52 58 59
41 52 59 53
42 53 59 60
43 53 60 54
44 54 60 61
45 54 61 55
46 55 61 62
47 55 62 21
48 21 62 20
49 10 11 56
50 56 11 63
51 56 63 57
52 57 63 64
53 57 64 58
54 58 64 65
55 58 65 59
56 59 65 66
57 59 66 60
58 60 66 67
59 60 67 61
60 61 67 68
61 61 68 62
62 62 68 69
63 62 69 20
64 20 69 19
65 11 2 63
66 63 2 12
67 63 12 64
68 64 12 13
69 64 13 65
70 65 13 14
71 65 14 66
72 66 14 15
73 66 15 67
74 67 15 16
75 67 16 68
76 68 16 17
77 68 17 69
78 69 17 18
79 69 18 19
80 19 18 3
2 2 2 64
81 5 29 48
82 48 29 70
83 48 70 47
84 47 70 71
85 47 71 46
86 46 71 72
87 46 72 45
88 45 72 73
89 45 73 44
90 44 73 74
91 44 74 43
92 43 74 75
93 43 75 42
94 42 75 76
95 42 76 8
96 8 76 41
97 29 30 70
98 70 30 77
99 70 77 71
100 71 77 78
101 71 78 72
102 72 78 79
103 72 79 73
104 73 79 80
105 73 80 74
106 74 80 81
107 74 81 75
108 75 81 82
109 75 82 76
110 76 82 83
111 76 83 41
112 41 83 40
113 30 31 77
114 77 31 84
115 77 84 78
116 78 84 85
117 78 85 79
118 79 85 86
119 79 86 80
120 80 86 87
121 80 87 81
122 81 87 88
123 81 88 82
124 82 88 89
125 82 89 83
126 83 89 90
127 83 90 40
128 40 90 39
129 31 6 84
130 84 6 32
131 84 32 85
132 85 32 33
133 85 33 86
134 86 33 34
135 86 34 87
136 87 34 35
137 87 35 88
138 88 35 36
139 88 36 89
140 89 36 37
141 89 37 90
142 90 37 38
143 90 38 39
144 39 38 7
$EndElements
Here is my Python code:
from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh("rectangles.msh", MPI.COMM_WORLD, gdim=2)
# Now I want to convert the mesh file to XDMF using meshio.
proc = MPI.COMM_WORLD.rank
import meshio
def create_mesh(mesh, cell_type, prune_z=False):
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("rectangles.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "triangle", prune_z=False)
line_mesh = create_mesh(msh, "line", prune_z=False)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mt.xdmf", line_mesh)
# Now read the mesh files produced above.
from dolfinx.io import XDMFFile
with XDMFFile(MPI.COMM_WORLD, "mesh.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, "mt.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
# Define function space
V = FunctionSpace(mesh, ("CG", 1))
volt = TrialFunction(V)
u = TestFunction(V)
rho_tensor = as_tensor([[1, 0], [0, 1]])
sigma_tensor = inv(rho_tensor)
# 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};
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
# We define the boundary conditions
from petsc4py.PETSc import ScalarType
u_bc = Function(V)
left_facets = ft.find(inj_current_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)
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))))
from ufl import dot
# Weak form.
F = inner(grad(u), sigma_tensor * grad(volt))*dx + u * J * ds(vanish_voltage_curve)
The very last line produces the error:
ERROR:UFL:Shapes do not match: <Grad id=140443392115072> and <ComponentTensor id=140443392775632>.
---------------------------------------------------------------------------
UFLException Traceback (most recent call last)
Cell In [67], line 3
1 from ufl import dot
2 # Weak form.
----> 3 F = inner(grad(u), sigma_tensor * grad(volt))*dx + u * J * ds(vanish_voltage_curve)
5 # Solve the PDE.
6 from dolfinx.fem.petsc import NonlinearProblem
File /usr/local/lib/python3.10/dist-packages/ufl/operators.py:158, in inner(a, b)
156 if a.ufl_shape == () and b.ufl_shape == ():
157 return a * Conj(b)
--> 158 return Inner(a, b)
File /usr/local/lib/python3.10/dist-packages/ufl/tensoralgebra.py:147, in Inner.__new__(cls, a, b)
145 ash, bsh = a.ufl_shape, b.ufl_shape
146 if ash != bsh:
--> 147 error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
149 # Simplification
150 if isinstance(a, Zero) or isinstance(b, Zero):
File /usr/local/lib/python3.10/dist-packages/ufl/log.py:158, in Logger.error(self, *message)
156 "Write error message and raise an exception."
157 self._log.error(*message)
--> 158 raise self._exception_type(self._format_raw(*message))
UFLException: Shapes do not match: <Grad id=140443392115072> and <ComponentTensor id=140443392775632>.
I tried a few things, such as changing TrialFunction to Function, to no avail. Or changing the sigma_tensor for a 3D version of it, to no avail. I also tried to use dot() instead of inner(), again, with a similar error. I suspect the problem has to do with the part “sigma_tensor * grad(volt)”, but I am not 100% sure. Any idea on how to solve the problem is very appreciated!