INdAM Workshop DREAMS

Oleg Davydov (University of Giessen, Germany)

Rida T. Farouki (University of California Davis, CA)

Bert Jüttler (Johannes Kepler University Linz, Austria)

Tom Lyche (University of Oslo, Norway)

Giancarlo Sangalli (University of Pavia, Italy)

Stefano Serra-Capizzano (University of Insubria, Italy)

Keynote presentations

Meshless finite difference methods
Oleg Davydov (University of Giessen, Germany)

Meshless finite difference methods for partial differential equations apply the methodology of the Finite Difference Method in the grid-free setting by using numerical differentiation formulas on scattered nodes. These formulas can be obtained by requiring polynomial consistency or via optimal recovery of differential operators with the help of kernel (radial basis) interpolation. Since no mesh has to be imposed on the nodes, they can be freely distributed following the exact geometry of the model and/or the features of the solution. After introducing the method I will illustrate its performance in a number of numerical experiments, in particular for time-dependent equations on evolving-in-time manifolds.
Joint work with Dang Thi Oanh, Hoang Xuan Phu, Robert Schaback, Andriy Sokolov, Ngo Manh Tuong, and Stefan Turek.

Precision CNC machining applications based on Pythagorean-hodograph curves
Rida T. Farouki (University of California Davis, CA)

Two novel applications for the Pythagorean-hodograph (PH) curves in CNC machining are presented: precision machining of rational swept surface forms, and high-speed cornering under specified tolerance and acceleration bounds. The range of possible rational swept surface forms can be greatly expanded to encompass transformations of a profile curve that depend on differential and integral properties of a PH sweep curve, and an algorithm has been developed to machine such surfaces directly from their high-level procedural definitions. In the high-speed cornering application, the sharp corners of a piecewise-linear path are rounded by G^2 PH quintic corner curves that satisfy a prescribed tolerance, and a family of feedrate functions for executing such curves with bounded acceleration is developed. Experimental implementation results for both applications are presented.

Algorithms for circular arc polygons: Medial axis computation, straight skeletons and generalized star-shaped domains
Bert Jüttler (Johannes Kepler University Linz, Austria)

Circular arcs are well known as highly useful geometric primitives for geometry processing, since they combine geometric flexibility and high approximation power with the simplicity of performing geometric computations. More precisely, circular arc spline curves are known to possess cubic approximation power, while the computation of intersections requires solely square roots, and bisector curves are quadratic implicit curves (conic sections). Based on these observations, an algorithm for the computation of medial axes and trimmed offsets of planar free-form shapes has been established by Aichholzer et al. [1]. We recall these results and describe two recent extensions of this approach. First we introduce straight skeletons of circular arc polygons. These objects are based on a procedural definition of offsets that generates the minimum number of segments and provides a local invertibility property. Second we generalize the notion of star-shaped domain, based on the visibility with respect to circular arcs.
Joint work with F. Aurenhammer, M. Haberleitner, M.-S. Kim, S. Maroscheck, and B. Weiss.

[1] O. Aichholzer, W. Aigner, F. Aurenhammer, T. Hackl, B. Jüttler, and M. Rablb. Medial axis computation for planar free-form shapes. Comput. Aided Des. 41 (2009), 339-349.

Tchebycheffian Splines and Tchebycheffian B-splines
Tom Lyche (University of Oslo, Norway)

We give a survey of Tchebycheffian Splines and Tchebycheffian B-splines. These piecewise functions are natural generalizations of polynomial splines and polynomial B-splines. Each segment of the spline belongs to a piecewise Extended Complete Tchebycheff (ECT)-space. An ECT-space can also be considered as the null space of a linear differential operator with variable coefficients defined on a suitable interval. ECT-spaces can either be spanned by a set of basis functions that are a natural generalization of the polynomial power basis or an alternative basis with similar properties to the Bernstein basis for polynomials. We define Tchebycheffian B-splines recursively, allowing different spaces on the segments. If time allows we also consider Tchebycheffian divided differences. Several examples will be given to illustrate the concepts.

How fast is the isogeometric k-method?
Giancarlo Sangalli (University of Pavia, Italy)

In this work we show the superiority, in terms of computational efficiency, of the isogeometric k-method with respect to low-degree isogeometric discretizations. With an innovative implementation, increasing the spline degree and regularity (i.e., k-method) significantly improves not only the accuracy, which is known, but also the accuracy to computational-time ratio. On the contrary, when the k-method is implemented as a classical finite element method, increasing the degree becomes soon prohibitive and, in fact, quadratic approximation is the most efficient. The new implementation of the isogeometric method combines three ingredients: the matrix-free approach (standard in high-order methods), weighted-quadrature (an ad-hoc strategy to compute the integrals of the Galerkin system) and a preconditioner based on the Fast Diagonalization method (an old idea to solve the Sylvester matrix equation).
Joint work with Mattia Tani.

The GLT class as a generalized Fourier analysis and applications
Stefano Serra-Capizzano (University of Insubria, Italy)

The class of Generalized Locally Toeplitz (GLT) sequences has been introduced in [1,2] as a generalization of both classical Toeplitz sequences and variable coefficient differential operators. For every sequence of the class, it is possible to give a rigorous description of the asymptotic spectrum in terms of a function (the symbol) that can be easily identified. The GLT class has nice algebraic properties; it has been proven that it is stable under linear combinations, products, and inversion when the sequence which is inverted shows a sparsely vanishing symbol (sparsely vanishing symbol = a symbol which vanishes at most in a set of zero Lebesgue measure). Furthermore, the GLT class virtually includes any approximation of partial differential equations (PDEs) by local methods (Finite Difference, Finite Element, Isogeometric Analysis etc) and, based on this, we demonstrate that our results on GLT sequences can be used in a PDE setting in various directions, including preconditioning, multigrid, spectral detection of branches, fast `matrix-less’ computation of eigenvalues, stability issues. We will discuss specifically the impact and the further potential of the theory with special attention to the IgA setting.

[1] S. Serra-Capizzano. Generalized Locally Toeplitz sequences: spectral analysis and applications to discretized partial differential equations. Linear Algebra Appl. 366 (2003), 371-402.
[2] S. Serra-Capizzano. The GLT class as a generalized Fourier analysis and applications. Linear Algebra Appl. 419 (2006), 180-233.

Invited presentations

An IgA approach to energetic BEM: Preliminary results
Alessandra Aimi (University of Parma, Italy)

The energetic Boundary Element Method (BEM) is a discretization technique for the numerical solution of wave propagation problems, introduced in [1] and applied in the last decade to wave propagation inside bounded domains or outside bounded obstacles, taking into account viscous and material damping too. The differential initial-boundary value problem at hand is converted into a space-time Boundary Integral Equation (BIE), which is then written in weak form through energy considerations and discretized by a Galerkin approach. Taking into account the model problem of 2D soft scattering of undamped waves by open arcs, the aim of the talk is to explore the introduction of the IgA approach into energetic BEM, for what concerns space discretization, in order to take the same benefits already observed in IgA-SGBEM applied to BIEs related to elliptic problems [2]. Numerical challenges to be faced for an efficient integration of singular kernels related to the fundamental solution of the wave operator will be outlined and first numerical results will be given and discussed.

[1] A. Aimi and M. Diligenti. A new space-time energetic formulation for wave propagation analysis in layered media by BEMs. Int. J. Numer. Meth. Engrg. 75 (2008), 1102-1132.
[2] A. Aimi, M. Diligenti, M.L. Sampoli, and A. Sestini. Isogeometric analysis and symmetric Galerkin BEM: a 2D numerical study. Appl. Math. Comput. 272 (2016), 173-186.

Scattered data fitting with THB-spline for industrial applications
Cesare Bracco (University of Florence, Italy)

We present an algorithm for the reconstruction of surfaces from scattered data based on THB-splines. The method is obtained by combining the structure of hierarchical quasi-interpolation with local polynomial approximations. The key idea is to exploit both the adaptivity features of the THB-splines and local polynomial approximations of variable degree. This approach leads to an algorithm which naturally adapts to the shape and density of the data, allowing to get competitive approximation performances and to avoid artifacts in the final results. We analyze the behaviour of the method on a wide selection of scattered data of high complexity, including point clouds generated in industrial applications.

A symmetry preserving fast assembling strategy for IGA
Andrea Bressan (University of Oslo, Norway)

The IsoGeometric Analysis (IGA) method was proposed as a link between analysis and design. It consists in using an isoparametric Galerkin discretisation based on spline spaces. Due to the smoothness of the basis functions, IGA is a promising setting for high order methods. Increasing the order requires new assembling strategies because of the computational cost. The standard element by element assembling has a cost of O(p^(3d) N ) where p is the polynomial degree, d is the domain dimension and N is the space dimension. The best known IGA strategy has a computational cost O(p^(d+1)), but due to the approximation of the involved integrals it breaks the symmetry of system matrix (without affecting the a priori error estimates). The proposed method costs O(p^(d+2) N), but it preserves the matrix symmetry providing an interesting tradeoff for symmetric problems.

A change of paradigm in the Galerkin isogeometric method formation and assembly: Row assembly, sum factorization and weighted quadrature
Francesco Calabrò (University of Cassino and Southern Lazio, Italy)

In this talk we present a new formation and assembly strategy which relies on two key ideas: assembly row by row, instead of element by element and an efficient formation strategy based on weighted quadrature and sum factorization that is applied to each specific row of the matrix. We detail the major conceptual change of paradigm with respect to the standard implementation: the idea of using weighted quadrature. We see that if the test function is incorporated in the integration weight while the trial function, the geometry parametrization and the PDEs coefficients are taken as the integrand function, then this approach is very effective in reducing the computational cost. Finally, we present some benchmark tests and show that the optimal order of approximation of the method and the accuracy of full gauss quadrature is maintained while the computational burden of forming the matrix equations is significantly reduced.
Joint work with G. Sangalli, M. Tani, R. Hiemstra, T.J.R. Hughes, P. Antolin, A. Buffa, and M. Martinelli.

Hierarchical quadrature schemes for adaptive IgA-SGBEM
Antonella Falini (INdAM c/o University of Florence, Italy)

The isogeometric formulation of the Symmetric Galerkin Boundary Element Method (SGBEM) is included in the adaptivity framework introducing hierarchical B-spline basis functions. In particular, dealing with integral equations, specific weighted quadrature rules will be suitably tailored for the use of the hierarchical basis. In fact, due to the singular nature of the integral kernels, ad hoc quadrature schemes need to be introduced. To this end, a new method based on a spline quasi-interpolant (QI) operator is investigated. The local nature of the QI perfectly fits within the hierarchical structure and allows to keep low the computational costs. The applicability of the proposed approach will be supported by some final examples.

Fast computation of Toeplitz eigenvalues through asymptotic expansions and extrapolation
Carlo Garoni (INdAM Marie-Curie fellow)

Extrapolation is known to be one of the most successful ways for accelerating the convergence of numerical methods [1]. In the words of Birkhoff and Rota, “its usefulness for practical computations can hardly be overestimated”. In the presence of an asymptotic expansion for the quantity to be computed/approximated, a “canonical'” extrapolation method arises; think, for example, to Romberg’s integration method, which arises from the Euler-Maclaurin expansion associated with the trapezoidal formula. In this presentation, we discuss a recently conjectured asymptotic expansion for the eigenvalues of banded symmetric Toeplitz matrices. We also describe the related extrapolation method, which allows the fast computation of the spectrum of such matrices.
Joint work with Sven-Erik Ekstrom and Stefano Serra-Capizzano.

[1] C. Brezinski and M. Redivo-Zaglia. Extrapolation Methods: Theory and Practice. North-Holland, 1991.

Adaptive spline technologies for aircraft engine design
David Großmann (MTU Aero Engines AG, Munich, Germany)

Industrial products are usually designed within Computer Aided Engineering (CAE) systems based on the B-spline technology and its non-uniform rational extension (NURBS). To overcome the limitations of their tensor-product structure, we invested in the industrial integration of recently developed generalizations: The truncated hierarchical B-splines (THB-splines) and the patchwork B-splines (PB-splines). The talk will give an overview about the aircraft engine design process at MTU Aero Engines AG and its use of the new spline technologies. First, for an adaptive surface fitting framework to reconstruct CAD models from (optical) measured point data where they lead to significant improvements with respect to the quality of the resulting geometric shape compared to the existing tensor-product spline technology. This enables us to transfer automatically the shapes of manufactured and operated parts back into the CAE systems. Second, for the exact lofting of B-spline curves by using the highly efficient structures of patchwork B-splines in the definition of the blade geometries. Third, we use the adaptive technology for the simulation-based deformation of CAD models which solves the the fundamental engine design problem of optimizing the parts for performing in hot working conditions while a cold, unloaded CAD model is required for manufacturing.

Hierarchical box splines and the immersed boundary method in isogeometric analysis
Tadej Kanduč (INdAM c/o University of Florence, Italy)

One of the important features of isogeometric analysis in numerical simulations is an accurate representation of computational domains. Box splines together with the immersed boundary method allow us to develop more complex single patch geometries than with a classical approach using tensor-product B-splines. Box splines are defined on regular triangulations. Its mesh structure has an important computational advantage over splines on general triangulations. Adaptivity is achieved by considering hierarchically nested sequence of box splines spaces. The derived simulation method is free of user defined penalties and stabilization parameters, and the final system of equations is symmetric for symmetric problems. A small number of additional degrees of freedom (DOF) is introduced along the domain boundary to impose weak boundary conditions. The amount of DOF can be additional reduced if we employ truncated variants of box splines. Several numerical examples demonstrate the optimal convergence of the adaptive scheme.

Rigid-body motion interpolation using cubic PH biarcs
Marjeta Knez (University of Ljubljana, Slovenia)

Polynomial Pythagorean-hodograph (PH) curves in space, which are characterized by the property that the unit tangent is rational, have many important features for practical applications. One of them is that these curves can be equipped with rational orthonormal frames called Euler-Rodriques (ER) frames, where the first frame vector coincides with the unit tangent. The second important property is that the arc-length function is a polynomial. Joining these two properties we can construct motions of a rigid-body that interpolate some given positions and have a prescribed length of the center trajectory. In the talk an interpolation scheme for G^1 Hermite motion data, i.e., interpolation of data points and rotations at the points, with cubic PH biarcs is presented, where the rotational part of the motion is determined by the ER frame. Further, the length of the biarc is interpolated too. It is shown that the solution exists for any data and any length greater than the difference between the interpolation points. Moreover, the interpolant depends on some free parameters, which can be chosen so that the center trajectory is of a nice shape and the rotation of vectors in a normal plane around the tangent is minimized. The derived theoretical results are illustrated with numerical examples.

Gaussian quadrature for C^1 cubic Clough-Tocher macro-triangles
Jirí Kosinka (University of Groningen, The Netherlands)

A numerical integration rule for multivariate cubic polynomials over n-dimensional simplices was designed by Hammer and Stroud in 1956. The quadrature rule requires n+2 quadrature points: the barycentre of the simplex and n+1 points that lie on the connecting lines between the barycentre and the vertices of the simplex. In the planar case, this particular rule belongs to a two-parameter family of quadrature rules that admit exact integration of bivariate polynomials of total degree three over triangles. We prove that this rule is exact for a larger space, namely the C^1 cubic Clough-Tocher spline space over macro-triangles if and only if the split-point is the barycentre.

Recognizing algebraic surfaces from special classes
Miroslav Lávička (University of West Bohemia, Plzen, Czech Republic)

In computer aided geometric design, basic modelling surfaces, with the property being simple and widely used, are applied to construct complex models. Typical examples are ruled surfaces, rotational surfaces, canal surfaces, swept surfaces, translation surfaces, etc. Recognition of these surfaces from their equations, investigation of suitable parameterization methods and other related topics became a frequent research area in the past. One of the challenging tasks is to determine the type of a given algebraic surface which is described only by its implicit representation (or by a parameterization not reflecting directly the type). In this contribution we focus on surfaces of revolution and canal surfaces, which are often used in geometric modelling, computer-aided design and technical practice (e.g. as blending surfaces smoothly joining two parts with circular ends). Our goal is to formulate simple and efficient algorithms whose input is a polynomial with the coefficients from some subfield of R and the output is the answer whether the surface is a rotational surface or a rational canal surface. In the affirmative case we also compute a rational parameterization.

Low-rank approximation for isogeometric analysis
Angelos Mantzaflaris (Johannes Kepler University Linz, Austria)

Low-rank approximation techniques have proven to be the solution of choice for reducing the complexity and restoring computational tractability in areas such as optimization, sensitivity analysis, inverse problems, biology and chemistry. In essence, these methods provide ways to recover structure in the parameters or carefully select a few, so that the evaluation of the, so called, reduced model, on large data grids becomes tractable. This is done by a suitable projection of the function on a lower-dimensional manifold of tensors, whose dimension is called the rank of the tensor. Different notions of ranks and the corresponding low-rank approximation formats have been introduced, having different approximation and computational complexity properties. In this talk we discuss applications of this approach to isogeometric analysis.

Isogeometric analysis for a 2D and 3D MHD subproblem: Spectral symbol and fast iterative solvers
Mariarosa Mazza (Max Planck Institute for Plasma Physics, Munich, Germany)

In this talk, we focus on large and highly ill-conditioned linear systems arising from a B-spline discretization of a parameter-dependent magnetohydrodynamics (MHD) subproblem encountered when macroscopically describing the behavior of the plasma. We show that the sequence of the involved coefficient matrices belongs to the Generalized Locally Toeplitz class and we compute the spectral symbol describing its asymptotic eigenvalue distribution. The study of the symbol and of its eigenvalue functions shows that the coefficient matrices are affected by various sources of ill-conditioning related to both physical and approximation parameters. We use then the retrieved spectral information to design a strategy made up of spectrally complementary iterative solvers able to satisfactory deal with the sources of ill-conditioning. In detail, we use a 2D and 3D vector extension of the multi-iterative approach already proposed in the literature for the scalar Laplacian operator as preconditioner for a Krylov-type method. Several numerical examples shows that the resulting solver is computationally attractive and robust with respect to the relevant parameters.

C^2 continuous time dependent feedrate scheduling with configurable kinematic constraints
Duccio Mugnaini (University of Insubria, Italy)

We present a configurable trajectory planning strategy on planar paths for offline definition of time-dependent C^2 piecewise quintic feedrates. The more conservative formulation ensures chord tolerance, as well as prescribed bounds on velocity, acceleration and jerk Cartesian components. Since the less restrictive formulations of our strategy can usually still ensure all the desired bounds while simultaneously producing faster motions, the configurability feature is useful not only when reduced motion control is desired but also when full kinematic control has to be guaranteed. Our approach can be applied to any planar path with a piecewise sufficiently smooth parametric representation. When Pythagorean-hodograph spline curves are considered, the corresponding accurate and efficient CNC interpolator algorithms can be exploited.

Surface patches with Pythagorean-hodograph isoparametric curves
Francesca Pelosi (University of Rome “Tor Vergata”, Italy)

Although the isoparametric curves have no intrinsic geometrical significance, they are nevertheless useful in practical applications, such as path planning for the machining or inspection of surfaces. The Pythagorean-hodograph (PH) curves have a distinct advantage over ordinary polynomial curves in this respect, since their arc lengths are simply polynomial functions of the curve parameter. In the present work we aim to investigate the feasibility of constructing surface patches with Pythagorean-hodograph isoparametric curves. The simplest non-trivial instances, interpolating four prescribed patch boundary curves, involve degree (5,4) tensor-product surface patches s(u,v), whose v = constant isoparametric curves are all spatial PH quintics. The construction can be reduced to solving a novel type of quadratic quaternion equation and a closed-form solution can be derived. Conditions for the existence of solutions are identified and the residual scalar freedoms are exploited to improve the interior shape of the patch. A selection of computed examples illustrate the performance of the method.
Joint work with Rida T. Farouki, Maria Lucia Sampoli, and Alessandra Sestini.

Advanced modeling and applications of isogeometric shells: From composites to fluid-structure interaction
Alessandro Reali (University of Pavia, Italy – Technical University of Munich, Germany)

This work deals with some recent advances on modeling and applications of shell structures allowed by the unique features of isogeometric analysis (IGA). In particular, we herein focus on three interesting problems. The first one is related to an inexpensive modeling strategy for composite structures. The proposed approach consists of an isogeometric discretization comprising a single element through the thickness and a post-processing technique able to recover an accurate out-of-plane stress state by direct integration of the equilibrium equations in strong form. We then present two novel isogeometric approaches for specific fluid-structure interaction (FSI) problems. The first one exploits a boundary integral formulation of Stokes equations to model the surrounding flow and a nonlinear Kirchhoff-Love shell theory to model the elastic behavior of the structure. The proposed method seems to be particularly attractive for the simulation of falling objects, since only the boundary representation (B-Rep) of the thin structure middle surface is needed to describe the entire studied problem. Finally, the goal of the last considered application is to show the high flexibility and potential of so-called “immersogeometric” methods to study complex problems as those typically found in Biomechanics. Within this context, FSI simulations of patient-specific aortic valve designs are successfully carried out and also compared with medical images.

New quadrature rules for IgA-BEM applications
Maria Lucia Sampoli (University of Siena, Italy)

Boundary element methods (BEMs) can be considered in some cases a valid alternative to classical domain methods, such as finite differences (FDMs) and finite elements (FEMs). Through the fundamental solution the differential problem can be reformulated by integral equations defined on the boundary of the given domain. These methods have two main advantages, the dimension reduction of the computational domain and the simplicity for treating external problems. As a major drawback, the resulting integrals can be singular and therefore robust and accurate quadrature formulas are necessary for their numerical computation. The solution is then obtained by collocation or Galerkin procedures. With the advent of Isogeometric Analysis (IGA) a new formulation of BEMs has been studied, where the discretization spaces are splines spaces represented in B-spline form. In order to take all the possible benefits from using B-splines instead of Lagrangian basis, an important issue is the development of specific new quadrature formulas for efficiently implementing the assembly phase of the method. In this talk the problem of constructing appropriate and accurate quadrature rules, tailored on B-splines, for Boundary Integral Equations is addressed.

Analysis-suitable C^1 multi-patch isogeometric spaces
Thomas Takacs (Johannes Kepler University Linz, Austria)

Multi-patch spline parametrizations are used in geometric design and isogeometric analysis to represent complex domains. We deal with a particular class of C^0 planar multi-patch spline parametrizations called analysis-suitable G^1 (AS-G^1) multi-patch parametrizations. This class of parametrizations has to satisfy specific geometric continuity constraints, and is of importance since it allows to construct, on the multi-patch domain, C^1 isogeometric spaces with optimal approximation properties.

Such AS-G^1 multi-patch parametrizations are suitable for modeling complex planar multi-patch domains. We present the theoretical foundations, and construct a basis, and an associated dual basis, for a sufficiently large subspace of the full C^1 space. The subspace maintains the reproduction properties of traces and normal derivatives along the interfaces. In contrast to the full space, its dimension does not depend on the domain parametrization, but only on the mesh topology. Moreover, we discuss constructions of AS-G^1 multi-patch parametrizations for arbitrary domains.

A fully adaptive method for the heat equation with hierarchical B-splines
Rafael Vázquez (EPFL, Lausanne, Switzerland)

We present an adaptive algorithm for the solution of the heat equation using hierarchical B-splines and the implicit Euler method for the spatial and time discretization, respectively. Our development follows closely the lines in [1], where fully adaptive schemes have been analyzed within the framework of classical finite elements and discontinuous Galerkin methods. Our adaptive method is based on an a posteriori error estimator that essentially consists of four indicators. On the one hand, we have a time error indicator and a consistency error indicator, that dictate the time-step size adaptation. On the other hand, we have a coarsening error indicator and a space error indicator that are used to obtain a suitably adapt the hierarchical mesh at each time step. In particular, we consider the function-based space error indicators introduced in [2]. The algorithm is guaranteed to reach the final time within a finite number of operations, and keeps the space-time error below a prescribe tolerance. We finally present some numerical tests to illustrate the performance of the proposed adaptive algorithm, using the implementation of hierarchical B-splines in GeoPDEs.
Joint work with Eduardo M. Garau, Fernando D. Gaspoz, and Pedro Morin.

[1] C. Kreuzer, C.A. M\”oller, A. Schmidt, and K.G. Siebert. Design and convergence analysis for an adaptive discretization of the heat equation. IMA J. Numer. Anal. 32 (2012), 1375-1403.
[2] A. Buffa and E.M. Garau. A posteriori error estimators for hierarchical B-spline discretizations. Preprint 2016.

To register for the INdAM Workshop DREAMS, please send an e-mail to the organizers by January 15, 2018 with your Name, Affiliation, e-mail address, arrival and departure dates.

Registration is free, but mandatory.

Istituto Nazionale di Alta Matematica “Francesco Severi” (INdAM)

Piazzale Aldo Moro 5, 00185 Rome, Italy