Equivalent for Expression in dolfinx

[tasora]

Hi, I am trying to define a time-dependent boundary condition in dolfinx. Previously, in dolfin, it was possible through the Expression object, for instance:

f = Expression(“sin(2.0*t)”, t=0.0, degree = 1)

What is the equivalent now in dolfinx? I couldn’t import from somewhere the Expression object, so I suppose this changed.

Also, in the for loop afterwards, it was possible to update the boundary condition by just updating the expression and reassembling the required matrix/vector, for instance:

while t<T:
f.t = t
b = assemble(L)
…

Is this different now in dolfinx?

Thank you in advance!

Tasos

[dokken]

The Expression class has been removed, as it caused a headache for developers and it is not clear what it does.
The new way of doing this is to interpolate an expression into a function:

from dolfinx import UnitSquareMesh, Function, FunctionSpace, MPI, Constant
from dolfinx.fem import assemble_vector
from ufl import TestFunction, dx, inner
import numpy as np

class MyExpression:
    def __init__(self):
        self.t = 0.0

    def eval(self, x):
        # Added some spatial variation here. Expression is sin(t)*x
        return np.full(x.shape[1], np.sin(self.t)*x[0])

# If your Expression is not spatially dependent, use Constant
# f = Constant(mesh, 0)
# L = inner(f, v)*dx
# f.value = np.sin(2)
mesh = UnitSquareMesh(MPI.comm_world, 10, 10)
V = FunctionSpace(mesh, ("CG", 1))
f = MyExpression()
f.t = 0
w = Function(V)
v = TestFunction(V)
L = inner(w, v)*dx
w.interpolate(f.eval)

print(assemble_vector(L).array)
f.t = 2
w.interpolate(f.eval)
print(assemble_vector(L).array)

[tasora]

Exactly what i needed ! Thanks !

[marchirschvogel]

Hello,

I would have a follow-up question to this topic about expression handlings in dolfin-x.

Assume a variational form with a time dependency, i.e. a traction/body force that is ramped up with time (no spatial dependency).

Please find this minimal example (adopted and modified from the dolfin-x elasticity demo):

#!/usr/bin/env python3

import time
import sys, os, subprocess, time
import math

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import *
from dolfinx.mesh import *
from dolfinx.cpp.mesh import *
from dolfinx.fem import locate_dofs_topological, locate_dofs_geometrical, assemble_matrix, assemble_vector, set_bc, apply_lifting
from dolfinx.io import *
from dolfinx.la import *
from ufl import *

# mesh = UnitCubeMesh(2, 2, 2)
mesh = BoxMesh(
    MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]),
                     np.array([2.0, 1.0, 1.0])], [12, 12, 12],
    CellType.tetrahedron, dolfinx.cpp.mesh.GhostMode.none)

x = SpatialCoordinate(mesh)
f = as_vector((1.0 * 1.0**2 * x[0], 1.0 * 1.0**2 * x[1], 0.0))

def sigma(v):
    return 2.0 * 1.0 * sym(grad(v)) + 1.0 * tr(sym(grad(v))) * Identity(
        len(v))

def trueload(t):
    return Constant(mesh, 10.0*t)
    #return np.float(10.0*t)

# Create function space
V = VectorFunctionSpace(mesh, ("Lagrange", 1))
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = inner(sigma(u), grad(v)) * dx
t=1.0
L = trueload(t) * inner(f, v) * dx

sys.exit()

In a more complex example, trueload(t) would be in a time loop and updated. Executing this example yields:
L = trueload(t) * inner(f, v) * dx
File “/home/mh/.local/lib/python3.6/site-packages/ufl/exproperators.py”, line 182, in _mul
return _mult(self, o)
File “/home/mh/.local/lib/python3.6/site-packages/ufl/exproperators.py”, line 114, in _mult
r1, r2 = len(s1), len(s2)
TypeError: object of type ‘float’ has no len()

while in dolfin, this works (using Constant(10.0*t)).

A “solution”: return return np.float(10.0*t) instead of return Constant(mesh, 10.0*t).

BUT: When the value is not declared as constant, the form compiler has to be called after each time step and the solution to a small problem takes excessively long (form compiler takes three times longer than the whole time step). This was the same for dolfin when returning a value not declared as constant…

So: Does anyone know how to return a constant and multiplying it to a form in dolfin-x?

Thanks!

Marc

[dokken]

Consider the following:

from dolfinx import UnitSquareMesh, Constant, FunctionSpace
from dolfinx.fem import assemble_matrix
from mpi4py import MPI
from ufl import TrialFunction, TestFunction, inner, dx

mesh = UnitSquareMesh(MPI.COMM_WORLD, 10,10)
V = FunctionSpace(mesh, ("CG", 1))
u, v = TrialFunction(V), TestFunction(V)
t = 0.1
c = Constant(mesh, 10*t)
a = c*inner(u, v)*dx
A0 = assemble_matrix(a)
A0.assemble()
t += 0.2
c.value = 10*t
A1 = assemble_matrix(a)
A1.assemble()
print(A0.norm(), A1.norm())

returning

0.05157249482255295 0.15471748446765723

[marchirschvogel]

Thanks for your reply!

I played around a little bit and indeed, also my approach/example works, if I use
from ufl import TrialFunction, TestFunction, inner, dx, SpatialCoordinate, sym, as_vector, grad, tr, Identity
instead of importing everything:
from ufl import *

With the specific import, my example does not return an error, while with the all import, it returns the reported one.

No idea why this is happening, something that is imported from ufl seems to cause a conflict (?).

Thx!
Marc

[dokken]

In general, one should never use wildcard (*) imports on more than one package, as there might be duplicates of namespaces. Even when it comes to dolfinx, dolfinx.cpp, dolfinx.fem, I would advise to not use wildcard imports, and rather import the functions you Need. This also matches the python styling convention

[marchirschvogel]

Ok thanks, I will consider this.

I tried to execute your example with the value re-assignment, but I’m getting a SIGSEGV segmentation error.

Actually, I cannot execute the code snippet

from mpi4py import MPI
from dolfinx import UnitSquareMesh, Constant

mesh = UnitSquareMesh(MPI.COMM_WORLD, 10,10)
b = Constant(mesh, 0)
b.value = 1

which yields
Loguru caught a signal: SIGSEGV
Segmentation fault (core dumped)

due to line b.value = 1. There seems to be some access violation when trying to change the Constant.
I tried it on two different machines, both with Ubuntu 18.04 LTS gcc-7.5.0. With a fresh installation of dolfin-x.
Have you had a similar problem already?

Thx!
Marc

P.S.: The dolfin equivalent works:

from dolfin import UnitSquareMesh, Constant

mesh = UnitSquareMesh(10,10)
b = Constant(0)
b.value = 1

[dokken]

I have not encountered this problem (I am using the docker image dolfinx/real, which uses: Ubuntu 20.04 and gcc 9.3.0).

[marchirschvogel]

Ok, finally I found a solution for a time dependency in a variational form that does not lead to form compiler overheads which I’d like to share.

Define an expression class for your time-dependent value:

class load_expr:
    def __init__(self):
        self.t = 0.0

    def eval(self, x):
        pr = 10.0
        # Added some spatial/temporal variation here
        return np.full(x.shape[1], pr*self.t/T)

Define a real function space (one global degree of freedom):
V_real = FunctionSpace(mesh, ("R", 0))

Interpolate the expression into the real function space:

load = load_expr()
load.t = 0.0
load_func = Function(V_real)
load_func.interpolate(load.eval)

In your variational form (can be defined outside of time loop):
varform = ... + load_func*J*dot(dot(inv(F).T,n0), var_u)*ds4

inside the time loop, call

load.t = t
load_func.interpolate(load.eval)

in order to update the expression.

Like this, there is no excessive repeated form compilation after each time step which was the case if I had a Constant(mesh, 10.*t/T) returned for every step.
(Maybe it could have been circumvented with updating the Constant with .value, however this produced the segfault as mentioned above.)

I hope this helps when someone has a similar problem.

Best,
Marc

[cvillette]

@dokken

print(assemble_vector(L).array)

Hello,
I am trying to run this minimal example in dolfinx 0.3.1.dev0. I had to change the module import calls a bit (see below). However, I run into the following error, any idea what might be causing it? Many thanks!

from dolfinx import mesh, fem
from dolfinx.fem import assemble_vector,FunctionSpace, Function,Constant
from ufl import TestFunction, dx, inner
import numpy as np
from mpi4py import MPI


class MyExpression:
    def __init__(self):
        self.t = 0.0

    def eval(self, x):
        # Added some spatial variation here. Expression is sin(t)*x
        return np.full(x.shape[1], np.sin(self.t)*x[0])

# If your Expression is not spatially dependent, use Constant
# f = Constant(mesh, 0)
# L = inner(f, v)*dx
# f.value = np.sin(2)
mesh = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(mesh, ("CG", 1))
f = MyExpression()
f.t = 0
w = Function(V)
v = TestFunction(V)
L = inner(w, v)*dx
w.interpolate(f.eval)
print(assemble_vector(L))

[dokken]

Here is a minimal updated example for the latest version of DOLFINx:

from dolfinx import mesh, fem
from dolfinx.fem import FunctionSpace, Function,Constant
from ufl import TestFunction, dx, inner
import numpy as np
from mpi4py import MPI


class MyExpression:
    def __init__(self):
        self.t = 0.0

    def eval(self, x):
        # Added some spatial variation here. Expression is sin(t)*x
        return np.full(x.shape[1], np.sin(self.t)*x[0])

# If your Expression is not spatially dependent, use Constant
# f = Constant(mesh, 0)
# L = inner(f, v)*dx
# f.value = np.sin(2)
mesh = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(mesh, ("CG", 1))
f = MyExpression()
f.t = 0
w = Function(V)
v = TestFunction(V)
L = inner(w, v)*dx
w.interpolate(f.eval)
print(fem.petsc.assemble_vector(fem.form(L)))

[cvillette]

Perfect, thank you so much!

[ampdes]

UFLValueError for vector valued expressions.

I am trying to interpolate into the vector-valued Function object, different expressions for different components.

Following the solution template for the scalar case, please find below a MWE. This throws UFLValueError.

I appreciate any directions on how to implement this correctly. Thanks.

I am using the docker image for the dolfinx-tutorial.

import numpy as np
import ufl
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.mesh import (CellType, create_rectangle)
from dolfinx.fem import (Function, FunctionSpace)

# ----------------------------------------------------
domain = create_rectangle(MPI.COMM_WORLD,
                      [np.array([0.0,0.0]), np.array([1.0,1.0])],
                      [2**5, 2**5],
                      CellType.triangle)

P2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2)
V = FunctionSpace(domain, P2)

# -----------------------------------------------------
class u_expression():
    def __init__(self, Ux=10.0, Uy=10.0, A=1.0, t=0.0):
        self.Ux = Ux
        self.Uy = Uy
        self.A = A
        self.t = t
        
    def eval(self, x):
        return ufl.as_vector((-self.Ux*np.exp(-self.A*self.t)*x[0]*(1-x[0])*x[1]*(1-x[1]),
                              -self.Uy*np.exp(-self.A*self.t)*x[0]*(1-x[0])*x[1]*(1-x[1]))
                            )

# ------------------------------------------------------
u_e = u_expression()

u_D = Function(V)
u_D.interpolate(u_e.eval)

Find below the full error message

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File /usr/local/lib/python3.9/dist-packages/dolfinx/fem/function.py:318, in Function.interpolate.<locals>._interpolate(u, cells)
    317 try:
--> 318     self._cpp_object.interpolate(u._cpp_object, cells)
    319 except AttributeError:

AttributeError: 'function' object has no attribute '_cpp_object'

During handling of the above exception, another exception occurred:

UFLValueError                             Traceback (most recent call last)
Input In [2], in <cell line: 35>()
     32 u_e = u_expression()
     34 u_D = Function(V)
---> 35 u_D.interpolate(u_e.eval)

File /usr/local/lib/python3.9/dist-packages/dolfinx/fem/function.py:336, in Function.interpolate(self, u, cells)
    333     map = mesh.topology.index_map(mesh.topology.dim)
    334     cells = np.arange(map.size_local + map.num_ghosts, dtype=np.int32)
--> 336 _interpolate(u, cells)

File /usr/lib/python3.9/functools.py:877, in singledispatch.<locals>.wrapper(*args, **kw)
    873 if not args:
    874     raise TypeError(f'{funcname} requires at least '
    875                     '1 positional argument')
--> 877 return dispatch(args[0].__class__)(*args, **kw)

File /usr/local/lib/python3.9/dist-packages/dolfinx/fem/function.py:320, in Function.interpolate.<locals>._interpolate(u, cells)
    318     self._cpp_object.interpolate(u._cpp_object, cells)
    319 except AttributeError:
--> 320     self._cpp_object.interpolate(u, cells)

Input In [2], in u_expression.eval(self, x)
     26 def eval(self, x):
---> 27     return ufl.as_vector((-self.Ux*ufl.exp(-self.A*self.t)*x[0]*(1-x[0])*x[1]*(1-x[1]),
     28                           -self.Uy*ufl.exp(-self.A*self.t)*x[0]*(1-x[0])*x[1]*(1-x[1]))
     29                         )

File /usr/local/lib/python3.9/dist-packages/ufl/tensors.py:313, in as_vector(expressions, index)
    310         error("Expecting a single Index object.")
    311     index = (index,)
--> 313 return as_tensor(expressions, index)

File /usr/local/lib/python3.9/dist-packages/ufl/tensors.py:243, in as_tensor(expressions, indices)
    239         error("Expecting nested list or tuple.")
    241     # Recursive conversion from nested lists to nested ListTensor
    242     # objects
--> 243     return _as_list_tensor(expressions)
    244 else:
    245     # Make sure we have a tuple of indices
    246     if isinstance(indices, list):

File /usr/local/lib/python3.9/dist-packages/ufl/tensors.py:191, in _as_list_tensor(expressions)
    189 def _as_list_tensor(expressions):
    190     if isinstance(expressions, (list, tuple)):
--> 191         expressions = [_as_list_tensor(e) for e in expressions]
    192         return ListTensor(*expressions)
    193     else:

File /usr/local/lib/python3.9/dist-packages/ufl/tensors.py:191, in <listcomp>(.0)
    189 def _as_list_tensor(expressions):
    190     if isinstance(expressions, (list, tuple)):
--> 191         expressions = [_as_list_tensor(e) for e in expressions]
    192         return ListTensor(*expressions)
    193     else:

File /usr/local/lib/python3.9/dist-packages/ufl/tensors.py:194, in _as_list_tensor(expressions)
    192     return ListTensor(*expressions)
    193 else:
--> 194     return as_ufl(expressions)

File /usr/local/lib/python3.9/dist-packages/ufl/constantvalue.py:437, in as_ufl(expression)
    435     return IntValue(expression)
    436 else:
--> 437     raise UFLValueError("Invalid type conversion: %s can not be converted"
    438                         " to any UFL type." % str(expression))

UFLValueError: Invalid type conversion: [-0.00000000e+00 -0.00000000e+00 -3.36102673e-17 ... -1.70761842e-17
 -2.36570835e-03 -0.00000000e+00] can not be converted to any UFL type.

[dokken]

@ampdes

    def eval(self, x):
        return ufl.as_vector((-self.Ux*np.exp(-self.A*self.t)*x[0]*(1-x[0])*x[1]*(1-x[1]),
                              -self.Uy*np.exp(-self.A*self.t)*x[0]*(1-x[0])*x[1]*(1-x[1]))
                            )

Simply replace this with:

class u_expression():
    def __init__(self, Ux=10.0, Uy=10.0, A=1.0, t=0.0):
        self.Ux = Ux
        self.Uy = Uy
        self.A = A
        self.t = t

    def eval(self, x):
        return (-self.Ux * np.exp(-self.A * self.t) * x[0] * (1 - x[0]) * x[1] * (1 - x[1]),
                -self.Uy * np.exp(-self.A * self.t) * x[0] * (1 - x[0]) * x[1] * (1 - x[1]))

The eval acts on numpy arrays, no ufl should be involved.

[Thomas996]

Hello,
I am following this code trying to implement it in a Vectorfunctionspace (Example linear elasticity and have following error:

Output exceeds the size limit. Open the full output data in a text editor
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/bcs.py:125, in DirichletBCMetaClass.__init__(self, value, dofs, V)
    124 try:
--> 125     super().__init__(_value, dofs, V)  # type: ignore
    126 except TypeError:

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.DirichletBC_float64(g: numpy.ndarray[numpy.float64], dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    2. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Constant<double>, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: numpy.ndarray[numpy.int32])
    4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: List[numpy.ndarray[numpy.int32][2]], V: dolfinx::fem::FunctionSpace)

Invoked with: <dolfinx.cpp.fem.Function_float64 object at 0x7fc437d0fdf0>, array([3839, 3840, 3841, 3842, 3844, 3845, 3847, 3849, 3853, 4081, 4082,
       4083, 4085, 4086, 4089, 4093, 4094, 4095, 4096, 4098, 4101, 4305,
       4306, 4307, 4309, 4310, 4313, 4317, 4318, 4319, 4321, 4325, 4326,
       4327, 4328, 4330, 4333, 4505, 4506, 4507, 4509, 4510, 4513, 4517,
       4518, 4519, 4521, 4525, 4526, 4527, 4529, 4533, 4534, 4535, 4536,
       4538, 4541, 4673, 4674, 4675, 4677, 4678, 4681, 4685, 4686, 4687,
       4689, 4693, 4694, 4695, 4697, 4701, 4702, 4703, 4705, 4709, 4710,
       4711, 4712, 4714, 4717, 4801, 4802, 4803, 4805, 4809, 4810, 4811,
       4813, 4817, 4818, 4819, 4821, 4825, 4826, 4827, 4829, 4881, 4882,
       4883, 4885, 4889, 4890, 4891, 4893, 4897, 4898, 4899, 4901, 4929,
       4930, 4931, 4933, 4937, 4938, 4939, 4941, 4953, 4954, 4955, 4957],
      dtype=int32), FunctionSpace(Mesh(VectorElement(Basix element (P, hexahedron, 1, gll_warped, unset, False), 3), 12), VectorElement(Basix element (P, hexahedron, 2, gll_warped, unset, False), 3))
...
       4711, 4712, 4714, 4717, 4801, 4802, 4803, 4805, 4809, 4810, 4811,
       4813, 4817, 4818, 4819, 4821, 4825, 4826, 4827, 4829, 4881, 4882,
       4883, 4885, 4889, 4890, 4891, 4893, 4897, 4898, 4899, 4901, 4929,
       4930, 4931, 4933, 4937, 4938, 4939, 4941, 4953, 4954, 4955, 4957],
      dtype=int32), <dolfinx.cpp.fem.FunctionSpace object at 0x7fc437838430>

And here is my code:

import numpy as np
import ufl

from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from dolfinx.fem import Function, FunctionSpace, Constant
from ufl import TestFunction, dx, inner


L = 20.0
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0],[L,1,1]],[20,5,5], mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 2))

def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], L)

fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with 2.as_integer_ratio
marked_facets = np.hstack([left_facets,right_facets])


marked_values = np.hstack([np.full_like(left_facets,1), np.full_like(right_facets,2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets],marked_values[sorted_facets])

# boundary condition on the left side: fixed
u_bc = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)
left_dofs = fem.locate_dofs_topological(V,facet_tag.dim, facet_tag.find(1))
bc_l = fem.dirichletbc(u_bc, left_dofs, V)


# B.C on the right side = a function to interpolate
import numpy as np
import ufl

from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from dolfinx.fem import Function, FunctionSpace, Constant
from ufl import TestFunction, dx, inner

from mpi4py import MPI
scale = 0.5
y0 = 0.5
z0 = 0.5

class MyExpression:
    def __init__(self):
        self.theta = np.pi/3

    def eval(self, x):
        return  (np.zeros_like(x[1]),
                 scale * (y0 + (x[1] - y0) * np.cos(self.theta) - (x[2] - z0) * np.sin(self.theta) - x[1]),
                 scale * (z0 + (x[1] - y0) * np.sin(self.theta) + (x[2] - z0) * np.cos(self.theta) - x[2]))

f = MyExpression()    
u_D = Function(V)
u_D.interpolate(f.eval)

# implement this function as DirichletBC at the right side of the beam
right_dofs = fem.locate_dofs_topological(V,facet_tag.dim, facet_tag.find(2))
bc_r = fem.dirichletbc(u_D,right_dofs, V)
bcs =  [bc_l, bc_r]   

It seems the problem is the u_D has a incompatible format for DirhicletBC? Can someone pls expain it?

[dokken]

@Thomas996

bc_r = fem.dirichletbc(u_D,right_dofs, V)

By just looking over this quickly, you do not need to send in V to the dirichletbc constructor, as u_D is in the same space as your variational form will be, i.e.
bc_r = fem.dirichletbc(u_D, right_dofs) should be sufficient

[Thomas996]

Thank you for your quick answer. Yes it is. The code works now.

[Thomas996]

Hi, dokken, when I implement such rotation as in code I posted before, i got following result:

the rotation was as expected at the boundary; but it seems only the surface has rotated the expected angle, and all the other cells except the of the last layer has deformed not enough. Which caused the strange torsion of last layer. What could be a solution here? I’d appreciate your help.

[dokken]

Please post a code that can reproduce the result.