Dolfinx : Block size not yet supported

Hello, trying to solve linear elasticity equation, but it looks there is a problem with the sparsity pattern (from what could find, but i don’t understand why). Can some one explain my error?
MWE:

import numpy as np
from dolfinx import fem, cpp, io, mesh
from mpi4py import MPI
from petsc4py import PETSc
import ufl

YOUNG=1.
NU=0.4
G=9.81
MASS=3.5
WEIGHT=G*MASS

L, W, H=0.303, 0.19, 0.064
EMP=0.06
AXE=0.025

NX, NY, NZ = 30, 19, 6

points=[np.array([0, 0, 0]), np.array([L, W, H])]
domain=mesh.create_box( 
    MPI.COMM_WORLD, 
    points,
    [NX,NY,NZ],
    cell_type=cpp.mesh.CellType.hexahedron,
    ghost_mode=mesh.GhostMode.shared_facet
    )
dim=domain.topology.dim

def axes(x):
    return np.logical_and(
               np.logical_or(np.sqrt( (x[0]-EMP)**2 + (x[2]-H/2)**2 ) <=AXE,
                             np.sqrt( (x[0]-L+EMP)**2 + (x[2]-H/2)**2 ) <=AXE),
               np.logical_or(np.isclose(x[1], 0, atol=1e-3), 
                             np.isclose(x[1], W, atol=1e-3))
           )

def bottom(x):
    return np.isclose(x[2], 0, atol=1e-3)

bottom_facets=mesh.locate_entities(domain, dim-1, bottom)

facet_indices=[bottom_facets]
facet_indices=np.array(np.hstack(facet_indices), dtype=np.int32)

facet_markers=[np.full(len(bottom_facets), 1)]
facet_markers=np.array(np.hstack(facet_markers), dtype=np.int32)

sorted_facets=np.argsort(facet_indices)
facet_tag=mesh.meshtags(
    domain, 
    dim-1, 
    facet_indices[sorted_facets],
    facet_markers[sorted_facets]
    )

U=fem.VectorFunctionSpace(domain, ("CG", 1))
u=fem.Function(U, name="Déplacement")

LAMBDA=YOUNG*NU/(1+NU)/(1-2*NU)
MU=YOUNG/(2*1+NU)

bcdofs=fem.locate_dofs_geometrical(U, bottom)
u_0=np.array((0,)*dim, dtype=PETSc.ScalarType)
bc=fem.dirichletbc(u_0, bcdofs, U)

ds= ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
BOTTOM_SURF=fem.assemble_scalar(fem.form(fem.Constant(domain, 1.)*ds(1)))

def eps(v):
    return ufl.sym(ufl.grad(v))
def sigma(v):
    return (LAMBDA*ufl.nabla_div(v)*ufl.Identity(dim)+2*MU*eps(v))
def dE(u, v):
    return ufl.inner(sigma(u), eps(v))

u_=ufl.TestFunction(U)
du=ufl.TrialFunction(U)
a=fem.form(ufl.inner(sigma(u_), eps(du))*ufl.dx)
b=fem.form(ufl.dot((WEIGHT/BOTTOM_SURF)*ufl.FacetNormal(domain), u_)*ds(1))

A=fem.assemble_matrix(a, bcs=[bc])
A.assemble()

The error message:

Traceback (most recent call last):
  File "/root/opti copy.py", line 81, in <module>
    A=fem.assemble_matrix(a, bcs=[bc])
  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/assemble.py", line 273, in _assemble_matrix_form
    A: la.MatrixCSRMetaClass = create_matrix(a)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/assemble.py", line 97, in create_matrix
    return la.matrix_csr(sp, dtype=a.dtype)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/la.py", line 81, in matrix_csr
    return matrixcls(sp)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/la.py", line 54, in __init__
    super().__init__(sp)
RuntimeError: Block size not yet supported

Replace this with
A=fem.petsc.assemble_matrix(a, bcs=[bc])

The other assemble matrix is for the built-in sparse matrix in DOLFINx, Which does not support block size (vector spaces) or Dirichlet BCs yet.

2 Likes