Hello,
I am trying to translate to Dolfinx this notebook, from the Numerical tours of continuum mechanics using FEniCS.
This is the code I got together so far, which is failing with a ValueError: too many values to unpack (expected 2) error on the last line, just before trying to solve.
import numpy as np
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from mpi4py import MPI
from dolfinx.fem import FunctionSpace, Function
from dolfinx.fem import VectorFunctionSpace
from dolfinx.fem import petsc
from dolfinx import fem, io
from dolfinx import mesh
from dolfinx.mesh import CellType, create_rectangle
from dolfinx.fem.forms import form as _create_form
import ufl
import pyvista
L, b, h = 1, 0.01, 0.03
Nx, Ny, Nz = 51, 5, 5
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, b,h])],
[Nx,Ny,Nz], cell_type=mesh.CellType.hexahedron)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
E , nu = 1e3, 0.
mu = fem.Constant(domain, E/2/(1+nu))
lmbda = fem.Constant(domain, E*nu/(1+nu)/(1-2*nu))
def eps(v):
return ufl.sym(ufl.grad(v))
def sigma(v):
return lmbda * ufl.nabla_div(du) * ufl.Identity(len(du)) + 2 * mu * eps(du)
# Compute the linearized unit preload
N0 = 1.
T = fem.Constant(domain, (-N0, 0, 0))
V = fem.VectorFunctionSpace(domain, ("CG", 2))
v = ufl.TestFunction(V)
du = ufl.TrialFunction(V)
def left(x):
return np.isclose(x[0], 0.)
def right(x):
return np.logical_and( np.logical_and(np.isclose(x[0], L), x[1] < b),x[2] < h)
boundary_facets = mesh.locate_entities_boundary(domain, fdim, left)
u_D = np.array([0,0,0], dtype=ScalarType)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
marked_facets = np.hstack([right_facets])
marked_values = np.hstack([ np.full_like(right_facets, 1)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
boundary_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, right)
boundary_dofs_y = fem.locate_dofs_topological(V.sub(1), domain.topology.dim - 1, boundary_facets)
boundary_dofs_z = fem.locate_dofs_topological(V.sub(2), domain.topology.dim - 1, boundary_facets)
zero_bc = fem.Constant(domain, ScalarType(0))
bcy = fem.dirichletbc(zero_bc, boundary_dofs_y, V.sub(1))
bcz = fem.dirichletbc(zero_bc, boundary_dofs_z, V.sub(2))
bcs = [bc, bcy,bcz]
# Linear Elasticity Bilinear Form
a = ufl.inner(sigma(du), eps(v))*ufl.dx
# External loading on the right end
l = ufl.dot(T, v)*ds(1)
A = fem.petsc.assemble_matrix(fem.form(a), bcs=bcs)
#################
########## NOTE 1
#### The original code goes
##########
#K = PETScMatrix()
#f = PETScVector()
#assemble_system(a, l, bcs, A_tensor=K, b_tensor=f)
#u = Function(V,name = "Displacement")
#solve(K, u.vector(), f, "mumps")
#################
# not so sure what it does really, as all is needed is the field solution u I thought to just solve the linear problem
#################
problem = fem.petsc.LinearProblem(a, l, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
u = problem.solve()
##########
########## NOTE 2
##########
#The original code goes
#kgform = -inner(sigma(u), grad(du).T*grad(v))*dx
#KG = PETScMatrix()
#assemble(kgform, KG)
# Zero out the rows to enforce boundary condition
#for bc in bcs:
# bc.zero(KG)
# not sure at all I understand any of it, especially the assemble(kgform, KG) command
#I attempt below my translation to the best of y understanding
#
##########
##########
kgform = -ufl.inner(sigma(u), ufl.grad(du).T*ufl.grad(v))*ufl.dx
KG = fem.petsc.assemble_matrix(fem.form(kgform), bcs=bcs)
getting the error
ValueError: too many values to unpack (expected 2)
Would be grateful for any hint please, thanks