A New Closure Strategy for Proper Orthogonal Decomposition Reduced-Order Models Imran Akhtar, Zhu Wang, Jeff Borggaard, Traian Iliescu Proper orthogonal decomposition (POD) is one of the most significant reduced-order modeling (ROM) techniques in fluid mechanics. However, the application of POD based reduced-order models (POD-ROMs) is primarily limited to laminar flows due to the decay of physical accuracy. A few nonlinear closure models have been developed for improving the accuracy and stability of the POD-ROMs, which are generally computationally expensive. In this paper we propose a new closure strategy for POD-ROMs that is both accurate and effective. In the new closure model, the Frobenius norm of the Jacobian of the POD-ROM is introduced as the eddy viscosity coefficient. As a first step, the new method has been tested on a one-dimensional Burgers equation with a small dissipation coefficient nu = 10^(-3). Numerical results show that the Jacobian based closure model greatly improves the physical accuracy of the POD-ROM, while maintaining a low computational cost. Remarks around 50 lines of Matlab: short finite element implementation Jochen Alberty, Carsten Carstensen and Stefan A. Funken A short Matlab implementation for P1-Q1 finite elements on triangles and parallelograms is provided for the numerical solution of elliptic problems with mixed boundary conditions on unstructured grids. According to the shortness of the program and the given documentation, any adaptation from simple model examples to more complex problems can easily be performed. Numerical examples prove the flexibility of the Matlab tool. On the Efficiency of Symbolic Computations Combined with Code Generation for Finite Element Methods MARTIN SANDVE ALNAES, KENT-ANDRE MARDAL Efficient and easy implementation of variational forms for finite element discretization can be accomplished with metaprogramming. Using a high-level language like Python and symbolic mathematics makes an abstract problem definition possible, but the use of a low-level compiled language is vital for run-time efficiency. By generating low-level C++ code based on symbolic expressions for the discrete weak form, it is possible to accomplish a high degree of abstraction in the problem definition while surpassing the run-time efficiency of traditional hand written C++ codes. We provide several examples where we demonstrate orders of magnitude in speedup. A stable finite element method for the Stokes equations, D N Arnold, F Brezzi, M Fortin, We present in this paper a new velocity-pressure finite element for the computation of Stokes flow. We discretized the velocity field with continuous piecewise linear functions enriched by bubble functions, and the pressure by piecewise linear functions. We show that this element statisfies the usual inf-sup condition and converges with first order for both velocities and pressure. Finally we relate this element to families of higher-order elements and to the popular Taylor-Hood element. deal.II—A General-Purpose Object-Oriented Finite Element Library W. BANGERTH, R. HARTMANN, G. KANSCHAT An overview of the software design and data abstraction decisions chosen for deal.II, a general purpose finite element library written in C++, is given. The library uses advanced object-oriented and data encapsulation techniques to break finite element implementations into smaller blocks that can be arranged to fit users requirements. Through this approach, deal.II supports a large number of different applications covering a wide range of scientific areas, programming methodologies, and application-specific algorithms, without imposing a rigid framework into which they have to fit. A judicious use of programming techniques allows us to avoid the computational costs frequently associated with abstract object-oriented class libraries. The paper presents a detailed description of the abstractions chosen for defining geometric information of meshes and the handling of degrees of freedom associated with finite element spaces, as well as of linear algebra, input/output capabilities and of interfaces to other software, such as visualization tools. Finally, some results obtained with applications built atop deal.II are shown to demonstrate the powerful capabilities of this toolbox. Adaptive Mesh Generation Marshall Bern These notes cover topics in mesh generation from a computational geometry perspective. This perspective means emphasis on difficiult domain geometry, unstructured triangular and tetrahedral meshes, and provable bounds on quality and complexity. We concentrate on new results, especially those results applicable to solution-adaptive methods in computational fluid dynamics. On the Finite Element Solution of the Pure Neumann Problem Pavel Bochev, R. B. Lehoucq This paper considers the finite element approximation and algebraic solution of the pure Neumann problem. Our goal is to present a concise variational framework for the finite element solution of the Neumann problem that focuses on the interplay between the alge braic and variational problems. While many of the results that stem from our analysis are known by some experts, they are seldom derived in a rigorous fashion and remain part of numerical folklore. As a result, this knowledge is not accessible (or appreciated) by many practitioners-both novices and experts-in one source. Our paper contributes a simple, yet insightful link between the continuous and algebraic variational forms that will prove useful. Solution to PDEs using radial basis function finite-differences (RBF-FD) on multiple GPUs Evan F. Bollig, Natasha Flyer, Gordon Erlebacher This paper presents parallelization strategies for the radial basis function-finite difference (RBF-FD) method. As a generalized finite differencing scheme, the RBF-FD method functions without the need for underlying meshes to structure nodes. It offers high-order accuracy approximation and scales as OðNÞ per time step, with N being with the total number of nodes. To our knowledge, this is the first implementation of the RBF-FD method to leverage GPU accelerators for the solution of PDEs. Additionally, this implementation is the first to span both multiple CPUs and multiple GPUs. OpenCL kernels target the GPUs and inter-processor communication and synchronization is managed by the Message Passing Interface (MPI). We verify our implementation of the RBF-FD method with two hyperbolic PDEs on the sphere, and demonstrate up to 9x speedup on a commodity GPU with unoptimized kernel implementations. On a high performance cluster, the method achieves up to 7x speedup for the maximum problem size of 27,556 nodes. POD and CVT-based reduced-order modeling of Navier–Stokes flows John Burkardt, Max Gunzburger, Hyung-Chun Lee A discussion of reduced-order modeling for complex systems such as fluid flows is given to provide a context for the construction and application of reduced-order bases. Reviews of the POD (proper orthogonal decomposition) and CVT (centroidal Voronoi tessellation) approaches to reduced-order modeling are provided, including descriptions of POD and CVT reduced-order bases, their construction from snapshot sets, and their application to the low-cost simulation of the Navier–Stokes system. Some concrete incompressible flow examples are used to illustrate the construction and use of POD and CVT reduced-order bases and to compare and contrast the two approaches to reduced-order modeling. Assembly of Finite Element Methods on Graphics Processors Cris Cecka, Adrian J. Lew, and E. Darve Recently, graphics processing units (GPUs) have had great success in accelerating many numerical computations. We present their application to computations on unstructured meshes such as those in finite element methods. Multiple approaches in assembling and solving sparse linear systems with NVIDIA GPUs and the Compute Unified Device Architecture (CUDA) are created and analyzed. Multiple strategies for efficient use of global, shared, and local memory, methods to achieve memory coalescing, and optimal choice of parameters are introduced. We find that with appropriate preprocessing and arrangement of support data, the GPU coprocessor using single-precision arithmetic achieves speedups of 30 or more in comparison to a well optimized double-precision single core implementation. We also find that the optimal assembly strategy depends on the order of polynomials used in the finite element discretization. Finite Part Integrals and Hypersingular Kernels Youn-Sha Chan, Albert C. Fannjiang, Glaucio H. Paulino, and Bao-Feng Feng Singular integral equation method is one of the most effective numerical methods solving a plane crack problem in fracture mechanics. Depending on the choice of the density function, very often a higher order of sigularity appears in the equation, and we need to give a proper meaning of the integration. In this article we address the Hadamard finite part integral and how it is used to solve the plane crack problems. Properties of the Hadamard finite part integral will be summarized and compared with other type of integrals. Some numerical results for crack problems by using Hadamard finite part integral will be provided. An efficient way to perform the assembly of finite element matrices in Matlab and Octave François Cuvelier, Caroline Japhet, Gilles Scarella We describe different optimization techniques to perform the assembly of finite element matrices in Matlab and Octave, from the standard approach to recent vectorized ones, without any low level language used. We finally obtain a simple and efficient vectorized algorithm able to compete in performance with dedicated software such as FreeFEM++. The principle of this assembly algorithm is general, we present it for different matrices in the P1 finite elements case and in linear elasticity. We present numerical results which illustrate the computational costs of the different approaches. Algorithm 866: IFISS, A Matlab Toolbox for Modelling Incompressible Flow HOWARD C. ELMAN, ALISON RAMAGE, DAVID J. SILVESTER IFISS is a graphical Matlab package for the interactive numerical study of incompressible flow problems. It includes algorithms for discretization by mixed finite element methods and a posteriori error estimation of the computed solutions. The package can also be used as a computational laboratory for experimenting with state-of-the-art preconditioned iterative solvers for the discrete linear equation systems that arise in incompressible flow modelling. A unique feature of the package is its comprehensive nature; for each problem addressed, it enables the study of both discretization and iterative solution algorithms as well as the interaction between the two and the resulting effect on overall efficiency. EXACT FULLY 3D NAVIER-STOKES SOLUTIONS FOR BENCHMARKING C. ROSS ETHIER, D. A. STEINMAN Unsteady analytical solutions to the incompressible Navier-Stokes equations are presented. They are fully three-dimensional vector solutions involving all three Cartesian velocity components, each of which depends non-trivially on all three co-ordinate directions. Although unlikely to be physically realized, they are well suited for benchmarking, testing and validation of three-dimensional incompressible Navier-Stokes solvers. The use of such a solution for benchmarking purposes is described. Error Estimates for the Approximation of a Class of Variational Inequalities Richard S. Falk In this paper, we prove a general approximation theorem useful in obtaining order of convergence estimates for the approximation of the solutions of a class of variational inequalities. The theorem is then applied to obtain an "optimal" rate of convergence for the approximation of a second-order elliptic problem with convex set K. J P FINK, W RHEINBOLDT On the Error Behavior of the Reduced Basis Technique for Nonlinear Finite Element Approximations This paper provides a theoretical foundation for the reduced basis technique and gives mathematical reasons for its effectiveness. Some rather general results about the existence of solution curves for a large class of nonlinear problems are developed. These results are then used to derive estimates between an exact solution curve and its reduced basis approximation. On Weakly Imposted Boundary Conditions for Second Order Problems Jouni Freund, Rolf Stenberg We consider the finite element method in which essential boundary conditions are imposed in a weak sense using a technique introduced by Nitsche. This technique is applied in a convection-diffusion-reaction problem and for slip boundary conditions in Stokes flow. We give the error estimates for the methods and also some numerical results. Radial Basis Function-generated Finite Differences: A Mesh-free Method for Computational Geosciences Natasha Flyer, Grady B. Wright, Bengt Fornberg Radial basis function generated finite differences (RBF-FD) is a mesh-free method for numerically solving partial differential equations (PDEs) that emerged in the last decade and has shown rapid growth in the last few years. From a practical standpoint, RBF-FD sprouted out of global RBF methods, which have shown exceptional numerical qualities in terms of accuracy and time stability for numerically solving PDEs, but are not practical when scaled to very large problem sizes because of their computational cost and memory requirements. RBF-FD bypass these issues by using local approximations for derivatives instead of global ones. Matrices in the RBF-FD methodology go from being completely full to 99% empty. Of course, the sacrifice is the exchange of spectral accuracy from the global RBF methods for high-order algebraic convergence of RBF-FD, assuming smooth data. However, since natural processes are almost never infinitely differentiable, little is lost and much gained in terms of memory and runtime. This article provides a survey of a group of topics relevant to using RBF-FD for a variety of problems that arise in the geosciences. Particular emphasis is given to problems in spherical geometries, both on surfaces and within a volume. Applications discussed include non-linear shallow water equations on a sphere, reaction-diffusion diffusion equations, global electric circuit, and mantle convection in a spherical shell. The results from the last three of these applications are new and have not been presented before for RBF-FD. Scattered Data Interpolation: Tests of Some Methods Richard Franke This paper is concerned with the evaluation of methods for scattered data interpolation and some of the results of the tests when applied to a number of methods. The process involves evaluation of the methods in terms of timing, storage, accuracy, visual pleasantness of the surface, and ease of implementation. To indicate the flavor of the type of results obtained, we give a summary table and representative perspective plots of several surfaces. Bounds on the spectral and maximum norms of the finite element stiffness, flexibility and mass matrices, Isaac Fried Upper and lower bounds are established on the spectral and maximum norms (and hence on the corresponding condition numbers) of the stiffness, flexibility (the inverse of the stiffness) and mass matrices generated from regular and irregular meshes of finite elements. Explicit expressions for these bounds are derived in terms of the intrinsic and discretization parameters, for second and fourth order problems in one, two and three dimensions discretized with linear, triangular and tetrahedronal elements. Simple Finite Element Methods for Approximating Predator–Prey Dynamics in Two Dimensions Using Matlab Marcus R. Garvie, John Burkardt, Jeff Morgan We describe simple finite element schemes for approximating spatially extended predator–prey dynamics with the Holling type II functional response and logistic growth of the prey. The finite element schemes generalize Scheme 1 in the paper by Garvie (Bull Math Biol 69(3):931–956, 2007). We present user-friendly, open-source Matlab code for implementing the finite element methods on arbitrary shaped two-dimensional domains with Dirichlet, Neumann, Robin, mixed Robin– Neumann, mixed Dirichlet–Neumann, and Periodic boundary conditions. Users can download, edit, and run the codes from http://www.uoguelph.ca/~mgarvie/. In addition to discussing the well posedness of the model equations, the results of numerical experiments are presented and demonstrate the crucial role that habitat shape, initial data, and the boundary conditions play in determining the spatiotemporal dynamics of predator–prey interactions. As most previous works on this problem have focussed on square domains with standard boundary conditions, our paper makes a significant contribution to the area. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities Christophe Geuzaine and Jean-Francois Remacle2 Gmsh is an open-source 3-D finite element grid generator with a build-in CAD engine and post processor. Its design goal is to provide a fast, light and user-friendly meshing tool with parametric input and advanced visualization capabilities. This paper presents the overall philosophy, the main design choices and some of the original algorithms implemented in Gmsh. Transfinite element methods: blending function interpolation over arbitary curved element domains, William Gordon, Charles Hall In order to better conform to curved boundaries and material interfaces, curved finite elements have been widely applied in recent years by practicing engineering analysts. The most well known of such elements are the isoparametric elements. As Zienkiewicz points out, there has been a certain parallel between the development of element types as used in finite element analyses and the independent development of methods for the mathematical description of general free-form surfaces. One of the purposes of this paper is to show that the relationship between these two areas of recent mathematical activity is indeed quite intimate. In order to establish this relationship, we introduce the notion of a transfinite element, which, in brief, is an invertible mapping from a square parameter domain onto a closed, bounded and simply connected region in the XY plane together with a transfinite blending-function type interpolant to the dependent variable defined over R. The subparametric, isoparametric and superparametric element types discussed by Zienkiewicz can all be shown to be special cases obtainable by various discretizations of transfinite elements. Actual error bounds are derived for a wide class of semi-discretized transfinite elements (with the nature of the mapping remaining unspecified) as applied to two types of boundary value problems. These bounds for semi-discretized elements are then specialized to obtain bounds for the familiar isoparametric elements. While we consider only two dimensional elements, extensions to higher dimensions is straightforward. On conforming finite element methods for the inhomogeneous Navier-Stokes Equations, Max Gunzburger, Janet Peterson We consider the stationary Navier-Stokes equations, written in terms of the primitive variables, in the case where both the partial differential equations and boundary conditions are inhomogeneous. Under certain conditions on the data, the existence and uniqueness of the solution of a weak formulation of the equations can be guaranteed. A conforming finite element method is presented and optimal estimates for the error of the approximate solution are proved. In addition, the convergence properties of iterative methods for the solution of the discrete nonlinear algebraic systems resulting from the finite element algorithm are given. Numerical examples, using an efficient choice of finite element spaces, are also provided. New development in freefem++ F. HECHT This is a short presentation of the freefem++ software. In Section 1, we recall most of the characteristics of the software, In Section 2, we recall how to to build the weak form of a partial differential equation (PDE) from the strong form. In the 3 last sections, we present different examples and tools to illustrate the power of the software. First we deal with mesh adaptation for problems in two and three dimension, second, we solve numerically a problem with phase change and natural convection, and the finally to show the possibilities for HPC we solve a Laplace equation by a Schwarz domain decomposition problem on parallel computer An Overview of the Trilinos Project, MICHAEL A. HEROUX, ANDREW G. SALINGER, ROSCOE A. BARTLETT, HEIDI K. THORNQUIST, VICKI E. HOWLE, RAY S. TUMINARO, ROBERT J. HOEKSTRA, JAMES M. WILLENBRING, JONATHAN J. HU, and ALAN WILLIAMS, TAMARA G. KOLDA, RICHARD B. LEHOUCQ, KEVIN R. LONG, ROGER P. PAWLOWSKI, ERIC T. PHIPPS, KENDALL S. STANLEY The Trilinos Project is an effort to facilitate the design, development, integration, and ongoing support of mathematical software libraries within an object-oriented framework for the solution of large-scale, complex multiphysics engineering and scientific problems. Trilinos addresses two fundamental issues of developing software for these problems: (i) providing a streamlined process and set of tools for development of new algorithmic implementations and (ii) promoting interoperability of independently developed software. Trilinos uses a two-level software structure designed around collections of packages. A Trilinos package is an integral unit usually developed by a small team of experts in a particular algorithms area such as algebraic preconditioners, nonlinear solvers, etc. Packages exist underneath the Trilinos top level, which provides a common look-and-feel, including configuration, documentation, licensing, and bug-tracking. Here we present the overall Trilinos design, describing our use of abstract interfaces and default concrete implementations. We discuss the services that Trilinos provides to a prospective package and how these services are used by various packages. We also illustrate how packages can be combined to rapidly develop new algorithms. Finally, we discuss how Trilinos facilitates highquality software engineering practices that are increasingly required from simulation software. Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder Volker John This paper presents a numerical study of a two-dimensional time-dependent flow around a cylinder. Its main objective is to provide accurate reference values for the maximal drag and lift coefficient at the cylinder and for the pressure difference between the front and the back of the cylinder at the final time. In addition, the accuracy of these values obtained with different time stepping schemes and different finite element methods is studied. Streamline Diffusion Methods for the Incompressible Euler and Navier-Stokes Equations Claes Johnson and Jukka Saranen We present and analyze extensions of the streamline diffusion finite element method to the time-dependent two-dimensional Navier-Stokes equations for an incompressible fluid in the case of high Reynolds numbers. The limit case with zero viscosity, the Euler equations, is also considered. Motivation for using radial basis functions to solve PDEs, Edward Kansa The numerical solution of partial differential equations (PDEs) has been dominated by eitehr finite difference methods (FDM), finite element methods (FEM), and finite volume methods (FVM). These methods can be derived from the assumptions of the local interpolation schemes. These methods require a mesh to support the localized approximations; the construction of a mesh in three or more dimensions is a non-trivial problem. Typically, with these methods only the function is continuous across meshes, but not its partial derivatives. In practice, only low order approximations are used because of the notorious polynomial snaking problem. While higher order schemes are necessary for more accurate approximations of the spatial derivatives, they are not sufficient without monotonicity constraints. Because of the low order schemes typically employed, the spatial truncation errors can only be controlled by using progressively smaller meshes. The mesh spacing, h, must be sufficiently fine to capture the function's partial derivative behavior and to avoid unnecessarily large amounts of numerical artifacts contaminating the solution. Spectral methods, while offering very high order spatial schemes typically depend upon tensor grids in higher dimensions. ON A GALERKIN-LAGRANGE MULTIPLIER METHOD FOR THE STATIONARY NAVIER-STOKES EQUATIONS OHANNES A. KARAKASHIAN A Galerkin-Lagrange multiplier formulation is used for the numerical solution of the stationary Navier-Stokese equations, in order to avoid the construction of zero-divergence elements. The formulation is based on different approximating spaces for the velocity field and the pressure. Optimal rate of convergence estimates are derived. Moreover, a Galerkin-Newton scheme for the solution of the nonlinear equations is shown to be quadratically locally convergent. Another scheme is shown to be linearly globally convergent. Finite Element Integration on GPUs MATTHEW G. KNEPLEY, ANDY R. TERREL We present a novel finite element integration method for low-order elements on GPUs. We achieve more than 100GF for element integration on first order discretizations of both the Laplacian and Elasticity operators on an NVIDIA GTX285, which has a nominal single precision peak flop rate of 1 TF/s and bandwidth of 159 GB/s, corresponding to a bandwidth limited peak of 40 GF/s. Subject-specific finite-element modeling of normal aortic valve biomechanics from 3D+t TEE images Michel R. Labrosse, Carsten J. Beller, Munir Boodhwani, Christopher Hudson, Benjamin Sohmer In the past decades, developments in transesophageal echocardiography (TEE) have opened new horizons in reconstructive surgery of the aortic valve (AV), whereby corrections are made to normalize the geometry and function of the valve, and effectively treat leaks. To the best of our knowledge, we propose the first integrated framework to process subject-specific 3D+t TEE AV data, determine age-matched material properties for the aortic and leaflet tissues, build a finite element model of the unpressurized AV, and simulate the AV function throughout a cardiac cycle. For geometric reconstruction purposes, dedicated software was created to acquire the 3-D coordinates of 21 anatomical landmarks of the AV apparatus in a systematic fashion. Measurements from ten 3D+t TEE datasets of normal AVs were assessed for interand intra-observer variability. These tests demonstrated mean measurement errors well within the acceptable range. Simulation of a complete cardiac cycle was successful for all ten valves and validated the novel schemes introduced to evaluate age-matched material properties and iteratively scale the unpressurized dimensions of the valves such that, given the determined material properties, the dimensions measured in vivo closely matched those simulated in late diastole. The leaflet coaptation area, describing the quality of the sealing of the valve, was measured directly from the medical images and was also obtained from the simulations; both approaches correlated well. The mechanical stress values obtained from the simulations may be interpreted in a comparative sense whereby higher values are indicative of higher risk of tearing and/or development of calcification. A Two-Level Discretization Method for the Navier-Stokes Equations W. LAYTON We propose and analyze a two level method of discretizing the nonlinear Navier-Stokes equations. This method has optimal accuracy and requires neither iteration nor the solution of more than a very small number of nonlinear equations. A NONLINEAR, SUBGRIDSCALE MODEL FOR INCOMPRESSIBLE VISCOUS FLOW PROBLEMS WILLIAM J. LAYTON We consider a nonlinear subgridscale model of the Navier-Stokes equations resulting in a Ladyzhenskaya-type system. The difference is that the power "p" and scaling coefficient mu(h) do not arise from macroscopic fluid properties and can be picked to ensure both L-stability and yet be of the order of the basic discretization error in smooth regions. With a properly scaled p-Laplacian-type artificial viscosity one can construct a higher-order method which is just as stable as first-order upwind methods. A Supernodal Approach to Incomplete LU Factorization with Partial Pivoting XIAOYE S. LI, MEIYUE SHAO We present a new supernode-based incomplete LU factorization method to construct a preconditioner for solving sparse linear systems with iterative methods. The new algorithm is primarily based on the ILUTP approach by Saad, and we incorporate a number of techniques to improve the robustness and performance of the traditional ILUTP method. These include new dropping strategies that accommodate the use of supernodal structures in the factored matrix and an area-based fill control heuristic for the secondary dropping strategy. We present numerical experiments to demonstrate that our new method is competitive with the other ILU approaches and is well suited for modern architectures with memory hierarchy. DOLFIN: Automated Finite Element Computing ANDERS LOGG, GARTH N. WELLS We describe here a library aimed at automating the solution of partial differential equations using the finite element method. By employing novel techniques for automated code generation, the library combines a high level of expressiveness with efficient computation. Finite element variational forms may be expressed in near mathematical notation, from which low-level code is automatically generated, compiled, and seamlessly integrated with efficient implementations of computational meshes and high-performance linear algebra. Easy-to-use object-oriented interfaces to the library are provided in the form of a C++ library and a Python module. This article discusses the mathematical abstractions and methods used in the design of the library and its implementation. A number of examples are presented to demonstrate the use of the library in application code. LOCAL BISECTION REFINEMENT FOR N-SIMPLICIAL GRIDS GENERATED BY REFLECTION JOSEPH M. MAUBACH A simple local bisection refinement algorithm for the adaptive refinement of n-simplicial grids is presented. The algorithm requires that the vertices of each simplex be ordered in a special way relative to those in neighboring simplices. It is proven that certain regular simplicial grids on [0, 1]^n have this property, and the more general grids to which this method is applicable are discussed. The edges to be bisected are determined by an ordering of the simplex vertices, without local or global computation or communication. Further, the number of congruency classes in a locally refined grid turns out to be bounded above by n, independent of the level of refinement. Simplicial grids of higher dimension are frequently used to approximate solution manifolds of parametrized equations, for instance, as in [W. C. Rheinboldt, Numer. Math., 53 (1988), pp. 165-180] and [E. Allgower and K. Georg, Utilitas Math., 16 (1979), pp. 123-129]. They are also used for the determination of fixed points of functions from R to Rn, as described in [M. J. Todd, Lecture Notes in Economic and Mathematical Systems, 124, Springer-Verlag, Berlin,1976]. In two and three dimensions, such grids of triangles, respectively, tetrahedrons, are used for the computation of finite element solutions of partial differential equations, for example, as in [O. Axelsson and V. A. Barker, Finite Element Solution ofBoundary Value Problems, Academic Press, Orlando, 1984], [R. E. Bank and B. D. Welfert, SIAM J. Numer. Anal., 28 1991), pp. 591-623], [W. E Mitchell, SIAMJ. Sci. Statist. Comput., 13 (1992), pp. 146-147], and [M. C. Rivara, J. Comput. Appl. Math., 36 (1991), pp. 79-89]. The new method is applicable to any triangular grid and may possibly be applied to many tetrahedral grids using additional closure refinement to avoid incompatibilities. A Comparison of HP Adaptive Strategies for Elliptic Partial Differential Equations WILLIAM F. MITCHELL and MARJORIE A. MCCLAIN The hp version of the finite element method (hp-FEM) combined with adaptive mesh refinement is a particularly efficient method for solving PDEs because it can achieve an exponential convergence rate in the number of degrees of freedom. hp-FEM allows for refinement in both the element size, h, and the polynomial degree, p. Like adaptive refinement for the h version of the finite element method, a posteriori error estimates can be used to determine where the mesh needs to be refined, but a single error estimate cannot simultaneously determine whether it is better to do the refinement by h or p. Several strategies for making this determination have been proposed over the years. These strategies are summarized, and the results of a numerical experiment to study the performance of these strategies is presented. It was found that the reference-solution-based methods are very effective, but also considerably more expensive, in terms of computation time, than other approaches. The method based on a priori knowledge is very effective when there are known point singularities. The method based on the decay rate of the expansion coefficients appears to be the best choice as a general strategy across all categories of problems, whereas many of the other strategies perform well in particular situations and are reasonable in general. A collection of 2D elliptic problems for testing adaptive grid refinement algorithms, William F. Mitchell Adaptive grid refinement is a critical component of the improvements that have recently been made in algorithms for the numerical solution of partial differential equations (PDEs). The development of new algorithms and computer codes for the solution of PDEs usually involves the use of proof-of-concept test problems. 2D elliptic problems are often used as the first test bed for new algorithms and codes. This paper contains a set of twelve parametrized 2D elliptic test problems for adaptive grid refinement algorithms and codes. The problems exhibit a variety of types of singularities, near singularities, and other difficulties. DATA OSCILLATION AND CONVERGENCE OF ADAPTIVE FEM, PEDRO MORIN, RICARDO H. NOCHETTO, AND KUNIBERT G. SIEBERT Data oscillation is intrinsic information missed by the averaging process associated with finite element methods (FEM) regardless of quadrature. Ensuring a reduction rate of data oscillation, together with an error reduction based on a posteriori error estimators, we construct a simple and efficient adaptive FEM for elliptic partial differential equations (PDEs) with linear rate of convergence without any preliminary mesh adaptation nor explicit knowledge of constants. Any prescribed error tolerance is thus achieved in a finite number of steps. A number of numerical experiments in two and three dimensions yield quasi-optimal meshes along with a competitive performance. WHEN ZOMBIES ATTACK!: MATHEMATICAL MODELLING OF AN OUTBREAK OF ZOMBIE INFECTION Philip Munz, Ioan Hudea, Joe Imad, Robert J. Smith Zombies are a popular figure in pop culture/entertainment and they are usually portrayed as being brought about through an outbreak or epidemic. Consequently, we model a zombie attack, using biological assumptions based on popular zombie movies. We introduce a basic model for zombie infection, determine equilibria and their stability, and illustrate the outcome with numerical solutions. We then refine the model to introduce a latent period of zombification, whereby humans are infected, but not infectious, before becoming undead. We then modify the model to include the effects of possible quarantine or a cure. Finally, we examine the impact of regular, impulsive reductions in the number of zombies and derive conditions under which eradication can occur. We show that only quick, aggressive attacks can stave off the doomsday scenario: the collapse of society as zombies overtake us all. Adaptive anisotropic meshing for steady convection-dominated problems Hoa Nguyen, Max Gunzburger, Lili Ju, John Burkardt Obtaining accurate solutions for convection–diffusion equations is challenging due to the presence of layers when convection dominates the diffusion. To solve this problem, we design an adaptive meshing algorithm which optimizes the alignment of anisotropic meshes with the numerical solution. Three main ingredients are used. First, the streamline upwind Petrov–Galerkin method is used to produce a stabilized solution. Second, an adapted metric tensor is computed from the approximate solution. Third, optimized anisotropic meshes are generated from the computed metric tensor by an anisotropic centroidal Voronoi tessellation algorithm. Our algorithm is tested on a variety of two-dimensional examples and the results shows that the algorithm is robust in detecting layers and efficient in avoiding non-physical oscillations in the numerical approximation. Optimizations for Quadrature Representations of Finite Element Tensors through Automated Code Generation KRISTIAN B. ØLGAARD, GARTH N. WELLS We examine aspects of the computation of finite element matrices and vectors that are made possible by automated code generation. Given a variational form in a syntax that resembles standard mathematical notation, the low-level computer code for building finite element tensors, typically matrices, vectors and scalars, can be generated automatically via a form compiler. In particular, the generation of code for computing finite element matrices using a quadrature approach is addressed. For quadrature representations, a number of optimization strategies which are made possible by automated code generation are presented. The relative performance of two different automatically generated representations of finite element matrices is examined, with a particular emphasis on complicated variational forms. It is shown that approaches which perform best for simple forms are not tractable for more complicated problems in terms of run-time performance, the time required to generate the code or the size of the generated code. The approach and optimizations elaborated here are effective for a range of variational forms. Benchmark solutions for the incompressible Navier-Stokes equations in general coordinates on staggered grids, C W Oosterlee, P Wesseling, A Segal, E Brakkee Benchmark problems are solved with the steady incompressible Navier-Stokes equations discretized with a finite volume method in general curvilinear co-ordinates on a staggered grid. The problems solved are skewed driven cavity problems, recently proposed as non-orthogonal grid benchmark problems. The system of discretized equations is solved efficiently with a non-linear multigrid algorithm, in which a robust line smoother is implemented. Furthermore, another benchmark problem is introduced and solved in which a 90 degree change in grid line direction occurs. Finite-element methods in electronic-structure theory J.E. Pask, B.M. Klein, P.A. Sterne, C.Y. Fong We discuss the application of the finite-element (FE) method to ab initio solid-state electronic-structure calculations. In this method, the basis functions are strictly local, piecewise polynomials. Because the basis is composed of polynomials, the method is completely general and its convergence can be controlled systematically. Because the basis functions are strictly local in real space, the method allows for variable resolution in real space; produces sparse, structured matrices, enabling the effective use of iterative solution methods; and is well suited to parallel implementation. The method thus combines the significant advantages of both real-space-grid and basis-oriented approaches and so promises to be particularly well suited for large, accurate ab initio calculations. We discuss the construction and properties of the required FE bases and develop in detail their use in the solution of the Schrödinger and Poisson equations subject to boundary conditions appropriate for a periodic solid. We present results for the Schrödinger equation illustrating the rapid, variational convergence of the method in electronic band-structure calculations. We present results for the Poisson equation illustrating the rapid convergence of the method, both pointwise and in the L2 norm, and its linear scaling with system size in the context of a model charge-density and Si pseudo-charge-density. Finally, we discuss the application of the method to large-scale ab initio positron distribution and lifetime calculations in solids and present results for a host of systems within the range of a conventional LMTO based approach for comparison, as well as results for systems well beyond the range of the conventional approach. The largest such calculation, involving a unit cell of 4092 atoms, was shown to be well within the range of the FE approach on existing computational platforms. Finite element methods in ab initio electronic structure calculations J E Pask and P A Sterne We review the application of the finite element (FE) method to ab initio electronic structure calculations in solids. The FE method is a general approach for the solution of differential and integral equations which uses a strictly local, piecewise-polynomial basis. Because the basis is composed of polynomials, the method is completely general and its accuracy is systematically improvable. Because the basis is strictly local in real space, the method allows for variable resolution in real space; produces sparse, structured matrices, enabling the effective use of iterative solution methods; and is well suited for parallel implementation. The method thus combines significant advantages of both real-space-grid and basis-oriented approaches, and so is well suited for large, accurate ab initio calculations. We reviewthe construction and properties of the required FE bases and their use in the self-consistent solution of the Kohn–Sham equations of density functional theory. Free Tools and Strategies for the Generation of 3D Finite Element Meshes: Modeling of the Cardiac Structures E. Pavarino, L. A. Neves, J. M. Machado, M. F. de Godoy, Y. Shiyou, J. C. Momente, G. F. D. Zafalon, A. R. Pinto, and C. R. Valêncio The Finite Element Method is a well-known technique, being extensively applied in different areas. Studies using the Finite Element Method (FEM) are targeted to improve cardiac ablation procedures. For such simulations, the finite element meshes should consider the size and histological features of the target structures. However, it is possible to verify that some methods or tools used to generate meshes of human body structures are still limited, due to nondetailed models, nontrivial preprocessing, ormainly limitation in the use condition. In this paper, alternatives are demonstrated to solidmodeling and automatic generation of highly refined tetrahedral meshes, with quality compatible with other studies focused on mesh generation.The innovations presented here are strategies to integrate Open Source Software (OSS). The chosen techniques and strategies are presented and discussed, considering cardiac structures as a first application context. Quantum simulation of materials at micron scales and beyond Qing Peng, Xu Zhang, Linda Hung, Emily A. Carter, and Gang Lu1 We present a multiscale modeling approach that can simulate multimillion atoms effectively via density functional theory. The method is based on the framework of the quasicontinuum QC approach with orbital free density-functional theory OFDFT as its sole energetics formulation. The local QC part is formulated by the Cauchy-Born hypothesis with OFDFT calculations for strain energy and stress. The nonlocal QC part is treated by an OFDFT-based embedding approach, which couples OFDFT nonlocal atoms to local region atoms. The method — QCDFT — is applied to a nanoindentation study of an Al thin film, and the results are compared to a conventional QC approach. The results suggest that QCDFT represents a new direction for the quantum simulation of materials at length scales that are relevant to experiments. A Simple Mesh Generator in MATLAB Per-Olof Persson, Gilbert Strang Creating a mesh is the first step in a wide range of applications, including scientific computing and computer graphics. An unstructured simplex mesh requires a choice of meshpoints (vertex nodes) and a triangulation. We want to offer a short and simple MATLAB code, described in more detail than usual, so the reader can experiment (and add to the code) knowing the underlying principles. We find the node locations by solving for equilibrium in a truss structure (using piecewise linear force-displacement relations) and reset the topol ogy by the Delaunay algorithm. The geometry is described implicitly by its distance function. In addition to being much shorter and simpler than other meshing techniques, our algorithm typically produces meshes of very high quality. We discuss ways to improve the robustness and the performance, but our aim here is simplicity. Readers can download (and edit) the codes from http://math.mit.edu/-persson/mesh. THE REDUCED BASIS METHOD FOR INCOMPRESSIBLE VISCOUS FLOW CALCULATIONS* JANET S. PETERSON The reduced basis method is a type of reduction method that can be used to solve large systems of nonlinear equations involving a parameter. In this work, the method is used in conjunction with a standard continuation technique to approximate the solution curve for the nonlinear equations resulting from discretizing the Navier-Stokes equations by finite-element methods. This paper demonstrates that the reduced basis method can be implemented to approximate efficiently solutions to incompressible viscous flows. Choices of basis vectors, issues concerning the implementation of the method, and numerical calculations are discussed. Two fluid flow calculations are considered, the driven cavity problem and flow over a forward facing step. Finite-Element Neural Networks for Solving Differential Equations Pradeep Ramuhalli, Lalita Udpa, and Satish S. Udpa The solution of partial differential equations (PDE) arises in a wide variety of engineering problems. Solutions to most practical problems use numerical analysis techniques such as finite- element or finite-difference methods. The drawbacks of these approaches include computational costs associated with the modeling of complex geometries. This paper proposes a finite-element neural network (FENN) obtained by embedding a finite-element model in a neural network architecture that enables fast and accurate solution of the forward problem. Results of applying the FENN to several simple electromagnetic forward and inverse problems are presented. Initial results indicate that the FENN performance as a forward model is comparable to that of the conventional finite-element method (FEM). The FENN can also be used in an iterative approach to solve inverse problems associated with the PDE. Results showing the ability of the FENN to solve the inverse problem given the measured signal are also presented. The parallel nature of the FENN also makes it an attractive solution for parallel implementation in hardware and software. Algorithm 792: Accuracy Tests of ACM Algorithms for Interpolation of Scattered Data in the Plane ROBERT J. RENKA, RON BROWN We present results of accuracy tests on scattered-data fitting methods that have been published as ACM algorithms. The algorithms include seven triangulation-based methods and three modified Shepard methods, two of which are new algorithms. Our purpose is twofold: to guide potential users in the selection of an appropriate algorithm and to provide a test suite for assessing the accuracy of new methods (or existing methods that are not included in this survey). Our test suite consists of five sets of nodes, with node counts ranging from 25 to 100, and 10 test functions. These are made available in the form of three Fortran subroutines: TESTDT returns one of the node sets; TSTFN1 returns a value and, optionally, a gradient value, of one of the test functions; and TSTFN2 returns a value, first partials, and second partial derivatives of one of the test functions. AUTOMATED GOAL-ORIENTED ERROR CONTROL I: STATIONARY VARIATIONAL PROBLEMS MARIE E. ROGNES, ANDERS LOGG This article presents a general and novel approach to the automation of goal-oriented error control in the solution of nonlinear stationary finite element variational problems. The approach is based on automated linearization to obtain the linearized dual problem, automated derivation and evaluation of a posteriori error estimates, and automated adaptive mesh refinement to control the error in a given goal functional to within a given tolerance. Numerical examples representing a variety of different discretizations of linear and nonlinear partial differential equations are presented, including Poisson’s equation, a mixed formulation of linear elasticity, and the incompressible Navier–Stokes equations. A Radial Basis Function (RBF)-Finite Difference (FD) Method for Diffusion and Reaction-Diffusion Equations on Surfaces Varun Shankar, Grady B. Wright, Robert M. Kirby, Aaron L. Fogelson In this paper, we present a method based on Radial Basis Function (RBF)-generated Finite Differences (FD) for numerically solving diffusion and reaction-diffusion equations (PDEs) on closed surfaces embedded in Rd. Our method uses a method-of-lines formulation, in which surface derivatives that appear in the PDEs are approximated locally using RBF interpolation. The method requires only scattered nodes representing the surface and normal vec- tors at those scattered nodes. All computations use only extrinsic coordinates, thereby avoiding coordinate distortions and singularities. We also present an optimization procedure that allows for the stabilization of the discrete differen- tial operators generated by our RBF-FD method by selecting shape parameters for each stencil that correspond to a global target condition number. We show the convergence of our method on two surfaces for different stencil sizes, and present applications to nonlinear PDEs simulated both on implicit/parametric surfaces and more general surfaces represented by point clouds. Analysis of Synaptic Transmission in the Neuromuscular Junction Using a Continuum Finite Element Model Jason L. Smart and J. Andrew McCammon There is a steadily growing body of experimental data describing the diffusion of acetylcholine in the neuromuscular junction and the subsequent miniature endplate currents produced at the postsynaptic membrane. To gain further insights into the structural features governing synaptic transmission, we have performed calculations using a simplified finite element model of the neuromuscular junction. The diffusing acetylcholine molecules are modeled as a continuum, whose spatial and temporal distribution is governed by the force-free diffusion equation. The finite element method was adopted because of its flexibility in modeling irregular geometries and complex boundary conditions. The resulting simulations are shown to be in accord with experiment and other simulations. Numerical methods for convection-diffusion problems or The 30 years war Martin Stynes Convection-diffusion problems arise in the modelling of many physical processes. Their typical solutions exhibit boundary and/or interior layers. Despite the linear nature of the differential operator, these problems pose still-unanswered questions to the numerical analyst. This talk will give a selective overview of numerical methods for the solution of convection-diffusion problems, while placing them in a historical context. It examines the principles that underpin the competing numerical techniques in this area and presents some recent developments. A posteriori error estimates for finite element discretizations of the heat equation R. Verfuerth We consider discretizations of the heat equation by A-stable θ- schemes in time and conforming finite elements in space. For these discretizations we derive residual a posteriori error indicators. The indicators yield upper bounds on the error which are global in space and time and yield lower bounds that are global in space and local in time. The ratio between upper and lower bounds is uniformly bounded in time and does not depend on any step-size in space or time. Moreover, there is no restriction on the relation between the step-sizes in space and time. Enhanced FEM-based modeling of brain shift deformation in Image-Guided Neurosurgery Lara M. Vigneron, Romain C. Boman, Jean-Philippe Ponthot, Pierre A. Robe, Simon K. Warfield, Jacques G. Verly We consider the problem of improving outcomes for neurosurgery patients by enhancing intraoperative navigation and guidance. Current navigation systems do not accurately account for intraoperative brain deformation. We focus on the brain shift deformation that occurs just after the opening of the skull and dura. The heart of our system is a nonrigid registration technique using a biomechanical model. We specifically work on two axes: the representation of the structures in the biomechanical model and the evaluation of the surface landmark displacement fields between intraoperative MR images. Using the modified Hausdorff distance as an image similarity measure, we demonstrate that our approach significantly improves the alignment of the intraoperative images. Large Scale Finite Element Analysis Via Assembly-Free Deflated Conjugate Gradient Praveen Yadav, Krishnan Suresh Large-scale finite element analysis (FEA) with millions of degrees of freedom (DOF) is becoming commonplace in solid mechanics. The primary computational bottleneck in such problems is the solution of large linear systems of equations. In this paper, we propose an assembly-free version of the deflated conjugate gradient (DCG) for solving such equations, where neither the stiffness matrix nor the deflation matrix is assembled. While assembly-free FEA is a well-known concept, the novelty pursued in this paper is the use of assembly-free deflation. The resulting implementation is particularly well suited for large-scale problems and can be easily ported to multicore central processing unit (CPU) and graphics-programmable unit (GPU) architectures. For demonstration, we show that one can solve a 50106 degree of freedom system on a single GPU card, equipped with 3 GB of memory. The second contribution is an extension of the rigid body agglomeration concept used in DCG to a curvature-sensitive agglomeration. The latter exploits classic plate and beam theories for efficient deflation of highly ill-conditioned problems arising from thin structures.