X Xiao

P Alliez

L Busé

L Rineau

Delaunay Meshing and Repairing of NURBS Models

Keywords:

CAD models represented by NURBS surface patches are often hampered with defects due to inaccurate representations of trimming curves. Such defects make these models unsuitable to the direct generation of valid volume meshes, and often require trial-and-error processes to fix them. We propose a fully automated Delaunay-based meshing approach which can mesh and repair simultaneously, while being independent of the input NURBS patch layout. Our approach proceeds by Delaunay filtering and refinement, in which trimmed areas are repaired through implicit surfaces. Beyond repair, we demonstrate its capability to smooth out sharp features, defeature small details, and mesh multiple domains in contact.

Introduction

CAD models represented by NURBS surface patches often exhibit defects due to inaccurate representations of trimming curves, imperfect modeling processes or format migration. Such defects range from gaps to self-intersections through overlaps and islands. While small-size defects are troubleless for visualization or manufacturing purposes, they hamper the direct generation of volume meshes that are required for simulation. A common approach consists of (1) fixing the CAD model, (2) generating a surface mesh and (3) generating the final volume mesh. This approach may become laborintensive due to the trial-and-error nature of the fixing process and the iterations required between steps (1) and (2) until the volume mesh generator takes as input a valid surface mesh (both watertight and intersection-free). Streamlining modeling-and-simulation processes requires fully automated approaches capable of repairing while meshing.

Beyond repair-and-meshing, computational engineering practitioners require additional capabilities closer to modeling: e.g., smoothing sharp features, defeaturing or meshing multiple domains in contact.

Related work

We now review the related work focused on repair and meshing of NURBS models, and discuss the requirements of Delaunay meshing approaches.

Repair. Consistent and accurate CAD models have been a quest for decades [START_REF] Piegl | Knowledge-guided computation for robust cad[END_REF]. Defects in NURBS models originate from a combination of inaccurate modeling processes, numerical rounding and inaccuracies, migrations between models and data formats and imperfect automated NURBS generation algorithms. Common means to repair proceed through patching or reconstruction via design history or knowledge-guided NURBS [START_REF] Yang | Repairing cad model errors based on the design history[END_REF][START_REF] Rajab | Cad model repair using knowledge-guided nurbs[END_REF]. Commercial solutions such as CADfix or CADdoctor have been proposed to check, repair and translate multi-CAD models, where processes are as automated as possible. Attene et al. [START_REF] Attene | A lightweight approach to repairing digitized polygon meshes[END_REF][START_REF] Attene | Polygon mesh repairing: An application perspective[END_REF] attempt to heal defects locally and generate a watertight triangle mesh. Bischoff and Kobbelt [START_REF] Bischoff | Structure preserving cad model repair[END_REF] combine the local surface repair with volumetric reconstruction following the detection of artifacts in a CAD tessellation, such that a manifold output is always guaranteed. Software tools (e.g., MeshFix, Mesh-Works, PolyMender) have been proposed for repairing meshes and polygon soups, with strict robustness requirements for 3D printing [LEM * 17]. The quest for valid volume meshing with unconditional robustness culminates with the recent significant advances for tetrahedral meshing in the wild [HZG * 18] and its recent efficient variant [HSW * 20]. The robustness of the approach is confirmed, though the mesh quality is not unconditionally consistent when large defects are present in the model.

Meshing NURBS models. Mesh generation methods for NURBS models differ depending on whether they proceed in parameter space [BLG00, GBA * 17] or directly in 3D. The dependence on the input NURBS patch layout, particularly when meshing each patch in parameter space, can be removed via a process referred to as mesh quilting, which stitches at the mesh level and possibly offers additional virtues for repairing and defeaturing. Virtual topology has been adopted in commercial software tools (e.g. Ansys Workbench, Abaqus CAE, Siemens NX) to clean up the geometry with virtual topological merging or splitting operations [FCF * 08, TSRA17], which provides an alternative approach to removing dependency on input NURBS layout. Most volume mesh generation approaches generate the volume mesh once a valid surface mesh is available [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF], while the Delaunay filtering and refinement approach interleaves surface and volume meshing in 3D, with full independence to the input patch layout [START_REF] Cheng S.-W | Delaunay Mesh Generation[END_REF]. The Delaunay-based meshing approach thus removes the need for a valid surface mesh as input but requires consistent answers to intersection queries either between line queries and the input NURBS surfaces, or between point queries and the 3D domain bounded by the input NURBS surfaces. The resulting 3D mesh is always valid "by design" as results from the filtering of a 3D Delaunay triangulation, even for multiple domains in contact. Nevertheless, Delaunay refinement further requires proper initialization and conservative sizing of mesh elements with respect to the local feature size. The latter measures the distance between surface points and the medial axis of the domain [START_REF] Amenta | Surface reconstruction by voronoi filtering[END_REF]. For Delaunay refinement with weighted Voronoi diagrams, surface boundary curves are sampled with protecting balls which are treated as weighted points during Delaunay refinement. The local feature size is then extended to blend with the local gap size of polyhedral domains to determine the ball sizes [START_REF] Cheng S.-W | Rineau / Delaunay Meshing and Repairing of NURBS Models (a) NURBS model with small geometric details (b) Balls covering all features (c) Blending balls for defeaturing (d) Mesh preserving small features (e) Mesh defeaturing small features Figure 16: Small features in NURBS model are removed in the final mesh with blending balls[END_REF]. A challenge for Delaunay meshing of NURBS models is trimming defects due to the inconsistent representations of trimming curves in CAD models. Figure 1 depicts four common types of defects in a trimmed NURBS model. The presence of trimming defects can derail Delaunay filtering and refinement processes, as illustrated in Figure 2. The protecting-ball approach is further applied to repair defects in NURBS models. The first contribution to such a Delaunay-based meshing for NURBS model has been pioneered by Busaryev et al. [BDL09], by considering all surface boundaries as sharp features and by relying upon a merging process to determine the protecting balls. This approach, however, generates a Delaunay mesh confined by the input NURBS layout, and sharp features are preserved without taking into account their local geometric features.

Contributions

We seek a robust approach which offers the capability to repair while meshing, without requiring a valid surface mesh as input, and to generate volume meshes which are (1) valid by design and (2) independent of the input NURBS patch layout, without resorting to post-processing steps such as mesh quilting or remeshing. The mesh quality is guaranteed by the Delaunay refinement methodology we adopt, which provides relevant user-defined facet and cell criteria to control the meshing process. Departing from the protecting ball approach which preserves all surface boundaries during Delaunay meshing [START_REF] Busaryev | Repairing and meshing imperfect shapes with delaunay refinement[END_REF], our approach treats sharp features and smooth boundary connections differently such that the Delaunay mesh elements in smooth regions are fully independent of the NURBS layout. We also provide a thorough analysis of ball sampling and sizing in order to control the element sizes of the final mesh, especially near sharp features. We now summarize the technical contributions which form the framework of our approach for repairing and meshing NURBS models with trimming defects:

• A novel ball generation process is proposed to cover all trimmed regions, distinguishing smooth regions from sharp features. Ball sizes are adaptive to geometric features and defect sizes, and the user-defined criteria are used as prior knowledge when determining ball sizes protecting sharp features (Section 4). • Given a NURBS model as input, we analyze it in order to identify three types of areas: sharp features to preserve, trimmed areas to repair near trimming curves and defect-free areas away from sharp features and trimming curves. These areas are dealt with differently in order to provide consistent intersection queries during Delaunay refinement. In particular, the defects are healed via implicit surfaces interpolating the defect-free areas (Section 5). • We demonstrate the versatility of our framework not just for meshing, but also for smoothing out sharp features, removing small features, and meshing multiple domains. Such operations are especially relevant for computational engineering and simulation practice (Section 6).

Background on NURBS models and Delaunay meshing

In this section, we briefly review the design of CAD models with trimmed NURBS surfaces and the challenges of the inherent trimming defects to Delaunay meshing.

Trimmed NURBS models

In computer-aided design, a solid is often defined via a boundary representation, commonly referred to as B-Rep, i.e. a collection of connected surface elements which define the boundary between the interior and exterior of the solid. NURBS surfaces are commonly used to represent these surface elements for the mathematical soundness and flexibility for free-form modeling.

In a B-Rep CAD model, NURBS surfaces are trimmed in the parameter space to represent arbitrary surface boundaries which are intersection curves between surfaces. In general, the intersection computation between two NURBS surfaces leads to solving a system of nonlinear equations, and the 3D intersection curve is usually not a NURBS curve but a piecewise algebraic curve of a high degree. As a consequence, 3D intersection curves between NURBS surfaces are usually approximated via several methods, e.g. lattice evaluation, recursive subdivision or marching methods [START_REF] Hasle | Geometric modelling, numerical simulation, and optimization[END_REF]. A subsequent curve fitting scheme is often performed to generate a 3D intersection NURBS curve C which does not lie on either of the two NURBS surfaces in general, except in low-degree surface cases. This 3D intersection curve is also represented in each parameter space of the two intersecting NURBS surfaces as two low-degree B-spline (or NURBS) 2D trimming curves C trim , via curve fitting of points sampled in the parameter space [START_REF] Marussig | A review of trimming in isogeometric analysis: challenges, data exchange and simulation aspects[END_REF]. Theses 2D trimming curves can be embedded in 3D by evaluating surfaces along C trim . Consequently, the intersection curve between two NURBS surfaces has three distinct representations in the CAD model, i.e. a 3D curve and two 2D trimming curves. The fact that the images of these three curves do not coincide with each other leads to a variety of trimming defects, as already shown in Figure 1.

Delaunay meshing

Isotropic Delaunay meshing of an input 3D domain Ω hinges upon the concept of filtering and refining a 3D weighted Delaunay triangulation, whose vertices are located either inside Ω or on the domain boundary ∂Ω. The Delaunay triangulation and the Voronoi diagram of a 3D point set are dual to each other. Specifically, a Voronoi edge is the dual of a tetrahedron facet, and Voronoi vertices are located at tetrahedron circumsphere centers.

For an input 3D domain Ω, we retain the tetrahedra whose circumcenters are located inside Ω, and the boundary facets have dual Voronoi edges intersecting ∂Ω. Filtering thus requires computing intersections either between 3D points and Ω, or between Voronoi edges and ∂Ω. In addition, a restricted Delaunay triangulation derived from an initial coarse sampling of ∂Ω must be refined by inserting new vertices, referred to as Steiner points, to satisfy the mesh criteria. The Steiner points are inserted at circumcenters of poorly shaped facets and cells, which are computed as the intersections between Voronoi edges and ∂Ω.

The inevitable trimming defects in NURBS models can jeopardize the Delaunay meshing which relies on reliable and consistent intersection computation with the 3D domain. Due to the existence of defects, the intersection can be either a false negative (i.e., an intersection is missed due to gaps) or dubious points due to surface overlaps and intersections. Wrong intersections hamper the Delaunay filtering to identify valid cells and surface facets and to insert Steiner points during Delaunay refinement. As a result, the defects in a NURBS model may terminate the Delaunay refinement process early or even break the termination guarantees.

Overview

Overall procedure

The following pseudo-code Algorithm 1 outlines the proposed algorithm for meshing NURBS models with trimming defects. Each step is discussed in full detail in the following sections.

Algorithm 1 Isotropic meshing of NURBS models with trimming defects. Input: NURBS model, mesh quality parameters Output: Isotropic tetrahedral mesh 1: Ball covering of trimmed regions 2: Implicit surface approximation inside blending balls 3: Delaunay refinement in three types of subdomains NURBS models require preprocessing before Delaunay meshing, in which the trimmed regions (including the boundaries of surface patches) are singled out through ball covering for separate treatment. Sharp features (creases, corners, cusps) are covered with so-called protecting balls as in weighted Delaunay meshing approaches [DL09, CDS12], while other boundary curves are covered with blending balls inside which NURBS surfaces are substituted by approximate implicit surfaces. The connection between two adjacent NURBS surfaces is considered smooth when the surface normal directions along the common boundaries deviate by a small angle, where small refers to a user-defined angle threshold. As depicted in Figure 3, a NURBS model is divided into three subdomains for Delaunay meshing: vicinity of sharp creases, vicinity of smooth boundary curves and NURBS surfaces away from boundary curves. The trimmed regions processed by ball covering, i.e. the vicinities of boundary curves, are distinguished from defectfree areas.

The Delaunay refinement process differs depending on the three types of subdomains partitioned. In the defect-free subdomain, intersections between Voronoi dual edges and trimmed NURBS surfaces are computed efficiently via a hierarchical bounding volume

Cell criteria uniform tetrahedron circumradius Ca

Tetrahedral radius/edge ratio Cs hierarchy data structure storing Bézier surface patches, and efficient line-Bézier intersections. The subdomain covered with protecting balls is never queried for intersection during Delaunay refinement. The sizes of protecting balls depend on local feature sizes and user-defined discretization parameters, which later determine the discretization of sharp creases in the final mesh. Inside blending balls, NURBS surfaces are approximated by implicit surfaces defined through points sampled on the NURBS surfaces, such that the defects are repaired during Delaunay mesh refinement.

User parameters

The output mesh quality is controlled by user-defined parameters impacting both the ball covering and the Delaunay refinement steps. Table 1 records the main user-defined parameters required by our meshing algorithm.

In the ball covering process, a feature angle θ f between the two surface normals is prescribed to classify a sharp or smooth surface connection (set by default to 10 • ). The ball sizes, which take into account geometric features, are bounded by the prescribed minimum and maximum sizes r min and rmax. This provides the user with control over the discretization of sharp features as well as over the area of blending surfaces inside blending balls. Smaller protecting balls lead to a denser mesh around sharp features, and smaller blending balls give more faithful approximations of trimmed regions. In practice, for a better mesh quality around sharp features, protecting balls are expected to be reasonably large to cover curves with small curvatures and at the same time agree with the mesh parameters in the vicinity, and the largest ball size is restricted to rmax. The minimum ball size r min is mainly assigned to blending balls, in order to fully cover trimmed areas with small enough balls, see Section 4.2 for more details on the determination of ball sizes.

The mesh criteria in the Delaunay refinement are usually userdefined [START_REF] Jamin | Cgalmesh: a generic framework for delaunay mesh generation[END_REF], which determines the facet and cell properties in the mesh. Facet criteria control facet sizes and shapes. Facet sizes can be either uniform or varying with a sizing field, which are related to user-defined parameters on the triangle circumradius Fa and the surface approximation error F δ . The facet shape is controlled with a prescribed smallest angle Fs of the triangle. The property of tetrahedron cells is related to cell sizes and shapes. The cell size is controlled via the circumsphere radius Ca, and the cell shape is controlled via the radius-edge ratio Cs of the tetrahedron. The parameter for the cell shape should be chosen carefully to ensure the termination of the refinement process.

Intersection oracle

Computing intersections is central to the Delaunay meshing process, from the initial sampling to Delaunay filtering and refinement. The intersections are abstracted through an oracle which takes different surface representations as input [START_REF] Jamin | Cgalmesh: a generic framework for delaunay mesh generation[END_REF]. For a NURBS model, the oracle considers intersections with trimmed NURBS surfaces and blended implicit surfaces. Different intersection algorithms are activated based on whether the line segment query passes through blending balls.

Trimmed NURBS surfaces. Each intersection query is first computed with untrimmed NURBS surfaces defined over complete parametric domains. For the sake of computational efficiency, untrimmed NURBS surfaces are first decomposed into rational Bézier surfaces and stored into an AABB tree for accelerating intersection queries [START_REF] Alliez | 3D fast intersection and distance computation[END_REF]. Computing intersections between line primitives and Bézier surfaces is performed through recursive subdivision [START_REF] Sintef | The sintef spline library[END_REF]. Since the pre-images of an intersection point are obtained in the complete parametric domain, they are likely to lie outside the trimmed parametric domain bounded by trimming curves and thus should be filtered out. Only intersection points with pre-images inside the trimmed parametric domain are valid. Details on intersection computation with trimmed NURBS surfaces are provided by Shen et al. [START_REF] Shen | A line/trimmed nurbs surface intersection algorithm using matrix representations[END_REF].

Blending regions. When a line segment query intersects blending balls, it is first divided into sub-segments delineated by the boundaries of blending balls. For sub-segments inside blending balls, the intersection is computed with implicit surfaces that are substituted for trimmed regions. Details on solvers for implicit surface and intersections inside blending balls are provided in Section 5. Other sub-segments outside balls are processed by the above trimmed NURBS surface intersection. The complete set of resulting intersection points is obtained through computing intersections with all sub-segments.

Ball covering

All the trimmed regions of a NURBS model are covered with balls as a preprocessing before Delaunay meshing. Preprocess-ing via ball covering to repair trimming defects has been proposed by Busaryev et al. [START_REF] Busaryev | Repairing and meshing imperfect shapes with delaunay refinement[END_REF], in which all boundary curves are preserved with protecting balls such that the defects are never queried during Delaunay refinement. By contrast, our approach hinges upon a two-phase ball generation process, with protecting balls covering sharp features and blending balls covering smooth trimmed regions, such that the two different trimmed regions are dealt with differently later during Delaunay refinement.

The complete ball covering of a 3D intersection curve is required for a valid repair, that is, all defects must be contained inside balls. The completeness of ball covering depends on ball sizes and ball sampling processes on the curves. Busaryev et al. proposed to sample protecting balls uniformly on the images of trimming curves, and a subsequent merging operation is performed on the balls to cover the defects in trimmed regions [START_REF] Busaryev | Repairing and meshing imperfect shapes with delaunay refinement[END_REF]. The ball distribution and sizing heavily depends on this merging operation. It requires an appropriate pairing of balls to merge, which, however, cannot always be guaranteed. To avoid the side effect of the merging operation, we choose to sample balls directly along 3D intersection curves, taking into account local size of defects while determining the ball sizes. In our approach, all the 3D curves are first covered with protecting balls of adaptive sizes. After distinguishing sharp and smooth features, balls covering smooth trimmed regions are then replaced by uniformly sized balls considering defect sizes.

Sampling

For each intersection curve of two NURBS surfaces, balls are generated iteratively from one end to the other end. We follow the concept of using auxiliary balls introduced by Dey et al. [START_REF] Dey | Delaunay meshing of piecewise smooth complexes without expensive predicates[END_REF] to sample balls along a 3D curve. The ball generation begins with placing two end balls b 0 and bn, with centers at the two end points of the curve. After a ball b i with the center c i and the radius r i is determined, the next ball b i+1 is positioned with the aid of an auxiliary ball b a i . As depicted by Figure 4a, the auxiliary ball b a i is centered at the forward intersection between b i and the curve, with the radius of r i /3 to allow sufficient separation between centers of adjacent balls. The ball center c i+1 is at the forward intersection between b a i and the curve. There are two different scenarios when the ball generation reaches the other end of the curve, as depicted in Figure 4b and4c. In order to guarantee a full covering of balls over the curve and a sufficient overlapping of adjacent balls at the end, the following rules are adopted. In the first scenario when the last auxiliary ball intersects the end ball bn, the last ball b n-1 is positioned at the center of the auxiliary ball b a n-2 , with the radius of 2 3 r n-2 . In the second scenario when the ball b n-1 intersects the end ball bn, the ball is replaced with an enlarged size of 7 6 r n-1 . To ensure that the curve is completely covered by balls, adjacent balls must overlap sufficiently. In addition, a ball center should not be contained in another ball, which is required in the weighted Delaunay meshing to preserve surface features covered by protecting balls. According to [START_REF] Dey | Delaunay meshing of piecewise smooth complexes without expensive predicates[END_REF], the following sufficient overlap of two balls is enforced during ball generation, where c i and c i+1 are the two ball centers; the lower limit controls the separation between a ball center and the other ball, and the upper limit controls the minimum overlap between two balls.

max{r i , r i+1 } < ||c i+1 -c i || ≤ r i + r i+1 - 1 7 min{r i , r i+1 } , (1)

Adaptive sizing

The ball size is determined considering the ball connectivity and the desired mesh approximation. A ball should intersect sufficiently with its adjacent balls while not intersecting non-adjacent balls.

Departing from the ball refinement approach proposed by Dey and Levine [START_REF] Dey | Delaunay meshing of piecewise smooth complexes without expensive predicates[END_REF] to ensure the separation between non-adjacent balls, we achieve ball separation by measuring several feature distances at ball centers, including the distance to non-adjacent surfaces and the local feature size of the curve. Such an adaptive sizing also considers the approximation accuracy of sharp features. In the following, we use the ball curve to denote the curve covered by a ball.

Preprocessing of NURBS surfaces

To avoid the complexity of direct computations on NURBS surfaces, feature distances are computed on a tessellation of NURBS surfaces which is obtained by simply mapping 2D Delaunay meshes in parametric domains to 3D Euclidean space. Trimming curves are discretized with a good approximation relying on a faithful mapping to a set of 3D Bézier curves. Based on the fact that a 2D Bézier curve can be mapped into a 3D Bézier curve faithfully if the curve lies on a tensor-product Bézier surface, we follow the methodology presented by Lasser et al. [START_REF] Lasser | Bézier representation of trim curves[END_REF]. To summarize, the Bézier equivalents of trimming curves are divided by the knot spans of the decomposed Bézier surfaces of the NURBS surface such that each Bézier curve segment is restricted within one Bézier surface. Each trimming curve is thus exactly mapped into a set of 3D Bézier curves. The simple tessellation of NURBS surfaces with a controlled surface approximation can therefore be obtained straightforwardly, and the triangles are stored into an AABB tree [START_REF] Alliez | 3D fast intersection and distance computation[END_REF] for accelerating subsequent queries.

Surface separation distance d s sep

In order to ensure the separation between balls on different curves, the ball size considers the shortest distance between the ball center and surfaces non-adjacent to the ball curve. The shortest distance is computed by traversing the AABB tree created in the preprocessing. The triangles belonging to surfaces adjacent to the ball curve are tagged as ghost in the AABB tree, see Figure 5a, and these ghost triangles are intentionally ignored during tree search. In addition, for in-between balls on the curve, only triangles that do not intersect the two end balls are considered. However, it can be risky to get an undesired large separation distance if triangles intersecting the two end balls are excluded. We solve this issue with a hybrid scheme depicted in Figure 5b by considering dense sampling points inside these extra ghost triangles but outside end balls, including intersection points between ghost triangles and end balls. In general, the intersection points between the additional ghost triangles and balls are critical. The shortest distance between the ball center and the sampling points is then compared with the distance query with the AABB tree. The separation distance takes the smaller value.

Curve separation distance d c sep

Apart from the separation of any two balls on different curves, nonadjacent balls on the same curve should not intersect as well. This separation distance is estimated by increasing the ball size gradually until it intersects with a non-adjacent curve segment. The intersection is computed with discretized curves obtained in the NURBS preprocessing considering a prescribed approximation error. The process of gradually growing balls is illustrated in Figure 6a. The curve separation distance is therefore obtained by measuring the distance between the ball center and the non-adjacent line segment intersected.

Curve approximation d c approx

Sizes of protecting balls should also take into account the userdefined mesh quality for a good approximation of sharp features in the final mesh. As the centers of protecting balls are treated as weighted vertices in the Delaunay meshing process, the approximation of sharp features depends on the line segments connecting protecting ball centers. The deviation of the line segments from the curve is restricted within a prescribed edge approximation error δ c , which can be controlled by increasing the ball size r i gradually, as in the process of estimating the edge separation distance, until it reaches the edge approximation error δ c , see Figure 6b. The distance d c approx is the largest ball size that satisfies the edge approximation error. The ball size takes two thirds of the distance obtained, i.e. r i = 2 3 d c approx . This guarantees that the line segment connecting two adjacent balls deviates from the curve within δ c . The ball growing process terminates when either the curve separation or the curve approximation is satisfied. In addition, the curve curvature κ is considered in the curve approximation, that is, d c approx is also bounded by the curvature radius 1/κ to avoid a large ball in the region with a high curvature.

Trimming error ε

A 3D intersection curve has a trimming curve in each of its two adjacent surfaces. The images of the two trimming curves in 3D space are not identical in general, and such a discrepancy leads to trimming errors. In order to estimate the trimming error around a ball center, we compute at the ball center the normal plane to the curve and its intersection points with the images of two trimming curves. The trimming error is estimated as the largest distance among the ball center and the two intersection points. The trimming error indicates a lower bound for the ball size, and the sufficient overlapping between adjacent balls enforced with (1) in general ensures that the trimming error can be fully covered along a curve. The complete covering of trimming errors can always be guaranteed by checking if the intersections between a ball and the images of two trimming curves lie inside the adjacent balls.

Sizing formula

In addition to the feature distances described above, a ball size is often bounded by the user-defined maximum size rmax, which is necessary for the control of sharp feature discretization. Hence, the ball size is determined by the following formula:

r = min{ 1 3 min{d s sep , d c sep , 2d c approx } , rmax} , (2a)
r ≥ 2ε . ( 2b
)
Blending balls have uniform sizes given by the minimum size r min and the trimming error ε, that is,

r = max{r min , 2ε} . (3)
As a final step, the ball size obtained must satisfy the sufficient overlap between adjacent balls (1), i.e.

r ∈ {r|sufficient overlap} .

Algorithms

The ball generation for each curve is a two-phase process: adaptivesized balls for sharp feature curves and uniform-sized balls for smooth feature curves. Before the uniform-sized balls for smooth feature curves are generated, it is necessary to identify whether a curve delineates a sharp feature. The angle criterion is used for detecting sharp features, that is, the curve segment inside a ball is regarded as a sharp feature if the angle between surface normals inside the ball is greater than a user-defined angle threshold. The surface normals are evaluated at sampling points on boundary curves. All balls covering smooth features are then replaced with uniformsized balls. Finally, after all curves are covered with balls, balls at the curve ends should be merged. The overall algorithm for generating balls covering curves is outlined as follows in Algorithm 2 and the algorithm for generating balls on a single curve is provided in Appendix A.

Algorithm 2 Ball covering of trimmed NURBS models Input: Trimmed NURBS model Output: A set of balls {B i } covering all intersection curves {C i } for all curves do -Generate adaptive-sized balls (Appendix A) using sizing formula (2) -Distinguish balls covering sharp feature -Replace balls covering smooth feature curve with uniformsized ones (Appendix A) using sizing formula (3) end for Merge end balls and update ball connection information.

Repair of trimming defects

Trimming defects covered with protecting balls and blending balls are repaired in different ways. The defects are in fact neglected inside protecting balls in the Delaunay refinement, as they are in the regions that are never interrogated. The preservation of sharp features with protecting balls has been studied for the Delaunay meshing [START_REF] Busaryev | Repairing and meshing imperfect shapes with delaunay refinement[END_REF][START_REF] Dey | Delaunay meshing of piecewise smooth complexes without expensive predicates[END_REF]. In our implementation, the approximation of the sharp feature depends on the adaptive-sized protecting balls which are controlled with surface features and user-defined parameters, as described in Section 4.

In order to repair trimming defects inside blending balls, e.g. non-watertightness, non-manifold, self-intersection, etc., which jeopardize the Delaunay meshing process, we propose a new blending technique to approximate and blend the surfaces covered by blending balls in such a way that implicit surfaces are substituted for the original NURBS surfaces inside blending balls. A line segment query emitted during Delaunay refinement is divided into sub-segments after computing its intersection with boundaries of blending balls. The implicit intersection is then computed for subsegments inside blending balls.

Blending surfaces

Given discrete points in space, Duchon [START_REF] Duchon | Splines minimizing rotation-invariant semi-norms in sobolev spaces[END_REF] introduced an interpolating spline based on the minimization of a semi-norm generalizing the Fourier transform of the thin-plate bending energy. The 1D cubic spline and 2D thin-plane spline are two special cases of the Duchon's interpolating spline. We adopt its 3D Hermite interpolation variant which interpolates the values and gradients at points. This interpolating spline was used by Huang et al. [START_REF] Huang | Variational implicit point set surfaces[END_REF] for surface reconstruction from a set of sampling points, in which it is demonstrated that the interpolating scheme generates reasonable implicit surface with sparse and non-uniform sampling points, and is resilient to sampling imperfections. In our setting, Duchon's interpolating spline is a relevant choice as reliable sampling points and their surface normals can be obtained readily around blending balls. Using Duchon's implicit interpolant, the implicit surfaces inside blending balls interpolate sampling points and surface normals. In order to trade approximation accuracy for computational cost, a handful of sampling points are strategically selected on NURBS surfaces in four regions as depicted in Figure 7, which are blending ball boundaries, intersections between adjacent balls, truncated cone boundaries and 3D intersection curves, such that the implicit surface fits the blending region with satisfactory accuracy using a small number of sampling points. The polyharmonic spline with the triharmonic radial basis function (RBF) φ(x, y) = ||y -x|| 3 is considered to interpolate sampling points on surfaces. Given points p i ∈ R 3 (i = 1, 2, • • • , m) on a NURBS surface and their surface normals n i , an implicit surface interpolating the points and their normals is defined as follows,

f (x) = m ∑ j α j φ(x, p j ) + m ∑ j β T j Dyφ(x, p j ) + γ T x + η (5) with the unknown coefficient vectors α ∈ R m , β ∈ R 3m , γ ∈ R 3 and η ∈ R; Dyφ(•,
•) denotes the first derivative of triharmonic RBF w.r.t. the second variable.

The unknown coefficient vectors are obtained by considering the interpolation of point coordinates and normals as well as additional orthogonality conditions, which are

f (p i ) = 0, f (p i ) = n i (interpolation) , (6)
m ∑ i α i p i + ∑ i β i = 0, m ∑ i α i = 0 (orthogonality) . ( 7
)
The linear equation system obtained from the above conditions is expressed in the following matrix form,

Ku = F (8) with K = K 0 K 1 K T 1 0 , u =     α β γ η     and F =     0 n 0 0     . (9)
Details on entries of matrix K are provided in Appendix B.

Implicit surface interrogation

The implicit function (5) has a zero value at the intersection point between a line segment and the implicit surface defined. If the subsegment inside a blending ball intersects the implicit surface, it is well assumed that there is a single intersection since the subsegment clipped by the blending ball is very short and it becomes increasingly normal to the surface as the Delaunay refinement proceeds. We hence resort to bisection to compute the intersection. In the case of multiple intersections between the sub-segment and the implicit surface, it is enough to insert one of the intersection points in Delaunay refinement.

Since the implicit surfaces are computed individually inside each blending ball, they are not supposed to be continuous across adjacent balls, i.e. the implicit surfaces inside the overlapping domain of two adjacent balls are different. To overcome the hurdle of the non-continuity of the implicit surfaces across two adjacent balls, the following blending function is proposed to construct a continuous transition from one implicit surface to the other,

f b (x) = (1 -w) f 1 (x) + w f 2 (x) , (10)
where f 1 and f 2 are implicit functions in two adjacent balls (centered at c 1 and c 2 respectively, as illustrated in Figure 8a), w is the weight given by

w = ||x -v 1 || ||v 2 -v 1 || (11)
with v 1 and v 2 the intersections between the line c 2 -x and two planes P 1 and P 2 orthogonal to the direction c 2 -c 1 .

Hence, in each bisection iteration when computing the intersection with implicit surfaces, if the point lies in the overlapping domain of adjacent balls it is evaluated with the blending function (10), and the implicit function (5) is used for points in one blending ball exclusively.

Since NURBS surfaces are approximated with implicit surfaces inside blending balls, the discrepancy between the two representations is inevitable, though small enough to be negligible, near the balls. For enhancing the robustness of the intersection oracle, both the NURBS intersection and the implicit intersection are computed in the transition region near the ball boundary, and the implicit intersection takes priority if both intersections exist. In order to enforce the watertightness in the transition region, we introduced sufficient interpolating points located strategically around blending balls as illustrated in Figure 7. In practice, the extreme case of a dual Voronoi edge passing through the tiny gap around a blending ball does not hinder the Delaunay meshing process, because it may only occur in the first few steps of Delaunay refinement depending on the initial Delaunay sampling, and Voronoi edges tend to be normal to the surface as the refinement proceeds.

Results

A selection of examples are shown in the following to demonstrate the capability of our approach to simultaneously repair and mesh NURBS models with trimming defects and its added values to NURBS defeaturing. More examples are given in the appendix. The general Delaunay refinement framework of the CGAL library is used in our implementation [START_REF] Jamin | Cgalmesh: a generic framework for delaunay mesh generation[END_REF], and the SISL library [START_REF] Sintef | The sintef spline library[END_REF] is used for the spline intersection computation.

Parameterization-free meshing

The parameterization-free meshing is demonstrated with a mechanical part as shown in Figure 9a. Sharp features are covered with protecting balls, and other smooth trimmed regions are covered with blending balls. The edge approximation error is considered for protecting balls, and blending balls have uniform sizes. Figure 9b shows the Delaunay surface mesh with a prescribed surface approximation error. The distribution of mesh elements is not constrained by the layout of trimmed NURBS surfaces, which is a key feature of our approach. All the sharp boundary curves are preserved and discretized with protecting ball centers. However, the weighted Delaunay meshing using protecting balls has an inherent issue of generating excessively small elements between sharp features subtending small angles. We resort to local remeshing for these cases [START_REF] Faraj | Multi-material adaptive volume remesher[END_REF][START_REF] Tournois | Tetrahedral remeshing[END_REF] to replace the ill-shaped elements. Figure 9c shows the comparison of the mesh elements in the cusp region before and after remeshing. The mesh quality is reported in Appendix C.

Our method can be extended to multidomain meshing straightforwardly, and different mesh quality criteria can be prescribed in domains. The multidomain meshing of turbines is given in Appendix D.

Defect repair

As depicted in Figure 1, four common types of defects exist in a trimmed NURBS model, including gaps, overlaps, intersections and mismatched end points. The NURBS model in Figure 10 is considered for the repair of surface overlaps. The blending balls are adaptive to the overlap size to ensure the complete covering of the defects, see Figure 10b. Inside blending balls, implicit surfaces are computed to bridge the two adjacent surfaces. The overlap defect is repaired in the Delaunay mesh controlled by the surface approximation error criterion, as shown in Figure 10d, and the mesh elements cross the overlap region in the NURBS model.

Figure 11 shows that multiple types of defects are repaired systematically with our approach, see Figure 1b for the close-up view of the defects, which includes gaps, intersection and mismatched ends between adjacent surfaces. The defects in sharp features are covered with protecting balls such that the defect regions are not interrogated in the Delaunay refinement and thus repaired automatically in the Delaunay mesh. The defects in the smooth trimmed regions are covered with blending balls. The two adjacent surfaces are bridged with implicit surfaces inside blending balls as shown in Figure 11c, such that the implicit surfaces are considered when the defect regions are interrogated during Delaunay meshing. Figure 11d shows the Delaunay mesh around the defect regions repaired with protecting and blending balls. The sharp features are repaired and preserved in the final mesh by connecting protecting ball centers, and other defect regions are meshed with crossing elements.

A watertight surface mesh is usually required by most NURBS meshers to generate a volume mesh. However, this requirement may not be guaranteed for some immoderately bad CAD designs. Figure 12a shows a NURBS model with tapering surfaces penetrating between surfaces, with unobserved defects around the tips. Gmsh [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF] cannot mesh this non-watertight case with small element sizes around penetrating tips, leaving the tapering surfaces unmeshed (see Figure 12b) and thus preventing from generating a volume mesh. The Delaunay mesh in Figure 12c is generated with our approach, in which the narrow tips are covered with protecting balls such that the intricate details due to a bad design are neglected in the meshing. The protecting balls are obtained with the sizing formula presented in Section 4.2. It is feasible to replace the small sized protecting balls with larger ones, in order to improve the mesh quality in the defect regions.

It is worth mentioning that a robust tetrahedral meshing algorithm TetWild proposed recently by Hu et al. [HZG * 18] can generate a tetrahedral mesh from a defect-laden surface mesh. Though it is not devised specifically for NURBS meshing, it provides another means of obtaining the tetrahedral mesh from an initial NURBS tessellation. The algorithm has been demonstrated unconditionally robust for general defects. One needs to select proper parameters in order to get a volume mesh with consistent quality. We compared with TetWild and its efficient variant FTetWild [HSW * 20] to generate a volume mesh from a tessellation which maintains trimming defects of a NURBS model (see Section 4.2.1). Figure 13 shows the diamond tessellation that challenges TetWild to generate a mesh with consistent quality, though different parameters are selected. Gaps can be robustly filled but may yield artifacts such as bumps in the mesh. FTetWild is faster and generates more consistent mesh qualities than TetWild for the diamond model, though we can still observe missing tetrahedra in the mesh as depicted in Figure 13e. Figure 14 compares TetWild with our method on a NURBS model with multiple defects including large gaps, overlaps and intersections. We observe that the mesh quality given by TetWild with default parameters is not consistent in similar regions. By contrast, our method provides a more flexible and consistent control on the mesh quality while repairing defects.

Comparison with Busaryev et al.

Though protecting balls have also been utilized by Busaryev et al. [START_REF] Busaryev | Repairing and meshing imperfect shapes with delaunay refinement[END_REF] to repair defects in NURBS models, our approach differs substantially and offers additional benefits. Figure 15 compares the Delaunay mesh obtained with the two approaches. As discussed in Section 4, the sizes of protecting balls are determined adaptively based on geometric features, avoiding unnecessarily small or large ball sizes due to the shrinking and merging processes proposed by Busaryev et al. Moreover, the use of blending balls to repair defects in smoothly connected trimmed regions removes the dependence of mesh elements on the NURBS layout. In addition to the defect repair, blending balls can also be leveraged for removing small features in a NURBS model. Such an operation has been explored in [START_REF] Busaryev | Repairing and meshing imperfect shapes with delaunay refinement[END_REF] Section 8 but by using protecting balls, which leads to unnecessary constraints in the final mesh. The small geometric details in Figure 16a are covered with protecting balls and preserved with small element sizes, see Figure 16b and16d

Performance

The timings of generating the Delaunay meshes presented in the examples are reported in Table 2. The computation was performed on a laptop computer with Intel(R) Xeon(R) CPU E5-2630 v4

Limitations

Since ball sizes are related to trimming errors in order to fully cover all the defects (see sizing formula (2) and (3)), immoderate defect sizes require large balls which put the proposed approach at risk when defect sizes are comparable to surface dimensions. Large protecting balls can swallow entire surface features, leading to undesired mesh quality in the repaired region, and they may break the ball separation rule, which can derail or terminate the Delaunay meshing process. Such cases with visually large defects require resorting to CAD repair as a preprocessing step.

Conclusions

We have presented a method for simultaneously repairing and meshing NURBS models with trimming defects, through providing a consistent interface to a Delaunay-based isotropic simplicial mesh generation framework. Defects such as gaps and overlaps due to trimming are healed through substituting local NURBS surfaces by implicit surfaces defined within a set of healing balls covering the defects. The ball covering process is guaranteed to fully cover all trimmed regions. It echoes the ball sampling process that has been proved to satisfy the requirement for weighted Delaunay meshing [START_REF] Dey | Delaunay meshing of piecewise smooth complexes without expensive predicates[END_REF]. The implicit surface using the Duchon's interpolating spline provides a tight approximation to the original NURBS surface. The blending property of the above implicit surfaces also provides us with a means to smooth out sharp features and defeature small details. We also contribute an algorithm for estimating the local feature size of NURBS surfaces and sharp features. Such a local feature size is used for sizing the above healing balls. We demonstrate the robustness and versatility of our method: the output volume meshes are valid by design (without overlaps or inverted elements), independent of the input NURBS patch layout and parameterization, and conforming at interfaces between multiple domains in contact.

The proposed approach is mainly limited by exceedingly large defect sizes such that balls covering the defects are larger than surface dimensions. Consequently, the entire surface feature may be swallowed by large protecting balls, and the blending balls may fail to heal the defects as non-adjacent surfaces are covered. The robustness of the implicit surface interrogation near the transition region is currently enhanced with a workaround for extremely rare cases of dual Voronoi edges almost tangential to surfaces and passing through gaps between implicit surfaces and original NURBS surfaces. In addition, the current mesh-and-repair algorithm provides insufficient feedback to users about the set of possible parameters for automatic defeaturing in the repair process. There is also an intrinsic issue of Delaunay-based approaches in which mesh elements are distorted and ill-shaped near sharp features subtending small angles. Resorting to a remeshing post-processing is a solution [START_REF] Faraj | Multi-material adaptive volume remesher[END_REF], but having a single integrated algorithm would be more elegant.

Our Delaunay-based mesh-and-repair algorithm has guaranteed to generate valid isotropic simplicial meshes for NURBS models. As future work, we would like to extend the approach to anisotropic meshing and higher-order mesh elements [FAB * 18], which are relevant for numerical analysis. In addition, we wish to complement our meshing algorithm by a scale-space analysis stage, providing users with a preview on the whole spectrum of possible operations on NURBS models in terms of repair, defeaturing and smoothing.

  (a) NURBS models with defects in shaded regions (b) Close-up views of defects

Figure 1 :

 1 Figure 1: Four common types of defects in NURBS models. From left to right: gaps, overlaps, intersections and mismatched end points.

( a )Figure 2 :

 a2 Figure 2: Delaunay triangulation restricted to 2D multidomain Ω and its boundary ∂Ω. Delaunay triangulation can fail with wrong intersections between Voronoi dual edges and the boundary due to boundary defects (problematic Voronoi edges are depicted in red in (d)).

Figure 3 :

 3 Figure 3: A NURBS model is divided into three types of subdomains: sharp creases and corners (red), smooth boundary curves (blue) and NURBS surfaces away from boundary curves (gray). The intersection computation differs among these three types of subdomains.

( a )

 a Ball positions determined with auxiliary balls. (b) First end scenario.(c) Second end scenario.

Figure 4 :

 4 Figure 4: Sampling of balls along a curve.

X

  . Xiao, P. Alliez, L. Busé & L. Rineau / Delaunay Meshing and Repairing of NURBS Models (a) Ghost triangles (dashed gray) in adjacent surfaces. (b) Additional ghost triangles (dashed) overlapped with end spheres and extra sampling points (blue dots).

Figure 5 :

 5 Figure 5: Computing the surface separation distance at a ball center.

( a )Figure 6 :

 a6 Figure 6: Determination of ball sizes considering curve separation distance and curve approximation with ball growing method. The ball size increment depends on the curve discretization, and each dashed circle represents an increment step.

X

  . Xiao, P. Alliez, L. Busé & L. Rineau / Delaunay Meshing and Repairing of NURBS Models (a) Sampling points (b) Surface normals

Figure 7 :

 7 Figure 7: Sampling points and their surface normals for constructing the interpolating implicit surface. Red surface points are on the ball boundary; blue points at intersections between adjacent balls; green points on the truncated cone inside the ball; yellow points on the 3D intersection curve.

(a)

  Blending function (b) Implicit surface visualization

Figure 8 :

 8 Figure 8: Continuous implicit surfaces between adjacent surfaces with a blending function. The implicit surface is visualized by sampling, with blue surface points inside individual balls and red surface points inside two adjacent balls evaluated with blending function.

( a )

 a Protecting and blending balls (b) Mesh with surface approximation (c) Mesh elements in cusp regions

Figure 9 :

 9 Figure 9: Parameterization-free meshing with sharp feature preservation. The inherent issue of Delaunay mesh in cusp regions is solved with remeshing. In (c): Left -ill-shaped elements in the cusp region; right -remeshed elements with sharp features preserved.

  (a) NURBS model (b) Blending balls covering overlaps (c) Implicit surface (d) Delaunay mesh repairing overlaps

Figure 10 :

 10 Figure 10: Repair of overlap between NURBS surfaces. The boundary curves of the overlapping surfaces are depicted in (b). Implicit surfaces inside blending balls are visualized with point sampling in (c).

Figure 11 :

 11 Figure 11: Repair of multiple types of defects, including gaps, intersections and end mismatch.

( a )

 a NURBS model with a tapering surface (b) Gmsh result with missing elements (in red regions) (c) Delaunay mesh (d) Close-up view at the penetrating tip

Figure 12 :

 12 Figure 12: Repair of a non-watertight NURBS model (from the ABC library) near penetrating tips.

Figure 13 :

 13 Figure 13: Defects in Diamond model challenge TetWild and fTetWild to generate a mesh with consistent quality in some regions (highlighted in red in (c), (d) and (e) for visualization purpose). L and e are ratios of the bounding box, controlling the target edge length and the envelop size, respectively.

  . To remove the small geometric details, blending balls are used to cover the target regions containing small details (see Figure 16c). The geometric details are replaced with implicit surfaces inside blending balls and thus neglected automatically during Delaunay refinement and defeatured in the final mesh as shown in Figure 16e. In Appendix D, another application of blending spheres to smooth out sharp features is demonstrated.

( a )

 a Triangulated tessellation of a NURBS model (from the ABC library) that contains large gaps, overlaps and intersections. (b) Tetrahedral mesh obtained with TetWild (L and e are default values, 0.05 and 0.001 respectively), with either missing or extra tetrahedra in some regions (highlighted in red for visualization purpose).

( c )

 c Defects are repaired in the mesh and mesh quality is consistent with our method.

Figure 14 :

 14 Figure 14: Multiple defects of large gaps, overlaps and intersections pose a challenge to TetWild in some regions without tuning meshing parameters.

X

  . Xiao, P. Alliez, L. Busé & L. Rineau / Delaunay Meshing and Repairing of NURBS Models (a) Delaunay mesh with Busaryev's method (figure from [BDL09]).(b) Delaunay mesh with our method.

Figure 15 :

 15 Figure 15: Comparison with Busaryev's method. Sharp features are preserved considering geometric features and mesh elements near smooth connections, and are not constrained by the NURBS layout.

( a )

 a Delaunay mesh with coarse cells. (b) Delaunay mesh with dense cells.

( d )

 d Mesh quality of mesh (b). The ratio of 1 represents a regular tetrahderon cell or an equilateral triangle facet.

Figure 17 :

 17 Figure 17: Mesh qualities measured by cell and facet shapes.

( a )

 a NURBS turbine model (b) Spherical envelope (c) Multidomain mesh (d) Turbine volume meshes with different mesh criteria

Figure 18 :

 18 Figure 18: Multidomain meshing of turbines inside a spherical envelope.

Figure 19 :

 19 Figure 19: Sharp features are smoothed out with blending balls. Left to right: 1) sharp features preserved with protecting balls; 2) protecting balls can be switched to blending balls to cover sharp features; 3) smoothed sharp features; 4) visualization of implicit surfaces by sampling inside blending balls.

Figure 20 :

 20 Figure 20: Axle

Figure 21 :

 21 Figure 21: CoolingDown

Figure 22 :

 22 Figure 22: CoolingUp

Figure 25 :

 25 Figure 25: Stator

Figure 26 :

 26 Figure 26: A complex mechanical part (609 NURBS surfaces) and its Delaunay surface mesh. Apart from trimming defects, the model also suffers from topological defects of trimming curves, see Appendix D for more details. The locally dense mesh is due to short sharp feature curves, as shown in close-up views.

Table 1 :

 1 User-defined parameters.

		maximum ball size rmax
	Ball covering	minimum ball size r min
		feature angle θ f
		uniform triangle circumradius Fa
	Facet criteria	smallest triangle angle Fs
		surface approximation F δ

Table 2 :

 2 Running time (in seconds) of ball generation and Delaunay meshing with balls and mesh statistics.

	Model	ball cover	meshing	#balls	#facets	#tetrahedra
	Fig. 9	137.1	262.2	20228	34134	100049
	Fig. 10	23.7	44.3	3239	12936	23046
	Fig. 11	5.8	18.4	22312	12492	25410
	Fig. 12	102.6	426.6	27037	67956	138655
	Fig. 14	11.1	1355.1	13801	62920	120820
	Fig. 15	6.0	139.2	24932	35412	67339
	Fig. 16	4.7	124.8	5287	26928	52327

  Mesh quality of mesh (a). The ratio of 1 represents a regular tetrahderon cell or an equilateral triangle facet.

		8.0%			10.0%			10.0%
		7.0%
		6.0%			8.0%			8.0%
	Distribution	3.0% 4.0% 5.0%		Distribution	4.0% 6.0%		Distribution	4.0% 6.0%
		2.0%			2.0%			2.0%
		1.0%
		1.08 0.0%	177.78		0.0073 0.0%	0.9800		0.0056 0.0%	1.0000
		Dihedral angles of tetrahedral cells			smallest edge / facet circumradius ratio		smallest edge / facet circumradius ratio
	Distribution	0.0% 1.0% 2.0% 6.0% 7.0% 8.0% (c) 3.61 3.0% 4.0% 5.0%	173.93	Distribution	0.0061 0.0% 2.0% 8.0% 10.0% 4.0% 6.0%	0.9991	Distribution	0.0059 0.0% 2.0% 8.0% 10.0% 4.0% 6.0%	1.0000
		Dihedral angles of tetrahedral cells			smallest edge / facet circumradius ratio		smallest edge / facet circumradius ratio

else // the first ending scenario 9: The interpolation condition at each sampling point is expressed as

Hence, the matrix K in Section 4.1 is computed via

where

and

Appendix C: Mesh quality measures Figure 17 gives the mesh quality of the mechanical part with respect to cell and facet shapes, measured by dihedral angles of tetrahedra and the ratios of smallest edge and cicumradius of tetrahedra and triangles. The ratios for tetrahedra are scaled by √ 6/4 such that a regular tetrahedron (with four equilateral triangular faces) has a ratio of 1. Similarly, the ratios for triangles are scaled by √ 3/3 such that an equilateral triangle has a ratio of 1.

The two meshes in Figure 17a and 17b are generated with the same parameters except the cell sizes. The meshing parameters (see Table 1 in Section 3.1) are rmax = 5, r min = 0.2, θ f = 10 • , Fs = 22 • , F δ = 0.1, Cs = 2. The cell sizes are 5 and 1, respectively. The mesh quality has been optimized with techniques implemented in CGAL.

Appendix D: Extended applications

We extend our approach to Delaunay meshing of multiple domains. The multidomain meshing is a key preprocessing step for the numerical analysis of fluid dynamics and aerodynamics, which usually requires volume meshes of both the objects and the ambient environment (e.g. fluid and air). We consider the meshing of NURBS turbine models immersed in a spherical envelope as shown in Fig- ure 18b. The Delaunay mesh in the multiple domains is depicted with different colors in Figure 18c. The mesh quality can be usertuned in a flexible way to set different mesh criteria for each domain. In this example, cell sizes are uniform outside turbines, except the bounding sphere region of the turbine colored in yellow (see Figure 18c). The volume meshes of the turbines are shown in Figure 18d. Two turbines are meshed with different surface approximation errors and cell sizes. The meshes of different domains conform at their interfaces.

Another possible application of blending balls is to smooth out sharp features. The protecting balls can be switched to be blending ones such that the sharp features are approximated with smooth implicit surfaces bridging the two adjacent surfaces. Figure 19 shows the comparison of sharp features covered with protecting balls and blending balls, where the sharp features are preserved inside protecting balls and smoothed out inside blending balls.

Appendix E: Additional examples

Figures 20 to 25 show the Delaunay meshes of STEP models in Visionair repository. Figure 26 is a complex mechanical part containing a number of small geometric details. Apart from trimming defects, the model also suffers from topological defects of nonclosed or self-intersected trimming curves. We detected and solved these specific defects in the preprocessing of NURBS models by reconstructing piecewise linear trimming curves. Short sharp feature edges lead to locally dense mesh, and this is one thing that the algorithm is to be improved to consider a merging operation of short sharp edges.


Keywords: Damage detection, Subspace methods, Temperature effects, Model interpolation, Laboratory beam structure, Climate chamber 1

Structural health monitoring (SHM) of civil structures often is limited due to changing environmental conditions, as those changes affect the structural dynamical properties in a similar way like damages can do. In this article, an approach for damage detection under changing temperatures is presented and applied to a beam structure. The used stochastic subspace-based algorithm relies on a reference null space estimate, which is confronted to data from the testing state in a residual function. For damage detection the residual is evaluated by means of statistical hypothesis tests. Changes of the system due to temperature effects are handled with a model interpolation approach from linear parameter varying system theory. From vibration data measured in the undamaged state at some few reference temperatures, a model of the dynamic system valid for the current testing temperature is interpolated. The reference null space and the covariance matrix for the hypothesis test is computed from this interpolated model. This approach has been developed recently and was validated in an academic test case on simulations of a mass-spring-damper. In this paper, the approach is validated experimentally on a beam structure under varying temperature conditions in a climate chamber. Compared to other approaches, the interpolation approach leads to significantly less false positive alarms in the reference state when the structure is exposed to different temperatures, while faults can still be detected reliably.

INTRODUCTION

Damage detection in the context of automated vibration-based structural health monitoring (SHM) often is limited due to environmental effects [START_REF] Sohn | Effects of environmental and operational variability on structural health monitoring[END_REF], such as temperature. Those methods aim at detecting changes in the dynamical behavior of the structure, as it is affected by damages. However, the dynamical properties usually are also affected by changing temperatures, which may lead to false alarms in the undamaged state. Consequently, a damage detection method is required, which can differentiate between changes due to damages and those caused by variation in temperature.

In the literature some approaches for temperature rejection in vibration-based damage detection methods can be found. Features robust to temperature effects are used, such as in [START_REF] Peeters | One year monitoring of the Z24 bridge: environmental influences versus damage effects[END_REF], where a regression analysis of the natural frequency is done after a system identification step. Others, e.g. [START_REF] Worden | Inferential parameterization using principal curves[END_REF], use novelty detection approaches, defining normal reference conditions firstly and then the corresponding system properties.

This paper uses an interpolation approach to account for temperature changes in stochastic subspace-based damage detection methods. Those output-only algorithms are promising for the application to automated vibration-based SHM, as they do not require information about the excitation and are able to handle the data uncertainties. A damage detection method introduced in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF] relies on a reference null space estimate of the output covariance Hankel matrix, which is confronted to data from the testing state in a residual function. The residual is evaluated statistically, to decide if it is significantly different from 0. The test value of this evaluation can be used as damage indicator.

First temperature rejection approaches for the stochastic subspace-based damage detection are proposed. They use averaged system matrices, computed from data at several reference temperature states [START_REF] Balmes | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF]. In this algorithm no information about the current testing temperature is needed. In another approach [START_REF] Fritzen | Elimination of temperature effects on damage detection within a smart structure concept[END_REF] several reference matrices are considered, which are valid for a specific temperature range in the damage detection step, using the information of the testing temperature.

In [START_REF] Viefhues | Subspacebased damage detection with rejection of the temperature effect and uncertainty in the reference[END_REF], the authors have presented an interpolation approach to account for temperature changes with the stochastic subspace-based damage detection, based on a model interpolation method from [START_REF] Zhang | LPV system local model interpolation based on combined model reduction[END_REF]. For a system whose behavior depends on an external parameter p, but shows linear behavior at fixed working points of this parameter, an interpolation approach for linear parameter varying (LPV) systems can be used to find a model that is valid at the new parameter value. From this interpolated model a reference null space is computed, corresponding to the parameter in the testing state. In the present context the parameter is set to be the temperature.

Up to now, the interpolation-based temperature rejection approach was validated theoretically on numerical data of a mass-spring-damper with a very simple temperature model in [START_REF] Viefhues | Subspacebased damage detection with rejection of the temperature effect and uncertainty in the reference[END_REF] and [START_REF] Viefhues | Fault detection for linear parameter varying systems under changes in the process noise covariance[END_REF]. This motivates for a deeper look into the applicability of the algorithms to real structures in this study, and an application is presented to data from a reinforced concrete beam recorded in a climate chamber at defined temperatures and at three different system states. Furthermore, some conclusions regarding the choice of the weighting functions are made.

This paper is organized as follows: In section 2 the theoretical background of the stochastic subspace-based damage detection

Stochastic subspace-based damage detection of a temperature affected beam structure

Eva Viefhues 2 Univ. Gustave Eiffel, Inria, COSYS-SII, I4S, Rennes, France email: eva.viefhues@bam.de, michael.doehler@inria.fr, patrick.simon@bam.de, ralf.herrmann@bam.de, falk.hille@bam.de, laurent.mevel@inria.fr method and the model interpolation approach to account for temperature changes is recalled. Section 3 introduces the laboratory tests in the climate chamber and the application of the algorithms to the laboratory data. The section concludes with some suggestions for the choice of the weighting function.

THEORETICAL BACKGROUND

Stochastic subspace-based damage detection

The dynamical behavior of a civil structure at some temperature T = Tj can be described by a stochastic discrete time state-space model 𝑥 𝑘+1,𝑗 = 𝐴 𝑗 𝑥 𝑘,𝑗 + 𝑤 𝑘,𝑗 𝑦 𝑘,𝑗 = 𝐶 𝑗 𝑥 𝑘,𝑗 + 𝑣 𝑘,𝑗 ,

where 𝑥 ∈ ℝ 𝑛 is the state of a system with model order 𝑛, 𝑦 ∈ ℝ 𝑟 is the measured output at 𝑟 sensor positions, and 𝑘 denotes the time index within the sample. The system matrices 𝐴 𝑗 = 𝐴(𝑇 𝑗 ) and 𝐶 𝑗 = 𝐶(𝑇 𝑗 ) are the state transition matrix with dimension 𝑛 × 𝑛 and the observation matrix with dimension 𝑟 × 𝑛, respectively. For the state noise 𝑤 and the output noise 𝑣 white noise is assumed.

From a data sample of length 𝑁 a block Hankel matrix is built from the output covariances 𝑅 𝑖,𝑗 = 𝑬(𝑦 𝑘+𝑖,𝑗 𝑦 𝑘,𝑗 𝑇 ) = 𝐶 𝑗 𝐴 𝑗 𝑖-1 𝐺 𝑗 , where ⅈ denote the time shift and 𝐺 𝑗 = 𝑬(𝑥 𝑘+1,𝑗 𝑦 𝑘,𝑗 𝑇 ). The Hankel matrix writes as

ℋ j = [ 𝑅 1,𝑗 ⋯ 𝑅 𝑞,𝑗 ⋮ ⋱ ⋮ 𝑅 𝑝+1,𝑗 ⋯ 𝑅 𝑞+𝑝,𝑗 ]. (2)
Due to the relation between the system matrices and the Hankel matrix, changes in the dynamical behavior of the structure will lead to changes in the Hankel matrix.

For damage diagnosis the left null space 𝑆 𝑗 of the Hankel matrix is used, computed in the undamaged reference state. It holds 𝑆 𝑗 𝑇 ⋅ ℋ 𝑗 = 0 , and with the estimate ℋ ̂ of the Hankel matrix in the testing state a data-based residual 𝜁 can be formulated as

𝜁 = √𝑁 vec(𝑆 𝑗 𝑇 ℋ ̂), (3)
where ′vec′ denotes the stacking vectorization operator. If the system has changed, 𝑆 𝑗 is no longer the left null space of the tested Hankel matrix and the residual mean value will differ from zero. The residual is evaluated statistically in a hypothesis test (e.g. from [START_REF] Döhler | Subspace-based damage detection under changes in the ambient excitation statistics[END_REF]), which writes as

𝑡 = 𝜁 𝑇 Σ 𝜁 -1 𝜁. (4)
The asymptotic residual covariance Σ ζ can be computed from the asymptotic covariance Σ ℋ of the vectorized Hankel matrix with

Σ 𝜁 = 𝐽 𝜁,ℋ Σ ℋ 𝐽 𝜁,ℋ 𝑇 , (5)
where 𝐽 𝜁,ℋ = 𝐼 ⊗ 𝑆 𝑇 . An estimate of the covariance from a data sample 𝑦 𝑘 , 𝑘 = 1, … , 𝑁 separated into 𝑛 𝑏 blocks of length 𝑁 𝑏 , with 𝑛 𝑏 • 𝑁 𝑏 = 𝑁, writes as [10]

Σ ̂ℋ = 𝑁 𝑏 𝑛 𝑏 -1 ∑ vec(ℋ ̂(ℎ) -ℋ n𝑏 ℎ=1 ) vec(ℋ ̂(ℎ) -ℋ ̂)𝑇 . (6)
In a training phase a threshold for the asymptotically 𝜒 2distributed test statistic 𝑡 is computed with data from the healthy structure, allowing a given type I error (false positives). The test value in the monitoring state is compared to this threshold.

Model interpolation in the stochastic subspace-based damage detection for temperature rejection

When a mechanical structure is exposed to changing temperature, this will affect the dynamical behavior of the structure. Consequently, the Hankel matrix ℋ ̂ in the monitoring state under temperature T does not correspond to the null space 𝑆 𝑗 in the reference state if the reference temperature Tj varies from the current temperature T. Besides changes due to damage, also changes due to temperature variation will lead to a residual mean deviating from 0, which may result in an increasing number of false alarms.

For a robust and reliable damage detection with the stochastic subspace-based method the choice of an adequate reference null space and a proper covariance is of crucial importance. In this context the model interpolation approach for linear parameter varying systems from [START_REF] Zhang | LPV system local model interpolation based on combined model reduction[END_REF] is used in [START_REF] Viefhues | Subspacebased damage detection with rejection of the temperature effect and uncertainty in the reference[END_REF] for a damage detection algorithm robust to temperature variation.

For the model interpolation, local models of the state space models in (1) are formulated for 𝑗 = 1, . . . , 𝑚 reference temperatures 𝑇 𝑗 . Due to different state bases used for the estimation of the local state space models, the state space matrices cannot be interpolated directly. Instead the 𝑚 local models are combined to a global state-space model, without requiring formulation of coherent local models [START_REF] Zhang | LPV system local model interpolation based on combined model reduction[END_REF]:

𝑥 𝑘+1 = 𝐴 𝑥 𝑘 + 𝑤 𝑘 𝑦 𝑘 = 𝐶(𝑇) 𝑥 𝑘 + 𝑣 𝑘 , (7)
where

𝑥 𝑘 = [𝑥 𝑘,1 𝑇 … 𝑥 𝑘,𝑚 𝑇 ] 𝑇 , 𝐴 = [ 𝐴 1 ⋱ 𝐴 𝑚 ] , 𝑤 𝑘 = [𝑤 𝑘,1 𝑇 … 𝑤 𝑘,𝑚 𝑇 ] 𝑇 , 𝐶(𝑇) = [𝜌 1 (𝑇)𝐶 1 … 𝜌 𝑚 (𝑇)𝐶 𝑚 ] , 𝑣 𝑘 = ∑ 𝜌 𝑗 (𝑇)𝑣 𝑘,𝑗 (𝑡) 𝑚 𝑗=1
.

The local model for a testing temperature is interpolated with a weighting function 𝜌 𝑗 (𝑇). The weighting function can be chosen by taking the effect of the temperature to the dynamical behavior qualitatively into account. For section-wise approximated linear effects e.g. bell-shaped functions may work, centered at the testing temperature T:

∑ 𝜌 𝑗 (𝑇) = 1 𝑚 𝑗=1 . ( 8
)
Damage detection with the data-driven subspace-based residual (3) does not require the state space matrices of the reference state explicitly. It operates on the null space of the Hankel matrix of the output covariances [START_REF] Peeters | One year monitoring of the Z24 bridge: environmental influences versus damage effects[END_REF]. As the null spaces relating on an arbitrary state basis cannot be interpolated easily, the Hankel matrices in the reference state are interpolated instead. Consider the parameter dependent formulation of the output covariances as

𝑅 𝑖 (𝑇) = 𝐶(𝑇)𝐴 𝑖-1 𝐺 = ∑ 𝜌 𝑗 2 (𝑇)𝐶 𝑗 𝐴 𝑗 𝑖-1 𝐺 𝑗 𝑚 𝑗=1 = ∑ 𝜌 𝑗 2 (𝑇)𝑅 𝑖,𝑗 𝑚 𝑗=1 , ( 9
)
where the global 𝐺 writes as

𝐺 = [ 𝜌 1 (𝑇)𝐺 1 ⋮ 𝜌 𝑚 (𝑇)𝐺 𝑚 ] . ( 10
)
The interpolated Hankel matrix can be computed to

ℋ(𝑇) = ∑ 𝜌 𝑗 2 𝑚 𝑗=1 (𝑇)ℋ 𝑗 . ( 11
)
The interpolated null space matrix 𝑆(𝑇) is the left null space of [START_REF] Simon | Vibration-based structural health monitoring of a reinforced concrete beam subject to varying ambient temperatures using Bayesian methods[END_REF] and replaced in the residual

𝜁 = √𝑁 vec(𝑆(𝑇) 𝑇 ℋ ̂). ( 12
)
The covariance of the interpolated Hankel matrix follows as

𝛴 ℋ (𝑇) = ∑ 𝜌 𝑗 4 𝑚 𝑗=1 (𝑇)𝛴 ℋ 𝑗 , ( 13
)
and the residual covariance as

𝛴 𝜁 = 𝐽 𝜁,ℋ 𝛴 ℋ 𝐽 𝜁,ℋ 𝑇 (14)
where 𝐽 𝜁,ℋ = 𝐼 ⊗ 𝑆(𝑇) 𝑇 . Finally, the corresponding test writes

𝑡 = 𝜁 𝑇 𝛴 𝜁 -1 𝜁. ( 15
)

APPLICATION TO A BEAM STRUCTURE UNDER CHANGING TEMPERATURES

In former studies in [START_REF] Viefhues | Subspacebased damage detection with rejection of the temperature effect and uncertainty in the reference[END_REF] and [START_REF] Viefhues | Fault detection for linear parameter varying systems under changes in the process noise covariance[END_REF], the model interpolation approach for a robust subspace-based damage detection was evaluated on numerical data from a mass-spring-damper with very simple temperature models. The effect of temperature on the dynamical behavior of the structure was simulated by decreasing stiffness of the springs with increasing temperature. In this work, the method is evaluated under more realistic conditions on experimental data of a laboratory beam structure.

Laboratory tests of a beam structure in a climate chamber

A reinforced concrete beam is tested in a climate chamber at defined temperatures and at different system states. The test setup is shown in Figure 1. The beam is bore at two points and the span is 2.72 m. Three geophones measure vertical vibration velocities on the top of the beam.

The structure is excited by an external shaker with a white noise signal in the range of 30 Hz to 2000 Hz, to simulate ambient excitation as it will occur similarly in e.g. bridges under operational conditions. The shaker is installed upside down (Figure 2) and isolated against the changing temperature conditions. The data is recorded with 5 kHz. Three different system states were considered: the undamaged reference state, damage level 1 and damage level 2. For damage level 1 the beam was loaded in the midspan with 20 kN. In the second damage state the structure was loaded up to 28 kN, that was when obvious cracks occur in the middle of the beam. The load was removed after introducing the damage. The experimental temperatures were chosen to be -25°C, -5°C, 5°C, 25°C, and 40°C. The temperatures were measured at several points inside the chamber and at two locations inside the beam. In each system state and at all 5 testing temperatures the velocities were recorded, after reaching a constant temperature level in the beam. Further details on this experiment are described in [START_REF] Simon | Vibration-based structural health monitoring of a reinforced concrete beam subject to varying ambient temperatures using Bayesian methods[END_REF].

With an SSI-algorithm the eigenstructure of the system is identified. Due to high noise impact above 500 Hz, the data is sampled down to 400 Hz. Figure 3 shows the first three eigenfrequencies in dependence of the temperature, which can be clearly identified in every system state and at every temperature. The first and the second eigenfrequency are quite close. This is due to the construction of the shaker installation, which works as a tuned mass damper to the beam structure. Simulation results with a model of the corresponding structure with a mass damper have shown a good accordance to these first three eigenfrequencies.

Each frequency is plotted within its 5%-range to give an idea of the quantity of the temperature and the damage effect. The black full line marks the frequencies in the undamaged state. Effects due to temperature become more obvious for the second and the third frequency. While the damage level 2 leads to significant decrease of the eigenfrequencies of about 5%, the effect of damage level 1 does not become obvious in the first and the second eigenfrequency. The third frequency increases by about 3%. This was not expected and cannot be explained from the information gained from the experiments. However, this fact does not affect the usability of the data for the application study.

Damage detection of the concrete beam under changing temperatures

To illustrate the effect of the changing temperature to damage detection of the concrete beam, the proposed method based on interpolation will be compared to 1/ taking a fixed reference at one of the temperatures (without consideration of the temperature effect), and 2/ the averaging approach from [START_REF] Balmes | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF] that defines a fixed reference by mixing data from different reference temperatures.

Fixed reference

The tests are done firstly for the case where only data from one fixed reference temperature is considered in the reference state. The null space in (3) and the residual covariance in ( 5) are computed from a 5 minutes data sample at 5°C. Figure 4 shows the test values of the hypothesis test (4) at different testing temperatures in the undamaged state and at the two damage levels. For the setup of the testing Hankel matrices ℋ ̂ data samples of 100 seconds are used. From the test values in the reference state for the undamaged structure at 5°C (grey bars) a threshold is defined, which allows for false positive alarms of 1% maximum. It becomes obvious that changes in the temperature in the undamaged state (green bars) effect the test values in a similar way as the different damages do. This results in a high number of false alarms in the undamaged state, namely for all cases where the test temperature is different from the reference temperature. While all test values in the damaged cases also exceed the threshold, they cannot be distinguished from the test values in the reference state due to the temperature effect.

Averaging approach from [5]

In [START_REF] Balmes | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF] a temperature rejection approach has been proposed, which uses data from different temperatures in the reference state together. From this merged data set the averaged reference Hankel matrix and its left null space is computed, as well as the covariance of the resulting residual. In Figure 5 the test values from this method are presented, where 5 minute data samples from three reference temperatures, -25°C, 5°C and 40°C, are merged in one reference data set. The testing Hankel matrices ℋ ̂ are computed from a 100 s data sample. A threshold is defined from data in the reference state, thus of the undamaged structure at the three reference temperatures (grey bars), allowing for 1% false alarms.

The structure is correctly identified to be undamaged if the testing temperature corresponds to one of the used reference temperatures. Information about all reference states is contained in the merged reference data set and thus considered in the reference null space and the covariance matrix.

If the temperature differs from the reference temperatures, the changes in the structure due to temperature lead to higher test values even in the healthy states, exceeding the threshold. The structure is detected to be damaged. These changes due to new testing temperatures have quantitatively the same effect to the test values as damages do. While damages can be detected at all temperature levels, they cannot be distinguished from the healthy state at testing temperatures that are not contained in the reference data.

Proposed interpolation approach

In the interpolation approach, the information about the temperature depending behavior of the structure at testing temperatures different from the reference temperatures can qualitatively be carried by an adequate weighting function.

The test values in Figure 6 are computed with the interpolation approach using the parameter depending Hankel [START_REF] Simon | Vibration-based structural health monitoring of a reinforced concrete beam subject to varying ambient temperatures using Bayesian methods[END_REF] and covariance matrices (14) of the undamaged system at the three reference temperatures -25°C, 5°C, and 40°C. For the weighting of these matrices a Gaussian function with a variance of 5.5, centered at the testing temperature, is applied. Data samples of 5 minutes are used for the reference set up, and the testing Hankel matrices are computed from 100 s data samples. A threshold is computed in the reference state of the undamaged structure at the three reference temperatures (grey bars), allowing for 1% false alarms.

With this approach the undamaged state can be distinguished clearly from the damaged state for any testing temperature. There is no difference in the test values, whether data from the testing temperature was available for the reference set up or not. An adequate null space can be found in any case, and false alarms in the undamaged state are avoided. The damages can be detected reliably.

In general, the test values of the interpolation method are higher than in the other computation methods. This can be explained from the choice of the weighting of the local reference models. For any testing temperature, the Gaussian function with variance of 5.5 leads to a weighted mixture between all of the used reference models. This means that even at the reference temperatures, where data for a more precise reference null space exists, an interpolated reference null space is used, which is affected by the other reference null spaces. This leads to some bias and as a result to higher test values of the hypothesis test. However, in the damaged case this effect persists, and damages still can be detected reliably.

In the following, some conclusions are made for the choice of weighting function. It could be confirmed that the usage of bell-shaped functions centered at the testing temperature is adequate, if the temperature effects can be assumed to be approximately linear for each section. However, the choice of the variance of these weighting functions is quite important. If a function with a high variance is chosen, this leads to results very similar to the averaging approach in Figure 5. If the variance is very low, the test values are similar to those in Figure 4, when only one reference temperature is used. In the present case the interpolation worked well with a function with 𝜎 2 = 5.5 variance, but the effect of the width should be studied in further applications. In any case it is recommended not to use weighting functions that neglect nearby local models, even if a good local model is known from existing reference data. In this case the test values for testing temperatures different from the reference ones would be higher and could lead to false alarms.

CONCLUSIONS

In this paper the application of an interpolation approach is presented for the consideration of the temperature effect in stochastic subspace-based damage detection to experimental data. It could be shown that the interpolation approach is able to distinguish between damaged and undamaged states reliably, and false alarms can be avoided. Some suggestions are made regarding the choice of an adequate weighting function.

Further research should address the method's performance and the choice of weighting functions for more complex temperature effects, especially discontinuous behavior as shown in [START_REF] Peeters | One year monitoring of the Z24 bridge: environmental influences versus damage effects[END_REF]. Moreover, the effect of a more detailed consideration of the uncertainties of the reference null space as proposed in [START_REF] Viefhues | Asymptotic analysis of subspace-based data-driven residual for fault detection with uncertain reference[END_REF] should be investigated on the test performance of the interpolation method. In this context the question also arises about the minimum damage size that can be detected with an interpolated reference null space, and if the bias in the interpolated reference null space resulting from combining several local models affects the test sensitivity.

Figure 1 .

 1 Figure 1. Test setup of the concrete beam with locations of geophones, shaker, and loading for the damage introduction.

Figure 2 .

 2 Figure 2. Concrete beam in the climate chamber.

Figure 3 .

 3 Figure 3. 1 st to 3 rd eigenfrequencies at different system states and at different temperatures. The y-axis is scaled to about 5% of each frequency.

Figure 4 .

 4 Figure 4. Test values at different testing temperatures in the undamaged state and at two damage levels. Fixed reference null space and covariance matrix are computed with data from one data set recorded at 5°C.

Figure 5 .

 5 Figure 5. Test values at different testing temperatures in the undamaged state and at two damage levels. Averaged reference null space and covariance matrix are computed from a merged data set recorded at -25°C, 5°C, and 40°C.

Figure 6 .

 6 Figure 6. Test values at different testing temperatures in the undamaged state and at two damage levels. Interpolated reference null space and covariance matrix are computed from local matrices at -25°C, 5°C, and 40°C.

  1,2 , Michael Döhler 2 , Patrick Simon 1 , Ralf Herrmann 1 , Falk Hille 1 , Laurent Mevel 2 1 Federal Institute for Materials Research and Testing (BAM), Division 7.2 Buildings and Structures, Berlin, Germany