Anisotropic heterogeneous diffusion

variational formulation
1

Hello everyone,
anybody knows how to define a tensorial (constant or non-constant) function to insert into the weak form for an anisotropic heterogeneous diffusion term, e.g.,

ufl.dot(ufl.grad(v),ufl.dot(K,ufl.grad(u)))*ufl.dx

I tried to define K as a constant numpy matrix or as a function or using ufl.TensorConstant but without success and can’t find any example related to this

I love fenicsx compared to fences but the definition of constants and functions has become non intuitive IMHO

Thanks

2

Consider the following:

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

x = ufl.SpatialCoordinate(mesh)

constant_tensor = dolfinx.fem.Constant(
    mesh, np.array([[1, 0], [2, 1]], dtype=np.float64))
variable_tensor = ufl.as_tensor([[x[0], 0], [1, 3*x[1]]])
print(constant_tensor.ufl_shape)
print(variable_tensor.ufl_shape)

print(dolfinx.fem.assemble_scalar(
    dolfinx.fem.form(
        ufl.inner(constant_tensor, variable_tensor)*ufl.dx)))
6 Likes
3

Hi Dokken,

Would you mind giving a detailed example to show how to achieve anisotropic diffusion in dolfinx?

Here is my demo:

import dolfinx.fem.petsc
from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem import FunctionSpace
from dolfinx import fem
import numpy
import ufl
from petsc4py.PETSc import ScalarType

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = FunctionSpace(domain, ("CG", 1))

uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)

# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, ScalarType(-6))

ijk = ufl.SpatialCoordinate(domain)
K = ufl.as_tensor([[ijk[0], 0], [1, 3*ijk[1]]])
a = (K  * ufl.dot(ufl.grad(u), ufl.grad(v))) * ufl.dx
L = f * v * ufl.dx

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

There is an error: ERROR:UFL:Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2, 2) and free indices with labels ().

Would you mind helping me to solve the problem?

4

Your variational form is not correct.
You start with -div(K*grad(u))*v*dx which is integrated by part to be inner(K*grad(u), grad(v))*dx

5

Sorry for my stupid mistake and thank you Dokken.

The code works well now.