Help needed in analyzing a mindlin plate in fenicsx

dolfinx
1

Hello,
I am trying to convert this tutorial (Reissner-Mindlin plate with Quadrilaterals — Numerical tours of continuum mechanics using FEniCS master documentation) which is in fenics to fenicsx, but the final script I wrote is giving me erroneous values.
I feel like I am making a very stupid mistake but am unable to find it. Any help is greatly appreciated.

from mpi4py import MPI
import numpy as np
import dolfinx.cpp
import dolfinx.io
import dolfinx.la
import ufl

def border(x):
    return np.isclose(x[0], 0.0)

def strain2voigt(eps):
    return ufl.as_vector([eps[0, 0], eps[1, 1], 2*eps[0, 1]])
def voigt2stress(S):
    return ufl.as_tensor([[S[0], S[2]], [S[2], S[1]]])
def curv(u):
    (w, theta) = ufl.split(u)
    return ufl.sym(ufl.grad(theta))
def shear_strain(u):
    (w, theta) = ufl.split(u)
    return theta-ufl.grad(w)
def bending_moment(u):
    DD = ufl.as_tensor([[D, nu*D, 0], [nu*D, D, 0],[0, 0, D*(1-nu)/2.]])
    return voigt2stress(ufl.dot(DD,strain2voigt(curv(u))))
def shear_force(u):
    return F*shear_strain(u)

def clamped_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0],0),np.isclose(x[0],1)),np.logical_or(np.isclose(x[1],0),np.isclose(x[1],1)))
N = 4
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD,N, N,dolfinx.cpp.mesh.CellType.quadrilateral)

# Material parameters for isotropic linear elastic behavior are first defined::
E = 1e3
nu = 0.3


thick = 1e-3#ufl.Constant(mesh,1e-3)

D = E*thick**3/(1-nu**2)/12.
F = E/2/(1+nu)*thick*5./6.

f = -thick**3
deg = 1
We = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), deg)
Te = ufl.VectorElement("Lagrange", mesh.ufl_cell(), deg)
V = dolfinx.FunctionSpace(mesh, ufl.MixedElement([We,Te]))

locate_dofs_w = dolfinx.fem.locate_dofs_geometrical((V.sub(0), V.sub(0).collapse()), clamped_boundary)
locate_dofs_t1= dolfinx.fem.locate_dofs_geometrical((V.sub(1), V.sub(1).collapse()), clamped_boundary)

ubcw=  dolfinx.Function(V.sub(0).collapse())
with ubcw.vector.localForm() as uloc:
    uloc.set(0.)

ubct1 = dolfinx.Function(V.sub(1).collapse())
with ubct1.vector.localForm() as uloc:
    uloc.set(0.)


bcs = [dolfinx.DirichletBC(ubcw, locate_dofs_w, V.sub(0)),
       dolfinx.DirichletBC(ubct1, locate_dofs_t1, V.sub(1))]
       # dolfinx.DirichletBC(ubct2, locate_dofs_t2, V.sub(2))]

u = dolfinx.Function(V)
u_ = ufl.TestFunction(V)
du = ufl.TrialFunction(V)

dx_shear = ufl.dx(metadata={"quadrature_degree": 2*deg-2})

L = f*u_[0]*ufl.dx    #f = ufl.as_vector([0.0, 1.0 / 16])
                   #b1 = - ufl.inner(f, v) * ds(1)
a = ufl.inner(bending_moment(u_), curv(du))*ufl.dx + ufl.dot(shear_force(u_), shear_strain(du))*dx_shear

A = dolfinx.fem.assemble_matrix(a)
A.assemble()
problem = dolfinx.fem.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

(w, theta) = ufl.split(uh)

ww = uh.sub(0)
tt = uh.sub(1)
print("Kirchhoff deflection:", -1.265319087e-3*float(f/D))
print("Reissner-Mindlin FE deflection:", max(abs(ww.compute_point_values())))

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    ww.name = "Deflection"
    tt.name ="Rotation"
    xdmf.write_function(ww)
    xdmf.write_function(tt)
2

I don’t see any obvious errors in your code. Strangely enough the code fails to compile on my local machine on docker. Note that it runs without using reduced integration (naturally the result is incorrect but it at least runs!)

Which version of dolfinx are you running (and how are you running it) ? I’ll take a look at it a bit more over the weekend.

3

There is a bug in dolfinx, as I can reproduce this with:

import dolfinx
from mpi4py import MPI
import ufl

N = 4
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD,N, N,dolfinx.cpp.mesh.CellType.quadrilateral)

deg = 1
We = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), deg)
V = dolfinx.FunctionSpace(mesh, We)

u = ufl.TestFunction(V)
du = ufl.TrialFunction(V)

dx_shear = ufl.dx(metadata={"quadrature_degree": 2*deg-2})
a = ufl.inner(du, u) * ufl.dx + ufl.inner(du, u) * dx_shear
A = dolfinx.fem.assemble_matrix(a)

and error message

libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:897:50: error: use of undeclared identifier 'FE4_C0_D10_Qf81'; did you mean 'FE4_C0_D10_Qa2c'?
        const double J_c0 = coordinate_dofs[0] * FE4_C0_D10_Qf81[0][0][iq][0] + coordinate_dofs[3] * FE4_C0_D10_Qf81[0][0][iq][1] + coordinate_dofs[6] * FE4_C0_D10_Qf81[0][0][iq][2] + coordinate_dofs[9] * FE4_C0_D10_Qf81[0][0][iq][3];
                                                 ^~~~~~~~~~~~~~~
                                                 FE4_C0_D10_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:881:25: note: 'FE4_C0_D10_Qa2c' declared here
    static const double FE4_C0_D10_Qa2c[1][1][1][4] = { { { { -0.4999999999999998, 0.5, -0.4999999999999999, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:897:102: error: use of undeclared identifier 'FE4_C0_D10_Qf81'; did you mean 'FE4_C0_D10_Qa2c'?
        const double J_c0 = coordinate_dofs[0] * FE4_C0_D10_Qf81[0][0][iq][0] + coordinate_dofs[3] * FE4_C0_D10_Qf81[0][0][iq][1] + coordinate_dofs[6] * FE4_C0_D10_Qf81[0][0][iq][2] + coordinate_dofs[9] * FE4_C0_D10_Qf81[0][0][iq][3];
                                                                                                     ^~~~~~~~~~~~~~~
                                                                                                     FE4_C0_D10_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:881:25: note: 'FE4_C0_D10_Qa2c' declared here
    static const double FE4_C0_D10_Qa2c[1][1][1][4] = { { { { -0.4999999999999998, 0.5, -0.4999999999999999, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:897:154: error: use of undeclared identifier 'FE4_C0_D10_Qf81'; did you mean 'FE4_C0_D10_Qa2c'?
        const double J_c0 = coordinate_dofs[0] * FE4_C0_D10_Qf81[0][0][iq][0] + coordinate_dofs[3] * FE4_C0_D10_Qf81[0][0][iq][1] + coordinate_dofs[6] * FE4_C0_D10_Qf81[0][0][iq][2] + coordinate_dofs[9] * FE4_C0_D10_Qf81[0][0][iq][3];
                                                                                                                                                         ^~~~~~~~~~~~~~~
                                                                                                                                                         FE4_C0_D10_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:881:25: note: 'FE4_C0_D10_Qa2c' declared here
    static const double FE4_C0_D10_Qa2c[1][1][1][4] = { { { { -0.4999999999999998, 0.5, -0.4999999999999999, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:897:206: error: use of undeclared identifier 'FE4_C0_D10_Qf81'; did you mean 'FE4_C0_D10_Qa2c'?
        const double J_c0 = coordinate_dofs[0] * FE4_C0_D10_Qf81[0][0][iq][0] + coordinate_dofs[3] * FE4_C0_D10_Qf81[0][0][iq][1] + coordinate_dofs[6] * FE4_C0_D10_Qf81[0][0][iq][2] + coordinate_dofs[9] * FE4_C0_D10_Qf81[0][0][iq][3];
                                                                                                                                                                                                             ^~~~~~~~~~~~~~~
                                                                                                                                                                                                             FE4_C0_D10_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:881:25: note: 'FE4_C0_D10_Qa2c' declared here
    static const double FE4_C0_D10_Qa2c[1][1][1][4] = { { { { -0.4999999999999998, 0.5, -0.4999999999999999, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:898:50: error: use of undeclared identifier 'FE4_C1_D01_Qf81'; did you mean 'FE4_C1_D01_Qa2c'?
        const double J_c3 = coordinate_dofs[1] * FE4_C1_D01_Qf81[0][0][iq][0] + coordinate_dofs[4] * FE4_C1_D01_Qf81[0][0][iq][1] + coordinate_dofs[7] * FE4_C1_D01_Qf81[0][0][iq][2] + coordinate_dofs[10] * FE4_C1_D01_Qf81[0][0][iq][3];
                                                 ^~~~~~~~~~~~~~~
                                                 FE4_C1_D01_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:882:25: note: 'FE4_C1_D01_Qa2c' declared here
    static const double FE4_C1_D01_Qa2c[1][1][1][4] = { { { { -0.4999999999999999, -0.4999999999999999, 0.5, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:898:102: error: use of undeclared identifier 'FE4_C1_D01_Qf81'; did you mean 'FE4_C1_D01_Qa2c'?
        const double J_c3 = coordinate_dofs[1] * FE4_C1_D01_Qf81[0][0][iq][0] + coordinate_dofs[4] * FE4_C1_D01_Qf81[0][0][iq][1] + coordinate_dofs[7] * FE4_C1_D01_Qf81[0][0][iq][2] + coordinate_dofs[10] * FE4_C1_D01_Qf81[0][0][iq][3];
                                                                                                     ^~~~~~~~~~~~~~~
                                                                                                     FE4_C1_D01_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:882:25: note: 'FE4_C1_D01_Qa2c' declared here
    static const double FE4_C1_D01_Qa2c[1][1][1][4] = { { { { -0.4999999999999999, -0.4999999999999999, 0.5, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:898:154: error: use of undeclared identifier 'FE4_C1_D01_Qf81'; did you mean 'FE4_C1_D01_Qa2c'?
        const double J_c3 = coordinate_dofs[1] * FE4_C1_D01_Qf81[0][0][iq][0] + coordinate_dofs[4] * FE4_C1_D01_Qf81[0][0][iq][1] + coordinate_dofs[7] * FE4_C1_D01_Qf81[0][0][iq][2] + coordinate_dofs[10] * FE4_C1_D01_Qf81[0][0][iq][3];
                                                                                                                                                         ^~~~~~~~~~~~~~~
                                                                                                                                                         FE4_C1_D01_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:882:25: note: 'FE4_C1_D01_Qa2c' declared here
    static const double FE4_C1_D01_Qa2c[1][1][1][4] = { { { { -0.4999999999999999, -0.4999999999999999, 0.5, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:898:207: error: use of undeclared identifier 'FE4_C1_D01_Qf81'; did you mean 'FE4_C1_D01_Qa2c'?
        const double J_c3 = coordinate_dofs[1] * FE4_C1_D01_Qf81[0][0][iq][0] + coordinate_dofs[4] * FE4_C1_D01_Qf81[0][0][iq][1] + coordinate_dofs[7] * FE4_C1_D01_Qf81[0][0][iq][2] + coordinate_dofs[10] * FE4_C1_D01_Qf81[0][0][iq][3];
                                                                                                                                                                                                              ^~~~~~~~~~~~~~~
                                                                                                                                                                                                              FE4_C1_D01_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:882:25: note: 'FE4_C1_D01_Qa2c' declared here
    static const double FE4_C1_D01_Qa2c[1][1][1][4] = { { { { -0.4999999999999999, -0.4999999999999999, 0.5, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:899:50: error: use of undeclared identifier 'FE4_C1_D01_Qf81'; did you mean 'FE4_C1_D01_Qa2c'?
        const double J_c1 = coordinate_dofs[0] * FE4_C1_D01_Qf81[0][0][iq][0] + coordinate_dofs[3] * FE4_C1_D01_Qf81[0][0][iq][1] + coordinate_dofs[6] * FE4_C1_D01_Qf81[0][0][iq][2] + coordinate_dofs[9] * FE4_C1_D01_Qf81[0][0][iq][3];
                                                 ^~~~~~~~~~~~~~~
                                                 FE4_C1_D01_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:882:25: note: 'FE4_C1_D01_Qa2c' declared here
    static const double FE4_C1_D01_Qa2c[1][1][1][4] = { { { { -0.4999999999999999, -0.4999999999999999, 0.5, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:899:102: error: use of undeclared identifier 'FE4_C1_D01_Qf81'; did you mean 'FE4_C1_D01_Qa2c'?
        const double J_c1 = coordinate_dofs[0] * FE4_C1_D01_Qf81[0][0][iq][0] + coordinate_dofs[3] * FE4_C1_D01_Qf81[0][0][iq][1] + coordinate_dofs[6] * FE4_C1_D01_Qf81[0][0][iq][2] + coordinate_dofs[9] * FE4_C1_D01_Qf81[0][0][iq][3];
                                                                                                     ^~~~~~~~~~~~~~~
                                                                                                     FE4_C1_D01_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:882:25: note: 'FE4_C1_D01_Qa2c' declared here
    static const double FE4_C1_D01_Qa2c[1][1][1][4] = { { { { -0.4999999999999999, -0.4999999999999999, 0.5, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:899:154: error: use of undeclared identifier 'FE4_C1_D01_Qf81'; did you mean 'FE4_C1_D01_Qa2c'?
        const double J_c1 = coordinate_dofs[0] * FE4_C1_D01_Qf81[0][0][iq][0] + coordinate_dofs[3] * FE4_C1_D01_Qf81[0][0][iq][1] + coordinate_dofs[6] * FE4_C1_D01_Qf81[0][0][iq][2] + coordinate_dofs[9] * FE4_C1_D01_Qf81[0][0][iq][3];
                                                                                                                                                         ^~~~~~~~~~~~~~~
                                                                                                                                                         FE4_C1_D01_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:882:25: note: 'FE4_C1_D01_Qa2c' declared here
    static const double FE4_C1_D01_Qa2c[1][1][1][4] = { { { { -0.4999999999999999, -0.4999999999999999, 0.5, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:899:206: error: use of undeclared identifier 'FE4_C1_D01_Qf81'; did you mean 'FE4_C1_D01_Qa2c'?
        const double J_c1 = coordinate_dofs[0] * FE4_C1_D01_Qf81[0][0][iq][0] + coordinate_dofs[3] * FE4_C1_D01_Qf81[0][0][iq][1] + coordinate_dofs[6] * FE4_C1_D01_Qf81[0][0][iq][2] + coordinate_dofs[9] * FE4_C1_D01_Qf81[0][0][iq][3];
                                                                                                                                                                                                             ^~~~~~~~~~~~~~~
                                                                                                                                                                                                             FE4_C1_D01_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:882:25: note: 'FE4_C1_D01_Qa2c' declared here
    static const double FE4_C1_D01_Qa2c[1][1][1][4] = { { { { -0.4999999999999999, -0.4999999999999999, 0.5, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:900:50: error: use of undeclared identifier 'FE4_C0_D10_Qf81'; did you mean 'FE4_C0_D10_Qa2c'?
        const double J_c2 = coordinate_dofs[1] * FE4_C0_D10_Qf81[0][0][iq][0] + coordinate_dofs[4] * FE4_C0_D10_Qf81[0][0][iq][1] + coordinate_dofs[7] * FE4_C0_D10_Qf81[0][0][iq][2] + coordinate_dofs[10] * FE4_C0_D10_Qf81[0][0][iq][3];
                                                 ^~~~~~~~~~~~~~~
                                                 FE4_C0_D10_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:881:25: note: 'FE4_C0_D10_Qa2c' declared here
    static const double FE4_C0_D10_Qa2c[1][1][1][4] = { { { { -0.4999999999999998, 0.5, -0.4999999999999999, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:900:102: error: use of undeclared identifier 'FE4_C0_D10_Qf81'; did you mean 'FE4_C0_D10_Qa2c'?
        const double J_c2 = coordinate_dofs[1] * FE4_C0_D10_Qf81[0][0][iq][0] + coordinate_dofs[4] * FE4_C0_D10_Qf81[0][0][iq][1] + coordinate_dofs[7] * FE4_C0_D10_Qf81[0][0][iq][2] + coordinate_dofs[10] * FE4_C0_D10_Qf81[0][0][iq][3];
                                                                                                     ^~~~~~~~~~~~~~~
                                                                                                     FE4_C0_D10_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:881:25: note: 'FE4_C0_D10_Qa2c' declared here
    static const double FE4_C0_D10_Qa2c[1][1][1][4] = { { { { -0.4999999999999998, 0.5, -0.4999999999999999, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:900:154: error: use of undeclared identifier 'FE4_C0_D10_Qf81'; did you mean 'FE4_C0_D10_Qa2c'?
        const double J_c2 = coordinate_dofs[1] * FE4_C0_D10_Qf81[0][0][iq][0] + coordinate_dofs[4] * FE4_C0_D10_Qf81[0][0][iq][1] + coordinate_dofs[7] * FE4_C0_D10_Qf81[0][0][iq][2] + coordinate_dofs[10] * FE4_C0_D10_Qf81[0][0][iq][3];
                                                                                                                                                         ^~~~~~~~~~~~~~~
                                                                                                                                                         FE4_C0_D10_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:881:25: note: 'FE4_C0_D10_Qa2c' declared here
    static const double FE4_C0_D10_Qa2c[1][1][1][4] = { { { { -0.4999999999999998, 0.5, -0.4999999999999999, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:900:207: error: use of undeclared identifier 'FE4_C0_D10_Qf81'; did you mean 'FE4_C0_D10_Qa2c'?
        const double J_c2 = coordinate_dofs[1] * FE4_C0_D10_Qf81[0][0][iq][0] + coordinate_dofs[4] * FE4_C0_D10_Qf81[0][0][iq][1] + coordinate_dofs[7] * FE4_C0_D10_Qf81[0][0][iq][2] + coordinate_dofs[10] * FE4_C0_D10_Qf81[0][0][iq][3];
                                                                                                                                                                                                              ^~~~~~~~~~~~~~~
                                                                                                                                                                                                              FE4_C0_D10_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:881:25: note: 'FE4_C0_D10_Qa2c' declared here
    static const double FE4_C0_D10_Qa2c[1][1][1][4] = { { { { -0.4999999999999998, 0.5, -0.4999999999999999, 0.5 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:909:27: error: use of undeclared identifier 'FE3_C0_Qf81'; did you mean 'FE3_C0_Qa2c'?
            t0[i] = fw0 * FE3_C0_Qf81[0][0][iq][i];
                          ^~~~~~~~~~~
                          FE3_C0_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:880:25: note: 'FE3_C0_Qa2c' declared here
    static const double FE3_C0_Qa2c[1][1][1][4] = { { { { 0.25, 0.25, 0.25, 0.25 } } } };
                        ^
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:912:33: error: use of undeclared identifier 'FE3_C0_Qf81'; did you mean 'FE3_C0_Qa2c'?
                A[4 * i + j] += FE3_C0_Qf81[0][0][iq][j] * t0[i];
                                ^~~~~~~~~~~
                                FE3_C0_Qa2c
libffcx_forms_bbce8a268b76eee4dd2f90f2f7ed90996312136e.c:880:25: note: 'FE3_C0_Qa2c' declared here
    static const double FE3_C0_Qa2c[1][1][1][4] = { { { { 0.25, 0.25, 0.25, 0.25 } } } };
                        ^
18 errors generated.
Traceback (most recent call last):
  File "/usr/lib/python3.9/distutils/unixccompiler.py", line 117, in _compile
    self.spawn(compiler_so + cc_args + [src, '-o', obj] +
  File "/usr/lib/python3.9/distutils/ccompiler.py", line 910, in spawn
    spawn(cmd, dry_run=self.dry_run)
  File "/usr/lib/python3.9/distutils/spawn.py", line 87, in spawn
    raise DistutilsExecError(
distutils.errors.DistutilsExecError: command '/usr/bin/clang-12' failed with exit code 1

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/usr/local/lib/python3.9/dist-packages/cffi/ffiplatform.py", line 51, in _build
    dist.run_command('build_ext')
  File "/usr/lib/python3.9/distutils/dist.py", line 985, in run_command
    cmd_obj.run()
  File "/usr/lib/python3.9/distutils/command/build_ext.py", line 340, in run
    self.build_extensions()
  File "/usr/lib/python3.9/distutils/command/build_ext.py", line 449, in build_extensions
    self._build_extensions_serial()
  File "/usr/lib/python3.9/distutils/command/build_ext.py", line 474, in _build_extensions_serial
    self.build_extension(ext)
  File "/usr/lib/python3.9/distutils/command/build_ext.py", line 529, in build_extension
    objects = self.compiler.compile(sources,
  File "/usr/lib/python3.9/distutils/ccompiler.py", line 574, in compile
    self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts)
  File "/usr/lib/python3.9/distutils/unixccompiler.py", line 120, in _compile
    raise CompileError(msg)
distutils.errors.CompileError: command '/usr/bin/clang-12' failed with exit code 1

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/root/shared/debug/mwe.py", line 19, in <module>
    A = dolfinx.fem.assemble_matrix(a)
  File "/usr/lib/python3.9/functools.py", line 877, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/root/shared/dolfinx/python/dolfinx/fem/assemble.py", line 278, in assemble_matrix
    A = cpp.fem.create_matrix(_create_cpp_form(a))
  File "/root/shared/dolfinx/python/dolfinx/fem/assemble.py", line 27, in _create_cpp_form
    return Form(form)._cpp_object
  File "/root/shared/dolfinx/python/dolfinx/fem/form.py", line 42, in __init__
    self._ufc_form, module, self._code = jit.ffcx_jit(
  File "/root/shared/dolfinx/python/dolfinx/jit.py", line 61, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/root/shared/dolfinx/python/dolfinx/jit.py", line 215, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], parameters=p_ffcx, **p_jit)
  File "/root/shared/ffcx/ffcx/codegeneration/jit.py", line 167, in compile_forms
    impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
  File "/root/shared/ffcx/ffcx/codegeneration/jit.py", line 252, in _compile_objects
    ffibuilder.compile(tmpdir=cache_dir, verbose=True, debug=cffi_debug)
  File "/usr/local/lib/python3.9/dist-packages/cffi/api.py", line 725, in compile
    return recompile(self, module_name, source, tmpdir=tmpdir,
  File "/usr/local/lib/python3.9/dist-packages/cffi/recompiler.py", line 1564, in recompile
    outputfilename = ffiplatform.compile('.', ext,
  File "/usr/local/lib/python3.9/dist-packages/cffi/ffiplatform.py", line 22, in compile
    outputfilename = _build(tmpdir, ext, compiler_verbose, debug)
  File "/usr/local/lib/python3.9/dist-packages/cffi/ffiplatform.py", line 58, in _build
    raise VerificationError('%s: %s' % (e.__class__.__name__, e))
cffi.VerificationError: CompileError: command '/usr/bin/clang-12' failed with exit code 1

I will submit an issue and keep you posted.

1 Like
4

A workaround for now would be the following:

from mpi4py import MPI
import numpy as np
import dolfinx.cpp
import dolfinx.io
import dolfinx.la
import ufl
from petsc4py import PETSc

def border(x):
    return np.isclose(x[0], 0.0)

def strain2voigt(eps):
    return ufl.as_vector([eps[0, 0], eps[1, 1], 2*eps[0, 1]])
def voigt2stress(S):
    return ufl.as_tensor([[S[0], S[2]], [S[2], S[1]]])
def curv(u):
    (w, theta) = ufl.split(u)
    return ufl.sym(ufl.grad(theta))
def shear_strain(u):
    (w, theta) = ufl.split(u)
    return theta-ufl.grad(w)
def bending_moment(u):
    DD = ufl.as_tensor([[D, nu*D, 0], [nu*D, D, 0],[0, 0, D*(1-nu)/2.]])
    return voigt2stress(ufl.dot(DD,strain2voigt(curv(u))))
def shear_force(u):
    return F*shear_strain(u)

def clamped_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0],0),np.isclose(x[0],1)),np.logical_or(np.isclose(x[1],0),np.isclose(x[1],1)))
N = 4
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD,N, N,dolfinx.cpp.mesh.CellType.quadrilateral)

# Material parameters for isotropic linear elastic behavior are first defined::
E = 1e3
nu = 0.3


thick = 1e-3#ufl.Constant(mesh,1e-3)

D = E*thick**3/(1-nu**2)/12.
F = E/2/(1+nu)*thick*5./6.

f = -thick**3
deg = 1
We = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), deg)
Te = ufl.VectorElement("Lagrange", mesh.ufl_cell(), deg)
V = dolfinx.FunctionSpace(mesh, ufl.MixedElement([We,Te]))

locate_dofs_w = dolfinx.fem.locate_dofs_geometrical((V.sub(0), V.sub(0).collapse()), clamped_boundary)
locate_dofs_t1= dolfinx.fem.locate_dofs_geometrical((V.sub(1), V.sub(1).collapse()), clamped_boundary)

ubcw=  dolfinx.Function(V.sub(0).collapse())
with ubcw.vector.localForm() as uloc:
    uloc.set(0.)

ubct1 = dolfinx.Function(V.sub(1).collapse())
with ubct1.vector.localForm() as uloc:
    uloc.set(0.)


bcs = [dolfinx.DirichletBC(ubcw, locate_dofs_w, V.sub(0)),
       dolfinx.DirichletBC(ubct1, locate_dofs_t1, V.sub(1))]
       # dolfinx.DirichletBC(ubct2, locate_dofs_t2, V.sub(2))]

u = dolfinx.Function(V)
u_ = ufl.TestFunction(V)
du = ufl.TrialFunction(V)

dx_shear = ufl.dx(metadata={"quadrature_degree": 2*deg-2})

L = f*u_[0]*ufl.dx    #f = ufl.as_vector([0.0, 1.0 / 16])
                   #b1 = - ufl.inner(f, v) * ds(1)
a0 = ufl.inner(bending_moment(u_), curv(du))*ufl.dx 
a1 = ufl.dot(shear_force(u_), shear_strain(du))*dx_shear

A = dolfinx.fem.assemble_matrix(a0, bcs)
A.assemble()
dolfinx.fem.assemble_matrix(A, a1, bcs)
A.assemble()

b = dolfinx.fem.assemble_vector(L)
dolfinx.fem.apply_lifting(b, [a0, a1], [bcs, bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(b, bcs)

uh = dolfinx.Function(V)
# Setup petsc solver
solver = PETSc.KSP().create(mesh.mpi_comm())
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
pc = solver.getPC()
pc.setType(PETSc.PC.Type.LU)
solver.solve(b, uh.vector)
print(uh.vector.array)

(w, theta) = ufl.split(uh)

ww = uh.sub(0)
tt = uh.sub(1)
print("Kirchhoff deflection:", -1.265319087e-3*float(f/D))
print("Reissner-Mindlin FE deflection:", max(abs(ww.compute_point_values())))

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    ww.name = "Deflection"
    tt.name ="Rotation"
    xdmf.write_function(ww)
    xdmf.write_function(tt)

which yields results

Kirchhoff deflection: 1.381728443004e-05
Reissner-Mindlin FE deflection: [1.32270274e-05]
2 Likes
5

Thank you for the help!

6

Strange indeed. I am using the latest docker image as well running on ubuntu 20.04.

7

I am not sure what’s wrong, but when I run this code I get the following output:

Blockquote
Kirchhoff deflection: 1.381728443004e-05
Reissner-Mindlin FE deflection: [0.]

Could this be due to any version incompatibility?

8

Quickly checking, if you assemble the form that is integrated using reduced integration first, i.e.,

a0 = ufl.inner(bending_moment(u_), curv(du))*ufl.dx 
a1 = ufl.dot(shear_force(u_), shear_strain(du))*dx_shear

A = dolfinx.fem.assemble_matrix(a1, bcs)
A.assemble()
dolfinx.fem.assemble_matrix(A, a0, bcs)
A.assemble()

then you do get something more reasonable (corresponding to N = 50)

Kirchhoff deflection: 1.381728443004e-05
Reissner-Mindlin FE deflection: [1.38134448e-05]

However it isn’t obvious as to why that is the case. I would let Jorgen open an issue on the Github issue tracker and see it from there. Indeed with the order reversed I am getting a zero displacement on the latest docker image too. I’ll investigate later when I have more time, but for now there does seem to be a workaround (integrating and assembling the terms explicitly and separately)

2 Likes
9

The issue should be fixed now: FFCx can't handle forms with multiple quadrature degrees/rules in single form · Issue #394 · FEniCS/ffcx · GitHub

It would be great if you could check if your original code now runs.

1 Like
10

Hi, I tried running the program with dolfinx/dolfinx:latest docker image, but the problem still persists. I’m not sure if the docker image is still using the old version or not. Can you please help with this?

Thanks.

11

The dockerimage was updated about three hours ago. Try to pull again.

12

Yes, I can confirm that the fix works. Thank you so much for the help!