SAND2009-7358P
Computational mathematics and algorithms
Capability Couplings with Intrepid
Trilinos User’s Group Meeting
November 4, 2009
Kara Peterson
Pavel Bochev
David Hensinger
Denis Ridzal
Chris Siefert
Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company,
for the United States Department of Energy!s National Nuclear Security Administration
under contract DE-AC04-94AL85000.
SAND2009-7358P
Computational mathematics and algorithms
Solving PDEs with FEM in Trilinos
Pamgen: Mesh Shards: Ref. Cell Topology
" ˆ
"
ˆ
F :" # "
!
Intrepid: Transform, Integration,
! Discrete Spaces and! Operators
!
K e ij = % "# (x)"# (x)dx
i j
$
Lu = f in !
!
Epetra: Global Matrices
K uh = b
ML: Solve
uh = K-1 b
SAND2009-7358P
Computational mathematics and algorithms
Example Drivers
Poisson Equation with H(grad) Elements
%"u = f in #
& ! example_Poisson, example_03
'u = 0 on $# Functionality
• Intrepid for element matrix
LSFEM Div-Curl System with H(curl) Elements assembly
! • Epetra for global matrices
'" # u = f in $
) • PAMGEN for meshing
(" % u = g in $ ! example_CurlLSFEM, example_01
)n # u = 0 on &$ • ML to solve
*
• Intrepid for error estimates
• Runs on multiple processors
!
LSFEM Div-Curl System with H(div) Elements • Intrepid for element matrix
assembly
'" # u = f in $ • Epetra for global matrices
)
(" % u = g in $ ! example_DivLSFEM, example_02
)n % u = 0 on &$
*
Located in:
! • Trilinos/packages/trilinoscouplings/examples/scaling
• Trilinos/packages/intrepid/example/Drivers
SAND2009-7358P
Computational mathematics and algorithms
Current Status of Intrepid
Released with Trilinos 10.0: a suite of “mathematical” (stateless) tools for
• Cell topology
– Topology type (base vs. extended), cell type (standard = has reference cell vs. non-
standard), adjacency queries, permutations,…
• Cell geometry
– Jacobians, maps to and from reference cells (if applicable), surface normal, line
tangent, subcell measure,…
• Cell integration
– Cubature points on standard cells, direct integration on non-standard cells
• Discrete spaces on a cell workset
– Definition of finite dimensional approximation spaces, reconstruction (interpolation)
operators, transformation rules,…
• Discrete operators on a cell workset
– Basic differential operators acting on scalar, vector and tensor fields
• Discrete functionals on a cell workset
– Basic linear functionals acting on scalar, vector and tensor fields
• Utilities
– Everything else you may need but that does not fit in any of the above categories
SAND2009-7358P
Computational mathematics and algorithms
Mapping Functionality to Software
Cell Topology
CellTools src/Cell
Cell Geometry Shards::CellTopology
Cell Integration
Integration
Discrete Spaces
Basis src/Discretization
Discrete Operators
FunctionSpaceTools
Discrete Functionals
RealSpaceTools
ArrayTools
Utilities src/Shared
FieldContainer
Polylib
SAND2009-7358P
Computational mathematics and algorithms
Intrepid Example: Compatible LSFEM for
div-curl Systems
A model Div-Curl System
'" # u = f in $
) "#R3 ! bounded contractible domain
(" % u = g in $
)n # u = 0 on &$ $" ! Lipschitz continuous
*
A continuous least-squares principle
( J (u;f,g) = " # u $ f 2 + " % u $ g
*
2
! ) 0 0 min J (u;f,g)
X
* X = H 0 (curl,&) ' H ( div,&)
+
A semi-conforming LSFEM:
!
! ) J u ;f,g = " # u $ f 2 + " * % u $ g 2
+ ( h ) h 0 h h
* 0
X h = Ch is curl-conforming FE,
0
+ X h & H 0 (curl,') ( H ( div,')
, e.g., Nedelec
We gain some and loose some:
!
! • curl-conforming FE: natural for the tangential boundary condition
• curl-conforming FE: not in the domain of div - need discrete approximation!
SAND2009-7358P
Computational mathematics and algorithms
From Discretization to Linear System
For clarity we now focus on the lowest-order case
Ch " the lowest-order Nedelec space of the 1st kind
0
h
G0 " the lowest-order nodal C0 space
Discrete least-squares problem
! (" # uh ," # v h ) 0 + ("*h $ uh ,"*h $ v h ) 0 = (f," # v h ) 0 + ( g,"*h $ v h ) 0 %v h & Ch
0
(K + MC DVE M"1DT MC ) uh
G VE = b curl-curl + grad-div matrices
!
" MC curl-conforming mass matrix Note: Can use any O(h)
$
!
where # MG grad-conforming mass matrix approximation for M G
$ We will use mass lumping.
% DVE vertex-to-edge incidence matrix
!
!
! We already have an AMG solver for this system in Trilinos - a subsolver of the eddy
current
! ! solver developed in Bochev, Hu, Tuminaro and Siefert, SISC 2008.
! The ML solver requires 4 matrices: K , MC , MG , and DVE
! ! ! !
SAND2009-7358P
Computational mathematics and algorithms
Using Intrepid to Assemble the System
Procedure
1. Set up: define cell, integration method (cubature), and reference cell basis
a. Choose cell type: Hexahedron and generate grid (PAMGEN)
b. Define Cubature on reference cell and on a typical face (for B.C.)
c. Define & evaluate reference cell H(curl) basis: Basis_HCURL_HEX_I1_FEM
d. Define & evaluate reference cell H(grad) basis: Basis_HGRAD_HEX_C1_FEM
2. Set up: get transformation data
a. Define a physical cell workset
b. Define multi-dimensional arrays for cell data (FieldContainer)
c. Get pullback data at cubature points (Jacobians, determinants, face area…)
d. Get basis signs (=edge permutations) for H(curl) basis
3. Assemble
a. Assemble the H(grad) mass matrix
b. Assemble the H(curl) mass matrix
c. Assemble the H(curl) stiffness matrix
d. Generate the incidence matrix
e. Assemble right-hand side
4. Solve: using ML
SAND2009-7358P
Computational mathematics and algorithms
Step 1: Define cell, cubature, and basis
Define Cell Topology with Shards
typedef shards::CellTopology CellTopology;
CellTopology hex_8(shards::getCellTopologyData >() );
int numNodesPerCell = hex_8.getNodeCount();
int numEdgesPerCell = hex_8.getEdgeCount();
int spaceDim = hex_8.getDimension();
Get Cubature Points and Weights
DefaultCubatureFactory cubFactory;
int cubDegree = 2;
Teuchos::RCP > hexCub = cubFactory.create(hex_8, cubDegree);
int cubDim = hexCub->getDimension();
int numCubPoints = hexCub->getNumPoints();
FieldContainer cubPoints(numCubPoints, cubDim);
FieldContainer cubWeights(numCubPoints);
hexCub->getCubature(cubPoints, cubWeights);
SAND2009-7358P
Computational mathematics and algorithms
Step 1: Define cell, cubature, and basis
Define Basis on reference cell
Basis_HCURL_HEX_I1_FEM > hexHCurlBasis;
Basis_HGRAD_HEX_C1_FEM > hexHGradBasis;
int numFieldsC = hexHCurlBasis.getCardinality();
int numFieldsG = hexHGradBasis.getCardinality();
Evaluate Basis at Cubature Points
FieldContainer hexGVals(numFieldsG, numCubPoints);
FieldContainer hexCVals(numFieldsC, numCubPoints, spaceDim);
FieldContainer hexCurls(numFieldsC, numCubPoints, spaceDim);
hexHGradBasis.getValues(hexGVals, cubPoints, OPERATOR_VALUE);
hexHCurlBasis.getValues(hexCVals, cubPoints, OPERATOR_VALUE);
hexHCurlBasis.getValues(hexCurls, cubPoints, OPERATOR_CURL);
SAND2009-7358P
Computational mathematics and algorithms
Step 2: Pullback data: Jacobians
Define physical cell workset // Intialize hexWorkset from mesh data
int numCells = 100;
FieldContainer hexWorkset(numCells,
…
numNodesPerCell,
Physical coordinates of cell nodes
spaceDim);
Compute Jacobians of cell workset
FieldContainer hexJacobian(numCells, numCubPoints, spaceDim, spaceDim);
FieldContainer hexJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
FieldContainer hexJacobDet(numCells, numCubPoints);
CellTools::setJacobian(hexJacobian, cubPoints, hexWorkset, hex_8); DF
CellTools::setJacobianInv(hexJacobInv, hexJacobian ); DF "1
CellTools::setJacobianDet(hexJacobDet, hexJacobian ); det DF
!
!
!
ˆ
" "
ˆ
F :" # "
! !
SAND2009-7358P
Computational mathematics and algorithms
Step 2: Pullback data: H(curl) basis signs
FieldContainer hexBasisSigns(numCells, numFieldsC);
Local (reference) cell
ˆ
u0 ˆ
u1 ˆ
u2 ˆ
u3
2
3 1
0
! ! ! !
Orientations Basis functions
- u1 " = +#* (u1 )
ˆ u3 " = #$* (u3 )
ˆ
L R
- +
! !
+
ui " = # i(" )$* (ui(" ) )
ˆ
Global (physical) cell
SAND2009-7358P
Computational mathematics and algorithms
Step 3: Assembly
H(grad) mass matrix
numPt
(MG ) ij ˆ ˆ *
= $ " i (x)" j (x)dx = $ (% " i )(x)(%*" j )(x)dx = ˆ ˆ ˆ ˆ ˆ ˆ
$ "i ( x )" j ( x )J( x )dx & ˆ ˆ
(" ( x ˆ ˆ ˆ
)" j ( x p )J( x p )' p
i p
# # ˆ
# p=1
H(curl) mass matrix
!
(MC ) ij = $ ui (x) " u j (x)dx = $ (%*ui )(x) " (%*u j )(x)dx
ˆ ˆ
# #
numPt
= $ & & (DF
i j
'T
ui ( x )) " ( DF u j ( x )) J( x )dx (
ˆ ˆ ˆ ˆ'T
ˆ ˆ *& & (DF
i j
'T
ui ( x p )) " ( DF 'T u j ( x p )) J( x p )) p
ˆ ˆ ˆ ˆ ˆ
ˆ
# p=1
H(curl) stiffness matrix
!
ˆ ˆ ˆ ˆ
(K C ) ij = & " # ui (x) $ " # u j (x)dx = & ('* (" # ui ))(x) $ ('* (" # u j ))(x)dx
% %
& ( ( (J ˆ ˆ ˆ ˆ ˆ ˆ
=
ˆ
i j
)1
)(
DF(" # ui )( x ) $ J )1DF(" # u j )( x ) J( x )dx
ˆ ˆ )
%
numPt
* ,( ( ( J )1 ˆ ˆ ˆ )( ˆ ˆ ˆ
DF(" # ui )( x p ) $ J )1DF(" # u j )( x p ) J( x p )+ p
ˆ )
i j
p=1
SAND2009-7358P
Computational mathematics and algorithms
Step 3: Assembly
Right hand side
(b) i = & " # ui (x) $ f(x)dx + & "* $ ui (x)g(x)dx
% %
= & " # u (x) $ f(x)dx ' & u (x) $ "g(x)dx + & g(x)u (x) $ nd(
i i i
% % (
ˆ
= & ) (J
i
'1
DF(" # ui )(x)) $ f(x)J(x)dx ' & ) i DF 'T ui (x) $ "g(x)J(x)dx +
ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ & g(x)) DF
ˆ i
'T
ˆ ˆ
ui (x) $ nd(
ˆ
% ˆ
% ˆ
(
numPt numPt
" )# ( J i
$1
DF(% & ui )(x p )) ' (f(x p )) J(x p )( p $
ˆ ˆ ˆ ˆ )# (DF i
$T
(ui )(x p )) ' (%g(x p )) J(x p )( p
ˆ ˆ ˆ ˆ
p=1 p=1
nBndF numFPt
! + " "# (DF i
$T
ui (x fp )) % n f ( g(x fp )) J(x fp )& fp
ˆ ˆ ˆ ˆ
f =1 fp=1
!
!
SAND2009-7358P
Computational mathematics and algorithms
Step 3: Declare arrays
// Containers for element HGRAD mass matrix
FieldContainer massMatrixG(numCells, numFieldsG, numFieldsG);
FieldContainer weightedMeasure(numCells, numCubPoints);
FieldContainer weightedMeasureMuInv(numCells, numCubPoints);
FieldContainer hexGValsTrans(numCells, numFieldsG, numCubPoints);
FieldContainer hexGValsTransWeighted(numCells, numFieldsG, numCubPoints);
// Containers for element HCURL mass matrix
FieldContainer massMatrixC(numCells, numFieldsC, numFieldsC);
FieldContainer hexCValsTrans(numCells, numFieldsC, numCubPoints, spaceDim);
FieldContainer hexCValsTransWeighted(numCells, numFieldsC, numCubPoints, spaceDim);
// Containers for element HCURL stiffness matrix
FieldContainer stiffMatrixC(numCells, numFieldsC, numFieldsC);
FieldContainer weightedMeasureMu(numCells, numCubPoints);
FieldContainer hexCurlsTrans(numCells, numFieldsC, numCubPoints, spaceDim);
FieldContainer hexCurlsTransWeighted(numCells, numFieldsC, numCubPoints, spaceDim);
// Global arrays in Epetra format
Epetra_FECrsMatrix MassG(Copy, globalMapG, numFieldsG);
Epetra_FECrsMatrix MassC(Copy, globalMapC, numFieldsC);
Epetra_FECrsMatrix StiffC(Copy, globalMapC, numFieldsC);
Epetra_FEVector rhsC(globalMapC);
SAND2009-7358P
Computational mathematics and algorithms
Step 3: Assemble H(grad) mass matrix
fst::HGRADtransformVALUE(hexGValsTrans, ˆ ˆ
"i ( x p ) Transform to
Note that for H(Grad) basis values hexGVals); ˆ ˆ
" (x )
i p
physical
transformation is formal - simply coordinates
replicates values in rank-3 array
!
fst::computeCellMeasure(cellMeasure, ˆ
J( x p )" p
! Compute
hexJacobDet, ˆ
J( x p ) cell
cubWeights); "p measure
!
!
fst::multiplyMeasure(hexGValsTransWeighted, ˆ ˆ ˆ
" j ( x p )J( x p )# p
Multiply basis
cellMeasure, ! ˆ
J( x p )" p and cell
hexGValsTrans); ˆ ˆ
"i ( x p ) measure
!
! numPt
fst::integrate(massMatrixG, ˆ ˆ
$" ( x ˆ ˆ ˆ
)" j ( x p )J( x p )# p Integrate to
i p
!
hexGValsTrans,
p=1
ˆ ˆ compute
"i ( x p )
element mass
ˆ ˆ ˆ
" j ( x p )J( x p )# p
hexGValsTransWeighted, matrix
COMP_BLAS); !
numPt
!
ˆ ˆ ˆ ˆ
(MG ) ij " % # i ( x p )# j ( x!)J( x p )$ p
ˆ
p
p=1
SAND2009-7358P
Computational mathematics and algorithms
Step 3: Assemble H(curl) mass matrix
fst::HCURLtransformVALUE(hexCValsTrans, DF "T ( x p )ui ( x p )
ˆ ˆ ˆ
hexJacobInv, DF "T ( x p )
ˆ
hexCVals); ˆ ˆ
ui ( x p )
!
fst::multiplyMeasure(hexCValsTransWeighted,
! DF "T ( x p )ui ( x p )J( x p )# p
ˆ ˆ ˆ ˆ
reuse ! cellMeasure, ! ˆ
J( x p )" p
hexCValsTrans); DF "T ( x p )ui ( x p )
ˆ ˆ ˆ
!
numPt
fst::integrate(massMatrixC, ! % (DF "T
ui ( x p )) # ( DF "T u j ( x p )) J( x p )$ p
ˆ ˆ ˆ ˆ ˆ
p=1
!
hexCValsTrans, DF "T ( x p )ui ( x p )
ˆ ˆ ˆ
hexCValsTransWeighted, DF "T ( x p )ui ( x p )J( x p )# p
ˆ ˆ ˆ ˆ
COMP_BLAS); !
! numPt
fst::applyLeftFieldSigns(massMatrixC,
! &" (DF i
#T
ui ( x p )) $ ( DF #T u j ( x p )) J( x p )% p
ˆ ˆ ˆ ˆ ˆ
p=1
hexBasisSigns); "i
numPt
fst::applyRightFieldSigns(massMatrixC,
!
&" " (DFi j
#T
ui ( x p )) $ ( DF #T u j ( x p )) J( x p )% p
ˆ ˆ ˆ ˆ ˆ
p=1
hexBasisSigns); ! "j
numPt
(MC ) ij " '# i# j (DF $T ui ( x p )) % (DF! u j ( x p )) J( x p )& p
ˆ ˆ $T
ˆ ˆ ˆ
p=1
!
SAND2009-7358P
Computational mathematics and algorithms
Step 3: Assemble H(curl) stiffness matrix
fst::HCURLtransformCURL(hexCurlsTrans, ˆ ˆ ˆ ˆ
J "1 ( x p )DF( x p )(# $ u j )( x p )
ˆ
hexJacob, ˆ
DF( x p )
hexJacobianDet, ˆ
J( x p )
hexCVals); ! ˆ ˆ ˆ
(" # u j )( x p )
!
fst::multiplyMeasure(hexCurlsTransWght, ˆ ˆ ˆ ˆ
J "1 ( x p )DF( x p )(# $ u j )( x p ) (J( x p )% p )
ˆ ˆ
!
reuse ! cellMeasure, ! ˆ
J( x p )" p
hexCurlsTrans); ˆ ˆ ˆ ˆ
J "1 ( x p )DF( x p )(# $ u j ( x p ))
ˆ
!
numPt
fst::integrate(stiffMatrixC, ' (J ˆ ˆ ˆ ˆ ˆ ˆ
!
p=1
"1
)( )
DF(# $ ui )( x p ) % J "1DF(# $ u j )( x p ) J( x p )& p
ˆ
hexCurlsTrans,
! "1
ˆ ˆ ˆ ˆ ˆ
J ( x p )DF( x p )(# $ u j ( x p ))
hexCurlsTransWght, ˆ ˆ ˆ ˆ
J "1 ( x p )DF( x p )(# $ u j )( x p ) (J( x p )% p )
ˆ ˆ
COMP_BLAS); !
! numPt
(" ( J ˆ ˆ ˆ ˆ ˆ ˆ
fst::applyLeftFieldSigns(stiffMatrixC,
! p=1
i
#1
)( )
DF($ % ui )( x p ) & J #1DF($ % u j )( x p ) J( x p )' p
ˆ
hexBasisSigns); "i
numPt
(" " ( J ˆ ˆ ˆ ˆ ˆ ˆ
fst::applyRightFieldSigns(stiffMatrixC,!
p=1
i j
#1
)( )
DF($ % ui )( x p ) & J #1DF($ % u j )( x p ) J( x p )' p
ˆ
hexBasisSigns);
! "j
numPt
ˆ ˆ ! ˆ ˆ ˆ
(K C ) ij " )# i# j (J $1DF(% & ui )( x p )) ' (J $1DF(% & u j )( x p ))J( x p )( p
ˆ ˆ
p=1
!
SAND2009-7358P
Computational mathematics and algorithms
Step 3: Assemble right hand side
(b) i = & " # ui (x) $ f(x)dx ' & ui (x) $ "g(x)dx + & g(x)ui (x) $ nd(
% % (
numPt
fst::integrate(rhsf, ˆ ˆ ˆ
' (J "1
)
DF(# $ ui )( x p ) % (f( x p )) J( x p )& p
ˆ ˆ
fData, p=1
! ˆ
f(x p )
reuse ! hexCurlsTransWght,
ˆ ˆ ˆ ˆ
J "1 ( x p )DF( x p )(# $ u j )( x p ) (J( x p )% p )
ˆ ˆ
COMP_BLAS); !
!
numPt
fst::applyFieldSigns(rhsfData, ˆ ˆ ˆ
! (" ( J i
#1
)
DF($ % ui )( x p ) & (f( x p )) J( x p )' p
ˆ ˆ
p=1
hexBasisSigns);
"i
fst::integrate(rhsg, ! numPt
!
& (DF "T
ui ( x p )) # ($g( x p )) J( x p )% p
ˆ ˆ ˆ ˆ
gData, p=1
ˆ
"g(x p )
reuse ! hexCValsTransWght,
COMP_BLAS); DF "T ( x p )ui ( x p )J( x p )# p
ˆ ˆ ˆ ˆ
!
! numPt
fst::applyFieldSigns(rhsgData, ! '" (DF i
#T
ui ( x p )) $ (%g( x p )) J( x p )& p
ˆ ˆ ˆ ˆ
p=1
hexBasisSigns); "i
numPt numPt
(b) i " )# i ( J $1
DF(% & ui )(x p )) ' (f(x p )) J(x p )( p $
ˆ ˆ ˆ ˆ )# (DF
! i
$1
(ui )(x p )) ' (%g(x p )) J(x p )( p + B.T.
ˆ ˆ ˆ ˆ
p=1 p=1
!
SAND2009-7358P
Computational mathematics and algorithms
Step 3: Assemble right hand side
(b) i = & " # ui (x) $ f(x)dx ' & ui (x) $ "g(x)dx + & g(x)ui (x) $ nd(
% % (
For elements that have faces on the boundary:
!
CellTopology quad_4(shards::getCellTopologyData >() );
Teuchos::RCP > quadCub = cubFactory.create(quad_4, cubDegree);
quadCub->getCubature(quadFacePoints, faceWeights);
CellTools::mapToReferenceSubcell(refFacePoints,
quadFacePoints,
2, iface, hex_8);
hexHCurlBasis.getValues(faceCVals, refFacePoints, OPERATOR_VALUE); ˆ ˆ
u(x fp )
CellTools::setJacobian(faceJacobians,refFacePoints, ˆ
DF( x fp )
hexNodes, hex_8); !
CellTools::setJacobianInv(faceJacobInv, faceJacobians ); DF "1 ( x fp )
ˆ
!
fst::HCURLtransformVALUE(hexCValsTrans, DF "T (x fp )ui (x fp )
ˆ ˆ ˆ
hexJacobInv,hexCVals);
!
!
SAND2009-7358P
Computational mathematics and algorithms
Step 3: Assemble right hand side
(b) i = & " # ui (x) $ f(x)dx ' & ui (x) $ "g(x)dx + & g(x)ui (x) $ nd(
% % (
CellTools::getPhysicalFaceNormals(hexFaceN, nf
! faceJacobians, ˆ
DF( x fp )
iface, hex_8);
Compute the dot product and multiply by weights ! (DF "T ˆ ˆ
ui (x fp )) # n f J(x fp )$ fp
ˆ
!
numFPt
fst::integrate(rhsgBoundary,
% (DF "T
ui ( x fp )) # n f J( x fp )$ fp ( g( x fp ))
ˆ ˆ ˆ ˆ
gBoundaryData, !
fp=1
ˆ
g(x fp )
hexCValsDotNorm,
COMP_BLAS);
(DF "T
ui (x fp )) # n f J(x fp )$ fp
ˆ ˆ ˆ
!
! numFPt
fst::applyFieldSigns(rhsgBoundary, &" (DF #T
ui ( x fp )) $ n f J( x fp )% fp ( g( x fp ))
ˆ ˆ ˆ ˆ
! i
fp=1
hexBasisSigns); "i
numPt numPt
(b) i " )# i ( J DF(% & ui )(x p )) ' (f(x p ))J(x p )( p $ ! # i (DF $1(ui )(x p )) ' (%g(x p )) J(x p )( p
$1
ˆ ˆ ˆ ˆ ) ˆ ˆ ˆ ˆ
p=1 !
p=1
nBndF numFPt
+ " "# (DF i
$T
ui (x fp )) % n f ( g(x fp )) J(x fp )& fp
ˆ ˆ ˆ ˆ
f =1 fp=1
!
!
SAND2009-7358P
Computational mathematics and algorithms
Step 4: Solve
Mesh size
Error Order
10x10x10 20x20x20 40x40x40
L2 0.31318 0.13991 0.067749 1.1044
H(curl) 2.3036 1.1575 0.57905 0.9961
SAND2009-7358P
Computational mathematics and algorithms
Step 4: Solve
2D (x-y) projection of the 3D model
distorted mesh problem at the 3D Hex Mesh generated by PAMGEN (Trilinos package)
coarsest grid resolution.
Mesh size
Error Order
10x10x4 20x20x8 40x40x16
L2 0.89131 0.420695 0.212596 0.98466
H(curl) 4.89925 2.73187 1.55896 0.80930
L2 is optimal but H(curl) - suboptimal for non-affine Hex; see Falk et al. 2009