Fenicsx and anisotropic material properties

General
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.