Equivalent for Expression in dolfinx

[Thomas996]

Here is the code that reproduce the result:

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 = 10.0
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0],[L,1,1]],[10,5,5], mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 1))

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/4

    def eval(self, x):
        return (np.zeros_like(x[0]), # equal to x[0] - x[0],
                y0 + (x[1] - y0) * np.cos(self.theta) - (x[2] - z0) * np.sin(self.theta) - x[1],
                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)

right_dofs = fem.locate_dofs_topological(V,facet_tag.dim, facet_tag.find(2))
bc_r = fem.dirichletbc(u_D,right_dofs)
bcs =  [bc_l, bc_r]   

from petsc4py.PETSc import ScalarType

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

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(len(u)) + 2*mu*epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, ScalarType((0, 0, 0)))
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=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
#print(problem)
uh = problem.solve()

# visualization 
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))
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:
   pyvista.start_xvfb()
   figure_as_array = p.screenshot("deflection.png")

[Thomas996]

So basically what I am trying to achieve is a torsion of the whole beam in a continous and homogeneous way. (just like the body force to the beam). This boundary coundition method seems to be not like that. The change of the angle at the end is very abruptly and only happen at the last surface of the mesh. I wonder if there is way to do it?
Following is an image of such torsion:

[lclaire]

Hello, I just came across this thread, as I am trying to impose a not spatially dependent Dirichlet boundary condition that is varying with time. Here is an example of my code:

import dolfinx
from mpi4py import MPI
import numpy as np
from petsc4py.PETSc import ScalarType

L = 1. # total length
d = L/10. # thickness
h = d/6. # size of a cell

my_domain = dolfinx.mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, -0.5*d), (L, 0.5*d)), n=(int(L/h), int(d/h)),
                            cell_type=dolfinx.mesh.CellType.triangle)

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

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

V = dolfinx.fem.VectorFunctionSpace(my_domain, ("CG", 2))

left_dofs = dolfinx.fem.locate_dofs_geometrical(V, left)
right_dofs = dolfinx.fem.locate_dofs_geometrical(V, right)

u_right = dolfinx.fem.Function(V)
imposed_disp = dolfinx.fem.Constant(my_domain, ScalarType((0.0, 0.0)))
u_right.interpolate(imposed_disp)

It results in an infinite loop with the error message: “Couldn’t map ‘c_0’ to a float, returning ufl object without evaluation.” I am not sure how I nee to change the syntax…

Many thanks in advance for you help!
Claire

[dokken]

@lclaire

imposed_disp = dolfinx.fem.Constant(my_domain, ScalarType((0.0, 0.0)))
u_right.interpolate(imposed_disp)

You can either wrap it as

u_right = dolfinx.fem.Function(V)
imposed_disp = dolfinx.fem.Expression(dolfinx.fem.Constant(my_domain, ScalarType((0.0, 0.0))), V.element.interpolation_points())
u_right.interpolate(imposed_disp)

or do something along the lines of:

class TimeDepBC():
    def __init__(self, t=0):
        self.t = t
        self.values = None
    def __call__(self, x):
        if self.values is None:
            self.values = np.zeros((2, x.shape[1]), dtype=ScalarType)
        # Update values[0] and values[1] (the x and y component depending on the self.t parameter)
        # ...
        # For instance:
        self.values[0] = self.t
        return self.values


bc_expr = TimeDepBC()

u_right = dolfinx.fem.Function(V)
u_right.interpolate(bc_expr.__call__)
print(u_right.x.array)
bc_expr.t = 5
u_right.interpolate(bc_expr.__call__)
print(u_right.x.array)

[jkrokowski]

I would start by thinking a bit harder about what you are trying simulate with this script. This script has a “toy” material model that was implemented to visually inspect results, per page 55 of The legacy FEniCS tutorial. While it does model the equations of linear elasticity, it turns out that material behavior drives much of the actual output. For example, it’s plainly obvious that polyurethane and aluminum have extremely different responses to loads.

Consider these thoughts:

• You have applied a relatively large rotation to a relatively flexible object. What do you think would happen? (I would expect local failure, which is something you are seeing in your example) Try lowering that rotation angle to ~1-2 deg and see what you result looks like
• Additionally, try increasing the non-dimensional “elasticity” parameter \beta and see what happens (you should see something more like you expected)

For example, with \beta=50 and \theta=\frac{\pi}{4}, I get the following plot:

As a side note if you are interested in 1D torsional models and predicting torsional behaviour of slender objects, you may find this post helpful.

[yuri]

can someone please update the answer as for the latest version please and tag the new answer (eventually) as the SOLUTION?
please make the new solution using the following example:
r0 = Expression(‘Ru*pow(Rd/Ru, x[0]/L)’,
degree=2, Ru=Ru, Rd=Rd, L=L)
then I would like to access r0(x) or make arithmetic operation with r0 like y = r0x2+4… (for example)
thank you

[jkrokowski]

Hi @yuri. I’m not sure what you mean by “update the answer and tag as the solution”. If you read the thread above, there are multiple approaches for modeling torsion in a prismatic structure. I’m not sure which version you are referring to, but here are some ideas to help you get your code up and running:

• If you are looking at old code snippets and wanting to update to the latest version of dolfinx ( v0.8.0 last I checked), I would suggest starting with the changelog on the FEniCSx tutorial

• The user expression function Expression() no long exists in dolfinx (this is a legacy fenics syntax). The current syntax is to interpolate an expression. Take a look at this example

[LGerard]

@dokken

print(u_right.x.array)

Hi @dokken, I have been trying to implemented different materials with class, but I have the next error.

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

I add .any() or .all(), but the condition is not considered

from petsc4py import PETSc
from dolfinx.fem import Function, functionspace
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
from dolfinx.io import VTXWriter
import numpy as np

mesh = create_unit_square(MPI.COMM_WORLD, 20, 20)
Q = functionspace(mesh, ("Lagrange", 1))
H = 0.5
mu = 0.5
sigma = 0.1
class Porosity:
    def __call__(self, x):
        values = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
        phi0 = 0.6*np.exp(-1. / 2 * ((x[0] - mu) / sigma) ** 2)
        print(H)
        if x[1] <= H:
            values[0] = (1.0 - phi0) * np.sqrt(x[1] / H) + phi0
        else:
            values[0] = 2.0
        return values
phi = Function(Q)
phi.name = "phi"
porosity = Porosity()
phi.interpolate(porosity)
with VTXWriter(mesh.comm, "phi_test.bp", [phi], engine="BP4") as bp:
    bp.write(0.0)

[dokken]

This is because you get an array of size (3, num_points) as input to Porosity.__call__.
There are two ways you can do this;

1. Inefficient: Loop through all points and add to values[0,i] for each of them.
2. Efficient: Use vectorized evaluation, as shown below:

from dolfinx.fem import Function, functionspace
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
from dolfinx.io import VTXWriter
import numpy as np

mesh = create_unit_square(MPI.COMM_WORLD, 20, 20)
Q = functionspace(mesh, ("Lagrange", 1))
H = 0.5
mu = 0.5
sigma = 0.1


class Porosity:
    def __call__(self, x):
        phi0 = 0.6 * np.exp(-1.0 / 2 * ((x[0] - mu) / sigma) ** 2)
        less_marker = (x[1] <= H).astype(np.float64)
        less_val = (1.0 - phi0) * np.sqrt(x[1] / H + 1e-14) + phi0
        larger_val = np.full(x.shape[1], 2.0)
        return less_marker * less_val + (1 - less_marker) * larger_val


phi = Function(Q)
phi.name = "phi"
porosity = Porosity()
phi.interpolate(porosity)
with VTXWriter(mesh.comm, "phi_test.bp", [phi], engine="BP4") as bp:
    bp.write(0.0)

yielding

[emna_frikha]

EDITED BY FB due to incorrectly formatted code.

[francesco-ballarin]

@emna_frikha please open a new post, making sure you correctly format your code and read Read before posting: How do I get my question answered?