I am trying to port the NewFrac linear elasticity code to a modern version of Fenics. I get this error:
Traceback (most recent call last):
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/bcs.py", line 230, in dirichletbc
bc = bctype(_value, dofs, V)
TypeError: __init__(): incompatible function arguments. The following argument types are supported:
1. __init__(self, g: ndarray[dtype=float64, order='C', writable=False], dofs: ndarray[dtype=int32, shape=(*), order='C', writable=False], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
2. __init__(self, g: dolfinx.cpp.fem.Constant_float64, dofs: ndarray[dtype=int32, shape=(*), order='C', writable=False], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
3. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
4. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: collections.abc.Sequence[ndarray[dtype=int32, shape=(*), order='C', writable=False]], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
Invoked with types: dolfinx.cpp.fem.DirichletBC_float64, dolfinx.cpp.fem.Constant_float64, ndarray, dolfinx.fem.function.FunctionSpace
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/home/a/fenics/./test.py", line 180, in <module>
bc_bottom_y = dolfinx.fem.dirichletbc(zero_scalar, dofs_bottom_y, V)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/bcs.py", line 232, in dirichletbc
bc = bctype(_value, dofs, V._cpp_object)
TypeError: __init__(): incompatible function arguments. The following argument types are supported:
1. __init__(self, g: ndarray[dtype=float64, order='C', writable=False], dofs: ndarray[dtype=int32, shape=(*), order='C', writable=False], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
2. __init__(self, g: dolfinx.cpp.fem.Constant_float64, dofs: ndarray[dtype=int32, shape=(*), order='C', writable=False], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
3. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
4. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: collections.abc.Sequence[ndarray[dtype=int32, shape=(*), order='C', writable=False]], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
Invoked with types: dolfinx.cpp.fem.DirichletBC_float64, dolfinx.cpp.fem.Constant_float64, ndarray, dolfinx.cpp.fem.FunctionSpace_float64
I can’t find any documentation about this error. Here is the code I have (I cannot provide a minimal example, as I only run into the problem with this full code):
#!/usr/bin/env python3
import sys
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import dolfinx
import ufl
from mpi4py import MPI
from petsc4py import PETSc
import gmsh
from dolfinx.io.gmshio import extract_geometry,extract_topology_and_markers, ufl_mesh
from dolfinx.cpp.io import perm_gmsh
from dolfinx.cpp.mesh import to_type
from dolfinx.mesh import create_mesh
import basix
from dolfinx.fem import Constant
plt.rcParams["figure.figsize"] = (40,6)
def project(v, target_func, bcs=[]):
# Ensure we have a mesh and attach to measure
V = target_func.function_space
dx = ufl.dx(V.mesh)
# Define variational problem for projection
w = ufl.TestFunction(V)
Pv = ufl.TrialFunction(V)
a = ufl.inner(Pv, w) * dx
L = ufl.inner(v, w) * dx
# Assemble linear system
A = assemble_matrix(a, bcs)
A.assemble()
b = assemble_vector(L)
apply_lifting(b, [a], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)
solver = PETSc.KSP().create(A.getComm())
solver.setOperators(A)
solver.solve(b, target_func.vector)
def generate_mesh_with_crack(Lx=1,Ly=1,Lcrack=.3,lc=.015,refinement_ratio=10,dist_min=.05,dist_max=.2):
# For further documentation see
# - gmsh tutorials, e.g. see https://gitlab.onelab.info/gmsh/gmsh/-/blob/master/tutorial/python/t10.py
# - dolfinx-gmsh interface https://github.com/FEniCS/dolfinx/blob/master/python/demo/gmsh/demo_gmsh.py
#
gmsh.initialize()
gdim = 2
proc = MPI.COMM_WORLD.rank
if proc == 0:
model = gmsh.model()
model.add("Rectangle")
model.setCurrent("Rectangle")
p1 = model.geo.addPoint(0.0, 0.0, 0, lc)
p2 = model.geo.addPoint(Lcrack, 0.0, 0, lc)
p3 = model.geo.addPoint(Lx, 0, 0, lc)
p4 = model.geo.addPoint(Lx, Ly, 0, lc)
p5 = model.geo.addPoint(0, Ly, 0, lc)
l1 = model.geo.addLine(p1, p2)
l2 = model.geo.addLine(p2, p3)
l3 = model.geo.addLine(p3, p4)
l4 = model.geo.addLine(p4, p5)
l5 = model.geo.addLine(p5, p1)
cloop1 = model.geo.addCurveLoop([l1, l2, l3, l4, l5])
surface_1 = model.geo.addPlaneSurface([cloop1])
model.mesh.field.add("Distance", 1)
model.mesh.field.setNumbers(1, "NodesList", [p2])
#model.mesh.field.setNumber(1, "NNodesByEdge", 100)
#model.mesh.field.setNumbers(1, "EdgesList", [2])
#
# SizeMax - /------------------
# /
# /
# /
# SizeMin -o----------------/
# | | |
# Point DistMin DistMax
model.mesh.field.add("Threshold", 2)
model.mesh.field.setNumber(2, "IField", 1)
model.mesh.field.setNumber(2, "LcMin", lc / refinement_ratio)
model.mesh.field.setNumber(2, "LcMax", lc)
model.mesh.field.setNumber(2, "DistMin", dist_min)
model.mesh.field.setNumber(2, "DistMax", dist_max)
model.mesh.field.setAsBackgroundMesh(2)
model.geo.synchronize()
surface_entities = [model[1] for model in model.getEntities(2)]
model.addPhysicalGroup(2, surface_entities, tag=5)
model.setPhysicalName(2, 2, "Rectangle surface")
model.mesh.generate(gdim)
# Sort mesh nodes according to their index in gmsh
geometry_data = extract_geometry(model, name="Rectangle")[:,0:2]
topology_data = extract_topology_and_markers(model, "Rectangle")
# Broadcast cell type data and geometric dimension
gmsh_cell_id = MPI.COMM_WORLD.bcast(model.mesh.getElementType("triangle", 1), root=0)
# Extract the cell type and number of nodes per cell and broadcast
# it to the other processors
gmsh_cell_type = list(topology_data.keys())[0]
properties = gmsh.model.mesh.getElementProperties(gmsh_cell_type)
name, dim, order, num_nodes, local_coords, _ = properties
cells = topology_data[gmsh_cell_type]["topology"]
cell_id, num_nodes = MPI.COMM_WORLD.bcast([gmsh_cell_type, num_nodes], root=0)
else:
cell_id, num_nodes = MPI.COMM_WORLD.bcast([None, None], root=0)
cells, geometry_data = np.empty([0, num_nodes]), np.empty([0, gdim])
#Permute topology data from MSH-ordering to dolfinx-ordering
ufl_domain = ufl_mesh(cell_id, gdim, np.float64)
gmsh_cell_perm = perm_gmsh(to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm]
# Create distributed mesh
mesh = create_mesh(MPI.COMM_WORLD, cells, geometry_data[:, :gdim], ufl_domain)
return mesh
Lx = 1.
Ly = .5
Lcrack = 0.3
lc =.2
dist_min = .1
dist_max = .3
mesh = generate_mesh_with_crack(Lcrack=Lcrack,
Lx=Lx,
Ly=Ly,
lc=lc, # characteristic length of the mesh
refinement_ratio=20, # how much to refine at the tip zone
dist_min=dist_min, # radius of tip zone
dist_max=dist_max # radius of the transition zone
)
# element = basix.ufl.element('Lagrange',mesh.topology.cell_name(),degree=1) # ,shape=(2,))
# element = ufl.finiteelement.FiniteElement("Lagrange", mesh.ufl_cell(), degree=1, dim=(2,))
# V = dolfinx.fem.functionspace(mesh, element, (mesh.geometry.dim,))
# 2D Lagrange vector element of degree 1
element = basix.ufl.element(
"Lagrange", # element family
mesh.topology.cell_name(), # cell type as string
1, # polynomial degree
shape=(2,) # vector of 2 components
)
V = dolfinx.fem.functionspace(mesh, element, (mesh.geometry.dim,))
def bottom_no_crack(x):
return np.logical_and(np.isclose(x[1], 0.0),
x[0] > Lcrack)
def right(x):
return np.isclose(x[0], Lx)
V_x = V.sub(0)
V_y = V.sub(1)
dofs_bottom = dolfinx.fem.locate_dofs_geometrical(V, bottom_no_crack)
dofs_right = dolfinx.fem.locate_dofs_geometrical(V, right)
# --- Extract DOFs for each component ---
# Vector element interleaves: [u_x0, u_y0, u_x1, u_y1, ...]
dofs_bottom_x = dofs_bottom[::2] # x-component
dofs_bottom_y = dofs_bottom[1::2] # y-component
dofs_right_x = dofs_right[::2] # x-component
dofs_right_y = dofs_right[1::2] # y-component
# --- Zero BC values ---
zero_scalar = Constant(mesh, 0.0) # scalar, not a 2D vector
bc_bottom_y = dolfinx.fem.dirichletbc(zero_scalar, dofs_bottom_y, V)
bc_right_x = dolfinx.fem.dirichletbc(zero_scalar, dofs_right_x, V)
bcs = [bc_bottom_y, bc_right_x]
dx = ufl.Measure("dx",domain=mesh)
top_facets = dolfinx.mesh.locate_entities_boundary(mesh, 1, lambda x : np.isclose(x[1], Ly))
mt = dolfinx.mesh.MeshTags(mesh, 1, top_facets, 1)
ds = ufl.Measure("ds", subdomain_data=mt)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
E = 1.
nu = 0.3
mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
# this is for plane-stress
lmbda = 2*mu*lmbda/(lmbda+2*mu)
def eps(u):
"""Strain"""
return ufl.sym(ufl.grad(u))
def sigma(eps):
"""Stress"""
return 2.0 * mu * eps + lmbda * ufl.tr(eps) * ufl.Identity(2)
def a(u,v):
"""The bilinear form of the weak formulation"""
k = 1.e+6
return ufl.inner(sigma(eps(u)), eps(v)) * dx
def L(v):
"""The linear form of the weak formulation"""
# Volume force
b = dolfinx.fem.Constant(mesh,ufl.as_vector((0,0)))
# Surface force on the top
f = dolfinx.fem.Constant(mesh,ufl.as_vector((0,0.1)))
return ufl.dot(b, v) * dx + ufl.dot(f, v) * ds(1)
problem = dolfinx.fem.LinearProblem(a(u,v), L(v), bcs=bcs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
uh.name = "displacement"
Path("output").mkdir(parents=True, exist_ok=True)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "elastic.xdmf", "w") as file:
file.write_mesh(uh.function_space.mesh)
file.write_function(uh)
Thanks!
