High-order Mesh Research Articles

Measuring and improving the geometric accuracy of piece-wise polynomial boundary meshes

We present a new disparity functional to measure and improve the geometric accuracy of a curved high-order mesh that approximates a target geometry model. We have devised the disparity to account for compound models, be independent of the entity parameterization, and allow trimmed entities. The disparity depends on the physical mesh and the auxiliary parametric meshes. Since it is two times differentiable on all these variables, we can minimize it with a second-order method. Its minimization with the parametric meshes as design variables measures the geometric accuracy of a given mesh. Furthermore, the minimization with both the physical and parametric meshes as design variables improves the geometric accuracy of an initial mesh. We have numerical evidence that the obtained meshes converge to the target geometry (unitary normal) algebraically, in terms of the element size, with order 2p (2p−1, respectively), where p is the polynomial degree of the mesh. Although we obtain meshes with non-interpolative boundary nodes, we propose a post-process to enforce, if required by the application, meshes with interpolative boundary nodes and featuring the same order of geometric accuracy. In conclusion, we can obtain super-convergent orders, at least for sufficiently smooth parametric curve (surface) entities, for meshes of polynomial degrees up to 4 (3, respectively). In perspective, this super-convergence might enable using a lower polynomial degree to approximate the geometry than to approximate the solution without hampering the required geometric accuracy for high-order analysis.

Hr-Adaptivity for nonconforming high-order meshes with the target matrix optimization paradigm

We present an $hr$-adaptivity framework for optimization of high-order meshes. This work extends the $r$-adaptivity method for mesh optimization by Dobrev et al., where we utilized the Target-Matrix Optimization Paradigm (TMOP) to minimize a functional that depends on each element's current and target geometric parameters: element aspect-ratio, size, skew, and orientation. Since fixed mesh topology limits the ability to achieve the target size and aspect-ratio at each position, in this paper we augment the $r$-adaptivity framework with nonconforming adaptive mesh refinement to further reduce the error with respect to the target geometric parameters. The proposed formulation, referred to as $hr$-adaptivity, introduces TMOP-based quality estimators to satisfy the aspect-ratio-target via anisotropic refinements and size-target via isotropic refinements in each element of the mesh. The methodology presented is purely algebraic, extends to both simplices and hexahedra/quadrilaterals of any order, and supports nonconforming isotropic and anisotropic refinements in 2D and 3D. Using a problem with a known exact solution, we demonstrate the effectiveness of $hr$-adaptivity over both $r-$ and $h$-adaptivity in obtaining similar accuracy in the solution with significantly fewer degrees of freedom. We also present several examples that show that $hr$-adaptivity can help satisfy geometric targets even when $r$-adaptivity fails to do so, due to the topology of the initial mesh.

A new finite element method using MATLAB for Poisson’s Equation in axisymmetric domain

This paper demonstrates finite element procedure for two-dimensional axisymmetric domains. For many engineering applications like structural engineering, aerospace engineering, geo-mechanics etc., the solution domain and boundary conditions are axisymmetric. Henceforth, we can illuminate just the axisymmetric part of the solution domain that gives the data of the entire domain. This paper demonstrates the effectiveness of using MATLAB programming demonstrated by Persson et.al (2004) as the initial mesh for discretization of axisymmetric domains for higher order meshing. Further, solving some class of partial differential equations using finite element method with nodal relation given by subparametric transformations Rathod et.al (2008). In this paper a cubic order curved triangular meshing for some of the domains like ellipse and circle are demonstrated. These in turn finds its applications in the fields like stress analysis in mechanical engineering, torsion twist (shear strength) analysis in civil engineering, evaluation of stress intensity factor for quarter elliptical crack in pressure vessels in equipment industry etc,. The output data from the meshing scheme like meshing of the domain, nodal position, element connectivity and boundary edges are been used in the finite element procedures. The efficiency of the method is achieved by p-refinement scheme i.e., fixing the number of elements and increasing the polynomial order.

Implicit large Eddy simulation of the NASA CRM high-lift configuration near stall

Recent progress in adaptive high-order numerical methods, high-order mesh generation, and scalable implicit time integration approaches has made large eddy simulation (LES) more affordable than ever before. In the present study, we demonstrate the progress by performing implicit LES of a benchmark problem from the 3rd AIAA High-Lift Prediction Workshop: NASA's Common Research Model (CRM) high-lift configuration. The LES tool, hpMusic, is based on the flux reconstruction method capable of handling high-order mixed unstructured meshes. P-refinement studies are used to assess the mesh-dependence of the numerical solutions. Comparisons are also made with available experimental data, and a very good agreement was obtained in lift and drag coefficients between the present computation and wall-corrected experimental data.

2D Image-based higher-order meshing for 3D modelling in MATLAB

We present a higher-order meshing technique to 3D modelling applications from 2D images. This method uses 2D images that are converted to 3D stereolithography files on which tetrahedral meshes of higher order are developed. The resulting 3D meshed structures find several applications in Computational fluid dynamics and finite element analysis problems. The generateMesh is the initial mesh used by the presented mesh generator, which is a MATLAB function from the PDE Toolbox. The required inputs for the mesh generator are stereolithography file and mesh edge length upper bound of the 3D structure. A higher-order tetrahedral mesh, connectivity matrix, node coordinates, and nodes lying on the border are the results we get from the presented work. For many finite element analysis problems, the outcomes obtained from the presented mesh generator can be used very efficiently to enhance the numerical results. In medical simulation, photogrammetric fields, material engineering, computational materials science, computational electromagnetics, etc., this method can be used.

A NURBS-based discontinuous Galerkin method for conservation laws with high-order moving meshes

The objective of the present work is to develop a new numerical framework for simulations involving deformable domains, in the specific context of high-order meshes consistent with Computer-Aided Design (CAD) representations. Thus, the proposed approach combines ideas from isogeometric analysis, able to handle exactly CAD-based geometries, and Discontinuous Galerkin (DG) methods with an Arbitrary Lagrangian-Eulerian (ALE) formulation, able to solve complex problems with moving grids. The resulting approach is a DG method based on rational Bézier elements, that can be easily constructed from Non-Uniform Rational B-Splines (NURBS), formulated in a general ALE setting. We focus here on applications in two-dimensional compressible flows, but the method could be applied to other models as well. Two verification exercises are conducted, to assess rigorously the properties of the method and the convergence rates for representations up to sixth order. Finally, three problems are analysed in depth, involving compressible Euler and Navier-Stokes equations, for an oscillating cylinder and a pitching airfoil. In particular, the convergence of flow characteristics is investigated, as well as the impact of using curved boundaries in the context of deformable domains.

Curved p-version C1 finite elements for the finite deformation analysis of isotropic and composite laminated thin shells in irregular shape

In this work, we generalize our recently developed p-version curved C1 plate elements to the finite deformation analysis of general isotropic and laminated thin shells with arbitrary curvature and shape. The kinematic formulations are based on the geometrically-exact nonlinear thin shell model where the Kirchhoff-Love hypothesis is strictly implemented during the deformation. Since high accurate curved coarse meshes are required in the implementation of p-version finite element method, the quasi regional mapping technique is employed and modified to obtain the high-order meshes with fitted boundaries. Through this method, shells with irregular shapes can be easily modeled. To avoid dealing the difficult problem of constructing fully C1 conforming displacement fields on shell surface, the C1 conformity of present elements is implemented on the selected Gauss-Lobatto nodes, which had been found to be able to produce very fast convergence. Nodal variables including second-order directional derivatives on surfaces are derived and used for element assembly. In this manner, direct construction of displacement fields on shell mid-surface defined by the quasi regional mapping is allowed. Several benchmark problems of large deformation analysis of isotropic and laminated thin shells are presented to illustrate the accuracy as well as robustness of the present method.

HDGlab: An Open-Source Implementation of the Hybridisable Discontinuous Galerkin Method in MATLAB

This paper presents HDGlab, an open source MATLAB implementation of the hybridisable discontinuous Galerkin (HDG) method. The main goal is to provide a detailed description of both the HDG method for elliptic problems and its implementation available in HDGlab. Ultimately, this is expected to make this relatively new advanced discretisation method more accessible to the computational engineering community. HDGlab presents some features not available in other implementations of the HDG method that can be found in the free domain. First, it implements high-order polynomial shape functions up to degree nine, with both equally-spaced and Fekete nodal distributions. Second, it supports curved isoparametric simplicial elements in two and three dimensions. Third, it supports non-uniform degree polynomial approximations and it provides a flexible structure to devise degree adaptivity strategies. Finally, an interface with the open-source high-order mesh generator Gmsh is provided to facilitate its application to practical engineering problems.

On pixel-exact rendering for high-order mesh and solution

With the increasing use of high-order methods and high-order meshes, scientific visualization software need to adapt themselves to reliably render the associated meshes and numerical solutions. In this paper, a novel approach, based on OpenGL 4 framework, enables a GPU-based rendering of high-order meshes as well as an almost pixel-exact rendering of high-order solutions. Several aspects of the OpenGL Shading Language and in particular the use of dedicated shaders (GPU programs) allows to answer this visualization challenge. Fragment shaders are used to compute the exact solution for each pixel, made possible by the transfer of degrees of freedom and shape functions to the GPU with textures. Tessellation shaders, combined with geometric error estimates, allow us to render high-order curved meshes by providing an adaptive subdivision of elements on the GPU directly. A convenient way to compute bounds for high-order solutions is described. The interest of using Bézier basis instead of Lagrange functions lies in the existence of fast and robust evaluation of polynomial functions with de Casteljau algorithm. A technique to plot highly nonlinear isolines and wire frames with a desired thickness is derived. It is based on a finite difference scheme performed on GPU. In comparison with standard techniques, we remove the use of any linear interpolation step and the need to generate a priori a fixed subdivided mesh. This reduces the memory footprint, improves the accuracy and the speed of the rendering. Finally, the method is illustrated with various 3D examples.

Spectral Mesh-Free Quadrature for Planar Regions Bounded by Rational Parametric Curves

This work presents spectral, mesh-free, Green’s theorem-based numerical quadrature schemes for integrating functions over planar regions bounded by rational parametric curves. Our algorithm proceeds in two steps: (1) We first find intermediate quadrature rules for line integrals along the region’s boundary curves corresponding to Green’s theorem. (2) We then use a high-order quadrature rule to compute the numerical antiderivative of the integrand along a coordinate axis, which is used to evaluate the Green’s theorem line integral. We present two methods to compute the intermediate quadrature rule. The first is spectrally accurate (it converges faster than any algebraic order with respect to number of quadrature points) and is relatively easy to implement, but has no guarantee of polynomial exactness. The second guarantees exactness for polynomial integrands up to a pre-specified degree k with an a priori-known number of quadrature points and retains the convergence properties of the first, but is slightly more complicated. The quadrature schemes have applications to computation of geometric moments, immersogeometric analysis, conservative field transfer between high-order meshes, and initialization of multi-material simulations with rational geometry. We compare the quadrature schemes produced using our method to other methods in the literature and show that they are much more efficient both in terms of number of quadrature points and computational time. We provide an open-source implementation of the algorithm in MATLAB.

Higher-order mesh generation using CISAMR: A case study on bias in presentation and interpretation of results

A modified version of the Conforming to Interface Structured Adaptive Mesh Refinement (CISAMR) algorithm is presented for the construction of higher-order Lagrangian and NURBS-enhanced (NE) finite element meshes. CISAMR non-iteratively transforms a structured grid into a conforming mesh with an upper bound of three on elements aspect ratio. In this work, we introduce new algorithmic aspects for generating higher-order Lagrangian and NE meshes using CISAMR. For each type of element, a comprehensive study is then provided on the performance of first, second, and third-order meshes for solving linear elastic problems with smooth and oscillatory curvilinear edges. Local gradient recovery error and computational cost are used as performance metrics in these examples. Outcomes are then used as a case study to elucidate common sources of bias in presenting and interpreting results. In particular, we show how bias, whether purposeful or not, could lead to misleading conclusions regarding the performance of a method.

Simulation-driven optimization of high-order meshes in ALE hydrodynamics

In this paper we propose tools for high-order mesh optimization and demonstrate their benefits in the context of multi-material Arbitrary Lagrangian-Eulerian (ALE) compressible shock hydrodynamic applications. The mesh optimization process is driven by information provided by the simulation which uses the optimized mesh, such as shock positions, material regions, known error estimates, etc. These simulation features are usually represented discretely, for instance, as finite element functions on the Lagrangian mesh. The discrete nature of the input is critical for the practical applicability of the algorithms we propose and distinguishes this work from approaches that strictly require analytical information. Our methods are based on node movement through a high-order extension of the Target-Matrix Optimization Paradigm (TMOP) of [1]. The proposed formulation is fully algebraic and relies only on local Jacobian matrices, so it is applicable to all types of mesh elements, in 2D and 3D, and any order of the mesh. We discuss the notions of constructing adaptive target matrices and obtaining their derivatives, reconstructing discrete data in intermediate meshes, node limiting that enables improvement of global mesh quality while preserving space-dependent local mesh features, and appropriate normalization of the objective function. The adaptivity methods are combined with automatic ALE triggers that can provide robustness of the mesh evolution and avoid excessive remap procedures. The benefits of the new high-order TMOP technology are illustrated on several simulations performed in the high-order ALE application BLAST [2].

A matrix-free hyperviscosity formulation for high-order ALE hydrodynamics

The numerical approximation of compressible hydrodynamics is at the core of high-energy density (HED) multiphysics simulations as shocks are the driving force in experiments like inertial confinement fusion (ICF). In this work, we describe our extension of the hyperviscosity technique, originally developed for shock treatment in finite difference simulations, for use in arbitrarily high-order finite element methods for Lagrangian hydrodynamics. Hyperviscosity enables shock capturing while preserving the high-order properties of the underlying discretization away from the shock region. Specifically, we compute a high-order term based on a product of the mesh length scale to a high power scaled by a hyper-Laplacian operator applied to a scalar field. We then form the total artificial viscosity by taking a non-linear blend of this term and a traditional artificial viscosity term. We also present a matrix-free formulation for computing the finite element based hyper-Laplacian operator. Such matrix-free methods have superior performance characteristics compared to traditional full matrix assembly approaches and offer advantages for GPU based HPC hardware. We demonstrate the numerical convergence of our method and its application to complex, multi-material ALE simulations on high-order (curved) meshes.

Direct reconstruction method for discontinuous Galerkin methods on higher-order mixed-curved meshes II. Surface integration

As the second extension of the direct reconstruction method (DRM), DRM for the surface integration of the discontinuous Galerkin (DG) weak formulation is presented on multi-dimensional high-order mixed meshes. From the previous study, DRM was successfully extended to the volume integration of the DG method, and produced a significant reduction in both the computational cost and memory overhead, particularly for high-order solution approximations on 3-D high-order meshes. DRM has the potential to reduce the computational cost and memory overhead further by treating the surface integration term consistently. To realize this, we design a linear dependency detector and formulate an adaptive orthonormalization process from the modified Gram-Schmidt (MGS) process. DRM is then applied to the restricted function space of the surface integration in an adaptive and optimal manner. The resulting methods, namely the DRM surface integration with the brute force points (BFP) and shape function points (SFP) methods, perform very well on high-order mixed-curved meshes. In addition, the adaptive MGS process is robust and accurate in extracting linear faces/edges from a given set of mixed-curved meshes. Various benchmark tests from the compressible Euler and Navier-Stokes equations verify the impact of DRM for 3-D large-scale computations on mixed-curved meshes. As an example, the DRM framework applied to the 3-D surface and volume integrations with DG-P3 approximation on P3-mesh may yield O(10)× speed-up, compared to conventional quadrature-based methods.

Technologies for supporting high-order geodesic mesh frameworks for computational astrophysics and space sciences

Many important problems in astrophysics, space physics, and geophysics involve flows of (possibly ionized) gases in the vicinity of a spherical object, such as a star or planet. The geometry of such a system naturally favors numerical schemes based on a spherical mesh. Despite its orthogonality property, the polar (latitude-longitude) mesh is ill suited for computation because of the singularity on the polar axis, leading to a highly non-uniform distribution of zone sizes. The consequences are (a) loss of accuracy due to large variations in zone aspect ratios, and (b) poor computational efficiency from a severe limitations on the time stepping. Geodesic meshes, based on a central projection using a Platonic solid as a template, solve the anisotropy problem, but increase the complexity of the resulting computer code. We describe a new finite volume implementation of Euler and MHD systems of equations on a triangular geodesic mesh (TGM) that is accurate up to fourth order in space and time and conserves the divergence of magnetic field to machine precision. The paper discusses in detail the generation of a TGM, the domain decomposition techniques, three-dimensional conservative reconstruction, and time stepping.

Implicit shock tracking using an optimization-based high-order discontinuous Galerkin method

A novel framework for resolving discontinuous solutions of conservation laws, e.g., contact lines, shock waves, and interfaces, using implicit tracking and a high-order discontinuous Galerkin (DG) discretization was introduced in [39]. Central to the framework is an optimization problem whose solution is a discontinuity-aligned mesh and the corresponding high-order approximation to the flow that does not require explicit meshing of the unknown discontinuity surface. The method was shown to deliver highly accurate solutions on coarse, high-order discretizations without nonlinear stabilization and recover optimal convergence rates O(hp+1) even for problems with discontinuous solutions. This work extends the implicit tracking framework such that robustness is improved and convergence accelerated. In particular, we introduce an improved formulation of the central optimization problem and an associated sequential quadratic programming (SQP) solver. The new error-based objective function penalizes violation of the DG residual in an enriched test space and is shown to have excellent tracking properties. The SQP solver simultaneously converges the nodal coordinates of the mesh and DG solution to their optimal values and is equipped with a number of features to ensure robust, fast convergence: Levenberg-Marquardt approximation of the Hessian with weighted elliptic regularization, backtracking line search based on the ℓ1 merit function, and rigorous convergence criteria. We use the proposed method to solve a range of inviscid conservation laws of varying difficulty. We show the method is able to deliver accurate solutions on coarse, high-order meshes and the SQP solver is robust and usually able to drive the first-order optimality system to tight tolerances.

Discontinuous Galerkin via Interpolation: The Direct Flux Reconstruction Method

The discontinuous Galerkin (DG) method is based on the idea of projection using integration. The recent direct flux reconstruction (DFR) method by Romero et al. (J Sci Comput 67(1):351–374, 2016) is derived via interpolation and results in a scheme identical to DG (on hexahedral meshes). The DFR method is further studied and developed here. Two proofs for its equivalence with the DG scheme considerably simpler than the original proof are presented. The first proof employs the $$ 2K - 1 $$ degree of precision by a $$ K $$-point Gauss quadrature. The second shows the equivalence of DG, FR, and DFR by using the property that the derivative of the degree $$ K + 1 $$ Lobatto polynomial vanishes at the $$ K $$ Gauss points. Fourier analysis for these schemes are presented using an approach more geometric compared with existing analytic approaches. The effects of nonuniform mesh and those of high-order mesh transformation (a precursor for curved meshes in two and three spatial dimensions) on stability and accuracy are examined. These nonstandard analyses are obtained via an in-depth study of the behavior of eigenvalues and eigenvectors.

Automated higher-order mesh generation in MATLAB on biomaterials

We propose a method for automated higher-order triangle and tetrahedral mesh generation in MATLAB on biomaterials. The proposed mesh generator is based on the MATLAB function generate Mesh from the Partial Differential Equation Toolbox. Here, the inputs are STL file and an approximate upper bound on mesh edge length for the geometry under consideration using the higher-order triangle or tetrahedral elements. As, output we get a higher-order triangle or tetrahedral meshes, connectivity matrix, coordinates of the nodes, and boundary nodes. These outputs can be efficiently used for solving any partial differential equations by finite element method which arises in many applications like biomaterials, computational materials science, crystal plasticity, materials engineering, mechanics of materials etc.

On using geometry-based high-order finite element approximations for mid-frequency response of fluid-loaded elastic structures

Accurate solution of the exterior scattering and radiation problems in the mid-frequency regime remains of high interest and a computationally challenging task. For wave-dependent problems with large dispersion errors, high-order approximations provide several advantages, including increased accuracy and error convergence rate compared to low-order methods. The effective use of high-order finite element type approximations require special attention to geometry. In particular, one must use high-order (curved) meshes in order to preserve the accuracy and convergence advantages. We present recent developments in the use of geometry-based high-order approximations for computing scattering response of fluid-loaded elastic structures submerged in an acoustic medium. We present examples with validation against experimental data to demonstrate the effectiveness of these techniques as well as highlight future technical challenges to push to the higher ends of the mid-frequency regime.Accurate solution of the exterior scattering and radiation problems in the mid-frequency regime remains of high interest and a computationally challenging task. For wave-dependent problems with large dispersion errors, high-order approximations provide several advantages, including increased accuracy and error convergence rate compared to low-order methods. The effective use of high-order finite element type approximations require special attention to geometry. In particular, one must use high-order (curved) meshes in order to preserve the accuracy and convergence advantages. We present recent developments in the use of geometry-based high-order approximations for computing scattering response of fluid-loaded elastic structures submerged in an acoustic medium. We present examples with validation against experimental data to demonstrate the effectiveness of these techniques as well as highlight future technical challenges to push to the higher ends of the mid-frequency regime.

A high-order conservative remap for discontinuous Galerkin schemes on curvilinear polygonal meshes

A data transfer (called later remap) of physical fields between two meshes is an important step of arbitrary Lagrangian-Eulerian (ALE) simulations. This step is challenging for high-order discontinuous Galerkin schemes since the Lagrangian flow motion leads to high-order meshes with curved faces. It becomes even more challenging for unstructured polygonal meshes that do not have a polynomial map from the reference to a current cell. We propose and analyze a new framework to create remap schemes on curvilinear polygonal meshes based on the theory of virtual element projectors. We derive a conservative remap scheme that is high-order accurate in space and time. The properties of this scheme are studied numerically for smooth and discontinuous fields on unstructured quadrilateral and polygonal meshes.