Conformal and Non-conformal Adaptive Mesh Refinement with Hierarchical Array-based Half ...

13 Pages·2015·1.87 MB·English

Availableonlineatwww.sciencedirect.com
ProcediaEngineering00(2015)000–000
www.elsevier.com/locate/procedia
24thInternationalMeshingRoundtable(IMR24)
Conformal and Non-conformal Adaptive Mesh Refinement with
Hierarchical Array-based Half-Facet Data Structures
Xinglin Zhaoa, Rebecca Conleya, Navamita Rayb, Vijay S. Mahadevanb, Xiangmin Jiaoa,∗
aDepartmentofAppliedMathematicsandStatistics,StonyBrookUniversity,StonyBrook,NY11794,USA
bMathematicsandComputerScience,ArgonneNationalLaboratory,Argonne,IL60439,USA
Abstract
WepresentageneralizationoftheArray-basedHalf-Facet(AHF)meshdatastructure,calledHierarchicalAHF,forhierarchical
unstructuredmeshesgeneratedfromadaptivemeshrefinementforsolvingPDEs.ThisdatastructureextendstheAHFdatastructure
(V. Dyedov, et al. AHF: Array-based Half-Facet Data Structure for Mixed-Dimensional and Non-manifold Meshes) to support
meshes with hierarchical structure, which often are generated from adaptive mesh refinement (AMR). The design goals of our
datastructureincludegeneralitytosupportefficientneighborhoodqueries,refinementandderefinement,andhp-FEMwithmesh
smoothing. Ourdatastructureutilizesthesiblinghalf-facetsasacoreabstraction,coupledwithatreestructureforhierarchical
information. Tofacilitatetheinteroperabilityofmeshbasedapplications,auxiliarydatawillbedesignedontopofHierarchical
AHF.Wedescribethedatastructureandsoftwarerequirements,andpresentnumericalexperimentstodemonstrateitseffectiveness.
(cid:13)c 2015TheAuthors.PublishedbyElsevierLtd.
Peer-reviewunderresponsibilityoforganizingcommitteeofthe24thInternationalMeshingRoundtable(IMR24).
Keywords: adaptivemeshrefinement,hierarchicalmeshes,meshgeneration,datastructure,siblinghalf-facets
1. Introduction
Inlargescalesimulationproblems, meshgenerationandlinearsolversareoftenthetwomostexpensivestepsin
thesolutionsofpartialdifferentialequations(PDEs),usingfiniteelementmethodswithunstructuredmeshes. Amesh
withbillionsofelementswillberequired. Since,inmostcasesthecriteriaformeshresolutionisnotknownapriori,
inordertoavoidcomputationaloverhead,thegeneratedoriginalmeshisoftenrelativelycoarseoverallwithuniform
resolution. However, to obtain accuracy some regions need to be refined to reduce discretization errors while other
regions require finer models to approximate. Adaptive mesh refinement allow more efficient numerical simulations
byincreasingthecomputationaleffortnearinterestingfeaturesofthesolutions[4–6].
AMR has gradually become a vital step in large-scale numerical simulations since it optimizes the relationship
between accuracy and computational effort. One aspect of the refinement strategy is whether it requires the refined
meshtobeconformalornot. Ameshissaidtobeconformalifthepairwiseintersectionofanytwoentitiesiseither
alower-dimensionalentityorisempty. Otherwise,ameshisnon-conformal. Theconformalrequirementwillmake
∗ Correspondingauthor.Tel:[phone];fax:[phone].
E-mailaddress:[email protected]
1877-7058(cid:13)c 2015TheAuthors.PublishedbyElsevierLtd.
Peer-reviewunderresponsibilityoforganizingcommitteeofthe24thInternationalMeshingRoundtable(IMR24).
2 X.Zhao,R.Conley,N.Ray,V.S.Mahadevan,X.Jiao/ProcediaEngineering00(2015)000–000
nochangetotheunderlyingdatastructureforthemeshandtheformulationofnumericalalgorithms. Aconsiderable
amount of work has been done in this area [7,9,10,25]. However, to preserve conformity, some procedures need to
beappliedwhichwouldprobablydeliverafinermeshthanneeded,orevenpotentiallyaffecttheoverallmeshquality
which is crucial for the linear system in FEM [26]. This drives the research on refinement strategies allowing non-
conformalmeshes,i.e. hangingnodes,in[7,15,18,26]. Non-conformalrefinementwillincurextraworkinthePDE
solver part, but it will be much easier for the unification of hp-adaptivity for finite element method [7,8]. Here hp-
adaptivitymeansthatboththemeshsizehandthedegreepoftheapproximatingpiecewisepolynomialsareadapted.
Most of the implementations for mesh adaptation adopt a pointer-based mesh data structure, since they are rel-
atively easy to manipulate. In this paper, we develop an array-based mesh data structure to support adaptive mesh
refinementandderefinement.1 ItgeneralizesAHF[12],whichprovidesefficientmeshqueriesandmodification. The
array-based mesh data structures have many advantages in the context of numerical simulations, in terms of more
compact memory footprint, better interoperability with simulation codes, better efficiency on modern computer ar-
chitectures with deep memory hierarchy, and relative simplicity and higher efficiency for parallel implementations.
However, it is more challenging to support adaptive mesh refinement with array-based mesh data structures, which
requiredynamiccreationanddeletionofentities.
Thekeycontributionsofthispaperaremainlytwofold: First,weintroduceasimpledatamodelformesheswith
hierarchicalstructure. Ourdatamodeliseasytoimplementandisefficientinbothmemoryandcomputationalcost.
Thedatastructurefacilitatesbothstraightforwardrefinementandderefinementoperations,andalsoallowsbothcon-
formal and non-conformal meshes. In addition, a generic adaptive mesh refinement (AMR) framework on top of
Hierarchical AHF is developed and a prototype is implemented for both 2D triangular and 3D tetrahedral meshes.
Second, as an array-based data structure, AHF facilitates better interoperability across different application codes,
different programming languages (such as MATLAB, C/C++, FORTRAN, etc.), and different hardware platforms.
Meanwhile, sinceboththetreehierarchyandmeshdataarearraybased, bettermemorycompactnessandcomputa-
tionalefficiencycouldbeachieved. Moreover,thedatamodelcanbeeasilyintegratedwithmulti-levelmethodssuch
as multigrid solvers. Efficient intra- and inter-level mesh traversals are supported and the data structure is flexible
enoughtosupportbothuniformmeshrefinement(UMR)andAMR.WedevelopedUMRwithsurfacereconstruction
forthemultigridmethodofthefiniteelementmethod,whichwewillreportelsewhere. TheC++implementationon
topofMOAB[29]willbebasedontheworkofMOAB_AHF[12].
The remainder of the paper is organized as follows. Section 2 reviews some background knowledge and related
meshdatastructures. Section3describesourdatamodelforhierarchicalmeshes. Section4describesthealgorithms
for the construction, query, mesh modification operations, as well as their implementations in MATLAB. Section 5
presentssomenumericalresults. Section6concludesthepaperwithadiscussion.
2. BackgroundandRelatedWork
Inthissection,wefirstbrieflyreviewthemeshadaptationmethodologyfornumericalPDEs. Thensometerminol-
ogyformeshdatastructuresisexplainedandsomeexistingdatastructureswillbereferencedforcomparison,which
willestablishthefoundationofourproposeddatastructure.
2.1. MeshAdaptationforNumericalPDEs
AdaptivemethodsfornumericalPDEshavebeenanactiveresearchareasincethelate1970s[4,5]andarewidely
usedinpracticenowadays,tobalanceaccuracyandcomputationalefficiency. Inparticular,AdaptiveFEM(AFEM)
basedonthelocalmeshrefinementhasloopsofthefollowingform:
SOLVE→ESTIMATE→MARK→REFINE
toiterativelyimprovetheaccuracyofthenumericalapproximation. Generally,anaposteriorierrorestimatorisused
tomeasuretheaccuracyofobtainednumericalsolutions,whichisexactlytheESTIMATEmodulementionedabove.
1 Weusetheterm“derefinement”insteadof“coarsening”becausethealgorithmwouldonlyundotherefinementselectively,anditwouldnot
coarsenbeyondtheoriginalmesh.
X.Zhao,R.Conley,N.Ray,V.S.Mahadevan,X.Jiao/ProcediaEngineering00(2015)000–000 3
1 1
4 6 4 6
2 2
5 3 5 3
7
Figure1:Non-conformalrefinementwithhangingnodes. Figure2:Conformalrefinementswithtransientelements.
Elements with large errors are marked and adaptive mesh refinement/coarsening strategies are utilized to minimize
theerror. Someoftherefinementstrategieswilldeliverconformalmesheswhileotherswillnot. Ourgoalistodesign
adatastructurethatsupportsbothconformalandnon-conformaladaptivemeshrefinement.
Intheprocessofadaptivemeshrefinement,anaposteriorierrorestimatorwouldindicateelementswithlargeer-
ror. Theseelementswouldbemarkedandintheh-adaptivity,theREFINEmodulerefinesallthesemarkedelements.
Different refinement strategies have been considered. In 2D, during the middle of the 1980s Rivara introduced an
effectivemeshrefinementalgorithmbasedonlongestedgebisection[24]whileMitchelldevelopedarecursivealgo-
rithmforthenewestvertexbisection[21]andBankadoptedregularrefinementwithselectedtemporarybisections[2].
Bank’smethodisknownasred-greenrefinementandwasusedinthesoftwarepackagePLTMG[9]. Inthebeginning
ofthe1990s,RivaraandLevinextendedthelongestedgerefinementalgorithmtotetrahedralmeshes[25]. However,
it is not clear whether this algorithm degrades the mesh quality. Meanwhile, extensions of red-green refinement to
3D were considered in [11] and Bänsch generalized the newest vertex bisection method to 3D [10]. Both of them
preservethemeshqualityunderrefinement. SimilarapproachesweredevelopedbyLiuandJoe[19]andArnoldetal.
[3]. Moreover,Kossaczky[16]derivedarecursivevariantofBänsch’salgorithm,withabisectionrulefortetrahedra
usingthelocalorderofverticesandelementtype. Thisconceptisconvenientforimplementationandgeneralization
toanyspacedimension. Foramorecompletediscussionofmeshrefinement,see[13,22].
Generally,thenewestvertexrefinementwilldeliveraconformalmesh,buttheminimumanglesofthemeshwillbe
degraded. Ontheotherhand,regularrefinementwouldpreservetheanglesduringrefinement,butirregularnodes(or
hangingnodes)willbecreated. Forinstance,seeFigure1;triangle∆123ismarkedandrefinedinaregularway. The
refinementwouldresultinunbalancedvertices4,5,and6,a.k.a. “hangingnodes.” Generally,therearetwostrategies
to deal with hanging nodes in FEM: 1) associate degrees of freedom with the hanging nodes and eliminate them in
thelinearsystemaccordingtocontinuityconstraints[8,26];or2)converttheneighboringcellsofthehangingnodes
intotransientcells,asshowninFigure2. Thelatterisoftenreferredtoasred-greedrefinement[2].
Intermsofimplementation,1-irregularityruleisoftenappliedfornon-conformalAMR,whichrequiresthateach
edgehaveatmostonehangingnode. Ingeneral,k-irregularityrule[27]couldbeapplied,whichmeansthateachedge
couldhaveatmostkhangingnodes. Then,k = 0meansnohangingnodesareallowedandtherefinedmeshremains
conformal,andk=∞correspondstoadaptivitywitharbitrary-levelhangingnodes.
2.2. Terminology
Inoursetting,ameshisasimplicialcomplexdiscretelyrepresentingageometricortopologicalobject. Topolog-
ically,ad-dimensionalobjectisamanifoldwithboundaryifeverypointinithasaneighborhoodhomeomorphicto
either a d-dimensional ball or half-ball, where the points whose neighborhood is homeomorphic to a half-ball are
boundary(orborder)points. Wesayameshis1D,2D,or3Diftheobjectthatitrepresentsistopologically1D,2D,
or3D,respectively. Ameshiscomposedof0D,1D,2D,and3Dentities,whichwerefertoasvertices,edges,faces,
andcells,respectively.Typically,afaceiseitheratriangleorquadrilateral,andacellisatetrahedron,prism,pyramid,
orhexahedron,especiallyforfiniteelementmethods,althoughgeneralpolygonsandpolyhedraarealsooftenusedin
finitevolumemeshes.
Inad-dimensionalmesh, werefertothed-dimensionalentitiesaselementsandrefertothe(d−1)-dimensional
sub-entitiesasitsfacets. Morespecifically,thefacetsofacellareitsfaces,thefacetsofafaceareitsedges,andthe
facetsofanedgeareitsvertices. Eachfacethasanorientationwithrespecttothecontainingelement. Forexample,
4 X.Zhao,R.Conley,N.Ray,V.S.Mahadevan,X.Jiao/ProcediaEngineering00(2015)000–000
eachedgeofatrianglehasadirection,andalltheedgesformanorientedloop. Thusitmakessensetocallthefacets
ashalf-facets. Eachfacetmayhavemultipleincidentelements,especiallyfornon-manifoldentities. Werefertoall
suchhalf-facetsassiblinghalf-facets. Ahalf-facetwithoutanysiblingsisaborderhalf-facet,andwerefertoavertex
incidentonaborderhalf-facetasabordervertex.
2.3. Half-FacetDataStructure
The half-facet data structure is a generalization of the concept of the doubly-connected edge list (DCEL) for
surfaceandvolumemeshes[1,12,17]. Inanorientedmanifoldsurfacemesh,supposetheedgeswithineachfacecan
beorderedinacounter-clockwisedirectionwithrespecttotheoutwardnormal(orupwardnormalfor2Dmeshes).
Foravolumemesh, withineachcell, supposetheedgesofeachfaceareorientedinacounter-clockwiseorderwith
respecttotheoutwardnormalofthecell. In2D,theedgeswithineachfacearecalleddirectededgesorhalf-edges
whilewerefertotheorientedfacesashalf-facesin3D.Eachedgehastwoincidentfaces,andthetwohalf-edgeshave
opposite orientations and hence are said to be opposite or twin half-edges of each other. An edge on the boundary
doesnothaveatwinhalf-edge. Fortypicalmeshesinengineeringapplications,eachfaceintheinteriorofavolume
meshhastwocorrespondinghalf-faceswithoppositeorientations,whicharesaidtobeoppositeortwinhalf-facesof
eachother.
2.4. Pointer-BasedVersusArray-BasedImplementations
Ameshdatastructuremaybeimplementedusingeitherpointersorarrays. Thepointer-basedimplementationsare
morecommon,sincetheyarerelativelyeasytomanipulate. deal.II[7],aC++finiteelementlibrary,supportshp-FEM
in1D,2D(quadrilaterals)and3D(hexahedra),andallowshangingnodesintroducedinhp-adaptivity. Hangingnodes
areeliminatedaccordingtocontinuityconstraints[8]. Likewise,libMesh[15]isalsoaC++libraryforserial/parallel
adaptivealgorithms. Hermes2D[28]supportsadaptiveFEMin2Dbasedonthealgorithmin[26].
In our work, we choose to use an array-based, pointer-free implementation for a number of reasons. First, in an
array-based implementation, we can treat intermediate dimensional entities (such as half-facets) as implicit entities
and reference them without forming explicit objects. This can lead to significant savings in storage, especially on
computers with 64-bit pointers. Second, using arrays can also lead to faster memory access and hence better effi-
ciency. Inaddition,array-basedimplementationsalsoofferbetterinteroperabilityacrossapplicationcodes,different
programminglanguages,anddifferenthardwareplatforms(suchasbetweenGPUsandCPUs).
3. DataModelforHierarchicalMeshes
In this data model, we assume that each element has a standard numbering convention for its vertices and its
facets. Forstandardelements,wefollowtheconventionoftheCGNS(CFDGeneralNotationSystem)[23,30]. We
donotrequireexplicitrepresentationofintermediatedimensionalentitiesbetween1andd−1. Instead,wetreatthe
half-facetsasimplicitentities,andrefertoahalf-facetusingtheelementIDanditslocalIDwithintheelement.
In the process of AMR, to avoid the duplication of new vertices introduced by refinement, efficient adjacency
queries are critical. The AHF data structure provides efficient query operations with nice memory performance. A
hierarchicalstructureisgenerallynecessaryformulti-levelmethodsforthelinearsystemofnumericalPDEs. Inour
datamodelthehierarchyisstoredinanarray-basedtree-likestructure.WerefertothisdatamodelastheHierarchical
AHF.
3.1. HierarchicalStructure
Thedesignofthemeshdatastructureforadaptivemeshrefinementassumesthatwestartwithaconformalman-
ifold mesh. An initial conformal mesh is easy to generate and it is natural to form a hierarchical structure by mesh
adaptation. Therefinementandderefinementrequiresefficientadjacencyqueries,whichareprovidedbyAHF.
The initial mesh is adaptively refined and the results will be stored in a hierarchical structure. The resulting
elements from a subdivision of element K will be referred as the child elements of K, which in turn is called the
X.Zhao,R.Conley,N.Ray,V.S.Mahadevan,X.Jiao/ProcediaEngineering00(2015)000–000 5
Coordinatesforallvertices
Ω
verticesonL1 newverticesbyrefinement
e1 e1 e1
1 2 3
v2pe:newvertextooneofitsparents
v2hf: vertextohalffacets,anchorusedtolocatevertex,optional
e2 e2
HierarchicalMesh 1 5
conn connectivityoforiginalmesh e3
1
sibhfcs DatainAHFforadjacencyqueries
L1
e2pe NULL
e2ce elementtoitsfirstchildelement
conn connectivityofnewelementsinL2
sibhfcs updateAHFdataforrefinedelements,optional v2pe
L2
e2pe newelementtoitsparentelements,optional e2ce
e2pe
e2ce elementtoitsfirstchildelement
Figure3:Hierarchicalarray-basedhalf-facetdatastructureforamulti-levelmesh.
parent. IfelementK isrefined,thenitissaidtobeinactive. Elementsgeneratedfromsubsequentrefinementofthe
childrenofK willbecalledthedescendantsofK. Onthefirstlevel,theoriginalmeshisstored. Thensomeelements
ofthemesharemarkedforrefinement. Thesecondlevelwouldbethechildelementsoftheseelements;seeFigure3.
Oneachlevel,thechildelementsoftheupperlevelwillformaconformalmesh(whichmightnotbemanifold). This
isanalogoustoquad-tree. Forinstance,inFigure3,onlevel1,elemente1,e1,e1formtheinitialconformalmesh. On
1 2 3
thesecondlevel,thechildrenofe1 ande1 ,i.e. e2,...,e2 willalsoformaconformalmesh. Totraversethetree,we
1 2 1 8
store e2ce for each element, which is the mapping from the elements to the IDs of their child elements on the next
level. Oneachlevel, e2ceisrepresentedasanarray. Forregularrefinement, e2cewillonlystoretheIDofthefirst
childelement,sinceallchildrenarestoredinconsecutiveorderinanarray.
3.2. HierarchicalAHF
Inthehierarchicalmeshdatastructure,thetopologicalinformation,i.e.,theconnectivitytableofelementswillbe
storedforeachlevelofthemesh. Theoriginalmeshistreatedasthefirstlevelofthemesh. Duringtherefinement,
some elements of the mesh are marked to be refined. The second level would be the child elements; see Figure 3.
On each level, the connectivity will be stored in conn of the mesh data structure. Each level of the sub-mesh will
containverticesbothfromthecurrentlevelandpreviouslevels,thusstoringverticesforeachlevelwouldbeawaste
ofmemory.Thereforewestorethegeometricdata,i.e.coordinatesofallvertices,inaseparatearray.TheHierarchical
AHFrepresentationisillustratedinFigure3.
To support efficient intra- and inter-level queries, auxiliary information is necessary. For the intra-level queries,
the neighboring information, i.e. AHF data will be stored in an array for each level sub-mesh. Since the sub-mesh
is conformal on each level, AHF data can be built in a natural way and the data is represented as sibhfcs (sibling
half-facets)inthemeshdatastructureinFigure3. Intheprocessofmeshadaptation,theneighborinformationsibhfcs
6 X.Zhao,R.Conley,N.Ray,V.S.Mahadevan,X.Jiao/ProcediaEngineering00(2015)000–000
1 5 1 8 5
e1
e2ce e2pe
e2 4 6 7 4 elems child elems parent
2 e3 3 Reefi1ne 2 3 eee123 000 ee15−−ee48 12
e2 e4 1
elements newelements newsibhes
e1 1 2 5 e1 1 6 8 e1 0 h2,3i 0
e2 2 4 5 e2 6 7 8 e2 h3,3i h4,1i h1,2i
e3 2 3sibh4es ee34 68 27 75 ee34 h20,2i 00 h20,1i hh58,,11ii
eee123 eh23c00,e3i h200,3i hh120,,21ii eeee5678 117523 112732 117423 eeee5678 hhhh3764,,,,2112iiii hh8600,,22ii hh6500,,33ii
elems child
e1 1 Implicitrefinement
e2 5 Refine
e3 0 e4
1 11 5
8
10 13
v2pe 9
verts hlevel,eid,lidi 6 7 4 elemse2peparent
vvv678 hhh111,,,111,,,123iii 2 12 3 e1−e4 4
v9 h2,4,1i newelements newsibhes
vvvv11110123 hhhh2211,,,,4422,,,,2312iiii eeee1234 118911 19790 1115100 eeee1234 hh1200,,23ii hh2300,,13ii hh4200,,12ii
Figure4:ExampleofHierarchicalAHFunderrefinement.
willbeupdatedincrementally. Forinter-levelqueries,extrainformation(likee2peande2ce)isstoredinarrays. For
eachelementonacertainlevel,e2ceistheIDofthefirstchildelementonthenextlevel. Onthenextlevel,otherchild
elementsofthiselementwillbestorednexttothefirstchildelement. e2ceisstoredinanarrayforeachlevelanditis
necessaryforinter-leveltraversals. e2pe,elementtoparentelement,istheIDoftheparentelementonpreviouslevel
ofthiselementanditisoptional.
Tosupportefficientqueriestotheparentelementforeachnewvertex,aseparatemappingv2peisstored. Foreach
new vertex, v2pe is an array of tuples: level, eid, lid, where level is the level of its parent element, eid is the ID of
theparentelementinlevel,lid isthecanonicalIDofthisvertexinitsparentID.Ifthe1-irregularityruleisapplied,
the lid would be the same as the local ID of the refined edge. v2pe is generally necessary for multi-level methods.
Alsowecouldusev2petodeterminewhichvertexisahangingnodeonwhichlevel. Generally,ifthe1-irregularity
ruleisapplied,vertexvcouldonlybeahangingnodeonlevelv2pe(v).level+1. Thiscouldbefurtherdeterminedby
checkingifv2pe(v).eid’ssiblingelementsarerefined. Ifnot,thenvisahangingnode.
InFigure4,weillustratethedatastructurebyrefiningasimplemesh. First,auserspecifiesrefinementofe . The
1
newelementswillbecreatedande2ceinlevel1willbeupdated. Thentheuserspecifiese onthesecondleveltobe
4
refined. Herethe1-irregularityruleisappliedtokeepthemeshgraded. Thiswillintroduceanimplicitrefinementof
e onthefirstlevel. Correspondingly,dataonlevel2willbeupdated. v2pewillbestoredinaseparatearray.
2
4. ConstructionandModificationofHierarchicalAHF
In this section, we describe some detailed algorithms for the construction of Hierarchical AHF, as well as some
query and modification operations. Since AHF is array-based, these algorithms can be implemented in many pro-
gramming languages, including MATLAB, C/C++, FORTRAN, etc. We will also describe our implementation in
MATLAB.
X.Zhao,R.Conley,N.Ray,V.S.Mahadevan,X.Jiao/ProcediaEngineering00(2015)000–000 7
4.1. ConstructionofHierarchicalAHF
In the half-facet data structure, there are two components: sibhfcs (sibling half-facets) and v2hf (vertex to half-
facet). TheformeriscentraltoAHF,asnearlyalladjacencyqueriesrequireit. Thesesiblinghalf-facetsshouldmap
toeachotherandformacycle. Thelatterarray,v2hf,isoptionalformanyoperations;forHierarchicalAHF,itisnot
builtfornewvertices. Instead,v2pe(vertextoparentelement)isconstructedtostoreinformationofnewvertices.
Inahierarchicalmesh,theAHFfortheinitialmeshwillbeconstructedfirstandthenthesub-meshofeachlevel
isconstructedincrementally, takingadvantageofancestorinformation. Ingeneral, therefinementandderefinement
requiredifferentalgorithms. Inthefollowing,wedescribethesetwopartsinamannerindependentofthedimension
ofthemesh.
4.1.1. HierarchicalAHF:Refinement
Algorithm 1 outlines the steps for mesh refinement, which is applicable to half-facets in arbitrary dimensions,
and is particularly efficient in 1 to 3 dimensions. New elements are created and appended in corresponding levels.
Meanwhile,theadjacencyinformation,sibhfcs,isupdated. Thissteprequirestheinputofelementsthataremarked
forrefinement:
refTags: arraysstoredinahierarchicalstructure,whichstoretheelementstoberefinedoneachlevel.
Algorithm1UpdateSiblingHalf-FacetsforRefinement.
Input: hielems: hierarchicalmeshdata,refTags
Output: sibhfcs: cyclicmappingsofsiblinghalf-facetsforeachlevelofmesh
1: foreachlevelinhielemsdo
2: foreachelementeinrefTags(level)do
3: foreachedgeinelementedo
4: Loopthroughelementsinlevelincidenttoedgetocheckifedgeisrefined
5: ifedgeisnotrefinedthen
6: Refineedgebyinsertingvertexvinthemiddle;
7: Updatev2peforvertexv,v2pe(v)=(cid:104)level,e,edge(cid:105);
8: Refineelementebypredefinedstrategyandupdatee2ce(e)
9: Constructsibhfcsforchildrenofelemente;
{Updatesibhfcsforsubmeshonlevel+1:}
10: foreach facetinelementedo
11: Checkoppositeelementof facetonlevelofsubmesh
12: ifoppositeelementisrefinedthen
13: Updatesibhfcsforchildrenofelemente;
14: Updatesibhfcsforchildrenofoppositeelement;
15: else
16: Updatesibhfcsforchildrenofelemente;
The computational cost of Algorithm 1 is linear, assuming that the number of elements incident on an edge is
boundedbyasmallconstantc. Forthestoragerequirement, let|C |denotethenumberofelementstoberefinedin
r
a certain level of the mesh. The amortized memory storage increased by refinement will be approximately (2dv +
c
2df +2d)|C |integers,fortheconnectivity,theneighborinformation,andinter-levelmaps,withextraspacefornew
c r
verticescoordinatesandv2pemap. Herev and f arethenumbersofverticespercellandthenumberoffacesper
c c
cells,2d isthenumberofchildrenperelement.
4.1.2. HierarchicalAHF:Derefinement
During the second step, we update the sibling half-facets during derefinement. This step requires the input of
elementsthataremarkedforderefinement:
8 X.Zhao,R.Conley,N.Ray,V.S.Mahadevan,X.Jiao/ProcediaEngineering00(2015)000–000
derefTags: ahierarchicalstructurewhichstoreselementstobederefinedoneachlevel. Foreachelementeinderef-
Tags, weassumethateisrefinedandthederefinementoperationwillremoveallchildrenofeandsetetobe
active.
Algorithm2outlinestheprocedureforthisstage,whichisapplicabletohalf-facetsofarbitrarydimensions. Particu-
larly,avertexisactiveifandonlyifithasincidentcells.Ahangingnodewillbesetasinactiveifallincidentelements
areremovedduringderefinement.
Algorithm2UpdateSiblingHalf-FacetsforDerefinement.
Input: hielems: hierarchicalmeshdata,derefTags
Output: sibhfcs: cyclicmappingsofsiblinghalf-facetsforeachlevelofmesh
1: foreachlevelinhielemsdo
2: foreachelementeinderefTags(level)do
3: foreachedgeinelementedo
4: Loopthroughelementsinlevelincidenttoedge;
5: ifnoneofincidentelementsisrefinedthen
6: Derefineedgebyremovingvertexvwhichisinthemiddle;
7: Updatev2peforvertexv,v2pe(v)=(cid:104)0,0,0(cid:105);
8: else
9: edgecannotbederefined;
10: vertexvwhichisinthemiddleisstillactive
11: Derefineelementeandsete2ce(e)=0;
12: Setallitschildrenmute;
{Updatesibhfcsforsubmeshonlevel+1:}
13: foreach facetinelementedo
14: Checkoppositeelementof facetonlevelofsubmesh
15: ifoppositeelementisrefinedthen
16: Updatesibhfcsforchildrenofoppositeelement;
17: Setsibhfcsforchildrenofelementetozeros;
SimilartoAlgorithm1,thecomputationalcostofAlgorithm2isalsolinear,assumingthatthenumberofelements
incidentonanedgeisboundedbyasmallconstantc.Toanalyzethestoragerequirement,let|C |denotethenumberof
d
elementstobecoarsenedinacertainlevelofthemesh. Theupdatewillintroduceapproximately(2dv +2df +2d)|C |
c c d
“holes”intheelementconnectivity,sibhfcsarraysandthemape2pe.Dynamicmemorymanagementcouldbeutilized
to reuse such holes, for instance, by building a queue to record the holes in the corresponding array introduced by
deletionandhavinganynewinsertionreusethememory.
4.2. MeshAdaptation
MeshrefinementandderefinementcanbeimplementedrelativelyeasilyinAHF.Forhierarchicalmeshes,AHFis
particularlyattractivebecausetheadaptivitycouldbeperformedefficientlyandAHFcanbemodifiedincrementally.
Thisleadstoverymodularadaptivitystrategies. Toavoidexcessivememorycopying,weexpandthearraybyasmall
percentage(e.g. by20%)eachreallocation,sothattheamortizedcostforthelocalmodificationsisconstant.
The data structure could support a refinement strategy whether the mesh is required to be conformal or not. In
ourMATLABimplementation,wesupportbothregularrefinementandred-greenrefinement(Figure5). Particularly,
we enforce the 1-irregularity rule for the non-conformal refinement. The Kelly error indicator [14] is utilized for
estimating accuracy and marking elements. The AHF code [12] is used to generate sibling half facet data for the
initialmesh.
X.Zhao,R.Conley,N.Ray,V.S.Mahadevan,X.Jiao/ProcediaEngineering00(2015)000–000 9
Figure5:Red-GreenRefinement
(a)OriginalMesh (b)6thAMR (c)Zoomed-InMeshAroundCen-
ter
Figure6:AFEM:Meshadaptationforre-entrantcornerproblem
5. NumericalExperiments
5.1. AdaptiveFiniteElementMethod(AFEM)
WeimplementedAMRforthere-entrantcornerproblemfrom[20]withtheKellyerrorindicator. Theequationis
∂2u ∂2u
− − =0
∂x2 ∂y2
ondomain[−1,1]×[−1,1]\{0≤ x≤1,y=0}. Theboundaryconditionisu=gandtheexactsolutionis
(cid:18)θ(cid:19)
u=r12 sin
2
(cid:112)
wherer= x2+y2,θ=tan−1(y/x)∈[0,2π).
Weapply6adaptiverefinementsovertheoriginalmeshandcompareitwithFEMonameshwithuniformregular
refinement. TheresultsareshowninFigure6,7. Theoriginalmeshhasacrack{0 ≤ x ≤ 1,y = 0}alongwhichthe
solutionisnotsmooth. Therefinementiscentralizedalongthecrackduetothenon-smoothnessofthesolution,see
Figure6b. TheresultinFigure7ashowsthattheadaptiveFEM´approachdeliversabetterconvergenceratethanFEM
over a uniformly refined mesh. The L error is computed as |u−u |2dA. The numerical results indicate that the
2 Ω h
sameaccuracycouldbeachievedwithmanylessDOFsornumberofelements.
5.2. ComparisonwithPointer-BasedDataStructure
We will compare the storage for Hierarchical AHF with the pointer-based data structure in libMesh, as it is the
most closely related to our data structure. Let C and V represent the set of cells and vertices of the given mesh,
and let |·| denote the size of a set. For Hierarchical AHF, let C and V denote the cells and vertices, respectively,
1 1
10 X.Zhao,R.Conley,N.Ray,V.S.Mahadevan,X.Jiao/ProcediaEngineering00(2015)000–000
10−1 105
Uniform Refinement Uniform Refinement
Adaptive Refinement
Adaptive Refinement
s
ent104
m
Error10−2 of ele
L 2 mber 103
u
N
10−130 2 103 104 105 1021 2 3 4 5 6 7
Number of DOFs Mesh Adaptation
(a)Convergencerate (b)Numberofelements
Figure7:AFEM:Numericalresultsforre-entrantcornerproblem
intheoriginal, unrefinedmesh. WewillconsideranimplementationofHierarchicalAHFthatincludestheelement
connectivity (elements), vertex to parent element mapping (v2pe), sibling half-facet mapping (sibhfcs), element to
parentelement(e2pe)andelementtochildelementmapping(e2ce). Sincethevertextoincidenthalf-facetmapping
(v2hf)isoptional, wewillnotincludeitinthisanalysis. Therefore, wehavethefollowingfivemapswhichrequire
thefollowingnumberofentities:
elementconnectivity: n =v |C|
c c
vertextoparentelementmap: n =|V|−|V |
p 1
siblinghalf-facetmap: n = f |C|
s c
elementtoparentmap: n =|C|−|C |
ep 1
elementtochildmap: n =|C|
ec
where v and f are the numbers of vertices per cell and the number of faces per cells, respectively. In general,
c c
the entities are stored as 32-bit integers. For the half-facet ID (cid:104)eid, lfid(cid:105), we encode it in a 32-bit integer. For the
vertex to parent element map we store (cid:104)level, eid, lid(cid:105) as two 32-bit integers, one for level and one for (cid:104)eid, lfid(cid:105).
Thusthestorageinbytesis
S =4n +8n +4n +4n +4n
AHF32 c p s ep ec
=4(2+v + f )|C|−4|C |+8|V|−8|V |
c c 1 1
Ifweweretostoretheentitiesas64-bitintegers,thestoragewouldeffectivelydouble.
For each cell, libMesh stores the element connectivity and the “face neighbors” of the cells. Two cells are face
neighborsiftheyshareaside;in1Dasideisavertex,in2Dasideisanedge,andin3Dasideisaface. LikeHier-
archicalAHF,adaptivemeshrefinementandcoarseningiscentraltolibMeshandhencethecellsandtheirancestors
arestoredinatree. Specifically,apointertotheparentofanelementandanarrayofpointerstoitschildren(ifany)
arestored. Ingeneral,ad-dimensionalelementisrefinedinto2d childrenofthesametypeexceptwhendealingwith
pyramids, which are refined into pyramids and tetrahedra. For the sake of simplicity, we will use 2d as the number
ofchildrenofanelement. NotethatthelevelofanelementisnotstoredinlibMesh, sincethiscanbefoundrecur-
sively from the parents. To store nodal information, libMesh has a node class. Each object in the node class stores
thecoordinatesofthenode,auniqueglobalIDnumberandthedegreeoffreedomindices. Sincewearecomparing
the storage for the topological information of the mesh, we will consider the storage cost of the global ID number.
Thereforewehave4mapsrequiringthefollowingnumberofentities:
elementconnectivity: n =v |C|
c c
neighboringobjects: n = f |C|
n c
hierarchical: n =|C|+2d|C |
h r
nodal: n =|V|
v