Hello,
I’m struggling with the translation of the last part of this code (dolfin) to dolfinx. It calculates the reaction forces via the virtual work principle.
residual = action(a, u) - l
v_reac = Function(V)
bcRx = DirichletBC(V.sub(0), Constant(1.), facets, 1)
bcRx.apply(v_reac.vector())
Rx = assemble(action(residual, v_reac))
Can someone give me an idea where to find the action and apply functions in dolfinx? Is there an easier way to do this?
Thanks in advance!
Action is a ufl function, thus from ufl import action should suffice.
Here is a minimal example of how to use action and an equivalent to apply:
import ufl
import dolfinx
import dolfinx.io
import numpy as np
from mpi4py import MPI
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 10, 10)
V = dolfinx.VectorFunctionSpace(mesh, ("CG", 1))
u = dolfinx.Function(V)
# Define dirichlet bc on subspace
V0 = V.sub(0).collapse()
u_bc = dolfinx.Function(V0)
u_bc.x.array[:] = 1
fdim = mesh.topology.dim - 1
facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, lambda x: np.full(x.shape[1], True))
dofs = dolfinx.fem.locate_dofs_topological((V.sub(0), V0), fdim, facets)
bc = dolfinx.DirichletBC(u_bc, dofs, V.sub(0))
# Apply BC to Function
dolfinx.fem.set_bc(u.vector, [bc])
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "u.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(u)
# Use action on form
u_, v_ = ufl.TrialFunction(V), ufl.TestFunction(V)
a = ufl.inner(u_, v_) * ufl.dx
aV = ufl.action(a, u)
print(dolfinx.fem.assemble_vector(aV).array)
Thank you very much for your answer @dokken! The equivalent to apply works and I was able to use the example for a linear problem. Is it also possible to use the action (or something equivalent) on a nonlinear form?
Here is a (more or less) minimal example of my code. A bar is clamped on one end with a traction T on the other end. I want to calculate the reaction forces on the clamped end (as shown in the link I posted above). This is just an example code to see if the reaction forces are calculated correctly, in my actual code I have a complex geometry and Dirichlet boundary conditions for the displacement instead of a traction force. If there is another way to calculate the reaction forces (other than integrating the stress field over my boundary, this is not exact enough) please give me an advise on how to do so.
Thanks again!
import numpy as np
from dolfinx.cpp.io import perm_gmsh, extract_local_entities
import dolfinx
import ufl
import dolfinx.io
from petsc4py import PETSc
def get_DirichletBC(V, subspace, coordinate, location, disp):
"""
Returns dolfinx.Dirichlet() condition
Input:
str V: Functionspace variable
int subspace: V.sub(0,1,2) spatial direction of function space
int coordinate: spatial direction of boundary
double location: value of boundary location in respective direction
double disp: prescribed displacement for respective boundary
"""
dofs = dolfinx.fem.locate_dofs_geometrical((V.sub(subspace), V.sub(subspace).collapse()),
lambda x: np.isclose(x[coordinate], location))
u_d = dolfinx.Function(V.sub(subspace).collapse())
with u_d.vector.localForm() as loc:
loc.set(disp)
dirichlet = dolfinx.DirichletBC(u_d, dofs, V.sub(subspace))
return dirichlet
# Green-Lagrange strain
def Epsilon():
return 0.5 * (F.T * F - I)
# 2nd Piola-Kirchhoff stress
def Sigma():
return lambda_ * ufl.tr(Epsilon())*I + 2 * mu * Epsilon()
mesh = dolfinx.BoxMesh(MPI.COMM_WORLD,[[0.0,0.0,-15.0], [1, 1, 15.0]], [5, 5, 150], dolfinx.cpp.mesh.CellType.hexahedron)
# create face-tags
boundaries = [(1, lambda x: np.isclose(x[2], -15.0)),
(2, lambda x: np.isclose(x[2], 15.0))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = dolfinx.mesh.locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_tag = dolfinx.MeshTags(mesh, fdim, np.array(np.hstack(facet_indices), dtype=np.int32),
np.array(np.hstack(facet_markers), dtype=np.int32)
# define measure for surface integration
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
# create function space
V = dolfinx.VectorFunctionSpace(mesh, ("CG", 1))
#define material parameters
E = 44
nu = 0.4
lambda_ = nu * E / ((1 - 2 * nu) * (1 + nu))
mu = 0.5 * E / (1 + nu)
# define trial and test function
#u = ufl.TrialFunction(V) # for linear problem
u = dolfinx.Function(V) # for nonlinear problem
v = ufl.TestFunction(V)
# define volume forces f and traction forces T
f = dolfinx.Constant(mesh, (0.0, 0.0, 0.0))
T = dolfinx.Constant(mesh, (0.0, 0.0, 0.01))
#define kinematic quantities used in the problem
#spatial dimension
d = len(u)
#identity tensor
I = ufl.variable(ufl.Identity(d))
#Deformation gradient
F = ufl.variable(I + ufl.grad(u))
#right Cauchy-Green tensor
C = ufl.variable(F.T*F)
#Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
# Define form F (we want to find u such that F(u) = 0)
Form = ufl.inner(ufl.grad(v), Sigma())*ufl.dx - ufl.inner(v, f)*ufl.dx - ufl.inner(v, T)*ds(2)
# define Dirichlet bonudary conditions
clamped_x = get_DirichletBC(V, 0, 2, -15.0, 0.0)
clamped_y = get_DirichletBC(V, 1, 2, -15.0, 0.0)
clamped_z = get_DirichletBC(V, 2, 2, -15.0, 0.0)
top_x = get_DirichletBC(V, 0, 2, 15.0, 0.0)
top_y = get_DirichletBC(V, 1, 2, 15.0, 0.0)
bc = [clamped_x, clamped_y, clamped_z, top_x, top_y]
problem = dolfinx.fem.NonlinearProblem(Form, u, bc)
solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
num_its, converged = solver.solve(u)
u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
u.name = "Displacement"
Hy,
I try to use the solution provided by Dokken.
unfortunately, it does not work.
It seems that there have been changes in fenicsx since July 2021
I corrected a few lines but it still doesn’t work.
Does anyone have the solution.
Thanks
Xavier
Here is my script
import ufl
from dolfinx import mesh
from dolfinx import fem
from dolfinx import io
import numpy as np
from mpi4py import MPIOmega = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = fem.VectorFunctionSpace(Omega, (“CG”, 1))
u = fem.Function(V)Define dirichlet bc on subspace
V0 = V.sub(0).collapse()[0]
u_bc = fem.Function(V0)
u_bc.x.array[:] = 1
fdim = Omega.topology.dim - 1
facets = mesh.locate_entities_boundary(Omega, fdim, lambda x: np.full(x.shape[1], True))
dofs = fem.locate_dofs_topological((V.sub(0), V0), fdim, facets)
bc = fem.dirichletbc(u_bc, dofs, V.sub(0))Apply BC to Function
fem.set_bc(u.vector, [bc])
with io.XDMFFile(MPI.COMM_WORLD, “u.xdmf”, “w”) as xdmf:
xdmf.write_mesh(Omega)
xdmf.write_function(u)Use action on form
u_, v_ = ufl.TrialFunction(V), ufl.TestFunction(V)
a = fem.form(ufl.inner(u_, v_) * ufl.dx)
aV = ufl.action(a, u)
print(fem.assemble_vector(aV).array)
and the error message
Unable to convert object to a UFL form: <dolfinx.fem.forms.Form object at 0x7f3b766cde90>
ERROR:UFL:Unable to convert object to a UFL form: <dolfinx.fem.forms.Form object at 0x7f3b766cde90>
Traceback (most recent call last):
File “/home/fenicsx/shared/essai.py”, line 41, in
aV = ufl.action(a, u)
File “/usr/local/lib/python3.10/dist-packages/ufl/formoperators.py”, line 108, in action
form = as_form(form)
File “/usr/local/lib/python3.10/dist-packages/ufl/form.py”, line 500, in as_form
error(“Unable to convert object to a UFL form: %s” % ufl_err_str(form))
File “/usr/local/lib/python3.10/dist-packages/ufl/log.py”, line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Unable to convert object to a UFL form: <dolfinx.fem.forms.Form object at 0x7f3b766cde90>
You need to use ufl.action prior to wrapping it as a fem.form, i.e.
a = ufl.inner(u_, v_) * ufl.dx
aV = ufl.action(a, u)
L = fem.form(aV)
print(fem.petsc.assemble_vector(L).array)
Than you
It works once I replace the last line by
print(fem.petsc.assemble_vector(L).array)