Bildbasierte L osung von Partiellen Di …...Konstruktion von CFE-R aumen Assemblierung,...

EinleitungDiskretisierung mit CFEs

Implementierung und numerische Ergebnisse

Bildbasierte Losung vonPartiellen Differentialgleichungenmit Composite Finite Elements

Vortrag im Diplomandenseminar 2010 / 2011

Sebastian Westerheide

[email protected]

31.01.2011

Sebastian Westerheide Bildbasierte Losung von PDEs mit Composite F. Elements



Ubersicht

1 EinleitungMotivationAnsatze fur ein numerisches Losungsverfahren

2 Diskretisierung mit CFEsGrundlagenKonstruktion von CFE-RaumenAssemblierung, Lokalisierung

3 Implementierung und numerische ErgebnisseImplementierungNumerische Ergebnisse




MotivationAnsatze fur ein numerisches Losungsverfahren

Motivation (1/3)

Motivation:

Die mathematische Modellierung und Simulation vielerphysikalischer Prozesse basiert auf dem Losen von partiellenDifferentialgleichungen (PDEs)

Zunehmende Popularitat von Anwendungen aus dem Bereichder Biologie und Medizin

z.B. makroskopische Simulationen, die Radiologen undChirurgen mit Informationen zur Behandlungsplanungversorgen





Motivation (2/3)

Dabei:

Gewinnung des Rechengebiets haufig durch bildgebendeVerfahren

Typische Vorgehensweise:CT

Segmentierung

Simulation

Die untersuchten Objekte enthalten oftmals eine große Anzahlan geometrischen Details





Motivation (3/3)

Diese mussen auch in makroskopischen Simulationenberucksichtigt werden, da das Auflosen der Mikroskalazu verlasslicheren Ergebnissen fuhrt

=⇒ Das Rechengebiet besitzt haufig eine extremkomplexe Form

Also benotigt: Verfahren zur numerischen Losung von PDEs aufpotentiell komplizierten Gebieten beliebiger Dimension, die durchsegmentierte Bilddaten gegeben sind.





Ansatze fur ein numerisches Losungsverfahren (1/2)

Ansatz 1:

Grundsatzlich ist die Finite Elemente Methode (FEM) eingeeignetes Verfahren

In Verbindung mit einem auf Bilddaten zugeschnittenenGittergenerator

Aber: Klassische FEM basiert auf einem Gitter, welches dasRechengebiet in einfache geometrische Objekte zerlegt

Problem:

Auflosung des Randes ⇒ extrem feine Gitter beiRechengebieten mit einer komplexen Form

=⇒:

Hohe Anzahl von FreiheitsgradenDiskretisierung auf einer rein groben Skala nicht moglich





Ansatze fur ein numerisches Losungsverfahren (2/2)

Ansatz 2:

Konstruktion von FE-Raumen unter Benutzung eines von derGeometrie unabhangigen Gitters

Dazu: Konstruiere Basisfunktionen so, dass sie zum Rand desRechengebiets passen und gegebene Randbedingungen erfullen

Wir wollen uns mit sogenannten Composite Finite Elements(CFEs) bzw. CFE-Raumen beschaftigen




GrundlagenKonstruktion von CFE-RaumenAssemblierung, Lokalisierung

Ubersicht








Modellproblem

Als Modellproblem betrachten wir dieDifferentialgleichung

−∇ · (A∇u) + cu = f in Ω,

zusammen mit der Randbedingung

u = gD auf ΓD,

(A∇u) · n = gN auf ΓN .

A : Ω→ Rd×d

c : Ω→ R

f : Ω→ R

gD : ΓD → R

gN : ΓN → R

Ω ⊂ Rd beschranktes Gebiet mit Lipschitz-Rand ∂Ω, d ∈ N∂Ω = ΓD ∪ ΓN

∃a0 ∈ R+ mit A(x)− a0I positiv definit ∀x ∈ Ω

c > 0 zur Vereinfachung der Darstellung





Schwache Losung

Die Bestimmung einer schwachen Losung fuhrt zu einemVariationsproblem: Finde ein u ∈ X, so dass

B(u, v) = `(v) ∀v ∈ X

Mit einer Bilinearform B : X ×X → R und einem linearenFunktional ` : X → R, welches von der Randbedingungabhangt

Je nach Randbedingung:

X = H1(Ω) oder X = H10 (Ω) :=

v ∈ H1(Ω) | v|ΓD

= 0

Existenz und Eindeutigkeit nach dem Lax-Milgram Lemma





Numerische Approximation

Grundlage fur eine numerische Approximation der Losungdieses Variationsproblems bildet die Galerkin-Methode

Wahle einen endlichdimensionalen Unterraum Xh ⊂ X undbetrachte das zugeordnete diskrete Variationsproblem:Finde ein uh ∈ Xh, so dass

B(uh, vh) = `(vh) ∀vh ∈ Xh (1)

Existenz und Eindeutigkeit nach dem Lax-Milgram Lemma

Wahl einer Basis von Xh ⇒ (1) kann aquivalent als LGSaufgefasst werden

Zu konstruieren: endlichdim. Unterraume von H1(Ω) und H10 (Ω)





CFEs und CFE-Raume (1/2)

Idee [Hackbusch und Sauter, 1995]:

Wahle ein Gitter T , welches das Rechengebiet lediglichuberdecken braucht

Zerlege den in Ω enthaltenen Teil jedes Gitterelements K ∈ Tin viele feinere geometrische Objekte

Weise den feineren geometrischen Objekten klassische FiniteElemente zu

Setze diese unter geeigneten Nebenbedingungen zu einemComposite Finite Element auf K zusammen

Fuge die CFEs wie bei der klassischen FEM global zu einemRaum zusammen (→ CFE-Raum)





CFEs und CFE-Raume (2/2)

Hier:

CFE-Raume basierend auf simplizialen Gittern und CFEs, dieaus simplizialen Lagrange Elementen zusammengesetzt werden

Vorgestellter Konstruktionsprozess:

CFE-Raume ursprunglich als Grundlage eines schnellen Loserskonzipiert

Basierend auf der Vergroberung von Finite Elemente Raumen

Konstruktionsprozess wird eine Hierarchie von immer feinerwerdenden Gittern liefern (→ CFE-Triangulierungen)und eine zugehorige Kette von CFE-Raumen aufsteigenderDimensionalitat





Erzeugung von CFE-Triangulierungen

Wir generieren eine Sequenz T0, . . . , Tlmax von immer feinerwerdenden Gittern mit Gebieten Ω0 ⊃ . . . ⊃ Ωlmax im Rd.

Tlmax wird die Grundlage fur das Zusammensetzen derklassischen FE liefern ⇒ es muss gelten: Ωlmax = Ω

H10 (Ω)⇒ Tlmax muss zudem ΓD ⊂ ∂Ω auflosen

Sei ein polygonal berandetes Rechengebiet Ω gegeben, z.B.

Erzeugung in drei Schritten





Schritt 1: Referenz-Triangulierungen

Verfeinerung eines initialen Gitters, welches Ω uberdeckt, sodass jedes Gitterelement vollstandig in Gitterelemente auf dernachst feineren Ebene zerfallt

Naturliche Hierarchie durch Ineinanderschachtellung





Schritt 2: Zwischen-Triangulierungen

Anpassung an den Rand ∂Ω durch Verschiebung randnaherEckknoten auf den Rand

Logische Hierarchie wird geerbt

Ω polygonal berandet, Verfeinerung in Schritt 1 fein genug⇒ feinste Zwischen-Triangulierung besitzt geeignetes Tlmax

als Teilmenge





Schritt 3: CFE-Triangulierungen

Entfernen aller Gitterelemente, die nicht zur Konstruktionder CFE-Raume benotigt werden (→ Ω0 ⊃ . . . ⊃ Ωlmax)

Logische Hierarchie wird geerbt





Konstruktion von CFE-Raumen (1/4)

Definieren nun eine Kette von zugehorigen CFE-RaumenSCFE

0 ⊂ . . . ⊂ SCFElmax

Gemeinsames Framework fur H1(Ω) und H10 (Ω)

Sei l ∈ 0, . . . , lmax fixiert und k ∈ N ein fixierterPolynomgrad

Notation (klassische Raume simplizialer Lagrange Elemente):

Sl :=v ∈ C0(Ωl) | v|K ∈ Pk(K) ∀K ∈ Tl

SCFEl kann als eine Anpassung des Raumes Sl an das

Rechengebiet Ω interpretiert werden

H10 (Ω) ⇒ zudem eine Anpassung an Null-Randwerte auf

ΓD ⊂ ∂Ω






l = lmax:

Hochstens Anpassung an Null-Randwerte notig, da Ωlmax = Ω

SCFElmax

:= Slmax oder SCFElmax

:= Slmax ∩H10 (Ω)

l < lmax:

Anpassung mit einem linearen Operator Ilmax,l : Sl → SCFElmax

(→ Prolongationsoperator)

Dieser setzt fur jedes vl ∈ Sl eine angepasste Funktionvlmax ∈ SCFE

lmaxzusammen

Linearitat ⇒ die entstehenden Funktionen bilden einenUnterraum von SCFE

lmax






Definition (CFE-Raume)

Fur l < lmax definieren wir

SCFEl := im(Ilmax,l) = v ∈ C0(Ω) | ∃vl ∈ Sl : v = Ilmax,l[vl]

wobei Ilmax,l : Sl → SCFElmax

in Abhangigkeit von der gewunschtenApproximationseigenschaft gewahlt wird.

Ilmax,l kann auf die Definition von linearen 1-Schritt-Operatoren Im+1,m : Sm → Sm+1, m ∈ l, . . . , lmax − 1zuruckgefuhrt werden:

Ilmax,l := Ilmax,lmax−1 Ilmax−1,lmax−2 · · · Il+1,l






1-Schritt-Operator Im+1,m : Sm → Sm+1 H1(Ω)

Werte ein vm ∈ Sm in den Lagrange Punkten von Sm+1 aus undinterpoliere stuckweise polynomiell in Sm+1

1-Schritt-Operator Im+1,m : Sm → Sm+1 H10 (Ω)

Analog. Passe die Werte jedoch in der Nahe des Dirichlet-RandesΓD ⊂ ∂Ω derart an, dass auf ΓD Null-Randwerte entstehen

SlIl+1,l−→ Sl+1

Il+2,l+1−→ Sl+2 −→ . . . −→ Slmax

(∩H1

0 (Ω))

Ωl ⊃ Ωl+1 ⊃ Ωl+2 ⊃ . . . ⊃ Ωlmax = Ω





Basisfunktionen

ϕlii=1,...,nl

Knotenbasis von Sl

=⇒ Ilmax,l[ϕli]i=1,...,nl

kanonische Basis von SCFEl

Freiheitsgrade liegen in den Lagrange Punkten von Sl(⇒ auch außerhalb von Ω)

H1(Ω), Polynomgrad k = 1:





Basisfunktionen - Zusammensetzung





Approximationvermogen

Die so konstruierten CFE-Raume SCFE0 , . . . , SCFE

lmaxubernehmen das

Approximationsvermogen der klassischen Raume S0, . . . , Slmax :

hl := maxK∈Tl diam(K) Gitterweite von TlDann existiert eine von hl unabhangiges C > 0, so dassfur jedes v ∈ Hk+1(Ω) bzw. jedes v ∈ H1

0 (Ω) ∩Hk+1(Ω)ein vl ∈ SCFE

l existiert mit

‖v − vl‖Hm(Ω) ≤ Chk+1−ml |v|Hk+1(Ω), m ∈ 0, 1

Beweis (Ω ⊂ R2, Polynomgrad k = 1):

H1(Ω): [Hackbusch und Sauter, 1995]H1

0 (Ω): [Sauter, 1997]





Assemblierung

Galerkin-Methode mit SCFEl samt kanonischer Basis

Ilmax,l[ϕli]i=1,...,nl

liefert LGS:

KlUl = Fl

(Kl)ij := B(Ilmax,l[ϕlj ], Ilmax,l[ϕ

li]), (Fl)i := `(Ilmax,l[ϕ

li])

l = lmax: assembliere elementweise, wie bei simplizialenLagrange Elementen ublich

l < lmax: direkte Assemblierung — Strategie? Quadratur?

Losung: Assembliere fur l′ = lmax und benutze schrittweisedas Produkt (→ Galerkin-Produkt)

Km = (Pm+1,m)TKm+1Pm+1,m

Fm = (Pm+1,m)TFm+1

m = lmax − 1, . . . , l

Pm+1,m darstellende Matrix von Im+1,m





Lokalisierung (1/3)

Situation:

Nur Diskretisierung fur ein festes l < lmax gewunscht

Nachteil bisher: LGS fur l = lmax muss assembliert werden

Beobachtung: Fernab von ∂Ω bildet Ilmax,l identisch ab(→ Ineinanderschatellung der Referenz-Triangulierungen)

Lokalisierung:

Generiere statt Tl+1, . . . , Tlmax

Randbereichs-Triangulierungen T locl+1, . . . , T loc

lmax

Wahle Teilmengen T locl+1, . . . , T loc

lmaxund Teilmenge Tl ⊂ Tl mit:

Tl ∪lmax⋃

m=l+1

T locm

ist eine (nicht-konforme) Triangulierung von Ω





Lokalisierung (2/3)





Lokalisierung (3/3)

Lokalisierung (Fortsetzung):

Die Funktionen in SCFEl konnen dann uber das

Zusammensetzen von simplizialen Lagrange Elementen auf dernicht-konformen Triangulierung von Ω dargestellt werden

Konstruktion von SCFEl

durch lokale 1-Schritt-Prolongationsoperatoren

Ermoglicht:

Schnelle Assemblierung durch lokale Galerkin-Produkte

Speichereffizienter Gittergenerator




ImplementierungNumerische ErgebnisseEnde

Ubersicht








Implementierung

Implementierung in C++, basierend auf dune-core, dune-femund dune-subgrid (ca. 6500 Zeilen Code):

Problemschnittstelle, kapselt das Modellproblem

Modellproblem fur den aufwandigeren Fall c = 0

Gittergenerator fur CFE-Triangulierungen(Lokalisierung moglich)

AluGrid (Referenz-Triangulierungen), GeometryGrid(Zwischen-Triangulierungen), SubGrid (CFE-Triangulierungen)

Familie von linearen Operatoren (Prolongationsoperatoren,Verknupfter Operator, Adjungierter Operator fur dasGalerkin-Produkt, . . .)

Algorithmus-Klasse mit drei Diskretisierungsmodi(klassische FEM, CFE, CFE mit Lokalisierung)





Numerische Ergebnisse (1/4)

−∆u = f in Ω,

u = gD auf ΓD = ∂Ω.

f und gD durch exakte Losung vorgegeben:

u(x) :=

d∏i=1

sin(2πxi)

Ω = Kugel in [0, 1]d mit Aussparung

d = 2

Erst k = 1, dann k = 2






k = 1:

l hl size CPU-time L2-error EOC H1-error EOC

0 0.424264 30 26.7549 0.30164 — 4.64485 —1 0.212132 98 28.5757 0.11055 1.44813 2.4414 0.9279272 0.106066 320 31.0443 0.0300921 1.87724 1.14962 1.086553 0.0534015 1204 33.3669 0.00617673 2.30752 0.480956 1.269864 0.026885 4629 39.658 0.0015517 2.01296 0.232293 1.060485 0.0136267 17950 84.2672 0.00038807 2.03951 0.114381 1.042576 0.00701445 70842 130.201 9.70511e-05 2.08709 0.0571176 1.045737 0.00371958 281705 273.402 2.43664e-05 2.17861 0.0285588 1.092678 0.00207055 1.12232e+06 585.962 6.21781e-06 2.3315 0.0142859 1.18247

Wir sehen: Theoretische Konvergenzordnung experimentell bestatigt






k = 2:


0 0.424264 30 27.5888 0.0702578 — 1.62389 —1 0.212132 98 79.7119 0.00765644 3.19791 0.355109 2.193122 0.106066 320 77.5402 0.000997681 2.94002 0.096342 1.882033 0.0534015 1204 311.802 0.000111473 3.19378 0.0217315 2.170054 0.026885 4629 955.536 1.27089e-05 3.16417 0.00524056 2.072555 0.0136267 17950 2235.61 1.5749e-06 3.07285 0.00129421 2.058076 0.00701445 70842 5801.75 1.97479e-07 3.12671 0.000323882 2.086077 0.00371958 281705 12253.5 2.47687e-08 3.27265 8.10292e-05 2.184198 0.00207055 1.12232e+06 19326.1 3.33147e-09 3.42468 2.13428e-05 2.2774

Wir sehen: Theoretische Konvergenzordnung experimentell bestatigt






−∆u = f in Ω,

∇u · n = gN auf ΓN = ∂Ω.

f und gD durch exakte Losung vorgegeben:

u(x) := 4

d∑i=1

(xi − 0.5)2 − 1

Ω = Kugel in [0, 1]d mit kleinen Lochern

d = 2

k = 1






k = 1:


0 0.424264 30 97.4322 0.125698 — 0.748371 —1 0.212132 98 105.007 0.0278702 2.17316 0.41001 0.8680932 0.106066 320 114.388 0.0042994 2.69651 0.209402 0.9693823 0.0534015 1206 129.733 0.000794572 2.46046 0.106377 0.9869454 0.026885 4656 191.945 0.000160159 2.33383 0.0534542 1.002765 0.0136599 18288 420.567 7.78649e-05 1.06512 0.0267956 1.019926 0.00704207 72238 1116.14 2.11697e-05 1.9657 0.0134071 1.045117 0.00400001 286832 4191.69 4.99169e-06 2.55442 0.00671126 1.223468 0.00236311 1.1319e+06 9816.43 1.15054e-06 2.78832 0.00336364 1.312449 0.00152264 4.49369e+06 11863.9 4.85162e-07 1.9646 0.00168469 1.57313

Wir sehen: Theoretische Konvergenzordnung (im Großen und Ganzen) experimentell bestatigt






−∆u = 0 in Ω,

u = 0 auf Γ1D,

u(x) = 5 + 4 · x2 auf Γ2D.

Ω = Kreis in [0, 1]2 mit 13”Kuhlstaben“ (→ Γ1

D)

d = 2

k = 1






d = 3: . . . moglich, hier ein erstes Ergebnis:





Ende

Vielen Dank fur die Aufmerksamkeit!

Gibt es irgendwelche Fragen?


Bildbasierte L osung von Partiellen Di …...Konstruktion von CFE-R aumen Assemblierung,...

Documents

Transcript of Bildbasierte L osung von Partiellen Di …...Konstruktion von CFE-R aumen Assemblierung,...