from contextlib import ExitStack
import numpy as np
import ufl
from dolfinx import la
from dolfinx.fem import (Expression, Function, FunctionSpace,
VectorFunctionSpace, dirichletbc, form,
locate_dofs_topological, locate_dofs_geometrical)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
set_bc)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_box,
locate_entities_boundary)
from ufl import dx, grad, inner
from mpi4py import MPI
from petsc4py import PETSc
dtype = PETSc.ScalarType
from dolfinx.io import gmshio
import meshio
from petsc4py import PETSc
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
def build_nullspace(V):
index_map = V.dofmap.index_map
bs = V.dofmap.index_map_bs
ns = [la.create_petsc_vector(index_map, bs) for i in range(6)]
with ExitStack() as stack:
vec_local = [stack.enter_context(x.localForm()) for x in ns]
basis = [np.asarray(x) for x in vec_local]
dofs = [V.sub(i).dofmap.list.array for i in range(3)]
for i in range(3):
basis[i][dofs[i]] = 1.0
x = V.tabulate_dof_coordinates()
dofs_block = V.dofmap.list.array
x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
basis[3][dofs[0]] = -x1
basis[3][dofs[1]] = x0
basis[4][dofs[0]] = x2
basis[4][dofs[2]] = -x0
basis[5][dofs[2]] = x1
basis[5][dofs[1]] = -x2
la.orthonormalize(ns)
assert la.is_orthonormal(ns)
return PETSc.NullSpace().create(vectors=ns)
slabdomaingmshio, cell_markers, facet_markers = gmshio.read_from_msh("mesh3D.msh", MPI.COMM_WORLD, gdim=3)
slabmeshio = meshio.read("mesh3D.msh")
tetra_mesh = create_mesh(slabmeshio, "tetra", prune_z=False)
meshio.write("mesh.xdmf", tetra_mesh)
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
slabdomaingmshio = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(slabdomaingmshio, name="Grid")
Q = FunctionSpace(slabdomaingmshio, ("DG", 0))
Emod = Function(Q)
concrete = ct.find(11)
Emod.x.array[concrete] = np.full_like(concrete, 4.3e6, dtype=PETSc.ScalarType)
steel = ct.find(12)
Emod.x.array[steel] = np.full_like(steel, 30e6, dtype=PETSc.ScalarType)
ω, ρ = 300.0, 10.0
x = ufl.SpatialCoordinate(slabdomaingmshio)
f = ufl.as_vector((ρ * ω**2 * x[0], ρ * ω**2 * x[1], 0.0))
ν = 0.3
μ = Emod / (2.0 * (1.0 + ν))
λ = Emod * ν / ((1.0 + ν) * (1.0 - 2.0 * ν))
def σ(v):
"""Return an expression for the stress σ given a displacement field"""
return 2.0 * μ * ufl.sym(grad(v)) + λ * ufl.tr(ufl.sym(grad(v))) * ufl.Identity(len(v))
V = VectorFunctionSpace(slabdomaingmshio, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = form(inner(σ(u), grad(v)) * dx)
L = form(inner(f, v) * dx)
dofs_L = locate_dofs_geometrical(V, lambda x: np.isclose(x[2], -72))
u_L = Function(Q)
bc_L = dirichletbc(PETSc.ScalarType(0), dofs_L,Q)
dofs_R = locate_dofs_geometrical(V, lambda x: np.isclose(x[2], 72))
u_R = Function(Q)
bc_R = dirichletbc(PETSc.ScalarType(0), dofs_R,Q)
bcs = [bc_L,bc_R]
A = assemble_matrix(a, bcs=[bcs])
A.assemble()
b = assemble_vector(L)
apply_lifting(b, [a], bcs=[[bcs]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bcs])
null_space = build_nullspace(V)
A.setNearNullSpace(null_space)
opts = PETSc.Options()
opts["ksp_type"] = "cg"
opts["ksp_rtol"] = 1.0e-10
opts["pc_type"] = "gamg"
opts["mg_levels_ksp_type"] = "chebyshev"
opts["mg_levels_pc_type"] = "jacobi"
opts["mg_levels_esteig_ksp_type"] = "cg"
opts["mg_levels_ksp_chebyshev_esteig_steps"] = 20
solver = PETSc.KSP().create(slabmeshio.comm)
solver.setFromOptions()
solver.setOperators(A)
uh = Function(V)
solver.setMonitor(lambda _, its, rnorm: print(f"Iteration: {its}, rel. residual: {rnorm}"))
solver.solve(b, uh.vector)
solver.view()
uh.x.scatter_forward()
sigma_dev = σ(uh) - (1 / 3) * ufl.tr(σ(uh)) * ufl.Identity(len(uh))
sigma_vm = ufl.sqrt((3 / 2) * inner(sigma_dev, sigma_dev))
W = FunctionSpace(slabmeshio, ("Discontinuous Lagrange", 0))
sigma_vm_expr = Expression(sigma_vm, W.element.interpolation_points())
sigma_vm_h = Function(W)
sigma_vm_h.interpolate(sigma_vm_expr)
with XDMFFile(slabmeshio.comm, "out_elasticity/displacements.xdmf", "w") as file:
file.write_mesh(slabmeshio)
file.write_function(uh)
with XDMFFile(slabmeshio.comm, "out_elasticity/von_mises_stress.xdmf", "w") as file:
file.write_mesh(slabmeshio)
file.write_function(sigma_vm_h)
unorm = uh.x.norm()
if slabmeshio.comm.rank == 0:
print("Solution vector norm:", unorm)
print('done')
I have a traceback at line 123:::
root@80fe948106d6:~/dolfinx-r1# cd /root/dolfinx-r1 ; /usr/bin/env /bin/python3 /root/.vscode-server/extensions/ms-python.python-2023.16.0/pythonFiles/lib/python/debugpy/adapter/../../debugpy/launcher 60645 -- /root/dolfinx-r1/gmsh_importmsh.py
Info : Reading 'mesh3D.msh'...
Info : 42 entities
Info : 6075 nodes
Info : 19690 elements
Info : Done reading 'mesh3D.msh'
Traceback (most recent call last):
File "/root/dolfinx-r1/gmsh_importmsh.py", line 123, in <module>
A = assemble_matrix(a, bcs=[bcs])
File "/usr/lib/python3.10/functools.py", line 889, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py", line 352, in _assemble_matrix_form
_assemble_matrix_mat(A, a, bcs, diagonal, constants, coeffs)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py", line 365, 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 0x7f781a007470>, <dolfinx.fem.forms.Form object at 0x7f78118321b0>, array([], dtype=float64), {(<IntegralType.cell: 0>, -1): array([[ 4300000.],
[ 4300000.],
[ 4300000.],
...,
[30000000.],
[30000000.],
[30000000.]])}, [[<dolfinx.fem.bcs.DirichletBC object at 0x7f7811832b10>, <dolfinx.fem.bcs.DirichletBC object at 0x7f7811832ca0>]]
What I am trying to do is replicate as much as possible demo_elasticty.py only modifying the parts that I need for a slab. I have a best attemp at the creation of Emod for the two materials concrete and steel. Why am I getting this traceback incompatible function arguments. The following argument types are supported…?
P.S. I will also have to create a Function for the poissons ratio where it is different for the two materials that I am trying to make work… The dolfix environment is a docker container dolfinx v0.6.0 as far as I know so far… A = assemble_matrix(a, bcs=[bcs]) <<-- line 123