Equivalent for Expression in dolfinx

dolfinx
21

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")
22

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:

images-20

23

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

24

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)
2 Likes
25

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:

image
image950×569 95.3 KB

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.

1 Like