Error in solve command while analysing a cantilever beam under prescribed displacement

General
1

Hello,
I am trying to model a cantilever beam with prescribed displacement at its end. Here is my 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>]]

Also, how can I find the reaction force at the loading point?

2

This should be

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

i.e. do not wrap the list inside a list for the boundary conditions.

1 Like
3

Thank you,
Any pointers on finding the reaction force?

4

There are many posts regarding this on the forum, including

and

1 Like