Different material properties in subdomains, tensor instead of scalar

General
1

I am attempting to solve Ohm’s law in a 2D region with 2 subregions, next to each other. Each subregion represent a different material. I defined the resistivity tensors as

rho_tensor_A = as_tensor([[1, 0], [0, 1]])
sigma_tensor_A = inv(rho_tensor_A)

rho_tensor_B = as_tensor([[10, 0], [0, 10]])
sigma_tensor_B = inv(rho_tensor_B)

I am following the tutorial Defining subdomains for different materials — FEniCSx tutorial, but I get stuck at the part kappa.x.array[bottom_cells] = np.full_like(bottom_cells, 1, dtype=ScalarType)

In my case, I tried

rho_tensor.x.array[11] = np.full_like(11, rho_tensor_A, dtype=ScalarType)

where “11” represents the Physical Surface defined with Gmsh GUI for the material A. However, this returns

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
TypeError: ListTensor.__float__ returned non-float (type NotImplementedType)

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
Cell In [11], line 19
     16 sigma_tensor_B = inv(rho_tensor_B)
     18 from petsc4py.PETSc import ScalarType
---> 19 rho_tensor.x.array[11] = np.full_like(11, rho_tensor_A, dtype=ScalarType)
     21 Q = FunctionSpace(mesh, ("DG", 0))
     22 kappa = Function(Q)

File <__array_function__ internals>:180, in full_like(*args, **kwargs)

File /usr/local/lib/python3.10/dist-packages/numpy/core/numeric.py:417, in full_like(a, fill_value, dtype, order, subok, shape)
    359 """
    360 Return a full array with the same shape and type as a given array.
    361 
   (...)
    414 
    415 """
    416 res = empty_like(a, dtype=dtype, order=order, subok=subok, shape=shape)
--> 417 multiarray.copyto(res, fill_value, casting='unsafe')
    418 return res

File <__array_function__ internals>:180, in copyto(*args, **kwargs)

ValueError: setting an array element with a sequence.

I suppose the error arises because I am not giving a scalar, but a list of lists (or a matrix, say). I am not sure how to fix the problem, I’ve been trying to replace the ScalarType by other things, without success.

Edit: I am still trying. Here is my latest code:

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)

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

# Define function space
V = FunctionSpace(mesh, ("CG", 1))
volt = Function(V) # TrialFunction(V)
u = TestFunction(V)

rho_tensor_A = as_tensor([[1, 0], [0, 1]])
#rho_tensor_A = 1
sigma_tensor_A = inv(rho_tensor_A)

rho_tensor_B = as_tensor([[10, 0], [0, 10]])
#rho_tensor_B = 1
sigma_tensor_B = inv(rho_tensor_B)

from petsc4py.PETSc import ScalarType
rho_tensor = Function(V)
rho_tensor.x.array[11] = np.full_like(11, rho_tensor_A, dtype=ScalarType)

which yields the error above.

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

Edit 2 : Trying to debug, I think the error shows up because of some incompatibility with ufl’s “as_tensor” and numpy arrays. The part np.full_like() expects an array whereas I provide a ufl.as_tensor() instead. I don’t really know how to overcome this problem, nor what is the suggested way people do it.

2

How do you expect a tensor to be assigned to a scalar function space?
The space ("CG", 1) contains one degree of freedom at any point. It does not make sense to try to add a scalar to it.
If you want to do so, you need to use a TensorElement or a TensorFunctionSpace as shown for a dummy rho below:

import dolfinx
from mpi4py import MPI
import ufl
import numpy as np
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

el = ufl.TensorElement("Lagrange", mesh.ufl_cell(), 1, shape=(2, 2))

V = dolfinx.fem.FunctionSpace(mesh, el)
u = dolfinx.fem.Function(V)


def rho_B(x):
    tensor = np.array([[1, 2], [3, 4.]])

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


u.interpolate(rho_B)

for i in range(2):
    for j in range(2):
        print(i, j, dolfinx.fem.assemble_scalar(
            dolfinx.fem.form(u[i, j]*ufl.dx)))

if you want the rho to be different in different subdomains, you should use the cells kwarg of interpolate, as it will then only do interpolation on dofs associated to those entities. In your case that would be
u.interpolate(rho_B, cells=ct.find(11))

1 Like
3

I’m not quite sure, I didn’t know about “degree of freedom” in FEM.

Hmm wait, if it doesn’t make sense, I shouldn’t want to do so, right? I would rather do things in a correct, usual way.

I am not sure I understand what you say. But yes, I want rho to be like the kappa of the tutorial, i.e. I want it to be region dependent (except that it is a tensor rather than a scalar). Here is my attempt:

from ufl import TensorElement
from dolfinx.fem import TensorFunctionSpace
# Define a tensor element (also possible to define a tensor space)
el = TensorElement("Lagrange", mesh.ufl_cell(), 1, shape=(2, 2))
V = FunctionSpace(mesh, el)
#V = TensorFunctionSpace("Lagrange", mesh.ufl_cell(), 1, shape=(2, 2))

volt = Function(V)

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])


#volt.interpolate(rho_B)

for i in range(2):
    for j in range(2):
        print(i, j, assemble_scalar(form(volt[i, j]*dx)))
        
volt.interpolate(rho_A, cells=ct.find(11))
volt.interpolate(rho_B, cells=ct.find(12))
print(volt.x.array)
print(len(volt.x.array))

It looks like in my example I am setting “volt” to be equal to rho(x), which is not what I wish. In other words, I do not understand why I would have to interpolate “volt” by using rho_A, rho_B. Maybe I goofed and I shouldn’t replaced your dummy “u” by volt. But it wouldn’t make sense to re place it by “rho” in my case. I think I am lost more than anything, at this point…

4

Given the two different rho’s, you should create a single rho, which should not replace volt, i.e

from mpi4py import MPI
import ufl
import numpy as np

# Define mesh and ct
# ...
# ...

tensor_el = ufl.TensorElement("DG", mesh.ufl_cell(), 0, shape=(2, 2))

T = dolfinx.fem.FunctionSpace(mesh, el)
rho = dolfinx.fem.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])


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

V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
volt = dolfinx.fem.Function(V)

# Define some variational form with rho and volt
v = ufl.TestFunction(V)
a = ulf.inner(ufl.dot(rho, ufl.grad(volt), ufl.grad(v))*ufl.dx
1 Like
5

Thank you very much dokken, this is much clearer to me now. I still don’t understand everything, for example why the “1” becomes a “0” in tensor_el = ufl.TensorElement("DG", mesh.ufl_cell(), 0, shape=(2, 2)), I suppose it’s because rho is a piecewise constant valued tensor.

I had to do a little adjustment to continue with the rest of my code that works when rho is of “as_tensor” type. This is because I have “sigma_tensor”, which is the inverse of the rho_tensor, in my expression of the weak form. However, despite the fact that I do not get any error, I get a result that doesn’t make sense. I.e., when I solve for “volt”, I get a voltage that doesn’t make sense, it doesn’t even satisfy the Dirichlet boundary condition.

I’ll try to solve the problem myself, and possibly come back here in case I get stuck.

6

You can use DG-1, but it is more computationally expensive than DG-0, and doesn’t modify the result (as the input tensor is a piecewise constant).

1 Like
7

Unfortunately, I am still stuck, my code “runs” but the result is not the expected one. I’ll provide part of my full code, as I suspect this is the faulty one (the rest of my code is reused in some of my other codes, and seems to work properly).

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([[10, 0], [0, 10.]]))
    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)

# 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_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()

The solution (volt), looks like so:

2023-01-06-171002_1920x1080_scrot
2023-01-06-171002_1920x1080_scrot1001×793 43.3 KB

which is completely bogus. The Dirichlet’s boundary condition (voltage = 0 on a curve) is not even respected:

The voltages are 1.1527004340747914e+158 and 1.2049001140225463e+144 in the curves where the current enters and leaves. (one of them should be 0). The simulated measured voltage is -1.1527004340747794e+158, which is astronomical… a few volts is the expected answer.
Is there something obviously wrong with the code I have posted here?

EDIT : I think I figured it out. I defined twice V, this messed things up. I have now a much more reasonable solution. Do not spend time trying to help me anymore… I just leave this post just in case.

8

Hello dokken,

Thanks for your answer. I followed the same approach for a different problem and it works! However, the tensor depends upon the spatial coordinates in the original question. Is it possible to update this tensor value depending on the time step under a while loop? Thanks in advance.

9

See for instance how the inlet condition is done in Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial

1 Like