The divergence-conforming immersed boundary method

(1)

Delft University of Technology

The divergence-conforming immersed boundary method

Application to vesicle and capsule dynamics

Casquero, Hugo ; Bona-Casas, Carles ; Toshniwal, Deepesh; Hughes, Thomas J.R. ; Gomez, Hector ;

Zhang, Yongjie Jessica

DOI

10.1016/j.jcp.2020.109872

Publication date

2021

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

Casquero, H., Bona-Casas, C., Toshniwal, D., Hughes, T. J. R., Gomez, H., & Zhang, Y. J. (2021). The

divergence-conforming immersed boundary method: Application to vesicle and capsule dynamics. Journal

of Computational Physics, 425, 1-26. [109872]. https://doi.org/10.1016/j.jcp.2020.109872

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.

(2)

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

The

divergence-conforming

immersed

boundary

method:

Application

to

vesicle

and

capsule

dynamics

Hugo Casquero

a

,

∗

,

Carles Bona-Casas

b

,

Deepesh Toshniwal

c

,

Thomas J.R. Hughes

d

,

Hector Gomez

e

,

Yongjie

Jessica Zhang

a

a_Department_of_Mechanical_Engineering,_Carnegie_Mellon_University,_Pittsburgh,_PA_15213,_United_States b_Departament_de_Física_&_IAC3_,_Universitat_de_les_Illes_Balears,_Palma_de_Mallorca,_07122,_Spain

c_Delft_Institute_of_Applied_Mathematics,_Delft_University_of_Technology,_Van_Mourik_Broekmanweg_6,_XE_Delft_2628,_the_Netherlands d_Oden_Institute_for_{Computational}_Engineering_and_Sciences,₂₀₁_East_24th_Street,_C0200,_Austin,_TX_78712-1229,_United_States e_School_of_Mechanical_Engineering,_Purdue_University,_West_Lafayette,_IN_47907,_United_States

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received20December2019 Receivedinrevisedform19September 2020

Accepted23September2020 Availableonline6October2020

Keywords:

Fluid-structureinteraction Immersedboundarymethod Volumeconservation Isogeometricanalysis Vesicles

Capsules

We extend the recently introduced divergence-conforming immersed boundary (DCIB) method[1] toﬂuid-structureinteraction(FSI)problemsinvolvingclosedco-dimensionone solids.Wefocusoncapsulesandvesicles,whosediscretizationisparticularlychallenging duetothehigher-orderderivativesthatappearintheirformulations.Intwo-dimensional settings, weemploy cubic B-splineswith periodicknot vectorstoobtain discretizations of closed curves with C2 _{inter-element} _continuity. _In _{three-dimensional} _settings, _we use analysis-suitable bi-cubic T-splines toobtain discretizations ofclosed surfaces with at least C1 inter-element continuity. Largespurious changesof the ﬂuidvolume inside closedco-dimensiononesolidsarea well-knownissueforIBmethods.TheDCIBmethod resultsinvolumechangesordersofmagnitudelowerthanconventionalIBmethods.This is a byproductof discretizingthe velocity-pressure pairwith divergence-conforming B-splines,whichleadtonegligibleincompressibilityerrorsattheEulerianlevel.Thehigher inter-element continuity ofdivergence-conformingB-splines isalso crucialto avoidthe quadrature/interpolationerrorsofIBmethodsbecomingthedominantdiscretizationerror. Benchmarkandapplicationproblemsofvesicleandcapsuledynamicsaresolved,including mesh-independencestudiesandcomparisonswithothernumericalmethods.

1. Introduction

Immersedapproachesforfluid-structureinteraction(FSI)havebeenusedtotackleawiderangeofopenquestions involv-ingtheinteractionofincompressibleviscousfluidswithincompressibleelasticsolids.However, thereliabilityofimmersed approaches, such asthe immersedboundary (IB)method [2–9] andthe ficticious domain(FD) method[10–14], is often jeopardized bylarge errorsin imposingthe incompressibility constraintattheEulerianand/or Lagrangian levels[15–19]. Althoughthisissuewasfirstmentionedmorethantwodecadesago[15],solvingthisproblematitsrootandwithoutside effects thatlimit theapplicability oftheresultingimmersed methodhasprovento beextremelychallenging. Inthe case ofdiscretizations based onfinitedifferencesor finitevolumes, thedominanterror isdueto thefact that theLagrangian

*

Correspondingauthor.

E-mailaddress:hugocp@andrew.cmu.edu(H. Casquero).

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

(3)

velocity obtained fromthe Eulerianvelocity using discretizeddelta functionsis farfrom beinga solenoidal field even if theEulerianvelocity isproperlyenforcedtobe divergence-freewithrespecttothefinite-difference/finite-volume approx-imation ofthedivergenceoperator [15,20,21]. Inthecaseofdiscretizations basedon finiteelements,the dominanterror is due tothe fact that the weakly divergence-freeEulerian velocity is far frombeing a solenoidal field [22,23,18]. Most ofthe proposed solutionsto reducethe incompressibility errorsdonot tackletheaforementioned rootcausesandtryto mitigatetheireffectsinstead.Tonameafew,penaltytermsareoftenaddedtotrytodecreasethespuriouschangeoffluid volumeinsideclosedco-dimensiononesolids[24–26],apost-processingcorrectionusingaLagrangemultiplierisusedafter each timestep toslightlychangethenodalcoordinatesofclosedco-dimensiononesolidstopreserve their innervolume [27,21,8],andextremelylargegrad-divstabilizationisaddednearthefluid-solidinterfacetoobtaina velocityfieldthatis closertoasolenoidalfieldinthisregion[22,18].Therearesomerecentsolutionsthateffectivelytackletheaforementioned rootcauses,butcompromise theapplicabilityofthe resultantnumericalmethod.In[20],afinite-difference discretization isproposedthatresultsinLagrangianvelocityfieldsthataresolenoidal,butthemethodislimitedtoperiodicdomains.In [28,18],an extendedfinite-elementdiscretizationisproposed thatleads toweaklydivergence-freeEulerianvelocitiesthat are goodapproximations ofa solenoidalfield by capturing thepressure discontinuityatthe fluid-solidinterface, butthe methodhasonlybeendevelopedfortwo-dimensionalsettingsthusfar.

Sincetheadventofisogeometricanalysis(IGA)[29,30],spline-baseddiscretizationsofimmersedapproachesforFSI prob-lemshaveproliferated[31–42].ThisincludesNURBS-basedandT-spline-basedgeneralizationsoftheIBmethod[32,34] and the FDmethod[33,35,38], whichare two ofthe mostwidespreadimmersed approachesforchallenging FSI applications. Unfortunately, the issuefound in classical finite-element discretizations ofweakly divergence-free Eulerianvelocities not beingagoodapproximationtoasolenoidalfieldpersistsinspline-baseddiscretizations[33,43].Adefinitivesolutiontothis problemistouseEuleriandiscretizationsthatresultinpointwisesatisfactionoftheincompressibility constraint.However, suchdiscretizations arescarceandapriceisoftenpaidinanotheraspectofthenumericalmethodincomparisonwitha standarddiscretization.Scott-Vogeliuselements[44] wereinitiallydevelopedfortwo-dimensionalsettings.Althoughthere havebeenattemptstogeneralizeScott-Vogeliuselements tothree-dimensionalsettings[45,46],afullgeneralization isstill out ofreach [46]. In addition, the proof of inf-sup stabilityof Scott-Vogelius elements is not trivial anddifferent proofs areavailableintheliteratureunderdifferentassumptions[44,47,48].Raviart-Thomaselements[49] arenot H1_-conforming and therefore the Galerkinmethod cannot be used to discretize the Navier-Stokes equations in primal form. In[50,51], thehigherinter-elementcontinuityofsplineswasleveragedtodefineasmoothgeneralizationofRaviart-Thomaselements, theso-calleddivergence-conformingB-splines.Divergence-conformingB-splinesarepointwise divergence-free,inf-sup stable, H1-conforming,non-negative,andpressure-robust[52–54].Moreover,theirhigherinter-elementcontinuityishighly bene-ficialinIBandFDmethodsavoidingthenecessityofperformingnumericalintegrationofintegrandsthatincludeEulerian functionswith jumpsin the interiorofLagrangian elements [5,55,56,1].As a result, divergence-conforming B-splinesare an idealcandidateforthe EuleriandiscretizationofimmersedapproachesforFSI.1 _{Divergence-conforming}_B-splines_have

alreadybeenusedinageneralizationoftheFDmethodthatconsidersopenco-dimensiononesolids,theso-called immer-sogeometric method[58], andina generalizationofthe IBmethodthatconsiders co-dimensionzerosolids,the so-called divergence-conforming immersedboundary(DCIB)method[1]. In[58,1],negligibleincompressibilityerrorswere obtained atthe Eulerianlevel.In [1], itwas also shownthat thespurious changeof solidvolume was orders ofmagnitudelower thaninotherIBmethods.

The capabilityofdivergence-conformingB-splinestopreserve thefluid volumeinsidea closedco-dimensiononesolid has not beenstudied yet. Thisis neededin a variety ofFSI applications, such asstudying the behavior ofcapsules and vesicles underflow.Both capsulesandvesicles containaviscous fluidintheir interiors,buttheir wallsare verydifferent. A capsule wallis a thinpolymeric sheet,has shearresistance, andisextensible. Avesicle wallis a phospholipidbilayer, hasbendingresistance,andisinextensible.Bothcapsulesandvesiclescanbeengineeredinlaboratoriesandareoftenused in bioengineeringapplications,such asdrugdelivery [59], encapsulationofhemoglobinforartificial blood[60], and con-struction ofcell-like bioreactors [61]. Vesicles are also presentat the boundary ofbiological cells, playinga key role in intercellular communicationandin pathologicalprocesses suchas cancerandautoimmune diseases [62,63]. Capsuleand vesicleformulationsareoftencombinedtobeusedasnumericalproxiesforredbloodcells[64–66].Developing discretiza-tions of capsule and vesicle formulations as well as coupling theseformulations with a fluid solver are active fieldsof research[67,68].ThemainbenchmarkproblemconsistsinreplicatingthemotionsofvesiclesandcapsulesinCouetteflow. Thetwo principalmotionsaretanktreading(TT) andtumbling(TU).IntheTTmotion,thewallandtheinnerfluidrotate resemblingthetreadofatank.IntheTUmotion,thewallandtheinnerfluidundergoaflippingmotionresemblingarigid body.Asignificantchallengeofcapsuleandvesicleformulationsistoaccuratelycomputethehigh-orderderivativespresent intheirformulations.Insteadofusingdirectlytheformulaeofdifferentialgeometrytoevaluatethenormalvector,the cur-vature,orthesecondderivativesofthecurvature,variousworkaroundshavebeenusedsuchasspring-likediscretizations fortwo-dimensionalvesicles[69,70], C0_-continuous_triangular_meshes_using_{trigonometric}_formulae_for_{three-dimensional} vesicles [71,72],C0-continuous triangularmeshesusingquadraticinterpolationforthree-dimensionalvesicles [25,73],and

1 _In_[₅₇_],_{inf-sup stable,}_pointwise_{divergence-free,}_H1_-conforming,_and_{pressure-robust}_tetrahedral_elements_on_simplicial_{triangulations}_are_constructed. Thepressure spaceis simplythespace ofpiecewiseconstantsandthe velocityspace consistsofpiecewisecubic polynomialsenrichedwith rational functions.Althoughthesetetrahedralelementsdonothavethehigherinter-elementcontinuityofdivergence-conformingB-splines,theiruseinimmersed approachesforFSIisalsoworthconsideration.

(4)

C0-continuoustriangularmeshesforthree-dimensionalcapsules[74–76].However, thehigherinter-element continuityof splinesopensthedoortodirectevaluationoftheaforementionedquantitiesusingtheformulaeofdifferentialgeometry.

Inthiswork,theDCIBmethodisappliedtoFSIproblemsinvolvingclosedco-dimensiononesolids.Theaccuracyofthe DCIBmethodpreservingthefluidvolumewithinclosedco-dimensiononesolidsiscomparedwithrespecttoconventional IB methods [77,16] and IB methods tailoredto haveimproved volumeconservation [15,20].In orderto discretize closed curveswithupto Cp−1 _{inter-element} _continuity,_instead_of_working_with_open_knot _vectors_as_it_is_done _in_conventional IGA, wework withB-splinesof degree p defined on periodic knotvectors. In orderto discretizeclosedsurfaces withat least C1 inter-element continuity,we use analysis-suitableT-splines [78–82]. Thesespline-baseddiscretizations of closed co-dimensiononesolidsinconjunctionwithperformingintegrationbypartsenablesthecomputationoftheforcesexerted by capsules andvesicles on the fluid by directly evaluating formulae of differential geometry. We consider a Heaviside function andits associatedequation iscoupledtotheNavier-Stokesequationsinorder toconsiderinnerandouter fluids withdifferentviscosities,whichisthemoststraightforwardmannertotriggerthetransitionfromTTmotiontoTUmotion inCouetteflow.Severalquantitativecomparisonsareincludedwithrespecttotheboundaryintegralmethod(BIM)[83],an IBmethodbasedonthelattice-Boltzmannmethod(LBM)[26],andIBmethodsbasedonfinitedifferences[84,20].

Thepaperisoutlinedasfollows.Section2definesthekinematicconceptsthatareneededtohandleclosedco-dimension one solidsinthe DCIBmethod.Section 3setsforththe governingequationsforclosedco-dimensiononesolidsthat have inner and outer fluids with different viscosities. Section 4 describes the variational form of the FSI problem, which is the startingpoint toperform thediscretizationprocess. Section 5.1presentsthespline-basedspatial discretizationof the DCIBmethodforclosed co-dimensiononesolids. Section 5.2describesthefully-implicit time discretizationbasedonthe generalized-αmethod.Section5.3describestheblock-iterativesolutionstrategyusedtosolvethefinalsystemofnonlinear algebraic equations.Section6includesfivenumericalexamples.Thefirstexampleisatwo-dimensionalprobleminvolving a closed curve with active behavior. This isa benchmark problem forstudying how accurately IBmethods preserve the areaoffluidinsidetheclosedcurve.TheperformanceoftheDCIBmethodiscomparedwithotherIBmethods.Thesecond andthird examples considercommonbenchmark problemsforformulations ofvesicles andcapsules underflow, respec-tively.Quantitiesofinterest,such astheTT inclinationangle,theTUperiod,theTaylordeformation parameter,andshear stresses,arecomparedwithothernumericalmethods.TakingadvantageofthegeometricalflexibilityoftheDCIBmethod, the fourthandfifth examples study vesicle dynamicsin Taylor-Couette flowand capsule segregationin Hagen-Poiseuille flow,respectively.ConclusionsaredrawninSection7.

2. Kinematics

TheIBmethodsolvesthe(Eulerian)Navier-Stokesequationsinboththeﬂuidandsoliddomainswithadded(Lagrangian) source terms that depend onthe position ofthe solid. Inorder to havea closed systemof equations,the Navier-Stokes equationsarecoupledwiththekinematicequation thatrelatestheLagrangiandisplacementofthesolidwiththeEulerian velocityoftheNavier-Stokesequations.

Letd

= {

2

,

3

}

and

(

0

,

T

)

bethenumberofspatialdimensionsandthetimeintervalofinterest,respectively.Intherest of thissection, we deﬁne Lagrangian and Eulerianquantities that are neededto state the governingequations of theIB methodforclosedco-dimensiononesolids,withemphasisontwo-dimensionalvesiclesandthree-dimensionalcapsules.

2.1. Lagrangiandescription

We considersolidswhose kinematicsare described by closedcurvesandclosed surfacesind

=

2

,

3, respectively.The solid occupies,attime t, theregion

t

⊂ R

d,which can beobtained astheimage ofa referenceconﬁguration

R

⊂ R

d

throughthemapping

ϕ

:

R

× (

0

,

T

)

→

t.Let X

∈

R beamaterialparticle.TheLagrangianunknownofthemathematical

modelistheLagrangiandisplacementu

:

R

× (

0

,

T

)

→ R

d,verifying

ϕ

=

X

+

u.Thedeformedconﬁguration,thereference

conﬁguration,andtheLagrangiandisplacementcanbeparametrizedinspaceusingd

−

1 parametriccoordinates.

2.1.1. Closedcurve

Let

ξ

L be theparametriccoordinate.Whenworkingwithtwo-dimensionalvesicles,itiscommontoreparametrizethe closedcurveintermsofitsarclength s,whichisdonetakingintoaccountthat

ds

=

d

ϕ

d

ξ

L

d

ξ

L, (1)

where

||

· ||

denotesthelengthofavector.Usingthearclengthastheparametriccoordinate,theunittangentvectortothe closedcurveisobtainedby

t

=

d

ϕ

ds. (2)

(5)

C =

det

⎛

⎜

⎝

d

ϕ

x ds d

ϕ

y ds d2

ϕ

x ds2 d2

ϕ

y ds2

⎞

⎟

⎠

. (3)

Thesignof

C

dependsontheorientationofthecurveprovidedbytheparametrization.Wealwaysuseclosedcurveswhose parametriccoordinatemovescounterclockwise.Asaresult,acirclehaspositivecurvature.Theunit outwardnormalvector totheclosedcurveisobtainedby

n

=

d

ϕ

y ds

,

−

d

ϕ

x ds

. (4) 2.1.2. Closedsurface

Let

ξ

L

= {ξ

₁L

,

ξ

₂L

}

be theparametric coordinates.Inthefollowing,indicesinGreek letterstakethevalues1

,

2 and sum-mationoverrepeatedGreekindicesisimplied.Non-unittangentvectorstotheclosedsurfaceareobtainedby

aα

=

∂

ϕ

∂ξ

L

α

. (5)

Theunitoutwardnormalvectortotheclosedsurfaceisobtainedby

n

=

a1

×

a2

||

a1

×

a2

||

, (6)

wheretheparametriccoordinates

ξ

L

1 and

ξ

2L havebeenchoseninsuchawaythata1

×

a2 pointsintheoutwarddirection. Thecovariantmetriccoeﬃcientsoftheclosedsurfacearedeﬁnedas

aαβ

=

aα

·

aβ. (7)

The contravariant metric coeﬃcients can be computed as the inverse matrix of the covariant coeﬃcients, i.e.,

aαβ

₌

aαβ

₋1

. The contravariantmetric coeﬃcientsare used to obtain the contravariantbase vectorsfrom the covariant base vectors,namely,aα

=

aαβ_a

β.Thecovariantcurvaturecoeﬃcientsoftheclosedsurfacearedeﬁnedas

bαβ

=

∂

a_α

∂ξ

_βL

·

n. (8)

TheChristoffelsymbolsaredeﬁnedas

γ_α_β

=

∂

aα

∂ξ

_βL

·

a

γ_. ₍₉₎

ThequantitiesinEqs.(5)-(9) aredefinedwithrespecttothedeformedconfiguration,butanalogousquantitiesaredefined withrespecttothereferenceconfigurationanddenotedby ˚aα, ˚n, ˚aαβ, ˚aαβ, ˚aα , ˚bαβ,and

˚γαβ.

2.2. Euleriandescription

Let

⊂ R

d bethetime-independentregioncontainingboththeﬂuidandthesolidandletx

∈

bea spatialposition. Forthesakeofbrevityandsinceitholdstrueinalltheexamplesofthispaper,weassumethatthesolidisfullyimmersed in the ﬂuid, i.e.,

t

∩ ∂

= ∅

∀

t

∈ (

0

,

T

)

. The Eulerianvariables of the mathematical modelare the Eulerian velocity v

:

× (

0

,

T

)

→ R

d,thepressurep

:

× (

0

,

T

)

→ R

,andtheHeavisidefunctionH

:

× (

0

,

T

)

→ R

(theHeavisidefunctionis onlyneededwhenﬂuidswithdifferentinnerandouterviscositiesareconsidered).Thephysicaldomain

andtheEulerian variablescanbeparametrizedinspaceusingd parametriccoordinates.

3. Governingequations

Atthecontinuous level,themathematicalmodeloftheIBmethodisequivalenttotheboundary-ﬁttedFSIformulation [5,55]. The strong form of the IB method forclosed co-dimension one solids with different inner and outer viscosities is stated as follows: Given

ρ

∈ R

+_,

_μ

_o

_{∈ R}

+_,

_μ

_i

_{∈ R}

+_, _g

V

∈

× (

0

,

T

)

→ R

d, v0

:

→ R

d, u0

:

R

→ R

d, and vB

:

∂

× (

0

,

T

)

→ R

d_,_ﬁnd_v

_:

_{× (}

₀

_,

_T

₎

_{→ R}

d_, _p

_:

_{× (}

₀

_,

_T

₎

_{→ R}

_, _H

_:

_{× (}

₀

_,

_T

₎

_{→ R}

_,_and_u

_:

R

× (

0

,

T

)

→ R

d,suchthat,

ρ

∂

v

∂

t

x

+ ∇

x

· (

v

⊗

v

)

= ∇

x

· (

2

μ

∇

xsymv

)

− ∇

xp

+

gV

+

t f

δ(

x

−

ϕ

(

X

,

t

))

d

t in

× (

0

,

T

)

, (10)

∇

x

·

v

=

0 in

× (

0

,

T

)

, (11)

(6)

∂

u

∂

t

X

=

v

δ(

x

−

ϕ

(

X

,

t

))

d

in

t

× (

0

,

T

)

, (12)

xH

= −∇

x

·

t n

δ(

x

−

ϕ

(

X

,

t

))

d

t in

× (

0

,

T

)

, (13)

μ

=

μ

iH

+

μ

o

(

1

−

H

)

in

× (

0

,

T

)

, (14) v

=

v0 on

× {

0

}

, (15) u

=

u0 on

0

× {

0

}

, (16) v

=

vB on

∂

× (

0

,

T

)

, (17) H

=

0 on

∂

× (

0

,

T

)

, (18)

where

ρ

isthedensity,

μ

i and

μ

o are theinner andouter ﬂuidviscosities, respectively, gV isanexternal force perunit

of volumeacting onthe system,

∇

symx

(

·)

is thesymmetric gradient operatordeﬁned by

∇

sym

x v

= (∇x

v

+ ∇x

vT

)/

2, v0 is theinitialvelocity, u0 istheinitialdisplacement, vB istheboundaryvelocity,

δ(x

−

ϕ

(X,

t

))

=

di=1

δ(

xi

−

ϕ

i

(X,

t

))

isthe

d-dimensionalDiracdeltafunction, and f

:

R

× (

0

,

T

)

→ R

d is theforce exerted bythe solidon theﬂuid. Thenotation

∂

v

∂

t

x and

∂u

∂

t

X

indicatesthatthetimederivativeistakenholdingx andX ﬁxed,respectively.

Eqs. (10)-(12) represent the linear momentum balanceequation, the mass conservation equation, andthe kinematic equation that relatestheLagrangiandisplacement withthe Eulerianvelocity,respectively. Eqs.(13)-(14) areaddedto the standard IB formulation when

μ

i

=

μ

o. The derivation of Eqs. (13)-(14) can be found in [85]. For brevity and since it

holdstrueinall theexamplesofthispaper,weassume thatthe innerﬂuid,theouter ﬂuid,andthesolid havethe same density

ρ

.InnerandouterﬂuidswithdifferentdensitiescanbeconsideredusingtheHeavisidefunction H akintothecase

μ

i

=

μ

o.Aclosedco-dimensiononesolidwithdifferentdensitythantheﬂuid canbeconsideredinananalogousmanner

aswe considered co-dimension zerosolidswithdifferent densitythan thefluid in [1] provided that a solid formulation with a consistent time evolution ofits thicknessis used. Eqs. (15)-(16) define the initial condition for the velocity and the displacement, respectively. Eqs. (17)-(18) definethe boundary condition for the velocity andthe Heaviside function, respectively.ForbrevityDirichletboundaryconditionsforthevelocityareappliedonthewholeboundaryinSections3,4, and5,butNeumannandperiodicboundaryconditionscanbeappliedintheDCIBmethodfollowingthestandardprocedures ofvariationalmethodsasdoneinSection6.

Theexpressionsof f fortwo-dimensionalvesiclesandthree-dimensionalcapsulesaredetailedbelow.

3.1. Two-dimensionalvesicles

Asderivedin[86],theLagrangianforceperunitlengthexertedbythevesicleontheﬂuidis

f

=

κ

d 2

_C

ds2

+

κ

C

3 2

− Cζ

n

+

d

ζ

dst

,

(19)

where

κ

isthebendingrigidityofthevesicle,

ζ

isaLagrangemultiplierthatenforcesthevesicledeformationstobelocally inextensible,andanullspontaneouscurvatureisconsidered.NotethatEq.(19) dependsonthedeformedconﬁguration,but notthereferenceconﬁguration.Inordertocircumventnumericalinstabilities,wereplace

ζ

bythefollowingexpression

ζ

=

4CI

λ

2

−

1

, (20)

whereCI isthedilatationmodulusand

λ

= ||

d

ϕ

(ξ

L

,

t

)/

d

ξ

L

||/||

d

ϕ

(ξ

L

,

0

)/

d

ξ

L

||

isthestretchratio.Eq.(20) isequivalent to

consideringacapsulewithperimeterstrain-energyfunctionorHelmholtzfreeenergygivenby

WI

=

CI

λ

2

−

1

2

. (21)

Eq.(21) isatwo-dimensionalanalogtothedilatationalcomponentofthesurfacestrain-energyfunctionproposedin[87].

3.2. Three-dimensionalcapsules

Asderivedin[88],theLagrangianforceperunitareaexertedbythecapsuleontheﬂuidis

f

=

∂

Tαβ

∂ξ

L α

+

α αλTλβ

+

β αλTαλ

aβ

+

Tαβbαβn, (22) with

(7)

Tαβ

=

2 Js

∂

Ws

∂

I1 ˚aαβ

₊

_{2 J} s

∂

Ws

∂

I2 aαβ, (23) I1

=

˚aαβaαβ

−

2, (24)

I2

=

det

(

˚aαβ

)

det

(

aαβ

)

−

1, (25)

Js

=

I2

+

1, (26)

where T αβ _are _the _{contravariant} _coeﬃcients_of _the _membrane _forces _in _the _deformed _{conﬁguration,} _W

s is the surface

strain-energyfunction,I1 andI2 aretheinvariantsofthedeformation,and Jsistheratioofdeformedlocalsurfaceareato

referencelocalsurfacearea.Notethat Eq.(22) dependsonboth thedeformedandreferenceconfigurations.Eqs.(22)-(26) define a membrane formulation, i.e., transverse shear andbending are neglected. This membraneformulation takesinto account bothgeometric andmaterialnonlinearities. Examplesofconstitutiveequationsare theneo-Hookeanlawandthe Skalaklaw[87] definedas

Ws

=

Gs 2

I1

−

1

+

1 I2

+

1

, (27) Ws

=

Gs 4

(

I 2 1

+

2I1

−

2I2

)

+

CI 4

(

J 2 s

−

1

)

2, (28)

respectively, where Gs is the surface shear modulus. When CI

Gs,the Skalak law leads to locally inextensible

mem-branes, butitcanbeusedtomodelother typesofmembraneswhenCI

Gs.The neo-HookeanandSkalaklawsresultin

signiﬁcantly differentmaterialresponsesunderlarge deformations,e.g., underuniaxialextension,the neo-Hookeanlawis strain-softeningwhereastheSkalaklawisstrain-hardening[89].

4. Variationalformulation

Givensuitable trial solutionspaces (

S

v,

S

p,

S

u,and

S

H) andweighting function spaces(

V

v,

V

p,

V

u,and

V

H), the

variationalformoftheIBmethodisstatedasfollows:Findv

∈ S

v,p

∈ S

p,u

∈ S

u,andH

∈ S

H,suchthat,

B

((

w

,

q

,

s

,

M

), (

v

,

p

,

u

,

H

))

−

L

(

w

)

=

0

∀(w

,

q

,

s

,

M

)

∈ V

v

× V

p

× V

u

× V

H, (29) with B

((

w

,

q

,

s

,

M

), (

v

,

p

,

u

,

H

))

=

w

,

ρ

∂

v

∂

t

x

− (∇

xw

,

ρ

v

⊗

v

)

− (∇

x

·

w

,

p

)

+

∇

xw

,

2

μ

o

∇

xsymv

+

∇

xw

,

2H

(

μ

i

−

μ

o

)

∇

xsymv

+ (

q

,

∇

x

·

v

)

− (

w

◦

ϕ

,

f

)

_t

+

s

,

∂

u

∂

t

X

−

v

◦

ϕ

t

+ (∇

xM

,

∇

xH

)

+ (∇

xM

◦

ϕ

,

n

)

t, (30) L

(

w

)

=

w

,

g_V

, (31)

where

(

·,

·)

and

(

·,

·)

t denotetheL

2 _inner_product_over_the_domains

_and

t,respectively.Asmathematicallyshownin

[5,6],thevariationalformulationoftheIBmethodenablestheeliminationoftheDiracdeltafunctions.NotethatEq.(13) is notconsideredin[5,6],butitsDiracdeltafunctioniseliminatedinthesameway.

4.1. Two-dimensionalvesicles

In this case, f contains a term involvingone factorwith fourth-order derivatives (d2

_C/

_ds2₎ _and _another_factor _with second-order derivatives (n). Therefore, integration by parts in the variational form can be applied once to remove the fourth-orderderivativesfromthevariationalform,namely,

t

κ

d 2

_C

ds2n

· (

w

◦

ϕ

)

ds

= −

t

κ

d

C

ds dn ds

· (

w

◦

ϕ

)

ds

−

t

κ

d

C

dsn

·

d

(

w

◦

ϕ

)

ds ds, (32)

wherewehaveassumedthat

κ

isaconstant.Asaresult,adiscretizationoftheclosedcurvewithC2 _{inter-element} conti-nuityisneededtodirectlycomputetheterm

(w

◦

ϕ

,

f

)

t.

(8)

4.2. Three-dimensionalcapsules

Inthiscase, f containsa terminvolving onefactorwithsecond-order derivatives(bαβ) andanotherfactorwith ﬁrst-order derivatives (T αβ). Therefore, integration by parts cannot be usedto decrease the order of thesederivatives anda discretizationoftheclosedsurfacewithC1 _{inter-element}_continuity_is_needed_to_directly_compute_the_term

_(w

_◦

_ϕ

_,

_f

₎

t.

5. Discretization

ThissectiondescribestheDCIBmethodforclosedco-dimensiononesolids.TheDCIBmethoddiscretizesthevariational formulationdeﬁnedinEqs.(29)-(31).Thus,devisingsmeareddeltafunctionsisnotneededintheDCIBmethod.Bycontrast, IBmethodsthatdiscretizethestrongformdeﬁnedinEqs.(10)-(18) requiresmeareddeltafunctions.

5.1. Spatialdiscretization

InadditiontothediﬃcultyofaccuratelyimposingtheincompressibilityconstraintinIBmethods,anadditionalchallenge istoaccurately computeintegralsinwhichbothEulerianandLagrangianfunctionsareinvolved.Thespatial discretization oftheDCIBmethodhasbeenpurposefullydesignedtotacklethesetwochallenges.

5.1.1. Closedcurves

We discretizeclosedcurvesusing B-splinesofdegree p with periodic knotvectorsandglobalCp−1 _continuity._Let_us startdeﬁningthemapping FL

: [

0

,

1

]

→

R usingB-splinesasfollows

FL

ξ

L

=

nL

C=1 PL_C

N_CL

ξ

L

,

ξ

L

∈ [

0

,

1

]

, (33)

wherethesuperscript L standsfor“Lagrangian”,

{

PL C

}

n

L

C=1 are thecontrolpoints thatdeﬁnethereferenceconﬁgurationof thesolid,and

{

N_CL

}

_CnL₌₁ areuni-variateB-splinebasisfunctions.

Uni-variateB-splinebasisfunctionsarecomputedfromaknotvector

usingtheCox-deBoorrecursionformula[90,29]. Aknotvector isaﬁnitenon-decreasingsequenceofrealnumbers

= {ξ

1

,

ξ

2

,

...,

ξ

nL₊_p₊₁

}

,where

ξ

i isthe i-thknot, p is

thepolynomial degree,andnL is thenumberofuni-variateB-splinebasis functions.The parametricdomain ofthecurve is deﬁnedby the interval

[ξ

p+1

,

ξ

nL₊₁

]

(we impose

ξ

p+1

=

0 and

ξ

nL₊₁

=

1). A knotspan isthe difference betweentwo

consecutiveknots(

ξ

i

= ξ

i+1

− ξ

i).Insidenonzeroknotspans,B-splinebasisfunctionsarepolynomialsofdegree p.Since

theknotscan berepeated,wedeﬁne asequenceofknotswithoutrepetitions(calledbreakpoints)

η

= {

η

1

,

η

2

,

...,

η

m

}

and

anothersequencewiththeknotmultiplicitiesr

= {

r1

,

r2

,

...,

rm

}

.Knotmultiplicity controlsthecontinuityofB-splinebasis

functionsatbreakpoints,namely,B-splinebasisfunctionshave

α

i

=

p

−

ricontinuousderivativesat

η

i.

Anopen (orclamped)knotvectorisobtainedwhenr1

=

rm

=

p

+

1.Openknotvectorsfacilitatethestrongimposition

ofDirichletboundaryconditions,whichhasmadeopenknotvectorsthestandardchoiceinIGA.However,whenusingopen knotvectors, imposingperiodicboundaryconditionsrequiresbuildingasystemofconstraints.Theseconstraintsgetmore complicatedasthecontinuityattheperiodicboundaryisincreased(see[91] foradescriptionofhowtoimposeC1 period-icitywithopenknotvectors).AnalternativeistobuildtheperiodicityintotheB-splinespace.Thisisaccomplishedusing periodic(oruniform)knotvectorsinsteadofopenknotvectors.Periodicknotvectorshavenorepeatedknotsandtheknots are evenly distributed. Reference[92] buildsthe geometry usingopen knotvectors. After that, the knot vectorsare un-clamped,i.e.,transformedintoperiodicknotvectors,whichrequiresmodifyingthecontrolpointstopreservethegeometry. ThisstrategyisappropriatetoimposeCp−1 _periodicity_in_geometries_that_are_not_closed,_e.g.,_imposing_Cp−1_periodicity_to thevelocity andthepressureina straight channel.However, whenthegeometryisclosedand Cp−1 _periodicity_must_be imposedtoboththegeometryandtheunknownsoftheproblem, periodicknotvectorsmustbe usedfromthebeginning. This isthe strategy thatwe follow hereforthe closedcurve andthe Lagrangian displacement.Toobtain Cp−1 periodic-ity, thelast p control points mustbe wrappedwiththe ﬁrst p controlpoints, namely, PL_nL₋_p₊₁

=

P

L 1, P L nL₋_p₊₂

=

P L 2,..., P_nLL

=

P L

p.AnexampleofaB-splineclosedcurveisshowninFigs.1and2.

Invokingtheisoparametricconcept,theLagrangiandisplacementisparameterizedasfollows

uhL

ξ

L

,

t

=

nL

C=1 uC

(

t

)

NLC

ξ

L

,

ξ

L

∈ [

0

,

1

]

, (34)

whereuC

(

t

)

arethecontrolvariablesoftheLagrangiandisplacement.Asaresult,thedeformedconﬁgurationis

parameter-izedasfollows

ϕ

hL

ξ

L

,

t

=

nL

C=1 P_CL

N_CL

ξ

L

+

nL

C=1 uC

(

t

)

NCL

ξ

L

,

ξ

L

∈ [

0

,

1

]

. (35)

(9)

Fig. 1. Giventheperiodicknotvector= {−0.6,−0.4,−0.2,0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6}andp=3,theB-splinebasisfunctions{NL C}8C=1 are plotted.

Fig. 2. TwoB-splinecurvesareplottedusing= {−0.6,−0.4,−0.2,0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6}andp=3.(a)Anopencurveinwhichallthe controlpointshavedifferentcoordinates.(b)AC2_-continuous_closed_curve_obtained_by_wrapping_the_last₃_control_points_with_the_ﬁrst₃_control_points, i.e.,PL

6=PL1,PL7=PL2,andPL8=P3L.ThecontrolpointsaredenotedbyredsquaresandtheB-splinecurvesarecoloredinblue.(Forinterpretationofthe colorsintheﬁgure(s),thereaderisreferredtothewebversionofthisarticle.)

Thenonzeroknotspansdeﬁnetheelementsofamesh

M

L _in

_[

₀

_,

₁

_]

_,_this_mesh_is_called_the_Lagrangian_mesh_from_now_on.

TheLagrangianmeshcanbepushed forwardtothereferenceconﬁgurationandthedeformedconﬁgurationusingEq.(33) andEq.(35),respectively.Alltheintegralsposedonthedomain

h_tL

=

ϕ

hL

₍

_[

₀

_,

₁

_],

_t

₎

_are_computed_using_Gauss_quadrature

intheelementsoftheLagrangianmesh,namely,p

+

1 quadraturepointsareusedperelement.UsingtheBubnov-Galerkin method,thediscretetrialandweightingfunctionspacesfortheLagrangiandisplacementaredeﬁnedasfollows

S

hL u

=

uhL

|

uhL

◦

FL

=

uhL

,

uhL

(

·,

t

)

∈

span

{

NL_C

(ξ

L

)

}

n_CL₌₁

, (36)

V

hL u

=

shL

|

shL

◦

FL

=

shL

,

shL

(

·,

t

)

∈

span

{

N_CL

(ξ

L

)

}

n_CL₌₁

. (37) 5.1.2. Closedsurfaces

Wediscretizeclosedsurfacesusingbi-cubicanalysis-suitableT-splines(AST-splines)withatleastC1 inter-element con-tinuity. AST-splinesareunstructuredandeightextraordinary pointswithvalence threeareenough tobuildsimpleclosed surfacessuchasasphereinthe“soccerball” parameterization.AdetailedexplanationofhowtoconstructtheAST-spline surfacesthatweuseinthisworkcanbefoundin[81].

5.1.3. Navier-Stokesequations

The spatial discretization of the Eulerian velocity, the pressure, and their weighting functions is performed using divergence-conformingB-splines,whicharesmoothgeneralizationsofRaviart-Thomasmixedﬁniteelements.

(10)

Letusstartbydeﬁningthemapping FE

:

→

hE usingnon-uniformrationalB-splines(NURBS)asfollows FE

ξ

E

=

nE

i=1 P_iE wi

NiE

ξ

E

nE j=1wj

NEj

ξ

E

,

ξ

E

∈

, (38)

where the superscript E stands for “Eulerian”,

= [

0

,

1

]

d_,

_{

_PE i

}

n

E

i=1 are the control points that deﬁne the geometry of the physicaldomain,

{

wi

}

n

E

i=1 arethe weightsassociated withthecontrol points,and

{

NiE

}

n

E

i=1 area setofd-variatebasis functions. Thed-variateB-splinebasisfunctionsareobtainedasthetensorproductofuni-variateB-splinebasis functions [90,29].NURBSbasisfunctionsreducetoB-splinebasisfunctionswhenwi

=

1

∀

i

∈

1

,

2

,

...,

nE sinceB-splinebasisfunctions

formapartitionofunity.Themainreasontoconsiderweights wi

=

1 istoexactlyrepresentcertaingeometriessuchasa

circularcylinder.

Thenonzeroknotspansofthed knotvectorsusedtoobtainthebasisfunctions

{

N_iE

}

n_i₌E₁ deﬁnetheelementsofamesh

M

E _in

_; _this_mesh_is_called_the_Eulerian_mesh_from_now_on._The_Eulerian_mesh_can _be_pushed _forward_to _the_physical

domain

hE _using_Eq.₍₃₈_)._All_the_integrals_posed_on_the_domain

hE _are_computed_using_Gauss_quadrature_in_the_elements

oftheEulerianmesh.

ThediscretetrialandweightingfunctionspacesforthevelocityandthepressurearedeﬁnedusingthePiolaand integral-preservingtransformations,respectively,asfollows

S

hE v

=

vhE

|

vhE

◦

FE

=

D F E

_vhE det

(

D FE

)

,

v hE

₍

_·,

_t

₎

_∈

_{V E L}hE

_,

_vhE

_·

_nhE

₌

_v B

·

nh E on

∂

× (

0

,

T

)

, (39)

S

hE p

=

⎧

⎪

⎨

⎪

⎩

p hE

_|

_phE

_◦

_FE

₌

ph E det

(

D FE

)

,

p hE

₍

_·,

_t

₎

_∈

_{P R E}hE

_,

phEd

=

0

⎫

⎪

⎬

⎪

⎭

, (40)

V

hE v

=

whE

|

whE

◦

FE

=

D F E

_whE det

(

D FE

)

,

w hE

(

·) ∈

V E Lh E

,

whE

·

nhE

=

0 on

∂

× (

0

,

T

)

, (41)

V

hE p

=

qhE

|

qhE

◦

FE

=

q hE det

(

D FE

)

,

q hE

(

·) ∈

P R Eh E

, (42) with

V E Lh E

=

S

k+1,k k,k−1

(

M

E

)

× S

k,k+1 k−1,k

(

M

E

)

if d

=

2,

S

k+1,k,k k,k−1,k−1

(

M

E

)

× S

k,k+1,k k−1,k,k−1

(

M

E

)

× S

k,k,k+1 k−1,k−1,k

(

M

E

)

if d

=

3, (43)

P R Eh E

=

S

k,k k−1,k−1

(

M

E

)

if d

=

2,

S

k,k,k k−1,k−1,k−1

(

M

E

)

if d

=

3, (44)

where D FE isthegradientofthemapping FE,

nhE istheunitoutwardnormalto

,and

S

p1,p2,...,pd

α1,α2,...,αd

(

M

E

₎

_is_the_d-variate

B-splinespaceofbasisfunctionsdeﬁnedon

M

E _that,_in_direction_i,_have_degree_p

i andC αi continuityatallinteriorknots.

InEqs. (39)-(44), wedeﬁneddivergence-conforming B-splinespaces ofmaximalcontinuitysince thesewill bethespaces usedthroughouttheexamplesofthispaper.NotethatonlythenormalDirichletboundaryconditionhasbeenimposedon the velocity discrete space, the tangential Dirichlet boundaryconditions will be imposed weakly usingNitsche’s method inthe semi-discrete form.Letusdeﬁne thebasis functions

{

N_vE

l,Al

}

nvl Al=1 and

{

N E p,B

}

n p

B=1 such thatspan

{

NvE1,A1

(ξ

E

₎

_}

nv1 A1=1

×

...

×

span

{

NE vd,Ad

(ξ

E

)

}

n_Avd_d₌₁

=

V E Lh E andspan

{

NE p,B

(ξ

E

)

}

n_Bp₌₁

=

P R Eh E

,wherenvl _and_np _are_the_total_number_of_degrees_of

freedomthatthel-thcomponentofthevelocityandthepressurehave,respectively.

The above choicesof spaces areinf-sup stableandthe discrete function spacesfor thevelocity and thepressure are

H1_-conforming_and_L2_-conforming,_{respectively,}_for_k

_≥

_1._In_addition,_imposing_the_discrete_Eulerian_velocity_to_be_weakly divergence-freeresultsinasolenoidalﬁeld,i.e.,

qhE

,

∇

x

·

vh E

hE

=

0

∀

q hE

_{∈ V}

hE p

=⇒ ∇

x

·

vh E

=

0

∀

x

∈

hE

,

∀

t

∈ [

0

,

T

]

, (45) as mathematically shownin [52,54]. Thatis, weak incompressibility impliesstrong (i.e.,pointwise) incompressibility for divergence-conformingB-splines.

(11)

5.1.4. Heavisidefunction

ThespatialdiscretizationoftheHeavisidefunctionanditsweightingfunctionisperformedusingB-splines.Thediscrete trialandweightingfunctionspacesfortheHeavisidefunctionaredeﬁnedas

S

hE H

=

HhE

|

HhE

◦

FE

=

HhE

,

HhE

(

·,

t

)

∈

H E Ah E

,

HhE

=

0 on

∂

× (

0

,

T

)

, (46)

V

hE H

=

MhE

|

MhE

◦

FE

=

MhE

,

MhE

(

·,

t

)

∈

H E Ah E

,

MhE

=

0 on

∂

× (

0

,

T

)

, (47) with

H E Ah E

=

S

k,k k−1,k−1

(

M

E

)

if d

=

2,

S

k,k,k k−1,k−1,k−1

(

M

E

)

if d

=

3. (48)

Letusdeﬁne thebasis functions

{

N_HE_,_D

}

n_DH₌₁ such that span

{

NE_H_,_D

(ξ

E

)

}

n_DH₌₁

=

H E Ah

E

,wherenH _is_the _number_of_degrees

of freedom that the discrete Heaviside function has.As mentioned in [85], some overshootsand undershoots may take placeneartheinterface. Whenthevalue of HhE

ata quadraturepointisgreater thanoneorlessthanzero,thisvalueis substitutedbyoneandzero,respectively.

5.1.5. Semi-discreteform

The semi-discrete form ofthe DCIB method for closedco-dimension one solids is stated asfollows: Find vhE

∈ S

hvE,

phE

∈ S

hE p ,uh L

∈ S

hL u ,andHh E

∈ S

hE

H,suchthatforall wh

E

∈ V

hE v ,qh E

∈ V

hE p ,sh L

∈ V

hL u ,andMh E

∈ V

hE H B

(

whE

,

qhE

,

shL

,

MhE

), (

vhE

,

phE

,

uhL

,

HhE

)

−

b

whE

,

vhE

−

L

whE

+

l

whE

=

0, (49) with b

whE

,

vhE

=

F∈∂ hE

F 2

μ

o

(

wh E

)

_||

· (∇

xsymvh E nhE

)

d

∂

−

F∈∂ hE

F

ρ

(

whE

)

_||

·

vhE

(

vB

·

nh E

)

₊d

∂

+

F∈∂ hE

F 2

μ

o

(

∇

xsymwh E nhE

)

· (v

hE

)

_||d

∂

−

F∈∂ hE

F 2

μ

o Cpen hF

(

whE

)

_||

· (v

hE

)

_||d

∂

, (50) l

whE

=

F∈∂ hE

F

ρ

(

whE

)

_||

·

vB

(

vB

·

nh E

)

₋d

∂

+

F∈∂ hE

F 2

μ

o

(

∇

xsymwh E nhE

)

· (

vB

)

||d

∂

−

F∈∂ hE

F 2

μ

o Cpen hF

(

whE

)

_||

· (v

B

)

||d

∂

, (51)

where theterms b

whE

,

vhE

andl

whE

havebeenadded to thesemi-discrete formto weakly imposethe tangential DirichletboundaryconditionsusingNitsche’smethod,nhE istheoutwardunitnormalvectorto

hE,

(

·)||

= (·)

−((·)

·

nhE

)n

hE isthevectortangentialcomponent,

(v

B

·

nh

E

)

₊

=

vB

·

nh E ifvB

·

nh E

>

0 and0 otherwise,

(v

B

·

nh E

)

₋

=

vB

·

nh E ifvB

·

nh E

≤

0 and0 otherwise,hF isthemeshsizeinthedirectionnormaltothefaceF ,andCpenistheNitsche’spenalizationparameter.

AllthenumericalresultsofthispaperusethevalueCpen

=

5

(

k

+

1

)

asproposedin[52].

Note that the integrals ofIB methods posed on thesolid domain includeEulerian functionssuch as vhE

, whE

, MhE_,

andtheirﬁrstderivatives.InordertoevaluatetheseEulerianfunctionsataquadraturepointwithgivenparametric coordi-nates

ξ

_GL intheLagrangianmesh,wefollowtwosteps.First,wecomputethephysicalcoordinatesxG

∈

h

E

associatedwith theparametric coordinates

ξ

L_G asxG

=

ϕ

h

L

(ξ

LG

,

t

)

.Then, wecompute theparametric coordinates

ξ

EG inthe Eulerianmesh

associated withxG as

ξ

EG

= (

FE

)

−1

(x

G

)

. The NURBS mapping FE can be invertedanalytically in manycases ofpractical

interest andwhen that isnot thecase, the mappingis invertedsolving a d

×

d system ofnonlinear algebraic equations. Once

ξ

E_G isknown,Eulerianfunctionscanbeevaluatedusingstandardprocedures.Thus,nointerpolationorapproximation hasbeenusedtoevaluatetheEulerianfunctionsatthequadraturepointsoftheLagrangianmesh.However,thequadrature errors oftheintegralsposedon

htL are largerthan instandardvariational methodssince theelementboundariesofthe

Eulerianmesh(wheretheEulerianfunctionshavereducedregularity)oftenintersecttheinterioroftheintegrationregions beingused(the elementsofthe Lagrangianmesh).As opposedtostandardﬁniteelements, wheretheinter-element con-tinuity of vhE

, whE

, MhE

is C0 _(and_its _ﬁrst_derivatives_have_jumps),2 _{divergence-conforming} _B-splines_enable_us_to_raise

2 _In_[₅₆_],_to_avoid_performing_numerical_integration_of_integrals_that_include_Eulerian_functions_with_jumps_inside_the_integration_regions,_vhE andwhE areprojectedintotheLagrangianmeshbyevaluatingtheEulerianfunctionsatthenodesoftheLagrangianmeshtousethosevaluesasnodalvaluesof