Defining a piecewise constant parameter over all elements

variational formulation
1

Hello everyone,

I have a parameter (that I call gamma) in my pde for which I want to take a different constant value in each element of the mesh. To be precise, I want it to have a constant value of 0.8 in each element if z-coordinate of the centroid of the element is greater than 0.15 and a random constant value otherwise. I am not sure exactly how to achieve that. Currently I have this code:

randvalues=np.random.rand(numelem,1)   #numelem is the total number of elements
class Gam(UserExpression):
	def __init__(self, mesh, **kwargs):
		self.mesh = mesh
		super().__init__(**kwargs)
	def eval_cell(self, value, x, ufc_cell):
		cell = Cell(self.mesh, ufc_cell.index)
		if (cell.midpoint().z() > 0.15):
			value[0]=0.8
		else:
			value[0] = 0.8+(1-0.8)*randvalues[ufc_cell.index] 

elem = FiniteElement('DG', triangle, 0) #Actual finite element for PDE is CG1
gamma = Gam(mesh,element=elem)

Y= FunctionSpace(mesh, "CG", 1)
y  = TestFunction(Y)             # Test function
z  = Function(Y)
R_z = gamma*y*2*z*(psi1)*dx

Does this look correct? I am not sure about whether the expression is being evaluated at the nodes or the quadrature points and if there is any interpolation involved.

Thanks!

2

During assembly gamma is evaluated at the quadrature points given by the quadrature rule and the polynomial degree of R_z. Based on the rule some of the quadrature points may coincide with nodes. This is easy to test

from dolfin import *
from random import random

class Gam(UserExpression):
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super().__init__(**kwargs)
    def eval_cell(self, value, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        if (cell.midpoint().x() > 0.5):
            value[0]=0.8
        else:
            value[0] = random()

        print('Eval at %r of cell %d -> %g' % (x, ufc_cell.index, value[0]))


mesh = UnitIntervalMesh(10)

gamma = Gam(mesh, degree=3)  # Change degree here to change degree of R_z

V = FunctionSpace(mesh, 'CG', 1)
v  = TestFunction(V)
L = gamma*v*dx   # Polyn degree of L is 1 + degree

assemble(L)

You are on a right track with your solution.

3

Thank you for your reply, Miroslav. I just want to clarify one thing with you. I ran this same code in 2D by changing to UnitSquareMesh. The print command:

print('Eval at %r of cell %d -> %g' % (x, ufc_cell.index, value[0]))

prints x as the nodal coordinates (and not the coordinates of the Gauss quadrature points) corresponding to the degree for degree>0. So does that mean the Expression is first being interpolated into a finite element space?

For my problem, I just want a constant value for the parameter gamma, so I guess I can put degree=0. But I am just curious about how exactly UserExpressions work.

Btw do you know how to get the actual ufc_cell.index in parallel? It shows the index for the cell allocated in each process.

Thanks again!