Convert Function to UFL type

dolfinx
1

Hello all,

I am trying to compute the reaction force on a boundary in a linear elasticity problem, but have run into some trouble. Here is the MWE.

# Scaled variable
L = 1
W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

import numpy as np
import ufl

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

from dolfinx import mesh, fem, plot, io

domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],
                  [20,6,6], cell_type=mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 1))

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

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0,0,0], dtype=ScalarType)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, ScalarType((0, 0, 0)))

ds = ufl.Measure("ds", domain=domain)

subdomain_clamp = mesh.meshtags(domain, fdim, np.sort(boundary_facets), np.full(boundary_facets.shape[0],1))
ds_clamp = ufl.Measure("ds", domain=domain, subdomain_data=subdomain_clamp)

def epsilon(u):
    return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu*epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, ScalarType((0, 0, -rho*g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

with io.XDMFFile(domain.comm, "deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)

nn = ufl.FacetNormal(domain)
flux_1 = fem.form(ufl.dot(ufl.dot(sigma,nn),nn))
TT_1 = fem.assemble_scalar(flux_1 * ds_clamp(1)) #Force resultant clamped boundary
print(TT_1)

When I run this, I get the following error.


flux_1 = fem.form(ufl.dot(ufl.dot(sigma,nn),nn))
  File "/usr/local/lib/python3.9/dist-packages/ufl/operators.py", line 172, in dot
    a = as_ufl(a)
  File "/usr/local/lib/python3.9/dist-packages/ufl/constantvalue.py", line 437, in as_ufl
    raise UFLValueError("Invalid type conversion: %s can not be converted"
ufl.log.UFLValueError: Invalid type conversion: <function sigma at 0x7f0e4d6005e0> can not be converted to any UFL type.

Can anyone please help me figure out what is the correct way to do this?

Thank you

2 
flux_1 = fem.form(ufl.dot(ufl.dot(sigma,nn),nn))

this isn’t a valid form since there’s no integration measure. I assume you’re interested in something like

\text{flux}_1 = \int_{\partial\Omega} \sigma \cdot \vec{n} \cdot \vec{n} \; \mathrm{d} s

which simply requires:

flux_1 = fem.form(ufl.dot(ufl.dot(sigma,nn),nn) * ufl.ds)

Edit: Reading your code in a little more detail, it looks like you’ve just got to move your measure ds_clamp(1) from the line

TT_1 = fem.assemble_scalar(flux_1 * ds_clamp(1))
3

Hello Nate,

Thank you for the reply. As per your suggestion, I tried

flux_1 = fem.form(ufl.dot(ufl.dot(sigma,nn),nn) * ds_clamp(1))
TT_1 = fem.assemble_scalar(flux_1)

But I still got the same error. It seems like the problem is that sigma in the first line is a function and not any UFL type?

Thank you

4

Ah. Well sigma(u) is a function. You should pass the appropriate argument. In this case I’m going to assume it should be sigma(uh).

5

Ok. I tried with
flux_1 = fem.form(ufl.dot(ufl.dot(sigma(uh),nn),nn) * ds_clamp(1))

but looks like that’s not the right way either. Now, I am getting an error about the argument being of an incompatible type.

flux_1 = fem.form(ufl.dot(ufl.dot(sigma(uh),nn),nn) * ds_clamp(1))
  File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/forms.py", line 139, in form
    return _create_form(form)
  File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/forms.py", line 134, in _create_form
    return _form(form)
  File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/forms.py", line 128, in _form
    return formcls(ufcx_form, V, coeffs, constants, subdomains, mesh, code)
  File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/forms.py", line 48, in __init__
    super().__init__(ffi.cast("uintptr_t", ffi.addressof(self._ufcx_form)),
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.Form_float64(spaces: List[dolfinx::fem::FunctionSpace], integrals: Dict[dolfinx::fem::IntegralType, Tuple[List[Tuple[int, object]], dolfinx.cpp.mesh.MeshTags_int32]], coefficients: List[dolfinx.cpp.fem.Function_float64], constants: List[dolfinx.cpp.fem.Constant_float64], need_permutation_data: bool, mesh: dolfinx.cpp.mesh.Mesh = None)
    2. dolfinx.cpp.fem.Form_float64(arg0: int, arg1: List[dolfinx::fem::FunctionSpace], arg2: List[dolfinx.cpp.fem.Function_float64], arg3: List[dolfinx.cpp.fem.Constant_float64], arg4: Dict[dolfinx::fem::IntegralType, dolfinx.cpp.mesh.MeshTags_int32], arg5: dolfinx.cpp.mesh.Mesh)

Invoked with: <cdata 'uintptr_t' 140425070601248>, [], [<dolfinx.cpp.fem.Function_float64 object at 0x7fb74ce902f0>], [], {<IntegralType.cell: 0>: None, <IntegralType.exterior_facet: 1>: <dolfinx.mesh.MeshTagsMetaClass object at 0x7fb742c31220>, <IntegralType.interior_facet: 2>: None, <IntegralType.vertex: 3>: None}, <dolfinx.mesh.Mesh object at 0x7fb75431ce00>

Thank you

6

Following the implementation given in Implementation — FEniCSx tutorial, I added the first four lines given below to project sigma to an appropriate FunctionSpace before trying to compute the flux.

V_stress = fem.FunctionSpace(domain, ("DG",0))
sigma_expr = fem.Expression(sigma(uh), V_stress.element.interpolation_points)
stress = fem.Function(V_stress)
stress.interpolate(sigma_expr)

flux_1 = fem.form(ufl.dot(ufl.dot(stress,nn),nn)*ds_clamp(1))
TT_1 = fem.assemble_scalar(flux_1) #Force resultant clamped boundary
print(TT_1)

But now I am getting a runtime error while trying to interpolate.

    stress.interpolate(sigma_expr)
  File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/function.py", line 336, in interpolate
    _interpolate(u, cells)
  File "/usr/lib/python3.9/functools.py", line 877, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/function.py", line 329, in _
    self._cpp_object.interpolate(expr._cpp_object, cells)
RuntimeError: Function value size not equal to Expression value size

Can someone please tell me what the error means?

Thank you

7

Interested! I am facing the similar problem. I wonder whether you did find a way to solve this issue.

In the given implementation example, the first argument of fem.Expression is a scalar value for each cell which is Von Misses stress. Yet, you’re trying to project a function onto the finite element domain, that’s why it gives the error. Please correct me if I am wrong.

1 Like
8

I agree with @yucesoya

9

Hello,
I’m also facing this problem when I tried to plot the strains. So according to @yucesoya we should not pass the function in the Expression. Then what we could do with that function?

Here is my code:

epsilon = eps(uh)
V_eps = fem.FunctionSpace(domain,("DG",0))
strain_expr = fem.Expression(epsilon, V_eps.element.interpolation_points())
strains = fem.Function(V_eps)
strains.interpolate(strain_expr)
10

The point is that sigma is a tensor (same for eps(uh)) and thus the function space is wrong. You should either:

  1. extract a component of epsilon, say epsilon[0,0]
  2. Define V_eps as a TensorFunctionSpace, which would interpolate a tensor expression into a tensor
11

Thank you for such a quick answer. I changed the FunctionSpace to a TensorFunctionSpace and it worked. Now I am a little bit confused about how to plot this now. Because when I try to do following, it raised the Error: data length of (24741) != required length (8043)

warped.cell_data["epsilon"] = strains.vector.array
warped.set_active_scalars("Epsilon")
p = pyvista.Plotter()
p.add_mesh(warped)
p.show_axes()
if not pyvista.OFF_SCREEN:
   p.show()
else:
   stress_figure = p.screenshot(f"Strains.png")
12

Well, it is a tensor, how should that be visualized?

Are you sure you dont want to work on each individual component of the tensor and visualize it?

13

I think you are right. What I want to do is to visualize the epsilon of every point of the whole mesh. Just like in the example of lin.elasticity the visualization of Von Mieses stresses and Displacement. But I am not sure how to do it.

14

As I stated above, then you should follow

and interpolate that into a function space. You can do this for each component and visualize them individually.

15

Thank you again, dokken. I have managed to run the code and visualize it successfully. But I got this TraitError at the end. What reason could it be ?

Output exceeds the size limit. Open the full output data in a text editor---------------------------------------------------------------------------
TraitError                                Traceback (most recent call last)
File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/ipywidgets/widgets/widget.py:756, in Widget._handle_msg(self, msg)
    754         if 'buffer_paths' in data:
    755             _put_buffers(state, data['buffer_paths'], msg['buffers'])
--> 756         self.set_state(state)
    758 # Handle a state request.
    759 elif method == 'request_state':

File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/ipywidgets/widgets/widget.py:625, in Widget.set_state(self, sync_data)
    622 if name in self.keys:
    623     from_json = self.trait_metadata(name, 'from_json',
    624                                     self._trait_from_json)
--> 625     self.set_trait(name, from_json(sync_data[name], self))

File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/traitlets/traitlets.py:1738, in HasTraits.set_trait(self, name, value)
   1736     raise TraitError(f"Class {cls.__name__} does not have a trait named {name}")
   1737 else:
-> 1738     getattr(cls, name).set(self, value)

File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/traitlets/traitlets.py:703, in TraitType.set(self, obj, value)
    702 def set(self, obj, value):
--> 703     new_value = self._validate(obj, value)
    704     try:
    705         old_value = obj._trait_values[self.name]

File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/traitlets/traitlets.py:735, in TraitType._validate(self, obj, value)
    733     return value
    734 if hasattr(self, "validate"):
--> 735     value = self.validate(obj, value)
    736 if obj._cross_validation_lock is False:
    737     value = self._cross_validate(obj, value)

File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/traitlets/traitlets.py:2302, in Union.validate(self, obj, value)
   2300         except TraitError:
   2301             continue
-> 2302 self.error(obj, value)

File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/traitlets/traitlets.py:841, in TraitType.error(self, obj, value, error, info)
    835 else:
    836     e = "The '{}' trait expected {}, not {}.".format(
    837         self.name,
    838         self.info(),
...
    839         describe("the", value),
    840     )
--> 841 raise TraitError(e)

TraitError: The 'target' trait of a DirectionalLight instance expected an Uninitialized or an Object3D, not the str 'IPY_MODEL_32c2dd9f-09e4-43ab-973b-decbbf240b0b'.

Here is my code:

import pyvista

pyvista.set_jupyter_backend('pythreejs')
# Create plotter and pyvista grid
p = pyvista.Plotter()
topology, cell_types, geometry = plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)


# Attach vector values to grid and warp grid by vector
grid["u"] = uh.x.array.reshape((geometry.shape[0], 3))
geometry.shape
actor_0 = p.add_mesh(grid, style = "wireframe", color="k")
warped = grid.warp_by_vector("u", factor = 1.0)
actor_1 = p.add_mesh(warped, show_edges = True)
p.show_axes()
if not pyvista.OFF_SCREEN:
    p.show()
else:
    figure_as_array = p.screenshot("deflection.png")

# Interpolate the Component of Tensor
epsilon = eps(uh)
epsilon_comp = epsilon[1,2]
V_eps = fem.FunctionSpace(domain,("DG",0))
strain_expr = fem.Expression(epsilon_comp, V_eps.element.interpolation_points())
strains = fem.Function(V_eps)
strains.interpolate(strain_expr)

# VIsualize
warped.cell_data["epsilon_comp"] = strains.vector.array
warped.set_active_scalars("epsilon_comp")
p = pyvista.Plotter()
p.add_mesh(warped)
p.show_axes()
if not pyvista.OFF_SCREEN:
   p.show()
else:
   stress_figure = p.screenshot(f"Strains.png")
16

The error happens inside ipywidgets. This means that

is not happening (either you haven’t installed pythreejs, or you are not running this as a notebook).

1 Like
17

When I try to print Component[0,0] of epsilon, it has following error and it only happend when I try to visualize [0,0], not other components like e.g. [0,1] or [2,2] etc…

Output exceeds the size limit. Open the full output data in a text editor---------------------------------------------------------------------------
FileExistsError                           Traceback (most recent call last)
File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/jit.py:67, in get_cached_module(module_name, object_names, cache_dir, timeout)
     65 try:
     66     # Create C file with exclusive access
---> 67     open(c_filename, "x")
     68     return None, None

FileExistsError: [Errno 17] File exists: '/Users/Apple/.cache/fenics/libffcx_expressions_1095dd4cc6194e6548d03b64f1baa3314ede6698.c'

During handling of the above exception, another exception occurred:

TimeoutError                              Traceback (most recent call last)
Cell In[10], line 5
      3 print(epsilon_comp)
      4 V_eps = fem.FunctionSpace(domain,("DG",0))
----> 5 strain_expr = fem.Expression(epsilon_comp, V_eps.element.interpolation_points())
      6 strains = fem.Function(V_eps)
      7 strains.interpolate(strain_expr)

File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/function.py:134, in Expression.__init__(self, ufl_expression, X, form_compiler_options, jit_options, dtype)
    130 else:
    131     raise RuntimeError(
    132         f"Unsupported scalar type {dtype} for Expression.")
--> 134 self._ufcx_expression, module, self._code = jit.ffcx_jit(mesh.comm, (ufl_expression, _X),
    135                                                          form_compiler_options=form_compiler_options,
    136                                                          jit_options=jit_options)
    137 self._ufl_expression = ufl_expression
    139 # Prepare coefficients data. For every coefficient in form take
    140 # its C++ object.

File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/jit.py:56, in mpi_jit_decorator..mpi_jit(comm, *args, **kwargs)
     51 @functools.wraps(local_jit)
     52 def mpi_jit(comm, *args, **kwargs):
     53 
     54     # Just call JIT compiler when running in serial
     55     if comm.size == 1:
---> 56         return local_jit(*args, **kwargs)
     58     # Default status (0 == ok, 1 == fail)
     59     status = 0

File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/jit.py:211, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
    208     r = ffcx.codegeneration.jit.compile_coordinate_maps(
    209         [ufl_object], options=p_ffcx, **p_jit)
    210 elif isinstance(ufl_object, tuple) and isinstance(ufl_object[0], ufl.core.expr.Expr):
--> 211     r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)
    212 else:
    213     raise TypeError(type(ufl_object))

File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/jit.py:217, in compile_expressions(expressions, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
    215 if cache_dir is not None:
    216     cache_dir = Path(cache_dir)
--> 217     obj, mod = get_cached_module(module_name, expr_names, cache_dir, timeout)
    218     if obj is not None:
    219         return obj, mod, (None, None)

File ~/miniconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/jit.py:89, in get_cached_module(module_name, object_names, cache_dir, timeout)
     87     logger.info(f"Waiting for {ready_name} to appear.")
     88     time.sleep(1)
...
---> 89 raise TimeoutError(f"""JIT compilation timed out, probably due to a failed previous compile.
     90 Try cleaning cache (e.g. remove {c_filename}) or increase timeout option.""")

TimeoutError: JIT compilation timed out, probably due to a failed previous compile.
        Try cleaning cache (e.g. remove /Users/Apple/.cache/fenics/libffcx_expressions_1095dd4cc6194e6548d03b64f1baa3314ede6698.c) or increase timeout option.
18

Please add the code that can reproduce this error, as otherwise it is guesswork from my end to debug what you have coded up.

19

It’s probably a cache problem of my notebook. When I run the code at uni’s computer it seems to work fine.