How to define a MeshFunction in terms of spatial coordinates?

1

Hi everyone,

I am trying to solve the Stokes equations in 2D with a prescribed divergence and an evolving mesh. I use dolfin to solve the Stokes equations at each timestep and then use ALE to move the mesh and boundary. The divergence is a spatially varying function defined in relation to the original geometric configuration.

I would like to resolve the divergence on the level of the mesh. For this purpose, I tried to define a MeshFunction at t=0 and associate a value for the divergence (ie a scalar) with each cell. The idea is that this would remain constant throughout the simulation and automatically be advected when I move the mesh with ALE.move.

How do I do this? Below is a sketch of the relevant code section. The step that is not working is defining the mesh function in terms of the spatial coordinates. Unfortunately the spatial profile of the divergence is too complex to work with manually defined subdomains, otherwise I would have used that.

In case this approach isn’t possible I would be very grateful for workaround suggestions!

from dolfin import *
from mshr import *

# generate mesh
mesh = generate_mesh(Circle(Point(0, 0), 1), 10)

# Define divergence in terms of initial coordinates
divergence = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

# Define MeshFunction
divMF = MeshFunction('double', mesh, 2, divergence) # This doesn't work

# Define UserExpression to calculate divergence in UFL form
class DivUE(UserExpression):
    def __init__(self, divMF, **kwargs):
        super().__init__(**kwargs)
        self.divMF = divMF
    def eval_cell(self, values, x, cell):
        values[0] = self.divMF[cell.index]
    def value_shape(self):
        return () #scalar
    
gamma = DivUE(divMF, degree=0)

# Define function spaces
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
a = (inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx-(1E-10)*p*q*dx
L = gamma*q*dx

# etc

Many thanks for your help!

2

You can use a DG 0 function, i.e.

from dolfin import *
from mshr import *

# generate mesh
mesh = generate_mesh(Circle(Point(0, 0), 1), 10)

# Define divergence in terms of initial coordinates
divergence = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
V = FunctionSpace(mesh, "DG", 0)
divfunc = project(divergence, V)
File("divfunc.pvd") << divfunc
3

Thanks! That’s much easier than I thought :slight_smile: