A hybrid mimetic spectral element method for three-dimensional linear elasticity problems

Delft University of Technology

A hybrid mimetic spectral element method for three-dimensional linear elasticity problems

Zhang, Yi; Fisser, Joël; Gerritsma, Marc

DOI

10.1016/j.jcp.2021.110179

Publication date

2021

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

Zhang, Y., Fisser, J., & Gerritsma, M. (2021). A hybrid mimetic spectral element method for

three-dimensional linear elasticity problems. Journal of Computational Physics, 433, [110179].

https://doi.org/10.1016/j.jcp.2021.110179

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

A

hybrid

mimetic

spectral

element

method

for

three-dimensional

linear

elasticity

problems

Yi Zhang

∗

, Joël Fisser,

Marc Gerritsma

DelftUniversityofTechnology,FacultyofAerospaceEngineering,Kluyverweg1,2629HSDelft,theNetherlands

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Availableonline4February2021

Keywords:

Mimeticspectralelementmethod Hybridization

Domaindecomposition Variationalprinciple Lagrangemultiplier DeRhamcomplex

We introduce adomain decomposition structure-preserving methodbased on ahybrid mimetic spectral element method for three-dimensional linear elasticity problems in curvilinearconforming structuredmeshes. Themethodis anequilibrium methodwhich satisfiespointwiseequilibriumofforces.Thedomaindecompositionisestablishedthrough hybridizationwhichfirstallowsforaninter-elementnormalstressdiscontinuityandthen enforcesthenormalstress continuityusingaLagrangemultiplierwhichturnsouttobe the displacementin the tracespace. Dualbasisfunctions are employedto simplify the discretization andtoobtainahighersparsity.Numericaltestssupportingthemethodare presented.

©2021TheAuthor(s).PublishedbyElsevierInc.Thisisanopenaccessarticleunderthe CCBYlicense(http://creativecommons.org/licenses/by/4.0/).

1. Introduction

Structure-preserving ormimeticmethods arenumericalmethodsthataimtopreservefundamentalpropertiesofthe con-tinuousproblems,likeconservationlaws(forexampleequilibriumofforcesinelasticity),atthediscretelevel.Themimetic spectral elementmethod(MSEM) [1] is anovel arbitraryordermimeticmethodthat usually usesthemixedformulation. Ithasbeenappliedto,forinstance,Stokesflow[2],thePoissonequation[3],theGrad-Shafranovequation[4],theshallow waterequations[5],and,recently,linearelasticity[6].Themixedformulation,together withhighordermethods,generally leads tolargeanddense matrices.Thisisparticularlythecasewhenoneconsidersthethree-dimensionalmixedelasticity formulationwhichsolvesforthreephysicalquantities,namelydisplacement,rotationandstress,simultaneouslyinaglobal discretesystem.Displacementandrotationfieldsconsistofthreecomponentswhilethestresstensorfieldhasnine compo-nents.ThismakestheMSEManexpensivemethodforthree-dimensionallinearelasticity.Aneffectivewaytoovercomethis drawbackistousethehybridfiniteelementmethod[7,8],adomaindecompositionmethodwhichbreaksuptheproblem intosmallersub-problems.

Withintheworldofmimeticmethodswedistinguishbetweenvariousnumericalmethods.Aparticularbranchofthese methodsisthe mimeticfinitedifferencemethods,see[9] anditsreferences.Thevirtualelementmethod[10,11] isanother wayin whichmimetic ideas are represented.Inyet another branchthe mathematicallanguage of differentialforms and similarities between differential forms and algebraic topology is exploited, [12,13]. A pioneering work of implementing theseideasinnumericalanalysiswasfirstconductedbyBossavit[14] incomputationalelectromagnetism.Later,acommon framework of using differential forms and algebraic topology for mimetic discretizations was developed by Bochev and Hyman[15].SomeimportantcontributionsinthecontextoffiniteelementdiscretizationswerethenmadebyArnold,Falk

*

Correspondingauthor.

E-mailaddress:zhangyi_aero@hotmail.com(Y. Zhang).

https://doi.org/10.1016/j.jcp.2021.110179

0021-9991/©2021TheAuthor(s).PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/).

andWinther,[16,17].TheextensiontospectralelementmethodsisgivenbytheMSEM.Initially,theMSEMwasintroduced usingthemathematicallanguageofdifferentialforms,whichispreferablebut lesspopularthanconventionalvector/tensor calculusforthemathematicaldescriptionofphysics[18].Itispossibletotranslatethemathematicallanguageofdifferential formsintovector/tensorcalculustoreachalargeraudience.Theworkin,e.g.,[4,6] andthepresentpaperareexamples.

Insolidmechanics,finiteelementmethodsaredevelopedbasedonvariationalprinciples.Inclassicfiniteelement meth-ods,stiffnessmatricesrepresentingthevariationalprincipleinelementsareestablishedandassembled.Elementsarethen coupled intheglobalstiffness matrix throughthe stronginter-element continuity. Unlikeclassicfinite elementmethods, hybridfiniteelementmethodsfirstallowforaninter-element discontinuityandthenimposethecontinuitywiththehelp ofaLagrangemultiplier[7].Thisprocessisusuallycalledhybridization.Thefirsthybridfiniteelementmethod,theassumed stresshybridmethod[19],wasproposedbyPianinthe1960s.Otherhybridfiniteelementmethodsinsolidmechanicsare, for example,the assumeddisplacement hybrid method[20] andthe assumedstress-displacement hybrid mixedmethod [21].The advantageoftheadditionoftheextracontinuityequationwiththeassociatedLagrange multiplieristhatitmay be possibleto solve fortheinterface Lagrange multiplier first,afterwhich theglobalsystemdecouples into independent problems atelement level.Thisis particularlyefficient forspectral element methodswhere thenumber ofinterface un-knowns isrelatively smallcompared tothe globalnumberofunknowns. Themortar element method[22–24], a domain decompositionmethodwhichcouplesdifferentnon-overlappingsubdomainsusesasimilaridea.The ideaofhybridization alsoplays an importantrole inthefinite elementtearingandinterconnecting (FETI)method[25,26]. Forlinearelasticity, oneofthemainproblemsinhybridfiniteelementmethodsistheappearanceofso-calledspuriouskinematicmodes orzero energymodes,see,forinstance,[27–29].Thesespuriousmodesconsist ofnon-solidbody deformationswhichdonotaffect thestressfieldindicatingthatsuchhybridformulationsarenon-wellposed.In[30],thewell-posednessofproblemsarising fromthehybridvariationalprinciplesandtheerrorbehaviorofthehybridmethodarestudied.Hybridizingcertainexisting mixedfiniteelementmethodsisstudiedin[31].

In this work, we introduce a hybrid mimetic spectral element method (hMSEM) based on a new hybrid variational principle for three-dimensional linearelasticity problems.The method is an arbitraryorder structure-preserving method which satisfies pointwise equilibriumof forces and, to our knowledge, is the first method that manages to reduce the computational cost of the MSEM without the introductionof spurious kinematic modes. Besides the computational cost reductionasaresultofthehybridization,additionalreductionisobtainedthroughtheuseofdualbasis functions[32–34] which significantlyincreasesthe sparsity ofthe discretesystem. Dual basis functionscanalso beapplied to theoriginal, non-hybrid,MSEMtoimproveitsefficiency.

The layoutofthepaperisasfollows: InSection2,abriefintroductiontothethree-dimensionallinearelasticity prob-lem isgiven, whichisfollowed byan introduction ofahybrid mixedformulation andits weak formulationinSection 3. The proposed methodisthen explainedindetailin Section4.Numerical experimentsare presentedinSection 5.Finally, conclusionsaredrawninSection6.

2. Three-dimensionallinearelasticity

In

R

3_{withcoordinatesystem}

_{

_x

_,

_y

_,

_z

_}

_,letu

₌



_u

x uy uz



Tbethedisplacementvector.Therotationvector

ω

isgiven

by

ω

=

⎧

⎨

⎩

ω

x

ωy

ωz

⎫

⎬

⎭

=

1 2

⎧

⎨

⎩

0

−∂/∂

z

∂/∂

y

∂/∂

z 0

−∂/∂

x

−∂/∂

y

∂/∂

x 0

⎫

⎬

⎭

⎧

⎨

⎩

ux uy uz

⎫

⎬

⎭

.

(1)

WeuseD torepresentthedivergencematrix,

D

=

⎧

⎨

⎩

∂/∂

x

∂/∂

y

∂/∂

z 0 0 0 0 0 0 0 0 0

∂/∂

x

∂/∂

y

∂/∂

z 0 0 0 0 0 0 0 0 0

∂/∂

x

∂/∂

y

∂/∂

z

⎫

⎬

⎭

.

Itstranspose,DT_{,thenisthegradientmatrix.Thestrainvector}

_ε

₌



_ε

xx

ε

yx

ε

zx

ε

xy

ε

y y

ε

zy

ε

xz

ε

yz

ε

zz



T

canbe writtenas

ε

=

DT_u

₊

_RT

_ω

_{,whereR isamatrixgivenby}

R

=

⎧

⎨

⎩

0 0 0 0 0 1 0

−

1 0 0 0

−

1 0 0 0 1 0 0 0 1 0

−

1 0 0 0 0 0

⎫

⎬

⎭

.

The stress vector

σ

=



σ

xx

σ

yx

σ

zx

σ

xy

σ

y y

σ

zy

σ

xz

σ

yz

σ

zz



T

can becomputed usingtheconstitutive relation,

C

=

1 E

⎧

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎨

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎩

1 0 0 0

−

ν

0 0 0

−

ν

0 1

+

ν

0 0 0 0 0 0 0 0 0 1

+

ν

0 0 0 0 0 0 0 0 0 1

+

ν

0 0 0 0 0

−

ν

0 0 0 1 0 0 0

−

ν

0 0 0 0 0 1

+

ν

0 0 0 0 0 0 0 0 0 1

+

ν

0 0 0 0 0 0 0 0 0 1

+

ν

0

−

ν

0 0 0

−

ν

0 0 0 1

⎫

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎬

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎭

,

isthecompliancetensor. E and

ν

representYoung’smodulus andPoisson’sratioofthematerial,respectively.Equilibrium offorcesstatesthat D

σ

+

f

=

0,where f

=



fx fy fz



T

isthebodyforcevector, andequilibriumofmomentsimplies R

σ

=

0.Ifweputthestresscomponentsina3 by 3 tensor,equilibriumofmomentsthenimpliesthat thestresstensoris symmetric.

We now consider a bounded, connected domain



with boundary

∂

= 

u

∪ 

t, where

u

∩ 

t

= ∅

,

u

= ∅

. On

u

,the displacement u isprescribed; u

|

u

= 

u

=





ux



uy



uz



T

.On

t

,the surface tractiont is prescribed; t

|

t

=

t

=



tx



ty



tz



T

.Thethree-dimensionallinearelasticityproblemthencanbeformulatedas

C

σ

−

DT_u

₋

_RT

_ω

₌

_{0 in}

_,

_(2a)

D

σ

+

f

=

0 in

,

(2b)

−

R

σ

=

0 in

,

(2c)

u

= 

u on



u

,

(2d)

t

=

N

σ

=

t on



t

,

(2e)

wherethebodyforce f isknown,andN isamatrixgivenby

N

=

⎧

⎨

⎩

nx ny nz 0 0 0 0 0 0 0 0 0 nx ny nz 0 0 0 0 0 0 0 0 0 nx ny nz

⎫

⎬

⎭

,

wherenx,ny,nz arecomponentsoftheunitoutwardnormalvector,n

=



nx ny nz



T

.Itisstraightforwardtoprovethat thesolutionsu and

ω

ofproblem(2) satisfyrelation(1).

3. Ahybridmixedformulation 3.1. Notations

Throughoutthepaper,we restrictourselvestotheHilbertspaces L2

()

, H1

()

, H

(

div

;

)

,trace spaces H1/2

_(∂

_)

_and



H1/2

(∂)

[35],andtheirextensionsto[

·

]3 in

R

3_{.Thedualspacesareexpressedwiththenotation}

₍

_·)

_.Forexample,



H1/2

(∂)



_

=

H−1/2

(∂).

Thesubspace H0

(

div

;



;

)

isdefinedas

H0

(

div

; ; ) := {

ϕ

|

ϕ

∈

H

(

div

; ),

tr

ϕ

=

N

ϕ

=

0 on



} .

Withoutlossofgenerality,inthissection,wewillassumetheNeumannboundaryconditiontobe



t



_

t

=

0. 3.2. Amixedformulation

Given f

∈



L2

()



3 and



u

∈



H1/2

(u)



3,for

(

σ

,

u

,

ω

)

∈

[H0

(

div

; ; 

t)]3

×



L2

()



3

×



L2

()



3,avariational formula-tionbasedontheprincipleofminimizingthecomplementaryenergy,(i),subjecttoconstraintsofequilibriumofforces,(ii), andequilibriumofmoments,(iii),iswrittenas,[6],

L

(

σ

,

u

,

ω

;

f

,



u

)

=

1 2

(

σ

,

C

σ

)



− 

u

,

t

u







(i)

+ (

_

u

,

D

σ

_

+

f

)



_

(ii)

− (

_{  }

ω

,

R

σ

)

 (iii)

,

(3)

where

(

·, ·)

 denotes the L2-inner product and the duality pairing,

·, ·

u, between



u

∈



H1/2

(u)



3 and t

=

tr

σ

∈



H1/2

_(u)





3 _{stands for a boundary integral. Its stationary points weakly solve the problem (2), andthe} correspond-ing mixedweak formulationis: Given f

∈



L2

()



3 and



u

∈



H1/2

(u)



3,find

(

σ

,

u

,

ω

)

∈

[H0

(

div

; ; 

t)]3

×



L2

()



3

×



L2

_()



3_suchthat



σ

,

C

σ

_

+



u

,

D

σ

_

−



ω

,

R

σ

_

=





u

,

t



_ u

,

∀

σ

∈

[H0

(

div

; ; 

t

)

] 3

_,

_(4a)



u

,

D

σ

_

= −



u

,

f

_

,

∀

u

∈



L2

()



3

,

(4b)

−



ω

,

R

σ

_

=

0

,

∀

ω

∈



L2

()



3

.

(4c)

3.3. Thehybridformulation

Toconstructa hybridmixedformulation, we firstlet amesh,denoted by



h,partition thedomain



into M disjoint subdomains,



m,m

=

1

,

2

,

· · · ,

M,anduse



i j todenotetheinterfacebetweensubdomains



i and



j.

Given f

∈



L2

(

m

)



3 and



u

∈



H1/2

_(∂

m

∩ 

u)



3 ,for

(

σ

,

u

,

ω

, λ)

∈

[H0

(

div

; 

m

; ∂

m

∩ 

t)]3

×



L2

(

m

)



3

×



L2

(

m

)



3

×



H1/2

(

i j

)



3

,ahybridvariationalformulationisexpressedas

L

(

σ

,

u

,

ω

;

f

,



u

)

=



m



1 2

(

σ

,

C

σ

)

m

− 

u

,

t

∂m∩u

+ (

u

,

D

σ

+

f

)

m

− (

ω

,

R

σ

)

m



−



i j



λ,

ti

+

tj



i j

,

(5)

wheretiandtj aresurfacetractions ofsubdomains



i and



j(



i

∩ 

j

= ∅

),i.e.,ti

=

Ni

σ

i,tj

=

Nj

σ

j,and,on



i j,Ni

=

−

Nj.Acrosstheinterfaceofsubdomains,thesurfacetractionstiandtjcandifferfromeachother.Toenforcethecontinuity,

weintroducetheLagrangemultiplier

λ

∈



H1/2

_(

i j

)



3

andaddthesurfaceintegralconstrainttothevariationalformulation. Ifweperformvariationalanalysisof(5) withrespectto

σ

inaparticularsubdomain,forexample,thesubdomain



mwho

isaneighborofsubdomains



n (



mn

= ∅

),wehave



σ

,

C

σ

_ m

−





u

,

t



_∂ m∩u

+



u

,

D

σ

_ m

−



ω

,

R

σ

_ m

−



n



λ,

t



_ mn

=

0

.

Ifthesolutionissufficientlysmooth,wecanperformintegrationbypartstothethirdterm,



u

,

D

σ

_ m

= −



DT_u

_,

_σ

m

+



u

,

t



∂m

,

andusethefact



ω

,

R

σ

_ m

=



σ

,

RT

_ω

m,wewouldobtain



σ

,

C

σ

−

DTu

−

RT

ω

_ m

+



u

− 

u

,

t



_∂ m∩u

+



n



u

− λ,

t



_ mn

=

0

,

∀

σ

∈

[H0

(

div

; 

m

; ∂

m

∩ 

t

)

] 3

_.

Thisimplies C

σ

−

DT_u

₋

_RT

_ω

₌

_{0 in}

_

m

,

and u

= 

u on

∂



m

∩ 

u

,

u

= λ

on

∂



m

\ ∂.

So theLagrangemultiplier that enforcesthe continuityacrossthesubdomainshasa physicalinterpretation;itrepresents thedisplacementalongtheinterface.SimilardiscussionsaboutthattheLagrangemultiplierinthetracespaceofthehybrid methodhasaphysicalmeaningcanbefoundin[31].Ifwe performvariationalanalysiswithrespectto

σ

,u,and

ω

inall subdomains



mandwithrespectto

λ

onallinterfaces



i j,wewillfindthatthestationarypointsof(5) weaklysolvesthe

followingproblem,

C

σ

−

DTu

−

RT

ω

=

0 in



m

,

(6a)

D

σ

+

f

=

0 in



m

,

(6b)

−

R

σ

=

0 in



m

,

(6c)

u

= 

u on



u

,

(6e)

t

=

t on



t

,

(6f)

ti

+

tj

=

0 on



i j

.

(6g)

Hybrid problem (6) is equivalent to the linear elasticity problem(2). The weak formulation derived from (5) is written as: For all subdomains



m and all interfaces



i j, given f

∈



L2

(

m

)



3 and



u

∈



H1/2

_(∂

m

∩ 

u)



3 , find

(

σ

,

u

,

ω

, λ)

∈

[H0

(

div

; 

m

; ∂

m

∩ 

t)]3

×



L2

_(

m

)



3

×



L2

_(

m

)



3

×



H1/2

_(

i j

)



3 suchthat



σ

,

C

σ

_ m

+



u

,

D

σ

_ m

−



ω

,

R

σ

_ m

−



n



λ,

t



_ mn

=





u

,

t



_∂ m∩u

,

∀

σ

∈

[H0

(

div

; 

m

; ∂

m

∩ 

t

)

] 3

_,

_(7a)



u

,

D

σ

_ m

= −



u

,

f

_ m

,

∀

u

∈



L2

(

m

)



3

,

(7b)

−



ω

,

R

σ

_ m

=

0

,

∀

ω

∈



L2

(

m

)



3

,

(7c)

−



λ,

ti

+

tj



i j

=

0

,

∀ λ ∈



H1/2

(

i j

)



3

.

(7d)

Wecall(7) thehybridmixedweakformulation. 4. Numericalmethod

Theexactsequence- thedeRhamcomplex[16,17,34],

R

→

H1

()

−→

grad H

(

curl

; )

−→

curl H

(

div

; )

−→

div L2

()

→

0

,

(8) isofessentialimportanceforstructure-preservingmethods.Forexample,wehavechosen

σ

∈

[H

(

div

; )

]3and f

∈



L2

()



3 suchthattheequilibriumofforces,D

σ

+

f

=

0,canbeexactlysatisfiedfortheweakformulations.However,thisdoesnot guaranteethattheequilibriumofforcesissatisfiedatthediscretelevelunlessthefinitedimensionalfunctionspacesused forthediscretizationalsoformadeRhamcomplex.

Inthissection,wewillfirstintroducetheconstructionofthemimeticpolynomialspaces which, aswillbeseen,areable toformadiscretedeRhamcomplexineitherorthogonalorcurvilinearmeshes.ThesespaceshavebeenusedintheMSEM forvariousproblems[2–6].ForthehMSEM,byextendingtheoriginalmimeticpolynomialstotracespacesandintroducing the dualpolynomials, asecond deRham complexintermsofparticulardiscrete weakoperators willbe constructed.The twodiscretedeRhamcomplexesarethenappliedtothehybridlinearelasticityproblem.

4.1. Mimeticpolynomialspaces

The finite dimensional function spaces used in thispaper are the mimetic polynomial spaces spanned by either the primal polynomials or their dual representations (dual polynomials). The primal polynomials are constructed using the Lagrangepolynomialsandtheedgepolynomials[36].Thedualpolynomialsarebasically linearcombinationsoftheprimal polynomials.

4.1.1. Primalpolynomialsin

R

Forcompleteness,westartwiththewell-knownLagrangepolynomials.Given

(

N

+

1

)

nodes,

−

1

≤ ξ

0

< ξ

1

<

· · · < ξ

N

≤

1,

overintervalI

= [−

1

,

1

]

,the

(

N

+

1

)

Lagrangepolynomialsofdegree N aredefinedas li

(ξ )

=

N



j=0,j=i

ξ

− ξ

j

ξ

i

− ξ

j

,

i

∈ {

0

,

1

,

2

,

· · · ,

N

} ,

whichsatisfytheso-calledKroneckerdeltaproperty,

li



ξ

j

= δ

i,j

=



1 if i

=

j 0 else

.

(9)

ExamplesoftheLagrangepolynomialsareshowninFig.1(Left).Letph

(ξ )

beapolynomialofdegreeN,

ph

(ξ )

=

N



i=0

pili

(ξ ).

(10)

Fig. 1. Lagrangepolynomials(Left)andedgepolynomials(Right)derivedfromasetof5 nodes,−1= ξ0< ξ1<· · · < ξ4=1.Theverticalgraydashed

linesindicatetheinternalnodesξ1,ξ2,ξ3.ThenodalKroneckerdeltaproperty(9) isobvious.TheintegralKroneckerdeltaproperty(12) canbeseen,for

example,fromtheedgepolynomiale2(ξ )(orangesolidline).Directcalculationswillrevealthat ξ2

ξ1 e2(ξ )dξ=1 and

ξj

ξj−1e2(ξ )dξ=0, j∈ {1,3,4}.(For

interpretationofthecolorsinthefigure(s),thereaderisreferredtothewebversionofthisarticle.) dph

_{(ξ )}

d

ξ

=

N



i=0 pi dli

(ξ )

d

ξ

=

N



i=1

!

(

pi

−

pi−1

)

N



k=i dlk

(ξ )

d

ξ

"

=

N



i=1

!

(

pi−1

−

pi

)

i−1



k=0 dlk

(ξ )

d

ξ

"

,

(11)

wheretheN polynomialsofdegree

(

N

−

1

)

arethecorrespondingedgepolynomials[36] denotedbyei

(ξ )

,

ei

(ξ )

=

N



k=i dlk

(ξ )

d

ξ

= −

i−1



k=0 dlk

(ξ )

d

ξ

,

i

∈ {

1

,

2

,

· · · ,

N

} .

WithNewton-LeibnizintegralruleandthenodalKroneckerdeltaproperty(9),itiseasyto seethat theedgepolynomials satisfytheKroneckerdeltapropertyinanintegralsense,namely,

ξj

#

ξj−1 ei

(ξ )

d

ξ

= δ

i,j

=



1 if i

=

j 0 else

.

(12)

ExamplesofedgepolynomialsareshowninFig.1(Right).Now,wecanwrite(11) as qh

(ξ )

=

dp h

_{(ξ )}

d

ξ

=

N



i=1

(

pi

−

pi−1

)

ei

(ξ )

=

N



i=1 qiei

(ξ ).

(13)

Ifwecollecttheexpansioncoefficientsordegreesoffreedomofphandqh,weobtaintwovectors, p andq, p

=



p0 p1

· · ·

pN



T

,

q

=



q1 q2

· · ·

qN



T

.

(14)

Throughoutthepaper,underlinedquantitieswillrepresentthevectorsofexpansioncoefficients.From(13),wehave q

=

Edp

,

wherethelinearoperatorEd,

Ed

=

⎧

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎨

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎩

−

1 1 0

· · ·

0 0 0 0

−

1 1

· · ·

0 0 0 0 0

−

1

· · ·

0 0 0

..

.

..

.

..

.

. .

.

..

.

..

.

..

.

0 0 0

· · · −

1 1 0 0 0 0

· · ·

0

−

1 1

⎫

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎬

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎭

,

is calledthe incidencematrix whichis a discretecounterpart of thederivative operator. The incidencematrix isa sparse matrix(if N

>

1)andbecomessparserwhen N grows. Itisalso atopologicalmatrix.Forexample,assuming N doesnot change andthe degrees offreedom are labeled in a consistent way (the topology of the nodesdoesnot change), ifwe useadifferentsetofnodesormapthedomain I intoadifferentdomainusingacontinuousmapping,thebasisfunctions change,but theincidencematrix remains thesame.One point toemphasize isthat (13) isexact, whichmeans that the discretizationofthederivativeoperatorwiththeincidencematrixEd isexact[37,38].

TheLagrangepolynomialsandedgepolynomialsaretheprimalpolynomialsin

R

.Letthefinitedimensionalpolynomial spaces spannedby theLagrange polynomialsandthe edgepolynomialsbe denoted by

L

N and

E

(N−1) respectively.From

(10) and (13),we knowthat therangeofthederivativeoperator on

L

N isin

E

(N−1).Therefore,we canconcludethat

L

N and

E

(N−1)formthefollowingdeRhamcomplex,

L

N Ed

L

N ⊂ d H1

₍

_I

₎

d

E

(N−1)

E

(N−1) ⊂ L2

(

I

)

Here theunderlined spaces

L

N and

E

(N−1) denotethespacesof expansioncoefficient vectorsoftheelements in

L

N and

E

(N−1)respectively.Thisconvention,inadditiontotheconventionin(14),isalsousedthroughoutthepaper.

4.1.2. Primalpolynomialsin

R

3

The primal polynomialsin

R

3 _{are constructed withthe primal polynomials (the Lagrange polynomialsand theedge} polynomials)in

R

usingthetensorproduct.Weconsiderthereferenceelement

(ξ,

η

,

ς

)

∈ 

ref

= [−

1

,

1

]

3 andthreesetsof nodes,

−

1

≤ ξ

0

< ξ

1

<

· · · < ξ

Nξ

≤

1,

−

1

≤

η

0

<

η

1

<

· · · <

η

Nη

≤

1,and

−

1

≤

ς

0

<

ς

1

<

· · · <

ς

Nς

≤

1.The tensorproduct

oftheprimalpolynomialsin

R

givesprimalpolynomialsin

R

3thatspanthefollowingprimalpolynomialspaces, P

:= L

Nξ

⊗ L

Nη

⊗ L

Nς

_,

E

:= E

(Nξ−1)

⊗ L

Nη

⊗ L

Nς

× L

Nξ

⊗ E

(Nη−1)

⊗ L

Nς

× L

Nξ

⊗ L

Nη

⊗ E

(Nς−1)

_,

S

:= L

Nξ

⊗ E

(Nη−1)

⊗ E

(Nς−1)

× E

(Nξ−1)

⊗ L

Nη

⊗ E

(Nς−1)

× E

(Nξ−1)

⊗ E

(Nη−1)

⊗ L

Nς

_,

V

:= E

(Nξ−1)

⊗ E

(Nη−1)

⊗ E

(Nς−1)

_,

wherethenotationsofthespaces, P ,E,S,andV ,standforpoints(nodes),edges,surfaces,andvolumes.Weusethisnotation because the degrees offreedom of corresponding elements representvalues of the elements evaluated at the points or integratedovertheedges,surfaces,orvolumes.Anelementin S iswrittenas

α

h

(ξ,

η

,

ς

)

=

⎧

⎪

⎪

⎨

⎪

⎪

⎩

α

h ξ

(ξ,

η

,

ς

)

α

h η

(ξ,

η

,

ς

)

α

h ς

(ξ,

η

,

ς

)

⎫

⎪

⎪

⎬

⎪

⎪

⎭

=

⎧

⎪

⎪

⎪

⎨

⎪

⎪

⎪

⎩

$

Nξ i=0

$

Nη j=1

$

Nς k=1a ξ i,j,kli

(ξ )

ej

(

η

)

ek

(

ς

)

$

Nξ i=1

$

Nη j=0

$

Nς k=1a η i,j,kei

(ξ )

lj

(

η

)

ek

(

ς

)

$

Nξ i=1

$

Nη j=1

$

Nς k=0a ς i,j,kei

(ξ )

ej

(

η

)

lk

(

ς

)

⎫

⎪

⎪

⎪

⎬

⎪

⎪

⎪

⎭

,

andanelementinV iswrittenas

β

h

(ξ,

η

,

ς

)

=

Nξ



i=1 Nη



j=1 Nς



k=1 bi,j,kei

(ξ )

ej

(

η

)

ek

(

ς

).

If

β

h

=

div

α

h_{,fromSection4.1.1,weknowthattheirexpansioncoefficientssatisfy}

bi,j,k

=

aξi,j,k

−

a ξ i−1,j,k

+

a η i,j,k

−

a η i,j−1,k

+

a ς i,j,k

−

a ς i,j,k−1

,

(15)

whichbasicallyimpliestheGauss’stheorem,seeRemark1.Ifwelabelallexpansioncoefficientsof

α

h_and

_β

h _andputthem

intovectors

α

and

β

,weget

β

=

Ediv

α

,

(16)

where Ediv

:

S

→

V is the incidencematrixrepresenting thediscrete divergence operator.For example,let Nξ

=

Nς

=

1,

Nη

=

2,andwelabelthecoefficientsof

α

and

β

asshowninFig.2.Wehave Ediv

=

%

−

1 0 1 0

−

1 1 0

−

1 0 1 0 0

−

1 0 1 0

−

1 1 0

−

1 0 1

&

.

Similarly, we can obtain Egrad

:

P

→

E and Ecurl

:

E

→

S. The fact that curl

·

grad

(

·)

≡

0 and div

·

curl

(

·)

≡

0 implies EcurlEgrad

≡

0 andEdivEcurl

≡

0.

Remark1.If

α

h _{is the discrete}

_ξ

_{-direction stress}

_σ

h

ξ, and

β

h is the discrete

ξ

-direction body force fξh, the

ξ

-direction

equilibriumof forcesforthe volume

[ξ

i−1

,

ξ

i

]

× [

η

j−1

,

η

j

]

× [

ς

k−1

,

ς

k

]

then impliesdiv

σ

hξ

+

f

h

ξ

=

0 or, inmatrixformat,

Ediv

σ

ξ

+

fξ

=

0 or,moredirectly,

Fig. 2. An example of labeling the expansion coefficients ofαandβfor Nξ=Nς=1 and Nη=2.

whereti,j,k,theexpansioncoefficientsof

σ

hξ,are the(integrated) surface tractions,f ξ

i,j,k, theexpansioncoefficientsof f h

ξ,

aretheintegralvaluesofthebodyforce.

Insummary,wehaveconstructedthefollowingdeRhamcomplex, P Egrad P ⊂ grad H1

(

ref

)

grad E Ecurl E curl H

(

curl

; 

ref

)

⊂ curl S Ediv S div H

(

div

; 

ref

)

⊂ div V V L2

_(

ref

)

⊂

Foracomprehensiveintroductiontomimeticpolynomialspaces,wereferto[2,3].Forsplinebasisfunctionspacesofsimilar structures,wereferto[39,40].

4.1.3. Primaltracepolynomialsin

R

3

We consider thediscrete vector valued function

α

h _{in S. The trace of}

_α

h _{on the face,forexample,}

_(ξ,

_η

_,

_ς

₎

_{∈ }

ξ−

=

−

1

× [−

1

,

1

]

× [−

1

,

1

]

is trξ−

α

h

=

⎧

⎪

⎪

⎨

⎪

⎪

⎩

α

h_ξ

(

−

1

,

η

,

ς

)

α

h η

(

−

1

,

η

,

ς

)

α

_ςh

(

−

1

,

η

,

ς

)

⎫

⎪

⎪

⎬

⎪

⎪

⎭

·

⎧

⎨

⎩

−

1 0 0

⎫

⎬

⎭

= −

Nξ



i=0 Nη



j=1 Nς



k=1 aξ_i,j,kli

(

−

1

)

ej

(

η

)

ek

(

ς

)

=

Nη



j=1 Nς



k=1 aξ_j,−_kej

(

η

)

ek

(

ς

).

Theprimal tracepolynomialsej

(

η

)

ek

(

ς

)

thenspanatrace spaceon



ξ−.Wedenotethistracespaceby

∂

Sξ−.Ifwehave

locallylabeledthecoefficientsaξ_j−_,k,wecancollectandputtheminavector,

α

ξ−.Itisclearthatthereisalinearoperator

N_ξ− whichmaps

α

into

α

_ξ−:

α

_ξ−

=

N_ξ−

α

.

ThematrixN_ξ−,calledthetracematrix[38],liketheincidencematrices,isalsoatopologicalmatrix.Forexample,forthe configurationinFig.2, N_ξ−

=

%

−

1 0 0 0 0 0 0 0 0 0 0 0

−

1 0 0 0 0 0 0 0 0 0

&

.

Ifweapplytheaboveprocesstoall6faces,wecanobtaintracespaces,

∂

Sξ−

,

∂

Sξ+

,

∂

Sη−

,

∂

Sη+

,

∂

Sς−

,

∂

Sς+

,

(17) andtracematrices,

Nξ−

,

Nξ+

,

Nη−

,

Nη+

,

Nς−

,

Nς+

.

(18)

4.1.4. Dualpolynomials

WeconsidertheprimalpolynomialspaceV ,andlet

ϕ

h

_,

_φ

h

_∈

_{V .The L}2_{-innerproductbetween}

_ϕ

h_and

_φ

h_is



ϕ

h

, φ

h



ref

=

ϕ

T_M

_φ,

₍₁₉₎

where M is the mass matrixwhich is symmetric, and, because theprimal polynomialsare linearly independent, isalso bijectiveandpositive-definite.Asaconsequence,itisalwaysinvertible.Thedualpolynomialscanthenbedefinedas



· · · , '

eeei,j,k

(ξ,

η

,

ς

),

· · ·



T

:=

M−1



· · · ,

ei

(ξ )

ej

(

η

)

ek

(

ς

),

· · ·



T

.

ThesedualpolynomialsformanotherbasisofthespaceV .Seealsoequation(28)in[41].Fromnowon,weusethenotation withatildetorepresentanelementexpandedwithdualpolynomials.Forexample,anelement,

φ

h,

φ

h

(ξ,

η

,

ς

)

=

Nξ



i=1 Nη



j=1 Nς



k=1

φ

i,j,kei

(ξ )

ej

(

η

)

ek

(

ς

),

inV hasauniquedualrepresentation,denotedby

(

φ

h,

(

φ

h

(ξ,

η

,

ς

)

=

Nξ



i=1 Nη



j=1 Nς



k=1

(

φ

i,j,keee

'

i,j,k

(ξ,

η

,

ς

),

whosedegreesoffreedomare

(

φ

=

M

φ.

(20)

Notethat

(

φ

hisexactlyequalto

φ

h,butonlytheirrepresentationsaredifferent.Now,theL2-innerproductbetween

ϕ

h _and

(

φ

h is



ϕ

h

, (

φ

h



ref

=

ϕ

T

(

_φ.

₍₂₁₎

It looks trivial;ifwe insert(20) into(21),we immediatelyretrieve(19),but, inpractice, itcan significantlysimplifythe discretization as the L2_{-inner product between them is equal to the vector inner product oftheir expansion coefficient} vectors. More discussions will be given inSection 4.3. We canapply thesame approach to other primal polynomialsto constructdualpolynomialsordualtracepolynomials.Foranexampleoftheusageofthedualtracespace,wereferto[24]. Examples.

1. If

β

h

₌

_div

_α

h

_∈

_{V isexpandedwithprimalpolynomialsand}

(

_φ

h

_∈

_{V isexpandedwithdual polynomials,from(16) and}

(21),wehave



(

φ

h

, β

h



ref

=



(

φ

h

,

div

α

h



ref

= (

φ

TEdiv

α

.

(22)

2. If

α

h

_∈

_{S isexpandedwithprimalpolynomialsand}



_φ

_{∈ ∂}

_{S isexpandedwithdualtracepolynomials,wehave}

)



φ,

tr

α

h

*

∂ref

= 

φ

TN

α

,

(23)

wherethetracespace

∂

S

:= ∂

S_ξ−

× ∂

S_ξ+

× ∂

_Sη−

× ∂

_Sη+

× ∂

_Sς−

× ∂

_Sς+,see(17),andNcanbeobtainedbyassembling thetracematricesin(18).



Remark2.Given

φ

∈

L2

_(

ref

)

,tocalculatethegradientof

φ

,sinceitsspacedoesnotadmitastronggradientoperation,we needtoemployintegrationbypartsanddoitinaweakerwaywithrespecttotheinnerproduct,

#

 grad

φ

·

α

d



:=

#

∂



φ

α

·

nd



−

#



φ

div

α

d

.

Atthediscretelevel,ifweexpand

φ

into

(

φ

h

∈

V withdualpolynomials,expand



φ

into



φ

h

∈ ∂

S withdualtracepolynomials, andexpand

α

into

α

h

_∈

_{S withprimalpolynomials,wehave}



grad

(

φ

h

,

α

h



ref

=

)



φ

h

,

tr

α

h

*

∂ref

−



(

φ

h

,

div

α

h



ref

.

Wethencandefineadiscreteweakgradient,



grad

:

V

× ∂

S

→

S.Let

(

ϑ

h

:= 

grad

(

φ

h

, 

φ

∈

S beexpandedwithdual polynomi-als,from(22) and(23),wecanfindthat

(

ϑ

=

NT



φ

−

ET_div

(

φ.

(24)

This discrete weak gradient essentially is an implementation of integration by parts using the introduced polynomials. Becausethe incidencematrix Ediv andthetrace matrixNareboth topological,such animplementationleads toasimple approachforperformingtheweakgradientatthediscretelevel.Similarly,onecandefinecurl and

+

div and

'

thusconstructa seconddiscretedeRhamcomplex[34]:

0

←

P

×

0

←−

div' E

× ∂

P

←−

curl+ S

× ∂

E

←−

gradV

× ∂

S

← R.

Inthesamewayasshownfor



grad,curl and

+

div also

'

havetopologicalmatrixrepresentationsconsistingofthecorresponding incidencematrixandtracematrix,whichisbeyondthescopeofthispaper.See[34] foracomprehensiveexplanation.

Note that the presented approach of constructing dual polynomials is the most straightforward one. For alternative approaches,wereferto,e.g.,[32,33].

4.2. Coordinatetransformation

The primal anddual polynomials introduced so far are just for the reference element



ref. Let



m be an arbitrary

elementand



m bea

C

1 diffeomorphism,



m:



ref

→ 

m,



x y z



T

= 

m

(ξ,

η

,

ς

),

whoseJacobianmatrixisdenotedby

J

.LetPm,Em, Sm,andVm representthecorresponding mimeticpolynomialspaces

in



m.The primalbasis functionsin



m are obtainedby transformingtheprimal polynomialsin



ref withthe following transformations,[6,37]:

1. Thetransformationbetween

ψ

h

(ξ,

η

,

ς

)

∈

P and

ψ

mh

(

x

,

y

,

z

)

∈

Pmisgivenby

ψ

_mh

(

x

,

y

,

z

)

=



ψ

h

◦ 

−_m1



(

x

,

y

,

z

),

ψ

h

(ξ,

η

,

ς

)

=



ψ

_mh

◦ 

m



(ξ,

η

,

ς

).

2. Thetransformationbetween

ϕ

h

_(ξ,

_η

_,

_ς

₎

_∈

_{E and}

_ϕ

h

m

(

x

,

y

,

z

)

∈

Emisgivenby

ϕ

_mh

(

x

,

y

,

z

)

=



J

T

−1



ϕ

h

◦ 

−1 m



(

x

,

y

,

z

),

ϕ

h

(ξ,

η

,

ς

)

=

J

T



ϕ

h_m

◦ 

m



(ξ,

η

,

ς

).

3. Thetransformationbetween

α

h

_(ξ,

_η

_,

_ς

₎

_∈

_{S and}

_α

h

m

(

x

,

y

,

z

)

∈

Smisgivenby

α

h_m

(

x

,

y

,

z

)

=

J

det

J



α

h

◦ 

−1 m



(

x

,

y

,

z

),

α

h

(ξ,

η

,

ς

)

=

J

−1det

J



α

h_m

◦ 

m



(ξ,

η

,

ς

).

(25) 4. Thetransformationbetween

β

h

(ξ,

η

,

ς

)

∈

V and

β

_mh

(

x

,

y

,

z

)

∈

Vmisgivenby

β

_mh

(

x

,

y

,

z

)

=

1 det

J



β

h

◦ 

−_m1



(

x

,

y

,

z

),

β

h

(ξ,

η

,

ς

)

=

det

J



β

_mh

◦ 

m



(ξ,

η

,

ς

).

For example,let

σ

h

_∈

_[S]3 _{be the Cauchy stress tensor,the Piola transformation(25) convertsbetween}

_σ

h _{and the first}

Piola-Kirchhoffstresstensor

σ

h

m

∈

[Sm]3.

Note that, although the mapping changes the primal polynomials and therefore changes the metric-dependent mass matrices, it doesnot affectthe metric-independenttopological incidencematrices. Thus we havethe followingde Rham complexfor



m, P_m Egrad Pm ⊂ grad H1

(

m

)

grad E_m Ecurl Em curl H

(

curl

; 

m

)

⊂ curl S_m Ediv Sm div H

(

div

; 

m

)

⊂ div V_m Vm ⊂ L2

(

m

)

Thewayofconstructingthedualpolynomialsremainsthesame,and, forexample,therelations(22)-(24) arestillvalidin



m.Therefore,weobtaintheseconddiscretedeRhamcomplexfor



m,

0

←

Pm

×

0 ' div

←−

Em

× ∂

Pm + curl

←−

Sm

× ∂

Em grad

←−

Vm

× ∂

Sm

← R.

And,becausethetracematricesarealsotopological,thematrixrepresentationsforthe



grad,curl and

+

div remain

'

unchanged underthemapping.

4.3. Discretization

We nowcan presentthe discretizationofthe hybridmixedweak formulation(7) withthemimetic polynomials con-structed inprevious subsection.Suppose



h isan orthogonalorcurvilinearconforminghexahedralmeshinthe computa-tional domain



and, foreach element,e.g.,



m,there existsa

C

1 diffeomorphism



m that mapsthereferenceelement



refintoit.WefirstusetheGauss-Lobatto-Legendre(GLL)nodesasthebasisnodestoconstructparticularGLLpolynomial spaces,and,fromnowon,theaforementionednotations,forexampleSmandVm,refertotheirtransformationsin



m.We

alsosetuptheGauss-Legendre(GL)polynomialspaceswhichareonedegreelower,andweusethenotation ˚Pm todenote

corresponding nodalspace. Thisparticularchoiceisforthediscreterotation,

ω

h

_∈



_˚P m



3

,whichistheLagrangemultiplier thatenforcesthesymmetryofthestresstensororequilibriumofmoments.Forcomparisonandcompleteness,wewillfirst brieflyintroducethediscretizationwiththemimeticspectralelementmethod(MSEM)[6].

4.3.1. Mimeticspectralelementmethod

With the MSEM, we discretize the mixedweak formulation(4) in a conventional continuous mesh. In each element



m,thespace[Sm]3 isselectedtoapproximate[H

(

∇·; 

m

)

]3 forthestress

σ

;thespace[Vm]3 isselectedtoapproximate



L2

(

m

)



3

forthe body force f ; the space



˚Pm



3

is selectedto approximate



L2

(

m

)



3

for thedisplacement u and the rotation

ω

.Alldiscretevariablesareexpandedwiththeprimalpolynomials.Suchadiscretizationwilleventuallyleadtothe followingdiscretesystem,

⎧

⎨

⎩

Mm

(

WE

)

T

−R

mT WE 0 0

−

Rm 0 0

⎫

⎬

⎭

⎧

⎨

⎩

σ

u

ω

⎫

⎬

⎭

=

⎧

⎨

⎩

Bm



u

−

Wf 0

⎫

⎬

⎭

.

(26)

Comparingthisdiscrete systemtothemixedweakformulation(4) revealswhateachentryrepresents.NotethatE repre-sentsthemetric-independenttopologicalincidencematrix,thematrixWisadensematrixfortheinnerproductbetween elementsfrom



˚Pm



3

and[Vm]3,andthematrixBmisaboundaryintegralmatrix.Formoreinsightsof(26),seeequations

(20)and(21)in[6].Wedenotethelefthandsidelocalmatrixforelement



m by

F

m.Oncethediscretesystemsforall

el-ementsareconstructedlocally,wecanassemblethem,whichensuresthecontinuityofthesurfacetractionacrosselements andleadstoaglobalsystem.Thisleadstoagloballinearsystemreadytobesolved,andwedenoteitslefthandsideglobal matrixby

F

.

4.3.2. Hybridmimeticspectralelementmethod

Withthehybridmimeticspectralelementmethod(hMSEM),wediscretizethehybridmixedweakformulation(7) ina discontinuousmeshwhereweconsidereachelementasaseparatesubdomain.Inelement



m,thesamefinitedimensional

spacesasmentionedinSection4.3.1areselectedtoapproximatethespacesfor

σ

,

ω

,and f .While,foru,thespace[Vm]3

isselected,and,toobtainahighersparsity,itisexpandedusingthedualpolynomials.Toapproximatethespaces



H1/2

(

·)



3 fortheLagrangemultiplier

λ

,thespaces



∂

S(·)



3

withadualpolynomial basisareselected.Withtheseselections,wecan obtainthefollowingdiscretehybridsystem:

⎧

⎪

⎪

⎨

⎪

⎪

⎩

Mm ET

−R

Tm

−N

T◦ E 0 0 0

−

Rm 0 0 0

−

N_◦ 0 0 0

⎫

⎪

⎪

⎬

⎪

⎪

⎭

⎧

⎪

⎪

⎨

⎪

⎪

⎩

σ

u

ω

λ

⎫

⎪

⎪

⎬

⎪

⎪

⎭

=

⎧

⎪

⎪

⎨

⎪

⎪

⎩

NT_u



u

−

f 0 0

⎫

⎪

⎪

⎬

⎪

⎪

⎭

.

(27)

For

σ

and

ω

,theexpansionsarethesameforMSEMandhMSEM.Wehave



σ

h

,

C

σ

h



m

=

σ

T Mm

σ

,

σ

h

,

σ

h

∈

[Sm]3

,

and



ω

h

,

R

σ

h



m

=

ω

T Rm

σ

,

ω

h

∈



˚Pm



3

,

σ

h

∈

[Sm]3

.