Error when implementing Subdomain Tutorial on CoLab

dolfinx
1

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,
                 dx, grad, inner)

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)
a = inner(kappa * grad(u), grad(v)) * dx
# 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)
---> 40 a = inner(kappa * grad(u), grad(v)) * dx
     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.
2

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

3

Thank you for replying!

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!

4

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

5

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!

6

Then just use kappa as done in the tutorial.

7

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)
---> 46 kappa * grad(u)
     47 a = inner(kappa * grad(u), grad(v)) * dx
     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.

8

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

9

I see. Thank you so much for this!

10

It works now on Colab! Thank you!