Composite Finite Elements for Trabecular Bone Microstructures

In many medical and technical applications, numerical simulations need to be performed for objects with interfaces of geometrically complex shape. We focus on the biomechanical problem of elasticity simulations for trabecular bone microstructures. The goal of this dissertation is to develop and implement an efficient simulation tool for finite element simulations on such structures, so-called composite finite elements. We will deal with both the case of material/void interfaces (complicated domains) and the case of interfaces between different materials (discontinuous coefficients). In classical finite element simulations, geometric complexity is encoded in tetrahedral and typically unstructured meshes. Composite finite elements, in contrast, encode geometric complexity in specialized basis functions on a uniform mesh of hexahedral structure. Other than alternative approaches (such as e.g. fictitious domain methods, generalized finite element methods, immersed interface methods, partition of unity methods, unfitted meshes, and extended finite element methods), the composite finite elements are tailored to geometry descriptions by 3D voxel image data and use the corresponding voxel grid as computational mesh, without introducing additional degrees of freedom, and thus making use of efficient data structures for uniformly structured meshes. The composite finite element method for complicated domains goes back to Wolfgang Hackbusch and Stefan Sauter and restricts standard affine finite element basis functions on the uniformly structured tetrahedral grid (obtained by subdivision of each cube in six tetrahedra) to an approximation of the interior. This can be implemented as a composition of standard finite element basis functions on a local auxiliary and purely virtual grid by which we approximate the interface. In case of discontinuous coefficients, the same local auxiliary composition approach is used. Composition weights are obtained by solving local interpolation problems for which coupling conditions across the interface need to be determined. These depend both on the local interface geometry and on the (scalar or tensor-valued) material coefficients on both sides of the interface. We consider heat diffusion as a scalar model problem and linear elasticity as a vector-valued model problem to develop and implement the composite finite elements. Uniform cubic meshes contain a natural hierarchy of coarsened grids, which allows us to implement a multigrid solver for the case of complicated domains. Besides simulations of single loading cases, we also apply the composite finite element method to the problem of determining effective material properties, e.g. for multiscale simulations. For periodic microstructures, this is achieved by solving corrector problems on the fundamental cells using affine-periodic boundary conditions corresponding to uniaxial compression and shearing. For statistically periodic trabecular structures, representative fundamental cells can be identified but do not permit the periodic approach. Instead, macroscopic displacements are imposed using the same set as before of affine-periodic Dirichlet boundary conditions on all faces. The stress response of the material is subsequently computed only on an interior subdomain to prevent artificial stiffening near the boundary. We finally check for orthotropy of the macroscopic elasticity tensor and identify its axes.Zusammengesetzte finite Elemente für trabekuläre Mikrostrukturen in Knochen In vielen medizinischen und technischen Anwendungen werden numerische Simulationen für Objekte mit geometrisch komplizierter Form durchgeführt. Gegenstand dieser Dissertation ist die Simulation der Elastizität trabekulärer Mikrostrukturen von Knochen, einem biomechanischen Problem. Ziel ist es, ein effizientes Simulationswerkzeug für solche Strukturen zu entwickeln, die sogenannten zusammengesetzten finiten Elemente. Wir betrachten dabei sowohl den Fall von Interfaces zwischen Material und Hohlraum (komplizierte Gebiete) als auch zwischen verschiedenen Materialien (unstetige Koeffizienten). In klassischen Finite-Element-Simulationen wird geometrische Komplexität typischerweise in unstrukturierten Tetraeder-Gittern kodiert. Zusammengesetzte finite Elemente dagegen kodieren geometrische Komplexität in speziellen Basisfunktionen auf einem gleichförmigen Würfelgitter. Anders als alternative Ansätze (wie zum Beispiel fictitious domain methods, generalized finite element methods, immersed interface methods, partition of unity methods, unfitted meshes und extended finite element methods) sind die zusammengesetzten finiten Elemente zugeschnitten auf die Geometriebeschreibung durch dreidimensionale Bilddaten und benutzen das zugehörige Voxelgitter als Rechengitter, ohne zusätzliche Freiheitsgrade einzuführen. Somit können sie effiziente Datenstrukturen für gleichförmig strukturierte Gitter ausnutzen. Die Methode der zusammengesetzten finiten Elemente geht zurück auf Wolfgang Hackbusch und Stefan Sauter. Man schränkt dabei übliche affine Finite-Element-Basisfunktionen auf gleichförmig strukturierten Tetraedergittern (die man durch Unterteilung jedes Würfels in sechs Tetraeder erhält) auf das approximierte Innere ein. Dies kann implementiert werden durch das Zusammensetzen von Standard-Basisfunktionen auf einem lokalen und rein virtuellen Hilfsgitter, durch das das Interface approximiert wird. Im Falle unstetiger Koeffizienten wird die gleiche lokale Hilfskonstruktion verwendet. Gewichte für das Zusammensetzen erhält man hier, indem lokale Interpolationsprobleme gelöst werden, wozu zunächst Kopplungsbedingungen über das Interface hinweg bestimmt werden. Diese hängen ab sowohl von der lokalen Geometrie des Interface als auch von den (skalaren oder tensorwertigen) Material-Koeffizienten auf beiden Seiten des Interface. Wir betrachten Wärmeleitung als skalares und lineare Elastizität als vektorwertiges Modellproblem, um die zusammengesetzten finiten Elemente zu entwickeln und zu implementieren. Gleichförmige Würfelgitter enthalten eine natürliche Hierarchie vergröberter Gitter, was es uns erlaubt, im Falle komplizierter Gebiete einen Mehrgitterlöser zu implementieren. Neben Simulationen einzelner Lastfälle wenden wir die zusammengesetzten finiten Elemente auch auf das Problem an, effektive Materialeigenschaften zu bestimmen, etwa für mehrskalige Simulationen. Für periodische Mikrostrukturen wird dies erreicht, indem man Korrekturprobleme auf der Fundamentalzelle löst. Dafür nutzt man affin-periodische Randwerte, die zu uniaxialem Druck oder zu Scherung korrespondieren. In statistisch periodischen trabekulären Mikrostrukturen lassen sich ebenfalls Fundamentalzellen identifizieren, sie erlauben jedoch keinen periodischen Ansatz. Stattdessen werden makroskopische Verschiebungen zu denselben affin-periodischen Randbedingungen vorgegeben, allerdings durch Dirichlet-Randwerte auf allen Seitenflächen. Die Spannungsantwort des Materials wird anschließend nur auf einem inneren Teilbereich berechnet, um künstliche Versteifung am Rand zu verhindern. Schließlich prüfen wir den makroskopischen Elastizitätstensor auf Orthotropie und identifizieren deren Achsen

Towards a robust Terra

In this work mantle convection simulation with Terra is investigated from a numerical point of view, theoretical analysis as well as practical tests are performed. The stability criteria for the numerical formulation of the physical model will be made clear. For the incompressible case and the Terra specific treatment of the anelastic approximation, two inf-sup stable grid modifications are presented, which are both compatible with hanging nodes. For the Q1hQ12h element pair a simple numeric test is introduced to prove the stability for any given grid. For the Q1h Pdisc 12h element pair and 1-regular refinements with hangig nodes an existing general proof can be adopted. The influence of the slip boundary condition is found to be destabilizing. For the incompressible case a cure can be adopted from the literature. The necessary conditions for the expansion of the stability results to the anelastic approximation will be pointed out. A numerical framework is developed in order to measure the effect of different numerical approaches to improve the handling of strongly varying viscosity. The framework is applied to investigate how block smoothers with different block sizes, combination of different block smoothers, different prolongation schemes and semi coarsening influence the multigrid performance. A regression-test framework for Terra will be briefly introduced

Éléments finis hp adaptatifs avec contraction d’erreur garantie et solveurs multi-niveaux inexacts

We propose new practical adaptive refinement algorrithms for conforming hp-finite element approximations of elliptic problems. We consider the use of both exact and inexact solevsr within the established framework of adaptive methods consisting of four concatenated modules : SOLVE, ESTIMATE, MARK, REFINE. The strategies are driven by guaranteed equilibrated flux a posteriori error estimators. Namely, for an inexact approximation obtained by an (arbitrary) iterative algebraic solver, the bounds for the total, the algebraic, and the discretization errors are provided. Our hp-refinement criterion hinges on from solving two local residual problems posed on patches of elements around marked vertices selected by a bulk-chasing criterion. They respectively emulate h-refinement and p-refinement. One particular feature of energy error in the next adaptative loop step with respect to the present one. Numerical experiments are presented to turns out to be excellent, with effectivity indices close to the optimal value of one. In practice, we observe asymptomatic exponential convergence rates, in both the exact and inexact algebraic solver settings. Finally, we also provide a theoretical analysis of the proposed strategies.Nous proposons de nouveaux algorithmes de raffinement adaptatif pour l'approximation des problèmes elliptiques par la méthode des éléments finis hp. Nous considérons des solveurs algébriques exacts puis inexacts au sein du cadre générique des méthodes adaptatives consistant en quatre modules concaténés: RESOLUTION, ESTIMATION, MARQUAGE, RAFFINEMENT. Les stratégies reposent sur la construction d'estimateurs d'erreur a posteriori par flux équilibrés. Notamment, pour une approximation inexacte obtenue par un solveur algébrique itératif (arbitraire), nous prouvons une borne sur l'erreur totale ainsi que sur l'erreur algébrique et l'erreur de discrétisation. La structure hiérarchique des espaces d'éléments finis hp est cruciale pour obtenir la borne supérieure sur l'erreur algébrique, ce qui nous permet de formuler des critères d'arrêt précis pour le solveur algébrique. Notre critère de raffinement hp repose sur la résolution de deux problèmes résiduels locaux, posés sur les macro-éléments autour des sommets du maillage qui ont été marqués. Ces derniers sont sélectionnés par un critère de type bulk-chasing. Ceux deux problèmes résiduels imitent l'effet du raffinement h et p. Une caractéristique de notre approche est que nous obtenons une quantité calculable qui donne une borne garantie sur le rapport entre l'erreur d'énergie (inconnue) à la prochaine étape de la boucle adaptative et l'erreur actuelle (i.e. sur le facteur de réduction d'erreur). Des simulations numériques sont présentées afin de valider les stratégies adaptatives. Nous examinons la précision de notre borne sur le facteur de réduction d'erreur qui s'avère être excellente, avec des indices d'efficacité proches de 1. En pratique, nous observons des taux de convergence asymptotiquement exponentiels, aussi bien dans le cadre de la résolution algébrique exacte que dans celui de la résolution inexacte.Enfin, nous menons une analyse théorique des stratégies proposées

Multigrid Methods for Elliptic Optimal Control Problems

In this dissertation we study multigrid methods for linear-quadratic elliptic distributed optimal control problems. For optimal control problems constrained by general second order elliptic partial differential equations, we design and analyze a $P_1$ finite element method based on a saddle point formulation. We construct a $W$-cycle algorithm for the discrete problem and show that it is uniformly convergent in the energy norm for convex domains. Moreover, the contraction number decays at the optimal rate of $m^{-1}$, where $m$ is the number of smoothing steps. We also prove that the convergence is robust with respect to a regularization parameter. The robust convergence of $V$-cycle and $W$-cycle algorithms on general domains are demonstrated by numerical results. For optimal control problems constrained by symmetric second order elliptic partial differential equations together with pointwise constraints on the state variable, we design and analyze symmetric positive definite $P_1$ finite element methods based on a reformulation of the optimal control problem as a fourth order variational inequality. We develop a multigrid algorithm for the reduced systems that appear in a primal-dual active set method for the discrete variational inequalities. The performance of the algorithm is demonstrated by numerical results

The simulation of incompressible flow using the artificial compressibility method.

In this work an algorithm is developed for simulating incompressible steady flow on two and three-dimensional unstructured meshes. The Navier-Stokes equations are briefly reviewed as the basic governing equation for fluid flow. Using this set of equations, in the limit of incompressible flow, the problem of imposing the time independent continuity equation on the momentum equations arises. This difficulty can be removed by employing the Artificial Compressibility approach. This approach modifies the continuity equation by adding a pseudo pressure time derivative. This modification makes the set of equations well conditioned for numerical solution. If the set of modified equations is used for the solution of steady state problems, the added pressure derivative tends to zero and the set of equations reduces to the steady state incompressible Navier-Stokes equations. The cell-vertex finite volume method is employed for solving the modified equations on unstructured triangular and tetrahedral meshes. The principles of central difference space discretisation are described and the basic ideas behind adding artificial dissipation term are reviewed. A normalisation procedure for the computation of the artificial dissipation term is adopted. Two different formulations based upon a cell-vertex finite volume method and a Galerkin finite element method for the discretisation of the viscous terms on unstructured triangular meshes are employed in two and three dimensions. A modification to the finite volume formulation is introduced for improving accuracy on unstructured grids. The issues relating to multi-stage time stepping, boundary conditions and some techniques for increasing computational efficiency are described. A general review of several methods for generating regular and irregular unstructured triangular and tetrahedral meshes are presented. The proposed algorithm is validated by solving several inviscid and viscous two-dimensional test cases. Extension of the algorithm to three dimensions are studied by using some further benchmark examples. Some engineering applications are considered to present the ability of the developed flow solver to simulate more complicated real world problems. Finally, some conclusions are drawn and a few guidelines for further research work are suggested

FINITE ELEMENT METHODS FOR AXISYMMETRIC PDES AND DIVERGENCE FREE FINITE ELEMENT PAIRS ON PARTICULAR MESH REFINEMENTS

This dissertation discusses the following two main topics. 1) Finite element approximation for Partial Differential Equations (PDEs) defined on axisymmetric domains: We introduce the Darcy equations on axisymmetric domains and we show the stability of a low--order Raviart-Thomas element pair. We provide numerical experiments to support our theoretical results. Also, we introduce the Stokes equations on axisymmetric domains and show that the axisymmetric Stokes equations can fit within a commutative de Rham complex. 2) Connection between the grad-div stabilized and divergence-free Stokes finite element pairs and low--order divergence-free elements on particular mesh refinements: We introduce the most recent results that connect the grad-div stabilized Taylor--Hood (TH) finite element pair and divergence-free Scott--Vogelius (SV) finite element pairs, and we use these results to extend and generalize this connection to other Stokes finite element pairs. Finally, we provide numerical examples for low order divergence-free Stokes finite element pairs defined on particular mesh refinements. This research is focused on the numerical implementation aspects of these finite element pairs