Dirichlet boundary condition on a point and the consistent reaction

General
1

Hello,
I am trying to solve a simple elasticity problem where a cantilever beam is subjected to a constant displacement at one end. I am having trouble with defining the constant displacement at a particular node. I did find an old question for legacy FEniCS, but that does not work now.

Here is the code

l= 1
h = 0.1
E=130e9
nu=0.22
E=E/(1-nu*nu)
lambda_=nu*E/(1+nu)/(1-2*nu)
mu=E/2/(1+nu)
w= 1e6
I= h**3/12

import numpy as np
import ufl

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

from dolfinx import mesh, fem, plot, io
model='plane_strain'
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0,0]), np.array([l, h])],
                  [200,20], cell_type=mesh.CellType.quadrilateral)
V = fem.VectorFunctionSpace(domain, ("CG", 1))

def clamped_boundary(x):
    return np.isclose(x[0], 0)
def load_point(x):
    return (np.isclose(x[0],l) and np.isclose(x[1],h/2))


fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
load_point_= mesh.locate_entities_boundary(domain, fdim-1, load_point)

u_D1 = np.array([0,0], dtype=ScalarType)
bc1 = fem.dirichletbc(u_D1, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

u_D2= np.array([0,-1e-2], dtype=ScalarType)
bc2= fem.dirichletbc(u_D2, fem.locate_dofs_topological(V, fdim-1, load_point_), V)

bc=[bc1,bc2]

Here is the error message

Traceback (most recent call last):
  File "/Users/aaquib/Desktop/fenics/test2.py", line 32, in <module>
    load_point_= mesh.locate_entities_boundary(domain, fdim-1, load_point)
  File "/Users/aaquib/opt/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/mesh.py", line 110, in locate_entities_boundary
    return _cpp.mesh.locate_entities_boundary(mesh, dim, marker)
  File "/Users/aaquib/Desktop/fenics/test2.py", line 27, in load_point
    return (np.isclose(x[0],l) and np.isclose(x[1],h/2))
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

I am also wondering how would I find the reaction force at that loading point.
Please help!

2

This is a numpy related error message.
The function load_point(x) takes in a numpy array, and expects a numpy array in return.
numpy does not support array_a and array_b as an operation, and thus you either have to use:
bitwise_and numpy.bitwise_and — NumPy v1.26 Manual

def load_point(x):
    return np.isclose(x[0],l) & np.isclose(x[1],h/2)

logical_and: numpy.logical_and — NumPy v1.26 Manual

def load_point(x):
    return np.logical_and(np.isclose(x[0],l), np.isclose(x[1],h/2))
1 Like
3

Thank you, this resolves the issue in the BCs, but I am now getting an error with the solve command
Here is the code

l= 1
h = 0.1
E=130e9
nu=0.22
E=E/(1-nu*nu)
lambda_=nu*E/(1+nu)/(1-2*nu)
mu=E/2/(1+nu)
w= 1e6
I= h**3/12

import numpy as np
import ufl

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

from dolfinx import mesh, fem, plot, io
model='plane_strain'
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0,0]), np.array([l, h])],
                  [200,20], cell_type=mesh.CellType.quadrilateral)
V = fem.VectorFunctionSpace(domain, ("CG", 1))

def clamped_boundary(x):
    return np.isclose(x[0], 0)
def load_point(x):
    return np.logical_and(np.isclose(x[0],l), np.isclose(x[1],h/2))

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
load_point_= mesh.locate_entities_boundary(domain, fdim-1, load_point)

u_D1 = np.array([0,0], dtype=ScalarType)
bc1 = fem.dirichletbc(u_D1, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

u_D2= np.array([0,-1e-2], dtype=ScalarType)
bc2= fem.dirichletbc(u_D2, fem.locate_dofs_topological(V, fdim-1, load_point_), V)

bc=[bc1,bc2]

#T = fem.Constant(domain, ScalarType((0, 0)))
#ds = ufl.Measure("ds", domain=domain)

if model=='plane_stress':
    lambda_=2*mu*lambda_/(lambda_+2*mu)

def epsilon(u):
    return ufl.sym(ufl.grad(u)) 
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)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx 

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

and here is the error message

Traceback (most recent call last):
  File "/Users/aaquib/Desktop/fenics/test2.py", line 59, in <module>
    uh = problem.solve()
  File "/Users/aaquib/opt/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 604, in solve
    _assemble_matrix_mat(self._A, self._a, bcs=self.bcs)
  File "/Users/aaquib/opt/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 362, in _assemble_matrix_mat
    _cpp.fem.petsc.assemble_matrix(A, a, constants, coeffs, bcs)
TypeError: assemble_matrix(): incompatible function arguments. The following argument types are supported:
    1. (A: mat, a: dolfinx::fem::Form<double>, constants: numpy.ndarray[numpy.float64], coeffs: Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float64]], bcs: List[dolfinx::fem::DirichletBC<double>], unrolled: bool = False) -> None
    2. (A: mat, a: dolfinx::fem::Form<double>, constants: numpy.ndarray[numpy.float64], coeffs: Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float64]], rows0: numpy.ndarray[numpy.int8], rows1: numpy.ndarray[numpy.int8], unrolled: bool = False) -> None

Invoked with: <petsc4py.PETSc.Mat object at 0x1682c8680>, <dolfinx.fem.forms.Form object at 0x10f0b1df0>, array([], dtype=float64), {(<IntegralType.cell: 0>, -1): array([], shape=(0, 0), dtype=float64)}, [[<dolfinx.fem.bcs.DirichletBC object at 0x168293880>, <dolfinx.fem.bcs.DirichletBC object at 0x10f0b2f70>]]