Projecting regular functions over meshes to write XDMF files

Hello FEniCS community,

I have a dynamic PDEs system solved by FEniCS. The computed solution is storaged in Finite-Element (FE) functions (<class ‘dolfin.function.function.Function’>): u1, u2, u3, ...., un.

I can quickly write this Function objects in XDMF files for future visualization, (eg. with Paraview) as follows:

xdmffile_Un = XDMFFile('./Un.xdmf')
xdmffile_Un.write(Un, t)

However, I would also like to visualize other functions which depend on the computed solution:

Eg: " Arbitrary_function = f(u1, u2, u3, ...., un) "

I think that I cannot define an expression (<class ‘dolfin.function.expression.Expression’>) for each function since these do not depend directly on the spatial position but the computed solution. Basically, I need to project a bunch of numbers on the mesh, and these values depend on the solution calculated for each node and time during the simulation. Once I have the numbers projected over the mesh, I guess that I could write an XDMF file the same way as shown above.

Any idea about a possible & efficient way to do this?

Thank you for your time and I look forward to hearing from you,
Santiago

You actually can use the computed solution in an Expression. See the following example:

from dolfin import *
mesh = UnitIntervalMesh(100)
x = SpatialCoordinate(mesh)

# Obtain some non-trivial Functions:
V = FunctionSpace(mesh,"CG",1)
u1 = project(cos(7.0*pi*x[0]),V)
u2 = project(exp(x[0]),V)

# Use Functions in Expression:
expr = Expression("u1/u2",degree=1,u1=u1,u2=u2)

# Function ready to visualize:
expr_viz = interpolate(expr,V)

# Interactive plotting:
from matplotlib import pyplot as plt
plot(expr_viz)
plt.show()
3 Likes

Nice Kamensky and thank you for your fast reply!
In that case, I should change my question:

The mathematical expressions I have to deal with in this case are much more complex to be coded through Simple user-defined JIT-compiled expressions.
Then, How can I define User-defined expressions with classes in python to formulate long mathematical equations?

For instance, in this simple example (I know that I am making a mistake with the class definition):

from dolfin import UnitIntervalMesh, SpatialCoordinate, FunctionSpace, \
                   project, Expression, interpolate, UserExpression
from ufl.operators import cos, exp
from numpy import pi

mesh = UnitIntervalMesh(100)`
x = SpatialCoordinate(mesh)

# Obtain some non-trivial Functions:
V = FunctionSpace(mesh, "CG", 1)
u1 = project(cos(7.0*pi*x[0]), V)
u2 = project(exp(x[0]), V)

class MyExpression(UserExpression):
    def eval(self, u1, u2):
        value = u1/u2
        return value

f0 = MyExpression(degree=0)
Val = f0.eval(u1=u1, u2=u2)
expr_Val = interpolate(Val, V)
print(expr_Val(0))

Thank you in advance,

Santiago

You can pass the Functions to a custom constructor, make them attributes of the subclass, then evaluate the attributes in the eval method. See the following example:

from dolfin import UnitIntervalMesh, SpatialCoordinate, FunctionSpace, \
                   project, Expression, interpolate, UserExpression, plot
from ufl.operators import cos, exp
from numpy import pi

mesh = UnitIntervalMesh(100)
x = SpatialCoordinate(mesh)

# Obtain some non-trivial Functions:
V = FunctionSpace(mesh, "CG", 1)
u1 = project(cos(7.0*pi*x[0]), V)
u2 = project(exp(x[0]), V)

class MyExpression(UserExpression):
    def __init__(self,u1,u2,**kwargs):
        self.u1 = u1
        self.u2 = u2
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = self.u1(x)/self.u2(x)
    def value_shape(self):
        return ()
        
f0 = MyExpression(u1,u2,degree=1)

from matplotlib import pyplot as plt
plot(interpolate(f0,V))
plt.show()
1 Like