Search

The recently proposed soft finite element method (SoftFEM) reduces the stiffness (condition numbers), consequently improving the overall approximation accuracy. The method subtracts a least-square term that penalizes the gradient jumps across mesh interfaces from the FEM stiffness bilinear form while maintaining the system's coercivity. Herein, we present two generalizations for SoftFEM that aim to improve the approximation accuracy and further reduce the discrete systems' stiffness. Firstly and most naturally, we generalize SoftFEM by adding a least-square term to the mass bilinear form. Superconvergent results of rates $h^6$ and $h^8$ for eigenvalues are established for linear uniform elements; $h^8$ is the highest order of convergence known in the literature. Secondly, we generalize SoftFEM by applying the blended Gaussian-type quadratures. We demonstrate further reductions in stiffness compared to traditional FEM and SoftFEM. The coercivity and analysis of the optimal error convergences follow the work of SoftFEM. Thus, this paper focuses on the numerical study of these generalizations. For linear and uniform elements, analytical eigenpairs, exact eigenvalue errors, and superconvergent error analysis are established. Various numerical examples demonstrate the potential of generalized SoftFEMs for spectral approximation, particularly in high-frequency regimes.

The Schr\"{o}dinger equation with random potentials is a fundamental model for understanding the behaviour of particles in disordered systems. Disordered media are characterised by complex potentials that lead to the localisation of wavefunctions, also called Anderson localisation. These wavefunctions may have similar scales of eigenenergies which poses difficulty in their discovery. It has been a longstanding challenge due to the high computational cost and complexity of solving the Schr\"{o}dinger equation. Recently, machine-learning tools have been adopted to tackle these challenges. In this paper, based upon recent advances in machine learning, we present a novel approach for discovering localised eigenstates in disordered media using physics-informed neural networks (PINNs). We focus on the spectral approximation of Hamiltonians in one dimension with potentials that are randomly generated according to the Bernoulli, normal, and uniform distributions. We introduce a novel feature to the loss function that exploits known physical phenomena occurring in these regions to scan across the domain and successfully discover these eigenstates, regardless of the similarity of their eigenenergies. We present various examples to demonstrate the performance of the proposed approach and compare it with isogeometric analysis.

Sea ice profoundly influences the polar environment and the global climate. Traditionally, Sea ice has been modeled as a continuum under Eulerian coordinates to describe its large-scale features, using, for instance, viscous-plastic rheology. Recently, Lagrangian particle models, also known as the discrete element method (DEM) models, have been utilized for characterizing the motion of individual sea ice fragments (called floes) at scales of 10 km and smaller, especially in marginal ice zones. This paper develops a multiscale model that couples the particle and the continuum systems to facilitate an effective representation of the dynamical and statistical features of sea ice across different scales. The multiscale model exploits a Boltzmann-type system that links the particle movement with the continuum equations. For the small-scale dynamics, it describes the motion of each sea ice floe. Then, as the large-scale continuum component, it treats the statistical moments of mass density and linear and angular velocities. The evolution of these statistics affects the motion of individual floes, which in turn provides bulk feedback that adjusts the large-scale dynamics. Notably, the particle model characterizing the sea ice floes is localized and fully parallelized, in a framework that is sometimes called superparameterization, which significantly improves computation efficiency. Numerical examples demonstrate the effective performance of the multiscale model. Additionally, the study demonstrates that the multiscale model has a linear-order approximation to the truth model.

In this work, we investigate the spectral problem $Au = {\lambda}u$ where $A$ is a fractional elliptic operator involving left- and right-sided Riemann-Liouville derivatives. These operators are nonlocal and nonsymmetric, however, share certain classic elliptic properties. The eigenvalues correspond to the roots of a class of certain special functions. Compared with classic Sturm-Liouville problems, the most challenging part is to set up the framework for analyzing these nonlocal operators, which requires developing new tools. We prove the existence of the real eigenvalues, find the range for all possible complex eigenvalues, explore the graphs of eigenfunctions, and show numerical findings on the distribution of eigenvalues on the complex plane.

The study of quantum three-body problems has been centered on low-energy states that rely on accurate numerical approximation. Recently, isogeometric analysis (IGA) has been adopted to solve the problem as an alternative but more robust (with respect to atom mass ratios) method that outperforms the classical Born-Oppenheimer (BO) approximation. In this paper, we focus on the performance of IGA and apply the recently-developed softIGA to reduce the spectral errors of the low-energy bound states. The main idea is to add high-order derivative-jump terms with a penalty parameter to the IGA bilinear forms. With an optimal choice of the penalty parameter, we observe eigenvalue error superconvergence. We focus on linear (finite elements) and quadratic elements and demonstrate the outperformance of softIGA over IGA through a variety of examples including both two- and three-body problems in 1D.

We extend the softFEM idea to isogeometric analysis (IGA) to reduce the stiffness (consequently, the condition numbers) of the IGA discretized problem. We refer to the resulting approximation technique as softIGA. We obtain the resulting discretization by first removing the IGA spectral outliers to reduce the system's stiffness. We then add high-order derivative-jump penalization terms (with negative penalty parameters) to the standard IGA bilinear forms. The penalty parameter seeks to minimize spectral/dispersion errors while maintaining the coercivity of the bilinear form. We establish dispersion errors for both outlier-free IGA (OF-IGA) and softIGA elements. We also derive analytical eigenpairs for the resulting matrix eigenvalue problems and show that the stiffness and condition numbers of the IGA systems significantly improve (reduce). We prove a superconvergent result of order $h^{2p+4}$ for eigenvalues where $h$ characterizes the mesh size and $p$ specifies the order of the B-spline basis functions. To illustrate the main idea and derive the analytical results, we focus on uniform meshes in 1D and tensor-product meshes in multiple dimensions. For the eigenfunctions, softIGA delivers the same optimal convergence rates as the standard IGA approximation. Various numerical examples demonstrate the advantages of softIGA over IGA.

Elliptic interface problems whose solutions are $C^0$ continuous have been well studied over the past two decades. The well-known numerical methods include the strongly stable generalized finite element method (SGFEM) and immersed FEM (IFEM). In this paper, we study numerically a larger class of elliptic interface problems where their solutions are discontinuous. A direct application of these existing methods fails immediately as the approximate solution is in a larger space that covers discontinuous functions. We propose a class of high-order enriched unfitted FEMs to solve these problems with implicit or Robin-type interface jump conditions. We design new enrichment functions that capture the imposed discontinuity of the solution while keeping the condition number from fast growth. A linear enriched method in 1D was recently developed using one enrichment function and we generalized it to an arbitrary degree using two simple discontinuous one-sided enrichment functions. The natural tensor product extension to the 2D case is demonstrated. Optimal order convergence in the $L^2$ and broken $H^1$-norms are established. We also establish superconvergence at all discretization nodes (including exact nodal values in special cases). Numerical examples are provided to confirm the theory. Finally, to prove the efficiency of the method for practical problems, the enriched linear, quadratic, and cubic elements are applied to a multi-layer wall model for drug-eluting stents in which zero-flux jump conditions and implicit concentration interface conditions are both present.

Explicit time-marching schemes are popular for solving time-dependent partial differential equations; one of the biggest challenges these methods suffer is increasing the critical time-marching step size that guarantees numerical stability. In general, there are two ways to increase the critical step size. One is to reduce the stiffness of the spatially discretized system, while the other is to design time-marching schemes with larger stability regions. In this paper, we focus on the recently proposed explicit generalized-$\alpha$ method for second-order hyperbolic equations and increase the critical step size by reducing the stiffness of the isogeometric-discretized system. In particular, we apply boundary penalization to lessen the system's stiffness. For $p$-th order $C^{p-1}$ isogeometric elements, we show numerically that the critical step size increases by a factor of $\sqrt{\frac{p^2-3p+6}{4}}$, which indicates the advantages of using the proposed method, especially for high-order elements. Various examples in one, two, and three dimensions validate the performance of the proposed technique.

In this paper, we initiate the study of isogeometric analysis (IGA) of a quantum three-body problem that has been well-known to be difficult to solve. In the IGA setting, we represent the wavefunctions by linear combinations of B-spline basis functions and solve the problem as a matrix eigenvalue problem. The eigenvalue gives the eigenstate energy while the eigenvector gives the coefficients of the B-splines that lead to the eigenstate. The major difficulty of isogeometric or other finite-element-method-based analyses lies in the lack of boundary conditions and a large number of degrees of freedom for accuracy. For a typical many-body problem with attractive interaction, there are bound and scattering states where bound states have negative eigenvalues. We focus on bound states and start with the analysis for a two-body problem. We demonstrate through various numerical experiments that IGA provides a promising technique to solve the three-body problem.

The discrete element method (DEM) is providing a new modeling approach for describing sea ice dynamics. It exploits particle-based methods to characterize the physical quantities of each sea ice floe along its trajectory under Lagrangian coordinates. One major challenge in applying the DEM models is the heavy computational cost when the number of floes becomes large. In this paper, an efficient Lagrangian parameterization algorithm is developed, which aims at reducing the computational cost of simulating the DEM models while preserving the key features of the sea ice. The new parameterization takes advantage of a small number of artificial ice floes, named the superfloes, to effectively approximate a considerable number of the floes, where the parameterization scheme satisfies several important physics constraints. The physics constraints guarantee the superfloe parameterized system will have similar short-term dynamical behavior as the full system. These constraints also allow the superfloe parameterized system to accurately quantify the long-range uncertainty, especially the non-Gaussian statistical features, of the full system. In addition, the superfloe parameterization facilitates a systematic noise inflation strategy that significantly advances an ensemble-based data assimilation algorithm for recovering the unobserved ocean field underneath the sea ice. Such a new noise inflation method avoids ad hoc tunings as in many traditional algorithms and is computationally extremely efficient. Numerical experiments based on an idealized DEM model with multiscale features illustrate the success of the superfloe parameterization in quantifying the uncertainty and assimilating both the sea ice and the associated ocean field.

It is well-known that outliers appear in the high-frequency region in the approximate spectrum of isogeometric analysis of the second-order elliptic operator. Recently, the outliers have been eliminated by a boundary penalty technique. The essential idea is to impose extra conditions arising from the differential equation at the domain boundary. In this paper, we extend the idea to remove outliers in the superconvergent approximate spectrum of isogeometric analysis with optimally-blended quadrature rules. We show numerically that the eigenvalue errors are of superconvergence rate $h^{2p+2}$ and the overall spectrum is outlier-free. The condition number and stiffness of the resulting algebraic system are reduced significantly. Various numerical examples demonstrate the performance of the proposed method.

We propose a new class of high-order time-marching schemes with dissipation user-control and unconditional stability for parabolic equations. High-order time integrators can deliver the optimal performance of highly-accurate and robust spatial discretizations such as isogeometric analysis. The generalized-$\alpha$ method delivers unconditional stability and second-order accuracy in time and controls the numerical dissipation in the discrete spectrum's high-frequency region. Our goal is to extend the generalized-$alpha$ methodology to obtain a high-order time marching methods with high accuracy and dissipation in the discrete high-frequency range. Furthermore, we maintain the stability region of the original, second-order generalized-$alpha$ method foe the new higher-order methods. That is, we increase the accuracy of the generalized-$\alpha$ method while keeping the unconditional stability and user-control features on the high-frequency numerical dissipation. The methodology solve $k>1, k\in \mathbb{N}$ matrix problems and updates the system unknowns, which correspond to higher-order terms in Taylor expansions to obtain $(3/2k)^{th}$-order method for even $k$ and $(3/2k+1/2)^{th}$-order for odd $k$. A single parameter $\rho^\infty$ controls the dissipation, and the update procedure follows the formulation of the original second-order method. Additionally, we show that our method is A-stable and setting $\rho^\infty=0$ allows us to obtain an L-stable method. Lastly, we extend this strategy to analyze the accuracy order of a generic method.

In this paper, we discuss the application of the Generalized Finite Element Method (GFEM) to approximate the solutions of quasilinear elliptic equations with multiple interfaces in one dimensional space. The problem is characterized by spatial discontinuity of the elliptic coefficient that depends on the unknown solution. It is known that unless the partition of the domain matches the discontinuity configuration, accuracy of standard finite element techniques significantly deteriorates and standard refinement of the partition may not suffice. The GFEM is a viable alternative to overcome this predicament. It is based on the construction of certain enrichment functions supplied to the standard space that capture effects of the discontinuity. This approach is called stable (SGFEM) if it maintains an optimal rate of convergence and the conditioning of GFEM is not worse than that of the standard FEM. A convergence analysis is derived and performance of the method is illustrated by several numerical examples. Furthermore, it is known that typical global formulations such as FEMs do not enjoy the numerical local conservation property that is crucial in many conservation law-based applications. To remedy this issue, a Lagrange multiplier technique is adopted to enforce the local conservation. A numerical example is given to demonstrate the performance of proposed technique., Comment: 32 pages, 32 figures

We propose, analyze mathematically, and study numerically a novel approach for the finite element approximation of the spectrum of second-order elliptic operators. The main idea is to reduce the stiffness of the problem by subtracting a least-squares penalty on the gradient jumps across the mesh interfaces from the standard stiffness bilinear form. This penalty bilinear form is similar to the known technique used to stabilize finite element approximations in various contexts. The penalty term is designed to dampen the high frequencies in the spectrum and so it is weighted here by a negative coefficient. The resulting approximation technique is called softFEM since it reduces the stiffness of the problem. The two key advantages of softFEM over the standard Galerkin FEM are to improve the approximation of the eigenvalues in the upper part of the discrete spectrum and to reduce the condition number of the stiffness matrix. We derive a sharp upper bound on the softness parameter weighting the stabilization bilinear form so as to maintain coercivity for the softFEM bilinear form. Then we prove that softFEM delivers the same optimal convergence rates as the standard Galerkin FEM approximation for the eigenvalues and the eigenvectors. We next compare the discrete eigenvalues obtained when using Galerkin FEM and softFEM. Finally, a detailed analysis of linear softFEM for the 1D Laplace eigenvalue problem delivers a sensible choice for the softness parameter. With this choice, the stiffness reduction ratio scales linearly with the polynomial degree. Various numerical experiments illustrate the benefits of using softFEM over Galerkin FEM., Comment: 32 pages

We introduce a boundary penalization technique to improve the spectral approximation of isogeometric analysis (IGA). The technique removes the outliers appearing in the high-frequency region of the approximate spectrum when using the $C^{p-1}, p$-th ($p\ge3$) order isogeometric elements. We focus on the classical Laplacian (Dirichlet) eigenvalue problem in 1D to illustrate the idea and then use the tensor-product structure to generate the stiffness and mass matrices for multiple dimensional problems. To remove the outliers, we penalize the product of the higher-order derivatives from both the solution and test spaces at the domain boundary. Intuitively, we construct a better approximation by weakly imposing features of the exact solution. Effectively, we add terms to the variational formulation at the boundaries with minimal extra computational cost. We then generalize the idea to remove the outliers for the isogeometric analysis to the Neumann eigenvalue problem (for $p\ge2$). The boundary penalization does not change the test and solution spaces. In the limiting case when the penalty goes to infinity, we perform the dispersion analysis of $C^2$ cubic elements for Dirichlet eigenvalue problem and $C^1$ quadratic elements for Neumann eigenvalue problem. We obtain the analytical eigenpairs for the resulting matrix eigenvalue problems. Numerical experiments show optimal convergence rates for the eigenvalues and eigenfunctions of the discrete operator., Comment: 30 pages

It is well-known that the finite difference discretization of the Laplacian eigenvalue problem $-\Delta u = \lambda u$ leads to a matrix eigenvalue problem (EVP) $A x= \lambda x$ where the matrix $A$ is Toeplitz-plus-Hankel. Analytical solutions to tridiagonal matrices with various boundary conditions are given in Strang and MacNamara \cite{strang2014functions}. We generalize the results and develop analytical solutions to the generalized matrix eigenvalue problems (GEVPs) $A x= \lambda Bx$ which arise from the finite element method (FEM) and isogeometric analysis (IGA). The FEM matrices are corner-overlapped block-diagonal while the IGA matrices are almost Toeplitz-plus-Hankel. In fact, IGA with a correction that results in Toeplitz-plus-Hankel matrices gives a better numerical method. In this paper, we focus on finding the analytical eigenpairs to the GEVPs while developing better numerical methods is our motivation. Analytical solutions are also obtained for some polynomial eigenvalue problems (PEVPs). Lastly, we generalize the eigenvector-eigenvalue identity (rediscovered and coined recently for EVPs) for GEVPs and derive some trigonometric identities., Comment: 22 pages

We propose a variational splitting technique for the generalized-$\alpha$ method to solve hyperbolic partial differential equations. We use tensor-product meshes to develop the splitting method, which has a computational cost that grows linearly with respect to the total number of degrees of freedom for multi-dimensional problems. We consider standard $C^0$ finite elements as well as smoother B-splines in isogeometric analysis for the spatial discretization. We also study the spectrum of the amplification matrix to establish the unconditional stability of the method. We then show that the stability behavior affects the overall behavior of the integrator on the entire interval and not only at the limits $0$ and $\infty$. We use various examples to demonstrate the performance of the method and the optimal approximation accuracy. For the numerical tests, we compute the $L_2$ and $H^1$ norms to show the optimal convergence of the discrete method in space and second-order accuracy in time.

In this paper, we introduce the isoGeometric Residual Minimization (iGRM) method. The method solves stationary advection-dominated diffusion problems. We stabilize the method via residual minimization. We discretize the problem using B-spline basis functions. We then seek to minimize the isogeometric residual over a spline space built on a tensor product mesh. We construct the solution over a smooth subspace of the residual. We can specify the solution subspace by reducing the polynomial order, by increasing the continuity, or by a combination of these. The Gramm matrix for the residual minimization method is approximated by a weighted H1 norm, which we can express as Kronecker products, due to the tensor-product structure of the approximations. We use the Gramm matrix as a preconditional which can be applied in a computational cost proportional to the number of degrees of freedom in 2D and 3D. Building on these approximations, we construct an iterative algorithm. We test the residual minimization method on several numerical examples, and we compare it to the Discontinuous Petrov-Galerkin (DPG) and the Streamline Upwind Petrov-Galerkin (SUPG) stabilization methods. The iGRM method delivers similar quality solutions as the DPG method, it uses smaller grids, it does not require breaking of the spaces, but it is limited to tensor-product meshes. The computational cost of the iGRM is higher than for SUPG, but it does not require the determination of problem specific parameters., Comment: 27 pages, 4 figures

The generalized-$\alpha$ time-marching method provides second-order accuracy in time and controls the numerical dissipation in the high-frequency region of the discrete spectrum. This method includes a wide range of time integrators. We increase the order of accuracy of the method while keeping the unconditional stability and the user-control on the high-frequency numerical dissipation. The dissipation is controlled by a single parameter as in the original method. Our high-order schemes require simple modifications of the available implementations of the generalized-$\alpha$ method.

The generalized-$\alpha$ method encompasses a wide range of time integrators. The method possesses high-frequency dissipation while minimizing unwanted low-frequency dissipation and the numerical dissipation can be controlled by the user. The method is unconditionally stable and is of second-order accuracy in time. We extend the second-order generalized-$\alpha$ method to third-order in time while the numerical dissipation can be controlled in a similar fashion. We establish that the third-order method is unconditionally stable. We discuss a possible path to the generalization to higher order schemes. All these high-order schemes can be easily implemented into programs that already contain the second-order generalized-$\alpha$ method.

Most variational forms of isogeometric analysis use highly-continuous basis functions for both trial and test spaces. For a partial differential equation with a smooth solution, isogeometric analysis with highly-continuous basis functions for trial space results in excellent discrete approximations of the solution. However, we observe that high continuity for test spaces is not necessary. In this work, we present a framework which uses highly-continuous B-splines for the trial spaces and basis functions with minimal regularity and possibly lower order polynomials for the test spaces. To realize this goal, we adopt the residual minimization methodology. We pose the problem in a mixed formulation, which results in a system governing both the solution and a Riesz representation of the residual. We present various variational formulations which are variationally-stable and verify their equivalence numerically via numerical tests.

We present a variationally separable splitting technique for the generalized-$\alpha$ method for solving parabolic partial differential equations. We develop a technique for a tensor-product mesh which results in a solver with a linear cost with respect to the total number of degrees of freedom in the system for multi-dimensional problems. We consider finite elements and isogeometric analysis for the spatial discretization. The overall method maintains user-controlled high-frequency dissipation while minimizing unwanted low-frequency dissipation. The method has second-order accuracy in time and optimal rates ($h^{p+1}$ in $L^2$ norm and $h^p$ in $L^2$ norm of $\nabla u$) in space. We present the spectrum analysis on the amplification matrix to establish that the method is unconditionally stable. Various numerical examples illustrate the performance of the overall methodology and show the optimal approximation accuracy.

We study the generalized finite element methods (GFEMs) for the second-order elliptic eigenvalue problem with an interface in 1D. The linear stable generalized finite element methods (SGFEM) were recently developed for the elliptic source problem with interfaces. We first generalize SGFEM to arbitrary order elements and establish the optimal error convergence of the approximate solutions for the elliptic source problem with an interface. We then apply the abstract theory of spectral approximation of compact operators to establish the error estimation for the eigenvalue problem with an interface. The error estimations on eigenpairs strongly depend on the estimation of the discrete solution operator for the source problem. We verify our theoretical findings in various numerical examples including both source and eigenvalue problems., Comment: 22 pages

We approximate the spectra of a class of $2n$-order differential operators using isogeometric analysis in mixed formulations. This class includes a wide range of differential operators such as those arising in elliptic, biharmonic, Cahn-Hilliard, Swift-Hohenberg, and phase-field crystal equations. The spectra of the differential operators are approximated by solving differential eigenvalue problems in mixed formulations, which require auxiliary parameters. The mixed isogeometric formulation when applying classical quadrature rules leads to an eigenvalue error convergence of order $2p$ where $p$ is the order of the underlying B-spline space. We improve this order to be $2p+2$ by applying optimally-blended quadrature rules developed in \cite{20,52} and this order is an optimum in the view of dispersion error. We also compare these results with the mixed finite elements and show numerically that mixed isogeometric analysis leads to significantly better spectral approximations., Comment: 21 pages

We study the spectral approximation of a second-order elliptic differential eigenvalue problem that arises from structural vibration problems using isogeometric analysis. In this paper, we generalize recent work in this direction. We present optimally blended quadrature rules for the isogeometric spectral approximation of a diffusion-reaction operator with both Dirichlet and Neumann boundary conditions. The blended rules improve the accuracy and the robustness of the isogeometric approximation. In particular, the optimal blending rules minimize the dispersion error and lead to two extra orders of super-convergence in the eigenvalue error. Various numerical examples (including the Schr$\ddot{\text{o}}$dinger operator for quantum mechanics) in one and three spatial dimensions demonstrate the performance of the blended rules.

We introduce the dispersion-minimized mass for isogeometric analysis to approximate the structural vibration which we model as a second-order differential eigenvalue problem. The dispersion-minimized mass reduces the eigenvalue error significantly, from the optimum order of $2p$ to the superconvergence order of $2p+2$ for the $p$-th order isogeometric elements with maximum continuity, which in return leads to more robust of the isogeomectric analysis. We first establish the dispersion error for arbitrary polynomial order isogeometric elements. We derive the dispersion-minimized mass in one dimension by solving a $p$-dimensional local matrix problem for the $p$-th order approximation and then extend it to multiple dimensions on tensor-product grids. We show that the dispersion-minimized mass can also be obtained by approximating the mass matrix using optimally blended quadratures. We generalize the results of optimally blended quadratures from polynomial orders $p=1,\cdots, 7$ that were studied in \cite{calo2017dispersion} to arbitrary polynomial order isogeometric approximations. Various numerical examples validate the eigenvalue and eigenfunction error estimates we derive., Comment: 28 pages

We study the approximation of the spectrum of a second-order elliptic differential operator by the Hybrid High-Order (HHO) method. The HHO method is formulated using cell and face unknowns which are polynomials of some degree $k\geq0$. The key idea for the discrete eigenvalue problem is to introduce a discrete operator where the face unknowns have been eliminated. Using the abstract theory of spectral approximation of compact operators in Hilbert spaces, we prove that the eigenvalues converge as $h^{2t}$ and the eigenfunctions as $h^{t}$ in the $H^1$-seminorm, where $h$ is the mesh-size, $t\in [s,k+1]$ depends on the smoothness of the eigenfunctions, and $s>\frac12$ results from the elliptic regularity theory. The convergence rates for smooth eigenfunctions are thus $h^{2k+2}$ for the eigenvalues and $h^{k+1}$ for the eigenfunctions. Our theoretical findings, which improve recent error estimates for Hybridizable Discontinuous Galerkin (HDG) methods, are verified on various numerical examples including smooth and non-smooth eigenfunctions. Moreover, we observe numerically in one dimension for smooth eigenfunctions that the eigenvalues superconverge as $h^{2k+4}$ for a specific value of the stabilization parameter., Comment: 32 pages

Calabro et al. (2017) changed the paradigm of the mass and stiffness computation from the traditional element-wise assembly to a row-wise concept, showing that the latter one offers integration that may be orders of magnitude faster. Considering a B-spline basis function as a non-negative measure, each mass matrix row is integrated by its own quadrature rule with respect to that measure. Each rule is easy to compute as it leads to a linear system of equations, however, the quadrature rules are of the Newton-Cotes type, that is, they require a number of quadrature points that is equal to the dimension of the spline space. In this work, we propose weighted quadrature rules of Gaussian type which require the minimum number of quadrature points while guaranteeing exactness of integration with respect to the weight function. The weighted Gaussian rules arise as solutions of non-linear systems of equations. We derive rules for the mass and stiffness matrices for uniform $C^1$ quadratic and $C^2$ cubic isogeometric discretizations. Our rules further reduce the number of quadrature points by a factor of $(\frac{p+1}{2p+1})^d$ when compared to Calabro et al. (2017), $p$ being the polynomial degree and $d$ the dimension of the problem, and consequently reduce the computational cost of the mass and stiffness matrix assembly by a similar factor.

We study the spectral approximation properties of isogeometric analysis with local continuity reduction of the basis. Such continuity reduction results in a reduction in the interconnection between the degrees of freedom of the mesh, which allows for large savings in computational requirements during the solution of the resulting linear system. The continuity reduction results in extra degrees of freedom that modify the approximation properties of the method. The convergence rate of such refined isogeometric analysis is equivalent to that of the maximum continuity basis. We show how the breaks in continuity and inhomogeneity of the basis lead to artefacts in the frequency spectra, such as stopping bands and outliers, and present a unified description of these effects in finite element method, isogeometric analysis, and refined isogeometric analysis. Accuracy of the refined isogeometric analysis approximations can be improved by using non-standard quadrature rules. In particular, optimal quadrature rules lead to large reductions in the eigenvalue errors and yield two extra orders of convergence similar to classical isogeometric analysis.

We develop quadrature rules for the isogeometric analysis of wave propagation and structural vibrations that minimize the discrete dispersion error of the approximation. The rules are optimal in the sense that they only require two quadrature points per element to minimize the dispersion error [1], and they are equivalent to the optimized blending rules we recently described. Our approach further simplifies the numerical integration: instead of blending two three-point standard quadrature rules, we construct directly a single two-point quadrature rule that reduces the dispersion error to the same order for uniform meshes with periodic boundary conditions. Also, we present a 2.5-point rule for both uniform and non-uniform meshes with arbitrary boundary conditions. Consequently, we reduce the computational cost by using the proposed quadrature rules. Various numerical examples demonstrate the performance of these quadrature rules., Comment: 15pages

We develop and analyze quadrature blending schemes that minimize the dispersion error of isogeometric analysis up to polynomial order seven with maximum continuity in the span ($C^{p-1}$). The schemes yield two extra orders of convergence (superconvergence) on the eigenvalue errors, while the eigenfunction errors are of optimal convergence order. Both dispersion and spectrum analysis are unified in the form of a Taylor expansion for eigenvalue errors. As a consequence, the schemes increase the accuracy and robustness of isogeometric analysis for wave propagation as well as the differential eigenvalue problems. We analyze the methods' robustness and efficacy and utilize numerical examples to verify our analysis of the performance of the proposed schemes., Comment: 26 pages

This paper introduces optimally-blended quadrature rules for isogeometric analysis and analyzes the numerical dispersion of the resulting discretizations. To quantify the approximation errors when we modify the inner products, we generalize the Pythagorean eigenvalue theorem of Strang and Fix. The proposed blended quadrature rules have advantages over alternative integration rules for isogeometric analysis on uniform and non-uniform meshes as well as for different polynomial orders and continuity of the basis. The optimally-blended schemes improve the convergence rate of the method by two orders with respect to the fully-integrated Galerkin method. The proposed technique increases the accuracy and robustness of isogeometric analysis for wave propagation problems.

We propose a simple post-processing technique for linear and high order continuous Galerkin Finite Element Methods (CGFEMs) to obtain locally conservative flux field. The post-processing technique requires solving an auxiliary problem on each element independently which results in solving a linear algebra system whose size is low for any order CGFEM. The post-processing could have been done directly from the finite element solution that results in locally conservative flux on the element. However, the normal flux is not continuous at the elemental boundary. To construct locally conservative flux field whose normal component is also continuous, we propose to do the post-processing on the nodal-centered control volumes which are constructed from the original finite element mesh. We show that the post-processed solution converges in an optimal fashion to the true solution in an H1 semi-norm. We present various numerical examples to demonstrate the performance of the post-processing technique., Comment: 22 pages

A typical two-phase model for subsurface flow couples the Darcy equation for pressure and a transport equation for saturation in a nonlinear manner. In this paper, we study a combined method consisting of continuous Galerkin finite element methods (CGFEMs) followed by a post-processing technique for Darcy equation and finite volume method (FVM) with upwind schemes for the saturation transport equation, in which the coupled nonlinear problem is solved in the framework of operator decomposition. The postprocessing technique is applied to CGFEM solutions to obtain locally conservative fluxes which ensures accuracy and robustness of the FVM solver for the saturation transport equation. We applied both upwind scheme and upwind scheme with slope limiter for FVM on triangular meshes in order to eliminate the non-physical oscillations. Various numerical examples are presented to demonstrate the performance of the overall methodology.

We consider the construction of locally conservative fluxes by means of a simple post-processing technique obtained from the finite element solutions of advection diffusion equations. It is known that a naive calculation of fluxes from these solutions yields non-conservative fluxes. We consider two finite element methods: the usual continuous Galerkin finite element (CGFEM) for solving non dominating advection diffusion equations and the streamline upwind/Petrov-Galerkin (SUPG) for solving advection dominated problems. We then describe the post-processing technique for constructing conservative fluxes from the numerical solutions of the general variational formulation. The post-processing technique requires solving an auxiliary Neumann boundary value problem on each element independently and it produces a locally conservative flux on a vertex centered dual mesh relative to the finite element mesh. We provide a convergence analysis for the post-processing technique. Performance of the technique and the convergence behavior are demonstrated through numerical examples including a set of test problems for advection diffusion equations, advection dominated equations, and drift-diffusion equations., Comment: 20 pages, 10 figures