Fenicsx and anisotropic material properties

General
1

Dear Fenicsx community,

Looking at the tutorial dolfinx-tutorial/chapter3/subdomains.ipynb at v0.4.0 · jorgensd/dolfinx-tutorial · GitHub

In particular the lines

kappa.x.array[cells_0] = np.full(len(cells_0), 1)
kappa.x.array[cells_1] = np.full(len(cells_1), 0.1)

I was wondering how to define anisotropic material properties there?
I found solutions for legacy Fenics, but I struggle to find the correct syntax with Fenicsx,

2

You have for instance: Anisotropic heterogeneous diffusion - #2 by dokken
If you have a specific set of properties in mind, it would be helpful if you could share a minimal reproducible example that serves as a base-line for what you want to achieve

3

I tried to modify the example from the above mentionned tutorial with the information in the link you provided:

import gmsh
import numpy as np

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
from dolfinx.mesh import create_unit_cube, locate_entities
from dolfinx.plot import create_vtk_mesh

from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
                 dx, grad, inner)

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
Q = FunctionSpace(mesh, ("DG", 0))

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)

k_1 = Constant(mesh, np.array([[4,0,0],[0,5,0],[0,0,6]], dtype=np.float64)) # Anisotropic

kappa.x.array[cells_0] = np.full(len(cells_0), 1.) # Isotropic
kappa.x.array[cells_1] = np.full(len(cells_1), k_1)

V = FunctionSpace(mesh, ("CG", 1))
u, v = TrialFunction(V), TestFunction(V)
a = inner(kappa*grad(u), grad(v)) * dx
x = SpatialCoordinate(mesh)
L = Constant(mesh, ScalarType(1)) * v * dx
dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
bcs = [dirichletbc(ScalarType(1), dofs, V)]

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

but I get the error

    kappa.x.array[cells_1] = np.full(len(cells_1), k_1)
    ~~~~~~~~~~~~~^^^^^^^^^
ValueError: setting an array element with a sequence.
4

Consider the following:


import numpy as np

from dolfinx.fem import (Function, FunctionSpace)
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_cube, locate_entities

from ufl import (TensorElement)

from mpi4py import MPI
from petsc4py import PETSc

mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)

el = TensorElement("DG", mesh.ufl_cell(), 0)
Q = FunctionSpace(mesh, el)
bs = Q.dofmap.bs


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)


def isotropic(x):
    values = np.ones((bs, x.shape[1]), dtype=np.float64)
    return values


def anisotropic(x):
    values = np.zeros((bs, x.shape[1]), dtype=np.float64)
    values[0] = 4
    values[4] = 5
    values[8] = 6
    return values


kappa.interpolate(isotropic, cells=cells_0)
kappa.interpolate(anisotropic, cells=cells_1)
kappa.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                         mode=PETSc.ScatterMode.FORWARD)


with XDMFFile(mesh.comm, "kappa.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(kappa)

Note the usage of a tensor vector space (as your parameter is a tensor), and the usage of cell-wise interpolation.

5

Yes this is working as intended! Thank you Dokken.
Interesting to see the new logic behind the syntax.
@Dokken If I may ask a last question, why

    values[0] = 4
    values[4] = 5
    values[8] = 6

and not values[0], values[1], values[2]?

6

It’s a tensor element, so you’re assigning values to the diagonal of the row-major tensor array. You could also do the following, for example:

def anisotropic(x):
    values = np.zeros((3, 3, x.shape[1]), dtype=np.float64)
    values[0, 0] = 4
    values[1, 1] = 5
    values[2, 2] = 6
    return values.reshape((bs, x.shape[1]))
2 Likes