# Error when implementing Subdomain Tutorial on CoLab

Hi,

When I tried to implement the Subdomain on built-in meshes Tutorial demo on CoLab following Defining subdomains for different materials — FEniCSx tutorial, the following error occurs becasue of the failure of kappa * grad(u). kappa is a vector and grad(u) is a matrix here. Any idea here how to solve this problem? Thanks!

from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import create_unit_square, locate_entities
from dolfinx.plot import vtk_mesh

from ufl import (SpatialCoordinate, TestFunction, TrialFunction,

from mpi4py import MPI

import meshio
import gmsh
import numpy as np
import pyvista

pyvista.start_xvfb()
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
Q = dolfinx.fem.functionspace(mesh, ("DG", 0, (mesh.geometry.dim,)))

def Omega_0(x):
return x[1] <= 0.5

def Omega_1(x):
return x[1] >= 0.5

kappa = Function(Q)
cells_0 = locate_entities(mesh, mesh.topology.dim, Omega_0)
cells_1 = locate_entities(mesh, mesh.topology.dim, Omega_1)
kappa.x.array[cells_0] = np.full_like(cells_0, 1, dtype=default_scalar_type)
kappa.x.array[cells_1] = np.full_like(cells_1, 0.1, dtype=default_scalar_type)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))
u, v = TrialFunction(V), TestFunction(V)
# x = SpatialCoordinate(mesh)
# L = Constant(mesh, default_scalar_type(1)) * v * dx
# dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
# bcs = [dirichletbc(default_scalar_type(1), dofs, V)]


The error information is

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-54-b35bd3cd4c0d> in <cell line: 40>()
38 V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))
39 u, v = TrialFunction(V), TestFunction(V)
41 # x = SpatialCoordinate(mesh)
42 # L = Constant(mesh, default_scalar_type(1)) * v * dx

1 frames
/usr/local/lib/python3.10/dist-packages/ufl/exproperators.py in _mult(a, b)
157
158     else:
--> 159         raise ValueError(f"Invalid ranks {r1} and {r2} in product.")
160
161     # TODO: I think applying as_tensor after index sums results in

ValueError: Invalid ranks 1 and 2 in product.


Here you define Q as a vector space, then

and

Does not make sende as the first product is a vector, while the second is matrix

I am trying to learn how to set up subdomains with different kappa, so I followed the Subdomain Demo Defining subdomains for different materials — FEniCSx tutorial, and want to set up different kappa for different domains. Every code is the same with the Demo, except for

Q = dolfinx.fem.functionspace(mesh, ("DG", 0, (mesh.geometry.dim,)))


and

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))


In the Demo, they are

Q = FunctionSpace(mesh, ("DG", 0))


and

V = FunctionSpace(mesh, ("Lagrange", 1))


But another error would occur if I used the code shown in the Demo. I thought it was becasue of the update of FuncationSpace. Anything else is the same. Then I was not sure how to set up different kappa for different domains.

I saw another way https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html in FeniCs instead of FeniCSx is through Expression.

class K(Expression):
def __init__(self, materials, k_0, k_1, **kwargs):
self.materials = materials
self.k_0 = k_0
self.k_1 = k_1

def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1

kappa = K(materials, k_0, k_1, degree=0)


I can try this to see if this will work. But it would be great if you may have other suggestions.

Thank you!

Should kappa still be a constant per cell (But unique per cell) Or what does kappa signify in your code?

Thank you for your quick reply! kappa is a constant per cell, i.e. for one domain kappa is value_1 and for another domain kappa is value_2.

Let’s say the conductivity is different within different domains. So for one domain, the governing equation is
\partial T / \partial t = kappa_1 * diffusion term,
For another subdomain, the governing equation is
\partial T / \partial t = kappa_2 * diffusion term.

And FYI, the interfaces between two subdomains will have a Neumann boundary condition. But this will be implemented later, it is not related to the question here.

Thanks!

Then just use kappa as done in the tutorial.

That was what I did before.

If I copied and pasted the tutorial code, then I will get the error shown below,

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-7-1195952df25b> in <cell line: 22>()
20
21 mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
---> 22 Q = FunctionSpace(mesh, ("DG", 0))
23
24 def Omega_0(x):

TypeError: FunctionSpace.__init__() missing 1 required positional argument: 'cppV'


If I copied and pasted the tutorial code but used “Q = dolfinx.fem.functionspace(mesh, (“DG”, 0, (mesh.geometry.dim,)))” instead, then I will get the error shown in my first message,

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-8-6c2490af87d4> in <cell line: 46>()
44 V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))
45 u, v = TrialFunction(V), TestFunction(V)
48 # x = SpatialCoordinate(mesh)

1 frames
/usr/local/lib/python3.10/dist-packages/ufl/exproperators.py in _mult(a, b)
157
158     else:
--> 159         raise ValueError(f"Invalid ranks {r1} and {r2} in product.")
160
161     # TODO: I think applying as_tensor after index sums results in

ValueError: Invalid ranks 1 and 2 in product.


I guess my question is I cannot run the tutorial demo, but I thought it should work because I copied and pasted the tutorial demo. Is it possible it is because I am using FeniCSx on Colab?

I copied and pasted the demo in this one Defining subdomains for different materials — FEniCSx tutorial.

You are using the main branch of DOLFINx, that has some API changes. The tutorial reflects the latest stable release (0.7.x).

You should use

import dolfinx.fem
Q = fem.functionspace(mesh, ("DG", 0))


Most changes relating to API changes are listed at

I see. Thank you so much for this!

It works now on Colab! Thank you!