Uniaxial tension boundary conditions

dolfinx
1

I’m trying to model a simple uniaxial tension test by setting up a block and applying boundary conditions Ux(0) = 0 and Ux(L) = 0.1, leaving the other directions Uy and Uz free. However, I’m struggling to apply these BCs

import gmsh
import numpy as np
import ufl
import dolfinx
from dolfinx import fem, io, mesh, plot, nls, log
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import meshio
import dolfinx.io.gmshio as gmshio

from meshing import generate_block
from matplotlib import pyplot as plt

L = 10
W = 1
H = 1
lambda_ = 1
mu = 1
rho = 0.001
g = 0.001
seed =4

generate_block(L, W, H, seed)

gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 3
fdim = gdim - 1
msh, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)

# Find facets(boundaries)
facets_bottom = facet_markers.indices[facet_markers.values == 1]
facets_top = facet_markers.indices[facet_markers.values == 2]
facets_left = facet_markers.indices[facet_markers.values == 3]
facets_right = facet_markers.indices[facet_markers.values == 4]
facets_front = facet_markers.indices[facet_markers.values == 5]
facets_top = facet_markers.indices[facet_markers.values == 6]

celltype = ufl.Cell("tetrahedron", 3)
u_el = ufl.VectorElement("CG", celltype, 1)
V = fem.FunctionSpace(msh, u_el)

Q = V.sub(0).collapse()[0]
facet_left_dofs = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_left)
bc_left = fem.dirichletbc(fem.Constant(msh, 0.0),facet_left_dofs)

facet_left_dofs = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_right)
bc_right = fem.dirichletbc(fem.Constant(msh, 0.1),facet_right)

bcs = [bc_left,bc_right]

T = fem.Constant(msh, ScalarType((0, 0, 0.)))  #applied tractions
f = fem.Constant(msh, ScalarType((0, 0, 0))) # body forces

However, this throws the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In [22], line 51
     49 Q = V.sub(0).collapse()[0]
     50 facet_left_dofs = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_left)
---> 51 bc_left = fem.dirichletbc(fem.Constant(msh, 0.0),facet_left_dofs)
     53 facet_left_dofs = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_right)
     54 bc_right = fem.dirichletbc(fem.Constant(msh, 0.1),facet_right)

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/bcs.py:183, in dirichletbc(value, dofs, V)
    180     raise AttributeError("Boundary condition value must have a dtype attribute.")
    182 formcls = type("DirichletBC", (DirichletBCMetaClass, bctype), {})
--> 183 return formcls(value, dofs, V)

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/bcs.py:129, in DirichletBCMetaClass.__init__(self, value, dofs, V)
    127         super().__init__(_value, dofs, V._cpp_object)  # type: ignore
    128 else:
--> 129     super().__init__(_value, dofs)

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.DirichletBC_float64(g: numpy.ndarray[numpy.float64], dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    2. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Constant<double>, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: numpy.ndarray[numpy.int32])
    4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: List[numpy.ndarray[numpy.int32][2]], V: dolfinx::fem::FunctionSpace)

Invoked with: <dolfinx.cpp.fem.Constant_float64 object at 0x7f3c5e496f70>, [array([2847, 2880, 2892, 2895, 2949, 2952, 2988, 2991, 3003, 3027, 3030,
       3033, 3039, 3042, 3045, 3054, 3069, 3072, 3081, 3087, 3102, 3108,
       3117, 3120, 3126, 3132, 3135, 3144, 3147, 3150, 3153], dtype=int32), array([ 949,  960,  964,  965,  983,  984,  996,  997, 1001, 1009, 1010,
       1011, 1013, 1014, 1015, 1018, 1023, 1024, 1027, 1029, 1034, 1036,
       1039, 1040, 1042, 1044, 1045, 1048, 1049, 1050, 1051], dtype=int32)]

Any help would be appreciated, thank you!

2

This should be

facet_left_dofs = fem.locate_dofs_topological(V.sub(0), fdim, facets_left)
bc_left = fem.dirichletbc(fem.Constant(msh, 0.0),facet_left_dofs, V.sub(0))

See:
https://jsdokken.com/dolfinx-tutorial/chapter3/component_bc.html?highlight=constant#boundary-conditions

3

Thanks for the quick response, correcting that has gotten past creating the BCs, and running this as a .py from the terminal works fine, however, the linear solver now gets stuck with this error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In [3], line 80
     77 print('starting to solve the problem at time %g s' %(tic))
     79 problem = fem.petsc.LinearProblem(a, L, bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
---> 80 uh = problem.solve()
     82 toc = (time.time_ns()-epoch)*1e-9
     83 print('finished solving the problem at time %g s' %(toc))

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/petsc.py:609, in LinearProblem.solve(self)
    607 # Assemble lhs
    608 self._A.zeroEntries()
--> 609 _assemble_matrix_mat(self._A, self._a, bcs=self.bcs)
    610 self._A.assemble()
    612 # Assemble rhs

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/petsc.py:365, in _assemble_matrix_mat(A, a, bcs, diagonal, constants, coeffs)
    363 constants = _pack_constants(a) if constants is None else constants
    364 coeffs = _pack_coefficients(a) if coeffs is None else coeffs
--> 365 _cpp.fem.petsc.assemble_matrix(A, a, constants, coeffs, bcs)
    366 if a.function_spaces[0] is a.function_spaces[1]:
    367     A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)

TypeError: assemble_matrix(): incompatible function arguments. The following argument types are supported:
    1. (A: mat, a: dolfinx::fem::Form<std::complex<double> >, constants: numpy.ndarray[numpy.complex128], coeffs: Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.complex128]], bcs: List[dolfinx::fem::DirichletBC<std::complex<double> >], unrolled: bool = False) -> None
    2. (A: mat, a: dolfinx::fem::Form<std::complex<double> >, constants: numpy.ndarray[numpy.complex128], coeffs: Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.complex128]], rows0: numpy.ndarray[numpy.int8], rows1: numpy.ndarray[numpy.int8], unrolled: bool = False) -> None

Invoked with: <petsc4py.PETSc.Mat object at 0x7fcf7beb30b0>, <dolfinx.fem.forms.Form object at 0x7fcf7bd0bd80>, array([], dtype=complex128), {(<IntegralType.cell: 0>, -1): array([], shape=(0, 0), dtype=complex128)}, [<dolfinx.fem.bcs.DirichletBC object at 0x7fcfb4ab79c0>, <dolfinx.fem.bcs.DirichletBC object at 0x7fcf7bd0b7e0>]

Did you forget to `#include <pybind11/stl.h>`? Or <pybind11/complex.h>,
<pybind11/functional.h>, <pybind11/chrono.h>, etc. Some automatic
conversions are optional and require extra headers to be included
when compiling your pybind11 module.

…when trying to solve within a jupyter notebook. I think it has something to do with how ufl.inner() and ufl.dot() are performed in notebooks compared to running .py files, specifically here: L = ufl.inner(f, v) * dx + ufl.inner(T, v) * ds(2), is there a solution to this? As I find notebooks easier to debug.
Thanks again dokken, you’re such a great help!

4

You are running dolfinx with complex mode petsc. Then the boundary conditions also has to be complex valued,
Ie

Has to be changed to
bc_left = fem.dirichletbc(fem.Constant(msh, ScalarType(0.0)),facet_left_dofs, V.sub(0))

see
https://jsdokken.com/dolfinx-tutorial/chapter1/complex_mode.html
for this and other things to note when running in complex mode.