Galerkin BEM Research Articles

GPU-accelerated solutions to forward problem of TMS

Abstract A forward model characterizes the relationship between source and a measurement in computational form. In TMS, such models are used for planning and targeting the stimulation. The contribution of the head volume conductor to the forward model is typically solved using BEM or FEM. Here we consider BEM models in the context of real-time TMS navigation and the performance benefits of C++ and GPU (Cuda) code over MATLAB. MATLAB is highly efficient in matrix computations but rather slow in other computations, while the GPUs excel in massively parallel problems. We implemented linear Galerkin (LG) BEM solvers in C++ and Cuda languages and benchmarked them against our optimized MATLAB solver [1] in constructing and using a TMS forward model. The test model had realistic gyral structure and contained 21000 surface nodes and 10200 cortical dipole triplets. On a standard 2019 PC and entry-level GPU (Nvidia GTX 1060), we built the BEM model in 75 seconds and solved full surface potentials of all cortical dipole triplets in 15 seconds using C++/Cuda, compared to over 20 minutes and 104 seconds in MATLAB (not using GPU). When this model was used in TMS simulation with a 42-dipole coil, the electric field was solved for 49 coil positions per second (cps); GPU-accelerated MATLAB and C++ were as fast. With a 15000-dipole coil model, the C++/Cuda operated at the speed of 44 cps, while the MATLAB code was over 50 times slower (< 0.1 cps). In TMS navigation, the use of GPU and C++/Cuda allows to build the whole 4-compartment model and defining regions of interest during the stimulation session, instead of building the model offline and saving and loading 1.8-2.5 GB model files. Further, the use of C++/Cuda enables real-time performance with practically any coil model. [1] Stenroos and Koponen, NeuroImage 2019 Research Category and Technology and Methods Basic Research: 10. Transcranial Magnetic Stimulation (TMS) Keywords: TMS navigation, forward model, GPU, real-time

Directional $\mathcal{H}^2$ Compression Algorithm: Optimisations and Application to a Discontinuous Galerkin BEM for the Helmholtz Equation

Directional $\mathcal{H}^2$ Compression Algorithm: Optimisations and Application to a Discontinuous Galerkin BEM for the Helmholtz Equation

Adaptive quadrature rules for Galerkin BEM

The singular integrals in the Galerkin Boundary Element Method are usually treated using singularity removing transformations. The regularized integrals are approximated by tensor product Gaussian quadrature rules. Although these integrals are smooth, the convergence rate depends strongly on the type of singularity and the aspect ratio of the triangulation. Here an adaptive quadrature scheme is presented that ensures a convergence at a user-specified rate. The effectiveness of the resulting quadrature scheme is compared with the non-adaptive approach.

Integral Evaluation for a Closed-Form 2-D Potential Formulation of the Galerkin BEM

Closed-form solutions are derived for the regular and adjacent-singular integrals involving the two-dimensional Laplacian’s Green’s function and its normal derivative that arise in the Galerkin BEM. Motivation for their use is provided by comparing the accuracy and time required to compute the BEM system matrix relative to traditional numerical integration. Specifically, it is shown that the use of the closed-form solutions reduces computation time to produce the influence matrices by approximately 95% relative to Gauss quadrature with similar accuracy. All elements of these matrices are then expressed in a form convenient for implementation.

Hp non-conforming a priori error analysis of an Interior Penalty Discontinuous Galerkin BEM for the Helmholtz equation

This work is concerned with the construction and the hp non-conforming a priori error analysis of a Discontinuous Galerkin DG numerical scheme applied to the hypersingular integral equation related to the Helmholtz problem in 3D. The main results of this article are an error bound in a norm suited to the problem and in the L2-norm. Those bounds are quasi-optimal for the h-convergence and the p-convergence. Various formulation choices and penalty functions are theoretically discussed. In particular we show that a penalty function of the shape h2p leads to a quasi-optimal convergence of the scheme. Some numerical experiments confirm the expected rates of convergence and the effect of the penalty function.

Adaptive BEM for elliptic PDE systems, part I: abstract framework, for weakly-singular integral equations

In the present work, we consider weakly-singular integral equations arising from linear second-order elliptic PDE systems with constant coefficients, including, e.g. linear elasticity. We introduce a general framework for optimal convergence of adaptive Galerkin BEM. We identify certain abstract conditions for the underlying meshes, the corresponding mesh-refinement strategy, and the ansatz spaces that guarantee that the weighted-residual error estimator is reliable and converges at optimal algebraic rate if used within an adaptive algorithm. These conditions are satisfied, e.g. for discontinuous piecewise polynomials on simplicial meshes as well as certain ansatz spaces used for isogeometric analysis. Technical contributions include the localization of (non-local) fractional Sobolev norms and local inverse estimates for the (non-local) boundary integral operators associated to the PDE system.

Modeling multicrack propagation by the fast multipole symmetric Galerkin BEM

The Fast Multipole Method coupled with the Symmetric Galerkin BEM is employed in this work to simulate fatigue crack growth. The resulted crack propagation code is accelerated with a fast matrix update, a parallel implementation and a sparse matrix format. By using multiple nodes, this code accommodates also multiple surface-breaking cracks. The numerical tests presented herein allow the propagation of multiple cracks in single or multilayer domains.

An adaptive IgA‐BEM with hierarchical B‐splines based on quasi‐interpolation quadrature schemes

SummaryThe isogeometric formulation of the boundary element method (IgA‐BEM) is investigated within the adaptivity framework. Suitable weighted quadrature rules to evaluate integrals appearing in the Galerkin BEM formulation of 2D Laplace model problems are introduced. The proposed quadrature schemes are based on a spline quasi‐interpolation (QI) operator and properly framed in the hierarchical setting. The local nature of the QI perfectly fits with hierarchical spline constructions and leads to an efficient and accurate numerical scheme. An automatic adaptive refinement strategy is driven by a residual‐based error estimator. Numerical examples show that the optimal convergence rate of the Galerkin solution is recovered by the proposed adaptive method.

Quadrature for parabolic Galerkin BEM with moving surfaces

The successful implementation of the Galerkin Boundary Element Method hinges on the accurate and effective quadrature of the influence coefficients. For parabolic boundary integral operators quadrature must be performed in space and time where integrals have singularities when source- and evaluation points coincide. For problems where the surface is fixed, the time integration can be performed analytically, but for moving geometries numerical quadrature in space and time must be used. For this case a set of transformations is derived that render the singular space–time integrals into smooth integrals that can be treated with standard tensor product Gauss quadrature rules. This methodology can be applied to the heat equation and to transient Stokes flow.

The Wave-Matching Boundary Integral Equation — An energy approach to Galerkin BEM for acoustic wave propagation problems

In this paper, a new Boundary Integral Equation (BIE) is proposed for solution of the scalar Helmholtz equation. Applications include acoustic scattering problems, as occur in room acoustics and outdoor and underwater sound propagation. It draws together ideas from the study of time-harmonic and transient BIEs and spatial audio sensing and rendering, to produce an energy-inspired Galerkin BEM that is intended for use with oscillatory basis functions. Pivotal is the idea that waves at a boundary may be decomposed into incoming and outgoing components. When written in its admittance form, it can be thought of setting the Burton–Miller coupling parameter differently for each basis function based on its oscillation; this is a discrete form of the Dirichlet-to-Neumann map. It is also naturally expressed in a reflectance form, which can be solved by matrix inversion or by marching on in reflection order. Consideration of this leads to an orthogonality relation between the incoming and outgoing waves, which makes the scheme immune to interior cavity eigenmodes. Moreover, the scheme is seen to have two remarkable properties when solution is performed over an entire obstacle: (i) it has a condition number of 1 for all positive-real wavenumber k on any closed geometry when a specific choice of cylindrical basis functions are used; (ii) when modelling two domains separated by a barrier domain, the two problems are numerical uncoupled when plane wave basis functions are used — this is the case in reality but is not achieved by any other BIE representation that the authors are aware of. Normalisation and envelope functions, as would be required to build a Partition-of-Unity or Hybrid-Numerical-Asymptotic scheme, are introduced and the above properties are seen to become approximate. The modified scheme is applied successfully to a cylinder test case: accuracy of the solution is maintained and the BIE is still immune to interior cavity eigenmodes, gives similar conditioning to the Burton–Miller method and iterative solution is stable. It is seen that for this test case the majority of values in the interaction matrices are extremely small and may be set to zero without affecting conditioning or accuracy, thus the linear system become sparse - a property uncommon in BEM formulations.

A Galerkin BEM for high-frequency scattering problems based on frequency-dependent changes of variables

In this paper we develop a class of efficient Galerkin boundary element methods for the solution of two-dimensional exterior single-scattering problems. Our approach is based upon construction of Galerkin approximation spaces confined to the asymptotic behaviour of the solution through a certain direct sum of appropriate function spaces weighted by the oscillations in the incident field of radiation. Specifically, the function spaces in the illuminated/shadow regions and the shadow boundaries are simply algebraic polynomials whereas those in the transition regions are generated utilizing novel, yet simple, \emph{frequency dependent changes of variables perfectly matched with the boundary layers of the amplitude} in these regions. While, on the one hand, we rigorously verify for smooth convex obstacles that these methods require only an $\mathcal{O}\left( k^{\epsilon} \right)$ increase in the number of degrees of freedom to maintain any given accuracy independent of frequency, and on the other hand, remaining in the realm of smooth obstacles they are applicable in more general single-scattering configurations. The most distinctive property of our algorithms is their \emph{remarkable success} in approximating the solution in the shadow region when compared with the algorithms available in the literature.

Non-polynomial spline alternatives in Isogeometric Symmetric Galerkin BEM

The application of the Isogeometric Analysis (IgA) paradigm to Symmetric Galerkin Boundary Element Method (SGBEM) is investigated. In order to obtain a very flexible approach, the study is here developed by using non-polynomial spline functions to represent both the domain boundary and the approximate solution. The numerical comparison between IGA-SGBEM and both curvilinear and standard SGBEM approaches shows the general capability of the presented method to produce accurate approximate solutions with less degrees of freedom.

Symmetric Galerkin BEM for Non Linear Analysis of Historical Masonries

The preservation of the historical and monumental buildings, but also of the considerable heritage of old constructions made by traditional techniques, is one of the actual problems of the structural mechanics. The level of knowledge of their structural behavior in presence of external actions is made through calculus methods and simple procedures in order to allow a reading of the material suffering degree and as a consequence of the related safety. Unfortunately, often the masonry panels show openings located in an irregular way and cracks having small or big dimensions. In these cases the employment of strategies, as for example the identifying of the masonry piers and the transformation of the wall as a frame, prove to be inapplicable. Recently, a new method was introduced, called Boundary Element Method, which applied in its symmetric formulation, in virtue of some peculiarities, proves to be better usable in comparison with other analysis methods. In this paper an elastic analysis of walls, also in presence of geometrical nonlinearity consisting in the contact/detachment phenomenon among stone blocks. The wall having any shape and zone-wise variable physical characteristics is loaded in its plane. For these structures some interventions of structural strengthening have as aim to improve the wall behavior by reducing the stress concentration, so to have a better safety in comparison with its initial value.

Isogemetric analysis and symmetric Galerkin BEM: A 2D numerical study

Isogeometric approach applied to Boundary Element Methods is an emerging research area (see e.g. Simpson et al. (2012) [33]). In this context, the aim of the present contribution is that of investigating, from a numerical point of view, the Symmetric Galerkin Boundary Element Method (SGBEM) devoted to the solution of 2D boundary value problems for the Laplace equation, where the boundary and the unknowns on it are both represented by B-splines (de Boor (2001) [9]). We mainly compare this approach, which we call IGA-SGBEM, with a curvilinear SGBEM (Aimi et al. (1999) [2]), which operates on any boundary given by explicit parametric representation and where the approximate solution is obtained using Lagrangian basis. Both techniques are further compared with a standard (conventional) SGBEM approach (Aimi et al. (1997) [1]), where the boundary of the assigned problem is approximated by linear elements and the numerical solution is expressed in terms of Lagrangian basis. Several examples will be presented and discussed, underlying benefits and drawbacks of all the above-mentioned approaches.

Computation of Effective Material Properties in Smart Composite Materials by a Symmetric Galerkin BEM

In this paper, the symmetric Galerkin boundary element method (SGBEM) will be developed and applied for boundary value problems with layered and fiber reinforced piezoelectric representative volume elements (RVE) and real macroscopic structures. Mechanical and electric loadings are considered to determine the effective material properties. For this purpose, the resulting boundary value problem is formulated as boundary integral equations (BIEs). The Galerkin method is applied for the spatial discretization of the boundary to solve the BIEs numerically. The required surface derivatives of the generalized displacements are computed directly with a boundary integral equation. Numerical examples will be presented and discussed to show the efficiency of the present SGBEM and the influence of the fiber variation on the effective material properties.

216 2次元波動問題に対する演算子積分Galerkin境界要素法・有限要素法結合解法

216 2次元波動問題に対する演算子積分Galerkin境界要素法・有限要素法結合解法

Sparse tensor product spectral Galerkin BEM for elliptic problems with random input data on a spheroid

We introduce and analyze a sparse tensor product spectral Galerkin Boundary Element Method based on spherical harmonics for elliptic problems with random input data on a spheroid. Problems of this type appear in geophysical applications, in particular in data acquisition by satellites. Aiming at a deterministic computation of the k-th order statistical moments of the random solution, we establish convergence theorems showing that the sparse tensor product spectral Galerkin discretization is superior to the full tensor product spectral Galerkin discretization in the case of mixed regularity of the data's k-th order moments, naturally implying mixed regularity of the k-th order moments of the random solution. We prove that analytic regularity of the data's k-th order moments implies analytic regularity of the solution's k-th order moments. We illustrate performance of the sparse and full tensor product discretization schemes on several numerical examples.

Simple error estimators for the Galerkin BEM for some hypersingular integral equation in 2D

A posteriori error estimation is an important tool for reliable and efficient Galerkin boundary element computations. For hypersingular integral equations in 2D with a positive-order Sobolev space, we analyse the mathematical relation between the (h − h/2)-error estimator from [S. Ferraz-Leite and D. Praetorius, Simple a posteriori error estimators for the h-version of the boundary element method, Computing 83 (2008), pp. 135–162], the two-level error estimator from [M. Maischak, P. Mund, and E. Stephan, Adaptive multilevel BEM for acoustic scattering, 585 Comput. Methods Appl. Mech. Eng. 150 (1997), pp. 351–367], and the averaging error estimator from [C. Carstensen and D. Praetorius, Averaging techniques for the a posteriori bem error control for a hypersingular integral equation in two dimensions, SIAM J. Sci. Comput. 29 (2007), pp. 782–810]. All of these a posteriori error estimators are simple in the following sense: first, the numerical analysis can be done within the same mathematical framework, namely localization techniques for the energy norm. Second, there is almost no implementational overhead for the realization.

Numerical solution of exterior Maxwell problems by Galerkin BEM and Runge–Kutta convolution quadrature

In this paper we consider time-dependent electromagnetic scattering problems from conducting objects. We discretize the time-domain electric field integral equation using Runge–Kutta convolution quadrature in time and a Galerkin method in space. We analyze the involved operators in the Laplace domain and obtain convergence results for the fully discrete scheme. Numerical experiments indicate the sharpness of the theoretical estimates.

The nsBETI method: an extension of the FETI method to non‐symmetrical BEM‐FEM coupled problems

SUMMARYThe finite element tearing and interconnecting (FETI) method is recognized as an effective domain decomposition tool to achieve scalability in the solution of partitioned second‐order elasticity problems. In the boundary element tearing and interconnecting (BETI) method, a direct extension of the FETI algorithm to the BEM, the symmetric Galerkin BEM formulation, is used to obtain symmetric system matrices, making possible to apply the same FETI conjugate gradient solver. In this work, we propose a new BETI variant labeled nsBETI that allows to couple substructures modeled with the FEM and/or non‐symmetrical BEM formulations. The method connects non‐matching BEM and FEM subdomains using localized Lagrange multipliers and solves the associated non‐symmetrical flexibility equations with a Bi‐CGstab iterative algorithm. Scalability issues of nsBETI in BEM–BEM and combined BEM–FEM coupled problems are also investigated. Copyright © 2012 John Wiley &amp; Sons, Ltd.