UserExpression in dolfinx

dolfinx
1

Hi all,

I am trying to implement a function that depends on the spatial position. I can do it in fenics using the following code:

from dolfin import *

mesh = UnitSquareMesh(16,16)
V = FunctionSpace(mesh,'CG',1)

p = Function(V)
t = 0
class InitialCondition(UserExpression):
    def eval_cell(self, value, x, ufc_cell):
        if x[0] <= 0.01:
            value[0] = 1.0
        else:
            value[0] = 0.0

p.interpolate(InitialCondition())
xdmf = XDMFFile('p.xdmf')
xdmf.write(p, t)

However, I am having problems when I try to rewrite this code in dolfinx, since UserExpression is not defined.
Is there any alternative function in dolfinx?

Thanks!

2

You can use lambda-functions, together with dolfinx.fem.locate_dofs_geoemtrical as shown in: Combining Dirichlet and Neumann conditions — FEniCSx tutorial

3

Thank you @dokken for the speedy reply.

I’m using that approach to impose boundary conditions. However, what I want to do is to impose an initial condition. Therefore, instead of using bc = DirichletBC() and then apply this bc to the linear algebra structures, I need to store the values of the initial condition in a vector, say p_n.

Is then possible to define such kind of initial conditions in dolfinx? (That is, initial conditions over part of the domain as defined in the previous code).

Thanks again!

4

Simply use interpolate

def initial_condition(x):
     return 1.0*(x[0]<0.01)
u.interpolate(initial_condition)
5

Thank you very much!
It works.

6

Thanks for your help @dokken. I have the same problem with assigning initial values to Function. The difference between my question and the one you answered is that I need the real geometric coordinates for another package to calculate some quantities I need for initial values. For example, in the following code:

from fenics import *
import numpy as np
from scipy.spatial import KDTree

mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 0.1), 10, 10, 1)

PF = VectorFunctionSpace(mesh, 'CG', 1, dim=5)
phi   = Function(PF, name='order parameters')           # Order parameters

points = np.random.rand(10,2) 
tree = KDTree(points)

class InitialConditions(UserExpression):
    def eval(self, values, x):
        distance, grain_index = tree.query([x[0], x[1]])
        for i in range(5):
            values[i] = 0.0
        values[grain_index % 5] = 1.0*np.exp(-1.0*(distance/0.1)**2)
    def value_shape(self):
        return (5,)
    
phi.interpolate(InitialConditions())
xdmf = XDMFFile('p.xdmf')
xdmf.write(phi)

I am converting it from FEniCS to FEniCSx. Follow your strategy, I write the following code:

def initial_condition(x):
     distance, grain_index = tree.query(x[0], x[1])
     for i in range(num_phi):
          values[i] = 0.0
     values[grain_index % num_phi] = 1.0 * np.exp(-1.0*(distance/0.1)**2)
     return values

phi.interpolate(initial_condition)

However, I found that the parameter I transfer into the tree.query() function is not the geometric coordinates. Could you please tell me how I can get the node coordinates and transfer it into the tree.query() function?

7

I think the discussion here may help

8

You haven’t defined values in this code. Please also supply a minimal reproducible example