ReferenceGrad error with Quad elements

dolfinx
, ,
1

I am trying to implement a structural shell analysis formulation. It works with triangular elements but whenever I switch to quads, the code throws an error “Currently no support for ReferenceGrad in CoefficientDerivative.”

image
image890×602 20.4 KB

Here is the code I am suing:

from mpi4py import MPI
import numpy as np
from ufl import (
    Jacobian, CellNormal, CellDiameter, 
    cross, dot, sqrt, as_matrix,
    TrialFunction, TestFunction
)
from basix.ufl import element, mixed_element
from dolfinx.fem import (
    Function, dirichletbc,
    functionspace, locate_dofs_topological,
)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import gmshio
from ufl import derivative, dx, grad, split, dot, sym, as_vector


# Read mesh
msh, markers, facets = gmshio.read_from_msh("Ibeam/I_beam_quad.msh", MPI.COMM_WORLD)
num_nodes = msh.geometry.x.shape[0]
num_cell = msh.topology.original_cell_index.shape[0]

# Material and geometry
E = 210e9
nu = 0.3
h = 1e-3
A = E*h/(1-nu**2)
D = A*h**2/12.
Fs = E/2/(1+nu)*h*5./6.

# Element definition
u1 = element("Lagrange", msh.basix_cell(), 1, shape=(3,))   
u2 = element("Lagrange", msh.basix_cell(), 1, shape=(3,))            
V = functionspace(msh, mixed_element([u1, u2]))
w = Function(V)
w_t = TrialFunction(V)
dw = TestFunction(V)
q, theta = split(w)

# Constitutive matrices
A_mat = np.tile(np.array([A, nu*A, 0, nu*A, A, 0, 0, 0, A*(1-nu)/2.]), (num_cell, 1))
B_mat = np.tile(np.zeros((9, )), (num_cell, 1))
D_mat = np.tile(np.array([D, nu*D, 0, nu*D, D, 0, 0, 0, D*(1-nu)/2.]), (num_cell, 1))
As_mat = np.tile(np.array([Fs, 0., 0., Fs]), (num_cell, 1))

Vabd = functionspace(msh, ("DG", 0, (3, 3)))
Vs = functionspace(msh, ("DG", 0, (2, 2)))
D_max = np.amax(D_mat)

A_fun = Function(Vabd)
B_fun = Function(Vabd)
D_fun = Function(Vabd)
As_fun = Function(Vs)

A_fun.x.petsc_vec.setArray(A_mat)
B_fun.x.petsc_vec.setArray(B_mat)
D_fun.x.petsc_vec.setArray(D_mat)
As_fun.x.petsc_vec.setArray(As_mat)

# Loading (self-weight)
f_array = np.tile(np.array([0, 0.2, -1]), (num_nodes, 1))
VF = functionspace(msh, ("CG", 1, (3,)))
f = Function(VF)
f.x.petsc_vec.setArray(f_array)

# Local basis
J = Jacobian(msh)
n = CellNormal(msh)
e3 = n/sqrt(dot(n, n))
t = as_vector([J[0, 0], J[1, 0], J[2, 0]])
e1 = t/sqrt(dot(t, t))
e2 = cross(e3, e1)

# Strains
H = as_matrix([[e1[i] for i in range(0, 3)],
                [e2[i] for i in range(0, 3)]])
gradq = dot(H, dot(grad(q), H.T))
eps = sym(gradq)
kappa = sym(dot(H, dot(grad(cross(e3, theta)), H.T)))
gamma = dot(H, grad(dot(q, e3))) - dot(H, cross(e3, theta))
omega = (gradq[0, 1] - gradq[1, 0])/2 + dot(theta, e3)

eps_bar = as_vector([eps[0, 0], eps[1, 1], 2*eps[0, 1]])
kappa_bar = as_vector([kappa[0, 0], kappa[1, 1], 2*kappa[0, 1]])

# Stresses
N = A_fun*eps_bar + B_fun*kappa_bar
M = B_fun*eps_bar + D_fun*kappa_bar
Q = As_fun*gamma

# Strain energy
h_mesh = CellDiameter(msh)
dx_inplane = dx(metadata={"quadrature_degree":1})
dx_shear = dx(metadata={"quadrature_degree":0})
U = 0.5*dot(N, eps_bar)*dx_inplane + \
    0.5*dot(M, kappa_bar)*dx_inplane + \
    0.5*dot(Q, gamma)*dx_shear + \
    0.5*(D_max*12*omega/h_mesh**2)*omega*dx
    
# BCs
Vu, _ = V.sub(0).collapse()
left_dofs = locate_dofs_topological((V.sub(0), Vu), 1, facets.find(1))
right_dofs = locate_dofs_topological((V.sub(0), Vu), 1, facets.find(2))

uD = Function(Vu)
bcs = [
    dirichletbc(uD, left_dofs, V.sub(0)),
    dirichletbc(uD, right_dofs, V.sub(0)),
]

# Total energy
Wext = dot(f, q)*dx
Pi = U - Wext

F = derivative(Pi, w, dw)
K = derivative(F, w, w_t)

problem = LinearProblem(
    K,
    -F,
    bcs=bcs,
    petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"},
)
w = problem.solve()

The mesh file with triangular elements: I_beam.msh

$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
17 16 5 0
1 0.07499999999999998 -0.05 -0.1299038105676658 0 
2 0.07499999999999998 0 -0.1299038105676658 0 
3 0.07499999999999998 0.05 -0.1299038105676658 0 
4 -0.07499999999999998 -0.1 0.1299038105676658 0 
5 -0.07499999999999998 0 0.1299038105676658 0 
6 -0.07499999999999998 0.1 0.1299038105676658 0 
7 1.925 -0.05 -0.1299038105676658 0 
8 1.925 0 -0.1299038105676658 0 
12 1 -0.05 -1.732050807568877 0 
14 1 0 -1.732050807568877 0 
16 1.925 0.05 -0.1299038105676658 0 
22 1 0.05 -1.732050807568877 0 
24 2.075 0 0.1299038105676658 0 
25 2.075 -0.1 0.1299038105676658 0 
30 1 -0.1 -1.732050807568877 0 
32 2.075 0.1 0.1299038105676658 0 
38 1 0.1 -1.732050807568877 0 
1 0.07499999999999998 -0.05 -0.1299038105676658 0.07499999999999998 0 -0.1299038105676658 1 1 2 1 -2 
2 0.07499999999999998 0 -0.1299038105676658 0.07499999999999998 0.05 -0.1299038105676658 1 1 2 2 -3 
3 -0.07499999999999998 0 -0.1299038105676658 0.07499999999999998 0 0.1299038105676658 1 1 2 2 -5 
4 -0.07499999999999998 -0.1 0.1299038105676658 -0.07499999999999998 0 0.1299038105676658 1 1 2 4 -5 
5 -0.07499999999999998 0 0.1299038105676658 -0.07499999999999998 0.1 0.1299038105676658 1 1 2 5 -6 
6 1.925 -0.05 -0.1299038105676658 1.925 0 -0.1299038105676658 1 2 2 7 -8 
7 0.07499999999999973 -0.05 -0.1299038105676664 1.925000000000001 -0.05 0.1148192852329688 0 2 1 -7 
8 0.07499999999999973 0 -0.1299038105676664 1.925000000000001 0 0.1148192852329688 0 2 2 -8 
10 1.925 0 -0.1299038105676658 1.925 0.05 -0.1299038105676658 1 2 2 8 -16 
12 0.07499999999999973 0.05 -0.1299038105676664 1.925000000000001 0.05 0.1148192852329688 0 2 3 -16 
14 1.925 0 -0.1299038105676658 2.075 0 0.1299038105676658 1 2 2 8 -24 
16 -0.07500000000000018 0 0.129903810567666 2.075000000000001 0 0.4143117327143495 0 2 5 -24 
18 2.075 -0.1 0.1299038105676658 2.075 0 0.1299038105676658 1 2 2 25 -24 
19 -0.07500000000000018 -0.1 0.129903810567666 2.075000000000001 -0.1 0.4143117327143495 0 2 4 -25 
22 2.075 0 0.1299038105676658 2.075 0.1 0.1299038105676658 1 2 2 24 -32 
24 -0.07500000000000018 0.1 0.129903810567666 2.075000000000001 0.1 0.4143117327143495 0 2 6 -32 
9 0.07499999999999973 -0.05 -0.1299038105676664 1.925000000000001 0 0.1148192852329688 1 1 4 1 8 -6 -7 
13 0.07499999999999973 0 -0.1299038105676664 1.925000000000001 0.05 0.1148192852329688 1 1 4 2 12 -10 -8 
17 -0.07500000000000018 0 -0.1299038105676664 2.075000000000001 0 0.4143117327143495 1 1 4 3 16 -14 -8 
21 -0.07500000000000018 -0.1 0.1299038105676658 2.075000000000001 0 0.4143117327143495 1 1 4 4 16 -18 -19 
25 -0.07500000000000018 0 0.1299038105676658 2.075000000000001 0.1 0.4143117327143495 1 1 4 5 24 -22 -16 
$EndEntities
$Nodes
33 24 1 24
0 1 0 1
1
0.07499999999999998 -0.05 -0.1299038105676658
0 2 0 1
2
0.07499999999999998 0 -0.1299038105676658
0 3 0 1
3
0.07499999999999998 0.05 -0.1299038105676658
0 4 0 1
4
-0.07499999999999998 -0.1 0.1299038105676658
0 5 0 1
5
-0.07499999999999998 0 0.1299038105676658
0 6 0 1
6
-0.07499999999999998 0.1 0.1299038105676658
0 7 0 1
7
1.925 -0.05 -0.1299038105676658
0 8 0 1
8
1.925 0 -0.1299038105676658
0 16 0 1
9
1.925 0.05 -0.1299038105676658
0 24 0 1
10
2.075 0 0.1299038105676658
0 25 0 1
11
2.075 -0.1 0.1299038105676658
0 32 0 1
12
2.075 0.1 0.1299038105676658
1 1 0 0
1 2 0 0
1 3 0 2
13
14
0.02500000000011773 0 -0.04330127018942585
-0.02499999999987315 0 0.04330127018900226
1 4 0 0
1 5 0 0
1 6 0 0
1 7 0 1
15
0.9999999999999999 -0.05 0.1179491924311227
1 8 0 1
16
0.9999999999999999 0 0.1179491924311227
1 10 0 0
1 12 0 1
17
0.9999999999999999 0.05 0.1179491924311227
1 14 0 2
18
19
1.974999999999882 0 -0.04330127018942576
2.024999999999873 0 0.04330127018900232
1 16 0 1
20
0.9999999999999999 0 0.4179491924311229
1 18 0 0
1 19 0 1
21
0.9999999999999999 -0.1 0.4179491924311229
1 22 0 0
1 24 0 1
22
0.9999999999999999 0.1 0.4179491924311229
2 9 0 0
2 13 0 0
2 17 0 2
23
24
0.9999999999999999 0 0.2179491924308874
1 0 0.3179491924308693
2 21 0 0
2 25 0 0
$EndNodes
$Elements
15 42 1 42
1 1 1 1
1 1 2 
1 2 1 1
2 2 3 
1 3 1 3
3 2 13 
4 13 14 
5 14 5 
1 4 1 1
6 4 5 
1 5 1 1
7 5 6 
1 6 1 1
8 7 8 
1 10 1 1
9 8 9 
1 14 1 3
10 8 18 
11 18 19 
12 19 10 
1 18 1 1
13 11 10 
1 22 1 1
14 10 12 
2 9 2 4
15 1 2 16 
16 1 16 15 
17 15 16 8 
18 15 8 7 
2 13 2 4
19 2 3 17 
20 2 17 16 
21 16 17 9 
22 16 9 8 
2 17 2 12
23 2 13 23 
24 2 23 16 
25 16 23 18 
26 16 18 8 
27 13 14 24 
28 13 24 23 
29 23 24 19 
30 23 19 18 
31 14 5 20 
32 14 20 24 
33 24 20 10 
34 24 10 19 
2 21 2 4
35 4 5 20 
36 4 20 21 
37 21 20 10 
38 21 10 11 
2 25 2 4
39 5 6 22 
40 5 22 20 
41 20 22 12 
42 20 12 10 
$EndElements

Mesh file with quads: I_beam_quad.msh

$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
17 16 5 0
1 0.07499999999999998 -0.05 -0.1299038105676658 0 
2 0.07499999999999998 0 -0.1299038105676658 0 
3 0.07499999999999998 0.05 -0.1299038105676658 0 
4 -0.07499999999999998 -0.1 0.1299038105676658 0 
5 -0.07499999999999998 0 0.1299038105676658 0 
6 -0.07499999999999998 0.1 0.1299038105676658 0 
7 1.925 -0.05 -0.1299038105676658 0 
8 1.925 0 -0.1299038105676658 0 
12 1 -0.05 -1.732050807568877 0 
14 1 0 -1.732050807568877 0 
16 1.925 0.05 -0.1299038105676658 0 
22 1 0.05 -1.732050807568877 0 
24 2.075 0 0.1299038105676658 0 
25 2.075 -0.1 0.1299038105676658 0 
30 1 -0.1 -1.732050807568877 0 
32 2.075 0.1 0.1299038105676658 0 
38 1 0.1 -1.732050807568877 0 
1 0.07499999999999998 -0.05 -0.1299038105676658 0.07499999999999998 0 -0.1299038105676658 1 1 2 1 -2 
2 0.07499999999999998 0 -0.1299038105676658 0.07499999999999998 0.05 -0.1299038105676658 1 1 2 2 -3 
3 -0.07499999999999998 0 -0.1299038105676658 0.07499999999999998 0 0.1299038105676658 1 1 2 2 -5 
4 -0.07499999999999998 -0.1 0.1299038105676658 -0.07499999999999998 0 0.1299038105676658 1 1 2 4 -5 
5 -0.07499999999999998 0 0.1299038105676658 -0.07499999999999998 0.1 0.1299038105676658 1 1 2 5 -6 
6 1.925 -0.05 -0.1299038105676658 1.925 0 -0.1299038105676658 1 2 2 7 -8 
7 0.07499999999999973 -0.05 -0.1299038105676664 1.925000000000001 -0.05 0.1148192852329688 0 2 1 -7 
8 0.07499999999999973 0 -0.1299038105676664 1.925000000000001 0 0.1148192852329688 0 2 2 -8 
10 1.925 0 -0.1299038105676658 1.925 0.05 -0.1299038105676658 1 2 2 8 -16 
12 0.07499999999999973 0.05 -0.1299038105676664 1.925000000000001 0.05 0.1148192852329688 0 2 3 -16 
14 1.925 0 -0.1299038105676658 2.075 0 0.1299038105676658 1 2 2 8 -24 
16 -0.07500000000000018 0 0.129903810567666 2.075000000000001 0 0.4143117327143495 0 2 5 -24 
18 2.075 -0.1 0.1299038105676658 2.075 0 0.1299038105676658 1 2 2 25 -24 
19 -0.07500000000000018 -0.1 0.129903810567666 2.075000000000001 -0.1 0.4143117327143495 0 2 4 -25 
22 2.075 0 0.1299038105676658 2.075 0.1 0.1299038105676658 1 2 2 24 -32 
24 -0.07500000000000018 0.1 0.129903810567666 2.075000000000001 0.1 0.4143117327143495 0 2 6 -32 
9 0.07499999999999973 -0.05 -0.1299038105676664 1.925000000000001 0 0.1148192852329688 1 1 4 1 8 -6 -7 
13 0.07499999999999973 0 -0.1299038105676664 1.925000000000001 0.05 0.1148192852329688 1 1 4 2 12 -10 -8 
17 -0.07500000000000018 0 -0.1299038105676664 2.075000000000001 0 0.4143117327143495 1 1 4 3 16 -14 -8 
21 -0.07500000000000018 -0.1 0.1299038105676658 2.075000000000001 0 0.4143117327143495 1 1 4 4 16 -18 -19 
25 -0.07500000000000018 0 0.1299038105676658 2.075000000000001 0.1 0.4143117327143495 1 1 4 5 24 -22 -16 
$EndEntities
$Nodes
33 55 1 55
0 1 0 1
1
0.07499999999999998 -0.05 -0.1299038105676658
0 2 0 1
2
0.07499999999999998 0 -0.1299038105676658
0 3 0 1
3
0.07499999999999998 0.05 -0.1299038105676658
0 4 0 1
4
-0.07499999999999998 -0.1 0.1299038105676658
0 5 0 1
5
-0.07499999999999998 0 0.1299038105676658
0 6 0 1
6
-0.07499999999999998 0.1 0.1299038105676658
0 7 0 1
7
1.925 -0.05 -0.1299038105676658
0 8 0 1
8
1.925 0 -0.1299038105676658
0 16 0 1
9
1.925 0.05 -0.1299038105676658
0 24 0 1
10
2.075 0 0.1299038105676658
0 25 0 1
11
2.075 -0.1 0.1299038105676658
0 32 0 1
12
2.075 0.1 0.1299038105676658
1 1 0 1
13
0.07499999999999998 -0.02500000001223876 -0.1299038105676658
1 2 0 1
14
0.07499999999999998 0.02499999998776124 -0.1299038105676658
1 3 0 1
15
3.67162689141054e-11 0 -6.359443527337305e-11
1 4 0 1
16
-0.07499999999999998 -0.05000000002447752 0.1299038105676658
1 5 0 1
17
-0.07499999999999998 0.04999999997552249 0.1299038105676658
1 6 0 1
18
1.925 -0.02500000001223876 -0.1299038105676658
1 7 0 3
19
20
21
0.5211847671599293 -0.05 0.05491197122655911
1.478815232098933 -0.05 0.05491197142514648
0.9999999999999999 -0.05 0.1179491924311227
1 8 0 3
22
23
24
0.5211847671599293 0 0.05491197122655911
1.478815232098933 0 0.05491197142514648
0.9999999999999999 0 0.1179491924311227
1 10 0 1
25
1.925 0.02499999998776124 -0.1299038105676658
1 12 0 3
26
27
28
0.5211847671599293 0.05 0.05491197122655911
1.478815232098933 0.05 0.05491197142514648
0.9999999999999999 0.05 0.1179491924311227
1 14 0 1
29
1.999999999963284 0 -6.359443527337305e-11
1 16 0 3
30
31
32
0.443539053068413 0 0.3446897189630251
1.556460947134918 0 0.3446897189085425
0.9999999999999999 0 0.4179491924311229
1 18 0 1
33
2.075 -0.05000000002447752 0.1299038105676658
1 19 0 3
34
35
36
0.443539053068413 -0.1 0.3446897189630251
1.556460947134918 -0.1 0.3446897189085425
0.9999999999999999 -0.1 0.4179491924311229
1 22 0 1
37
2.075 0.04999999997552249 0.1299038105676658
1 24 0 3
38
39
40
0.443539053068413 0.1 0.3446897189630251
1.556460947134918 0.1 0.3446897189085425
0.9999999999999999 0.1 0.4179491924311229
2 9 0 3
41
42
43
0.9999999999999999 -0.02499999999914538 0.1179491924311227
0.521184767159951 -0.02500000000040423 0.0549119712265651
1.478815232101 -0.02500000000040423 0.0549119714245927
2 13 0 3
44
45
46
0.9999999999999999 0.02500000000085462 0.1179491924311227
0.521184767159951 0.02499999999959578 0.0549119712265651
1.478815232101 0.02499999999959578 0.0549119714245927
2 17 0 3
47
48
49
0.999999999988275 0 0.2622288582336462
0.4837854862407983 0 0.1944804172286039
1.516214513477559 0 0.1944804173040565
2 21 0 3
50
51
52
0.9999999999999999 -0.04999999999766757 0.4179491924311229
0.4435390530684149 -0.05000000000045631 0.3446897189630251
1.556460947134878 -0.05000000000045629 0.3446897189085532
2 25 0 3
53
54
55
0.9999999999999999 0.05000000000233244 0.4179491924311229
0.4435390530684373 0.04999999999954369 0.3446897189630309
1.556460947134928 0.04999999999954371 0.3446897189085394
$EndNodes
$Elements
15 60 1 60
1 1 1 2
1 1 13 
2 13 2 
1 2 1 2
3 2 14 
4 14 3 
1 3 1 2
5 2 15 
6 15 5 
1 4 1 2
7 4 16 
8 16 5 
1 5 1 2
9 5 17 
10 17 6 
1 6 1 2
11 7 18 
12 18 8 
1 10 1 2
13 8 25 
14 25 9 
1 14 1 2
15 8 29 
16 29 10 
1 18 1 2
17 11 33 
18 33 10 
1 22 1 2
19 10 37 
20 37 12 
2 9 3 8
21 1 13 42 19 
22 13 2 22 42 
23 42 22 24 41 
24 19 42 41 21 
25 21 41 43 20 
26 41 24 23 43 
27 43 23 8 18 
28 20 43 18 7 
2 13 3 8
29 2 14 45 22 
30 14 3 26 45 
31 45 26 28 44 
32 22 45 44 24 
33 24 44 46 23 
34 44 28 27 46 
35 46 27 9 25 
36 23 46 25 8 
2 17 3 8
37 2 15 48 22 
38 15 5 30 48 
39 48 30 32 47 
40 22 48 47 24 
41 24 47 49 23 
42 47 32 31 49 
43 49 31 10 29 
44 23 49 29 8 
2 21 3 8
45 4 16 51 34 
46 16 5 30 51 
47 51 30 32 50 
48 34 51 50 36 
49 36 50 52 35 
50 50 32 31 52 
51 52 31 10 33 
52 35 52 33 11 
2 25 3 8
53 5 17 54 30 
54 17 6 38 54 
55 54 38 40 53 
56 30 54 53 32 
57 32 53 55 31 
58 53 40 39 55 
59 55 39 12 37 
60 31 55 37 10 
$EndElements

I would love to get some insights into this, as I just started using Fenicsx a few months back. I am using the latest release (v0.9).