Issue with assemble_matrix & quadrature

dolfinx
1

All,
I’m new of fenicsx and after the basics, I’m trying to replicate this well known tutorial on plasticity: Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation .

I’ve slightly modified the test case to use a bar subject to traction.

I currently use dolfinx 0.4.2.0 and I get stuck when I try convert “assemble_ system” function.

I used assemble_matrix, and assemble_vector to create my bilinear and linear forms but I get the following error.

 File "error.py", line 65, in <module>
    A = fem.assemble_matrix(a_Newton)
  File "spack/var/spack/environments/fenics-real/.spack-env/._view/cccxafi6cfibntvlu62rgazanffzvvxu/lib/python3.9/functools.py", line 888, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "spack/var/spack/environments/fenics-real/.spack-env/view/lib/python3.9/site-packages/dolfinx/fem/assemble.py", line 229, in assemble_matrix
    A = create_matrix(a)
  File "spack/var/spack/environments/fenics-real/.spack-env/view/lib/python3.9/site-packages/dolfinx/fem/assemble.py", line 98, in create_matrix
    sp = dolfinx.fem.create_sparsity_pattern(a)
  File "/gpfs/rr/rrfea/rruk/home/acino/myApps/spack/var/spack/environments/fenics-real/.spack-env/view/lib/python3.9/site-packages/dolfinx/fem/__init__.py", line 34, in create_sparsity_pattern
    topology = a.mesh.topology
AttributeError: 'Form' object has no attribute 'mesh'

Here are the relevant part of the code up to when the error occur

import numpy as np
import os
import ufl
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import mesh, fem, plot, io


L = 100; W = 10
E =210000 
nu = 0.3
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
Et = E/100.  # tangent modulus
H = E*Et/(E-Et)


comm = MPI.COMM_WORLD
domain = mesh.create_box(comm, [np.array([0,0,0]), np.array([L, W, W])],
                  [2,2,4], cell_type=mesh.CellType.tetrahedron )

deg_u = 2
deg_stress = 2

V = fem.VectorFunctionSpace(domain, ("CG", deg_u ))

v = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)


We = ufl.VectorElement("Quadrature",domain.ufl_cell(),degree=deg_stress,dim=6,quad_scheme='default')
W = fem.FunctionSpace(domain, We)
n_elas = fem.Function(W)

W0e = ufl.FiniteElement("Quadrature",domain.ufl_cell(),degree=deg_stress,quad_scheme='default')
W0 = fem.FunctionSpace(domain, W0e)
beta = fem.Function(W0)

u = fem.Function(V, name="Total displacement")


def eps(v):
    e = ufl.sym(ufl.grad(v))
    return ufl.as_tensor([[e[0, 0], e[0, 1], e[0, 2]],
                          [e[0, 1], e[1, 1], e[0, 2]],
                          [e[0, 2], e[1, 2], e[2, 2]]])
def sigma(eps_el):
    return lmbda * ufl.tr(eps_el) * ufl.Identity(u.geometric_dimension()) + 2*mu*eps_el


def as_3D_tensor(X):
    return ufl.as_tensor([[X[0], X[3], X[4]],
                      [X[3], X[1], X[5]],
                      [X[4], X[5], X[2]]])

def sigma_tang(e):
    N_elas = as_3D_tensor(n_elas)
    return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*ufl.inner(N_elas, e)*N_elas-2*mu*beta*ufl.dev(e)

metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = ufl.dx(metadata=metadata)
a_Newton = ufl.inner(eps(v), sigma_tang(eps(u_)))*dxm

A = fem.assemble_matrix(a_Newton)
2

This should be changed to:

A = fem.petsc.assemble_matrix(fem.form(a_Newton))
1 Like