A 2D problem, faulty weak form

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!

If your mesh is 2D, you should use prune_z=True. If you do not prune the z-coordinate, your mesh is interpreted as a 2D object embedded in 3D space. Since in this case there are 3 spatial variables, the grad operator will add a dimension of length 3 to the shape of the operand, e.g. for your problem volt is a scalar (thus its shape is the empty tuple, ()), so grad(volt) will have shape ()+(3, ) or (3,).

Indeed, using

print(grad(volt).ufl_shape)
print(sigma_tensor.ufl_shape)

I find

(3,)
(2, 2)

which clearly do not match (you cannot multiply a 2x2 matrix with a 3-vector). If you use

    triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
    line_mesh = create_mesh(msh, "line", prune_z=True)

in the mesh conversion part of your script, it will be treated as a 2D object in 2D space, and the shape of grad(volt) will be (2,). Its shape will then match with sigma_tensor for the tensor-vector multiplication, and also with grad(u) for the subsequent inner product.

2 Likes

Thanks a lot. I actually changed more than these 2 lines, there was another part where I forgot to set prune_z to True. So far so good, however this lead to an Arity mismatch error for the line with problem = NonlinearProblem(F, volt, bcs=bcs).

I didn’t know how to fix the problem, but when I changed the line volt = TrialFunction(V) to volt = Function(V), there was no more error when I run the code above. But I do not know exactly if what I have done is correct, or whether I am doing things I don’t even understand and don’t correspond to what I want to do.

Then, when I continue my code with

from dolfinx import log
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(volt)
assert(converged)
print(f"Number of interations: {n:d}")

I get the error

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In [49], line 3
      1 from dolfinx import log
      2 log.set_log_level(log.LogLevel.WARNING)
----> 3 n, converged = solver.solve(volt)
      4 assert(converged)
      5 print(f"Number of interations: {n:d}")

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/nls/petsc.py:41, in NewtonSolver.solve(self, u)
     38 def solve(self, u: fem.Function):
     39     """Solve non-linear problem into function u. Returns the number
     40     of iterations and if the solver converged."""
---> 41     n, converged = super().solve(u.vector)
     42     u.x.scatter_forward()
     43     return n, converged

RuntimeError: Newton solver did not converge because maximum number of iterations reached

but while I will try to debug the solver, I would like to know whether what I have done so far is OK?

Edit: I have increased the max iterations number to 10 k, to no avail.

Without your amended code, it’s hard for me to be 100% sure what you have done, and therefore hard for me to be 100% sure that it is correct. But it sounds like it is correct.

The “arity” of a form is the number of Arguments that it has:
https://fenics.readthedocs.io/projects/ufl/en/latest/manual/form_language.html#form-arguments

TrialFunction and TestFunction are just shortcuts to create Arguments in UFL. If a form contains both a TrialFunction and a TestFunction, it has an arity of 2. (If it is linear in both the TrialFunction and the TestFunction, it is a bilinear form.) A linear problem expects a bilinear form for the left-hand side and a linear form for the right-hand side, i.e. forms with arities of 2 and 1, respectively. A nonlinear problem expects a form of arity 1, so it should contain a combination of Function objects and a TestFunction. Hence, what you have done by changing TrialFunction to Function is correct. (The form you have provided is a linear problem, so I’m not sure why you are using a NonlinearProblem.)

2 Likes

Thanks a lot for all this information, I wasn’t sure I could still trust fenics documentation for fenicsx, but I guess that what was written for UFL is still entirely valid for fenicsx.

OK.

I am using a nonlinear solver because I later plan to complicate the problem by using, say, a non linear function sigma(T), i.e. a non linear electrical resistance dependence with temperature. I would rather not have to change my code much, I would just deal with the sigma part instead of redefining the problem, solver, Test/Trial functions AND the sigma part.

My current code is:

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=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("rectangles.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("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 = Function(V) # 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 = dot(grad(u), sigma_tensor * grad(volt))*dx + u * J * ds(vanish_voltage_curve)

# Solve the PDE.
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
#solve(F == 0, volt, bcs, solver_parameters={"newton_solver":{"relative_tolerance": 1e-11, "absolute_tolerance": 2e-9, "maximum_iterations":1000}})
problem = NonlinearProblem(F, volt, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-4
solver.report = True

from petsc4py import PETSc
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
opts[f"{option_prefix}ksp_max_it"]= 10000
ksp.setFromOptions()

from dolfinx import log
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(volt)
assert(converged)
print(f"Number of interations: {n:d}")

The error returned is

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In [74], line 3
      1 from dolfinx import log
      2 log.set_log_level(log.LogLevel.WARNING)
----> 3 n, converged = solver.solve(volt)
      4 assert(converged)
      5 print(f"Number of interations: {n:d}")

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/nls/petsc.py:41, in NewtonSolver.solve(self, u)
     38 def solve(self, u: fem.Function):
     39     """Solve non-linear problem into function u. Returns the number
     40     of iterations and if the solver converged."""
---> 41     n, converged = super().solve(u.vector)
     42     u.x.scatter_forward()
     43     return n, converged

RuntimeError: Newton solver did not converge because maximum number of iterations reached

I probably wrongly implemented the problem I wish to solve, as there is no reason the solver shouldn’t converge.

You may want to check your mesh. It looks like the mesh vertices form a 9x9 grid, so your CG-1 function space should have 81 dofs. However, print(volt.vector.size) gives me 90, i.e. the mesh is probably disconnected between subdomains 11 and 12, adding an extra 9 dofs to the problem.

This would explain why the problem doesn’t converge; as you are prescribing a current into the left subdomain, but there is no way for the current to exit.

1 Like

I am extremely thankful to you, and the time you’ve spent debugging my code. Wow, I wouldn’t have thought to check the mesh, especially in such details. Visually, with Gmsh, I do not see anything “wrong”, nor with Paraview, but I suppose you’re right, this would make sense.

In my Gmsh code however, I think/thought I would avoid counting twice the line or curve between the 2 regions (subdomains 11 and 12). The code I used is:

SetFactory("OpenCASCADE");
Rectangle(1) = {0, 0, 0, 1, 2, 0};
Rectangle(2) = {-1, 0, 0, 1, 2, 0};
BooleanFragments{ Curve{4}; Delete; }{ }
Line(9) = {8, 5};
Line(10) = {3, 2};
BooleanFragments{ Curve{8}; Delete; }{ }
BooleanFragments{ Curve{2}; Delete; }{ }
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};
Transfinite Surface {2};
Transfinite Surface {1};

I thought the part BooleanFragments{ Curve{4}; Delete; }{ } would handle this. But if it doesn’t, then I don’t really know how to fix the mesh…

Edit: After having read gmsh’s documentation and an example, it looks like I should have specified two curves in the BooleanFragments’ argument.
My new try at a mesh “fails” in that half of the domain is non meshed… but I do not understand why. Here’s my code:

//+
SetFactory("OpenCASCADE");
//+
Rectangle(1) = {0, 0, 0, 1, 2, 0};
//+
Rectangle(2) = {-1, 0, 0, 1, 2, 0};
//+

//+
Line(9) = {8, 5};
//+
Line(10) = {3, 2};

//+

//+
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};
BooleanFragments{ Curve{4,6}; Delete; }{}
BooleanFragments{ Curve{8,9}; Delete; }{ }
//+
BooleanFragments{ Curve{2,10}; Delete; }{ }
//+
Transfinite Surface {2};
//+
Transfinite Surface {1};

@raboynics: I actually discovered the mesh issue by accident. I tried simplifying the problem with Dirichlet BCs on the left and right boundaries and no Neumann BC, i.e.:

u_bc = Function(V)
left_facets = ft.find(vanish_voltage_curve)
left_dofs = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
print(left_dofs)
bcs = [dirichletbc(ScalarType(V_current_contact_out), left_dofs, V)]
right_facets = ft.find(inj_current_curve)
right_dofs = locate_dofs_topological(V, mesh.topology.dim-1, right_facets)
print(right_dofs)
bcs = [
    dirichletbc(ScalarType(V_current_contact_out), left_dofs, V),
    dirichletbc(ScalarType(1.0), right_dofs, V),
]

which gave me this solution that suggested the disconnected mesh:

How about the following .geo file?

SetFactory("OpenCASCADE");
Rectangle(1) = {0, 0, 0, 1, 2, 0};
Rectangle(2) = {-1, 0, 0, 1, 2, 0};
BooleanFragments{ Surface{1}; Delete; }{ Surface{2}; Delete; }
Physical Surface("material A", 11) = {2};
Physical Surface("material B", 12) = {1};
Physical Curve("left_mat_A", 13) = {7};
Physical Curve("right_mat_B", 14) = {2};
Transfinite Surface {2};
Transfinite Surface {1};
1 Like

Thank you again! Ah, I see… yes, good idea.
I think I want to define 2 more curves than in your example (as in mine). The reason is that in the end, the “contacts” will be subsets of the left and right curves, i.e. the current will enter and leave in some part of the right and left ends of materials A and B, rather than the full curves. I will try to modify the file.geo on my own… and I will come back here if I am unable to do it. But your example should help me with it.

Edit: I could make it, in the end. I’m very happy, thank you so much! I will complicate this example though, I’m sure I’ll face many problems in a near future :slight_smile: