Direct structural analysis of domains defined by point clouds

TL;DR

This paper introduces a novel numerical analysis method for solids represented by point clouds, utilizing the Finite Cell Method to perform accurate, tessellation-free boundary analysis directly on point cloud data.

Contribution

It presents a new approach that enables direct structural analysis from point clouds without surface reconstruction, integrating high-order immersed boundary techniques.

Findings

Achieves accuracy comparable to traditional surface-based methods.

Demonstrates applicability on 3D historical structures.

Provides controlled error through point cloud density.

Abstract

This contribution presents a method that aims at the numerical analysis of solids represented by oriented point clouds. The proposed approach is based on the Finite Cell Method, a high-order immersed boundary technique that computes on a regular background grid of finite elements and requires only inside-outside information from the geometric model. It is shown that oriented point clouds provide sufficient information for these point-membership classifications. Further, we address a tessellation-free formulation of contour integrals that allows to apply Neumann boundary conditions on point clouds without having to recover the underlying surface. Two-dimensional linear elastic benchmark examples demonstrate that the method is able to provide the same accuracy as those computed with conventional, continuous surface descriptions, because the associated error can be controlled by the density of the cloud. Three-dimensional examples computed on point clouds of historical structures show how the method can be employed to establish seamless connections between digital shape measurement techniques and numerical analyses.

Equations

$$\delta W(\mathbf{u},\delta \mathbf{u})=\int {\text{displaylimits}}_{\Omega }\mathbf{\sigma }:\left({\nabla }_{\text{sym}}\delta \mathbf{u}\right)\text{d}V$$

$$\mathbf{u}$$

$$\mathbf{\sigma }=\alpha \mathbf{C}:\mathbf{\epsilon },$$

$$\alpha (\mathbf{x})=\{\begin{array}{ll}1& \forall \mathbf{x}\in {\Omega }_{phy}\\ 10^{-q}& \forall \mathbf{x}\in {\Omega }_{fict}.\end{array}$$

$$\mathbf{u}=\sum _i{\text{N}}_i{\mathbf{u}}_i\text{ ; }\delta \mathbf{u}=\sum _i{\text{N}}_i\delta {\mathbf{u}}_i.$$

$$\mathbf{Ku}=\mathbf{f},$$

$${\mathbf{k}}^e=\int {\text{displaylimits}}_{{\Omega }^{\text{e}}}{\left[\mathbf{L}{\mathbf{N}}^e\right]}^T\alpha (\mathbf{x})C\left[\mathbf{L}{\mathbf{N}}^e\right]\text{d}{\Omega }^{\text{e}},$$

$$\left({\mathbf{p}}_i-\mathbf{x}\right)\cdot {\mathbf{n}}_i\ge 0.$$

$$\int {\text{displaylimits}}_{\Gamma }f(\mathbf{x})d\Gamma =\int {\text{displaylimits}}_{\Omega }f(\mathbf{x})\delta (\mathbf{x})d\Omega ,$$

$$\delta (\mathbf{x})=\{\begin{array}{ll}\infty & \forall \mathbf{x}\in \Gamma \\ 0& \text{otherwise}.\end{array}$$

$$\delta (x)\approx {\delta }_{ϵ}(x)=\{\begin{array}{ll}\frac{1}{2ϵ}\left(1+\cos \left(\frac{\pi x}{ϵ}\right)\right)& \text{if}\left|x\right|\le ϵ,\\ 0& \text{otherwise},\end{array}$$

$${\delta }_{\Gamma }(\mathbf{x})={\delta }_{ϵ}(d_{\Gamma }(\mathbf{x})).$$

$$\int {\text{displaylimits}}_{{\Gamma }_{\text{N}}}\delta \mathbf{u}\cdot \mathbf{t}\text{d}A\approx \int {\text{displaylimits}}_{{\Omega }_{\cup }}{\delta }_{ϵ}\left(d_{S_{\Gamma }}\left(\mathbf{x}\right)\right)\left(\delta \mathbf{u}\cdot \mathbf{t}\right)\text{d}\Omega$$

$$e=\frac{\left|U_{\text{ref}}-U_{\text{num}}\right|}{U_{\text{ref}}},$$

Full text

Direct structural analysis of domains defined by point clouds

László Kudela

Chair for Computation in Engineering, Technical University of Munich,

Arcisstr. 21, 80333 München, Germany

Stefan Kollmannsberger

Chair for Computation in Engineering, Technical University of Munich,

Arcisstr. 21, 80333 München, Germany

Umut Almac

Istanbul Technical University, Faculty of Architecture, Turkey

Ernst Rank

Chair for Computation in Engineering, Technical University of Munich,

Arcisstr. 21, 80333 München, Germany

Institute for Advanced Study, Technical University of Munich,

Lichtenbergstr. 2a, 85748 Garching, Germany

Abstract

This contribution presents a method that aims at the numerical analysis of solids represented by oriented point clouds. The proposed approach is based on the Finite Cell Method, a high-order immersed boundary technique that computes on a regular background grid of finite elements and requires only inside-outside information from the geometric model. It is shown that oriented point clouds provide sufficient information for these point-membership classifications. Further, we address a tessellation-free formulation of contour integrals that allows to apply Neumann boundary conditions on point clouds without having to recover the underlying surface. Two-dimensional linear elastic benchmark examples demonstrate that the method is able to provide the same accuracy as those computed with conventional, continuous surface descriptions, because the associated error can be controlled by the density of the cloud. Three-dimensional examples computed on point clouds of historical structures show how the method can be employed to establish seamless connections between digital shape measurement techniques and numerical analyses.

Keywords: finite cell method, point clouds, image-based structural analysis

1 Introduction

It is well known in the computational mechanics community that transferring a CAD model into an analysis-suitable finite element mesh may account for as much as 80% of the entire analysis time [1]. Recent years’ research efforts aiming at circumventing this bottleneck have resulted in numerous alternative approaches, leading to significant progress in the quest of establishing seamless connections between geometric modeling and finite element analysis.

However, in some applications of the FEM, the geometries of interest are not directly available in form of CAD models. Typically, this situation arises in the context of biomechanical simulations where models are recorded by means of medical imaging techniques, such as CT scans. As standard finite elements require a CAD model to start with, these methods require special algorithms to recover a geometric model and eventually a finite element mesh from the imaging data – see e.g. [2] for a conceptual overview of these multi-step pipelines.

Volumetric imaging is not always the most feasible approach to record the shape of physical structures. It is especially large objects that do not allow for a cost-effective application of CT scanning. Nonetheless, it can be of importance to be able to compute the structural behavior of large objects – e.g. in the field of cultural heritage preservation, as there are often no digital CAD models available for historical structures. Moreover, even if there are schematic drawings, the shape of the object may differ from them, especially if the structure is exposed to damaging effects such as erosion, floods, earthquakes, or wars. In these cases, other shape measurement techniques need to be employed. The two most popular methods for this purpose are terrestrial laser scanning and close range photogrammetry-based reconstructions. Especially photogrammetry has gained a lot of attention recently, due to the inexpensiveness of the required equipment and because of the rapid development of the computational resources as well as the associated algorithms that allow for efficient, almost real-time reconstructions [3].

The methods of laser scanning and photogrammetry both reproduce the shape of the geometry of interest in the form of point clouds: a set of points representing the surface of the object. Such point clouds are not directly suited for numerical analysis. In order to transform the recorded data into an analysis-suitable model, it needs to pass through several stages, similar to the necessary procedure for models stemming from volumetric imaging.

Usually, these measurement-to-analysis procedures are characterized by the following main steps (Figure 1):

*Shape acquisition

*A 3D shape measurement technique is employed to capture the shape of the domain of interest, resulting in a point cloud representing the surface of the object. 2. 2.

*Surface reconstruction

*A geometric model is derived from the point cloud information using geometric segmentation and surface fitting methods. The resulting model is stored using standardized geometric representation techniques, such as STL, STEP, or IGES files. 3. 3.

*Mesh generation

*The CAD model from the previous step is discretized into a finite element mesh. 4. 4.

*Finite Element Analysis

*The mesh is handed over to a finite element solver together with the corresponding material properties and structural constraints.

Numerous applications implement the steps above – see e.g. [4, 5, 6, 7, 8] for examples in the preservation of historical structures, or [9] for an application in the context of biomechanical experiments.

Research in different fields of computational science and engineering has resulted in well-established approaches that allow to perform these steps one-by-one. Still, their deep integration into a seamless chain is not trivial, as it requires the interplay of various algorithms. Other than the problems inherent to the data transfer between different implementations, an even bigger challenge is posed by generating a finite element mesh from the geometric model reconstructed in the second step. The fine details recovered by modern surface reconstruction algorithms (e.g.[10, 11, 12]) are not necessarily the details that need to be carried over to a finite element mesh, where the process of refinement is usually governed by the physics of the problem rather than aesthetic aspects. While a geometric defeaturing step may be applied to remove physically uninteresting details, manipulating the geometry carries the danger of introducing flaws in the geometric model, resulting in an invalid, ”dirty” geometry that cannot be meshed directly [13, 14].

One method that aims to avoid the difficult task of mesh generation is the Finite Cell Method (FCM), introduced in [15]. The FCM is based on the combination of immersed boundary methods and high-order finite element basis functions used in p-FEM [16] or Isogeometric Analysis [1].

Instead of generating a boundary-conforming discretization, the FCM extends the physical domain of interest by a so-called fictitious domain in such a way that their union forms a simple bounding box that can be meshed easily. To stay consistent with respect to the original problem, the material parameters in the fictitious domain are penalized by a small factor $\alpha$. The introduction of $\alpha$ shifts the analysis effort from mesh generation to numerical integration. The most notable advantages of the FCM are the drastically reduced engineering efforts for preprocessing, the almost costless meshing, and the high accuracy and efficiency of the computation.

In its simplest implementation, the only information that the FCM needs from a geometric model is the inside-outside state – the question whether a given point in space lies in the physical or the fictitious part of the domain. There are numerous geometric representations that are suitable to provide such point membership tests and that have been shown to work well in combination with the FCM, ranging from simple shapes provided by constructive solid geometry [17] to models as complex as metal foams [18].

In this paper, we will demonstrate that a direct analysis of geometries described by oriented point clouds is possible. To this end, we will combine the finite cell method with geometries that are represented by oriented point clouds. The members of the point cloud and the vectors associated to them provide enough information for point membership tests, allowing for structural analyses of objects directly on their cloud representation. This way, the tedious tasks of recovering a geometric model and generating a boundary conforming mesh can be avoided, allowing for significant simplifications in the measurement-to-analysis pipeline.

2 The finite cell method combined with oriented point clouds

This section summarizes the basics of the finite cell method as well as the concept of performing inside-outside tests on oriented point clouds. The description is limited only to the necessary minimum for the context of this paper. For further details, refer to [19].

2.1 The finite cell method

Figure 2 illustrates the core idea of the FCM. The boundaries of the physical domain of interest $\Omega_{\text{phy}}$ are extended by a fictitious part $\Omega_{\text{fict}}$. Their union $\Omega_{\text{phy}}\cup\Omega_{\text{fict}}$ forms the embedding domain $\Omega_{\cup}$. As $\Omega_{\cup}$ possesses a simple, box-like geometry, it can easily be meshed into a structured grid of rectangular finite elements in 2D and cuboids in 3D.

Similar to standard finite elements, the FCM is derived from the principle of virtual work [20]:

[TABLE]

where $\boldsymbol{\sigma},\boldsymbol{b},\boldsymbol{u},\delta\boldsymbol{u}\text{ and }\nabla_{\text{sym}}$ denote the Cauchy stress tensor, the body forces, the displacement vector, the test function, and the symmetric part of the gradient, respectively. Prescribed displacements $\boldsymbol{u}_{0}$ are defined on the boundary $\Gamma_{\text{D}}$, while the traction vector $\boldsymbol{t}$ specifies the Neumann boundary conditions on $\Gamma_{\text{N}}$, such that $\Gamma_{\text{D}}\cap\Gamma_{\text{N}}=\varnothing$.

The stresses and strains are related through the constitutive tensor $\mathbf{C}$:

[TABLE]

where $\alpha$ is an indicator function defined as:

[TABLE]

To avoid ill-conditioning of the resulting equation system, the value of $q$ is chosen between $6$ and $12$ for practical applications.

Homogeneous Neumann boundary conditions are automatically satisfied. Nonhomogeneous Neumann boundary conditions can be realized by evaluating the contour integral over $\Gamma_{\text{N}}$ in Equation 1. However, in contrast to standard finite elements, Dirichlet constraints cannot be applied directly, as $\Gamma_{\text{D}}$ usually does not coincide with the element boundaries. Therefore, these boundary conditions are often formulated in the weak sense, e.g. using the penalty method or Nitsche’s method, e.g. [21].

The unknown quantities $\delta\mathbf{u}$ and $\mathbf{u}$ are discretized by a linear combination of shape functions $\text{N}_{i}$ with unknown coefficients $\mathbf{u}_{i}$:

[TABLE]

The finite cell method uses shape functions $\boldsymbol{N}_{i}$ of higher order. Popular choices are either the integrated Legendre polynomials known from p-FEM [19] or spline-based shape functions used in isogeometric analysis [1].

Following the standard Bubnov-Galerkin approach [22, 20], substituting Equation 4 into Equation 1 leads to the discrete finite cell representation

[TABLE]

where $\mathbf{K}$,$\mathbf{u}$,$\mathbf{f}$ denotes the stiffness matrix, the unknown displacement coefficients, and the load vector, respectively. The stiffness matrix results from a proper assembly of the element stiffness matrices:

[TABLE]

where $\mathbf{L}$ is the standard strain-displacement operator, $\mathbf{N}^{e}$ is the matrix of shape functions associated to the element. In the context of the FCM, the elements in the background mesh are often referred to as cells.

Due to the introduction of the scaling factor $\alpha(\boldsymbol{x})$ in Equation 6, the integrand becomes discontinuous. Standard quadrature schemes of the finite element method, such as the Gauss-Legendre rule, lose their precision in the presence of this discontinuity. In order to reduce the associated integration error, the above integral is usually evaluated by means of specially constructed quadrature schemes. The most popular method is based on a composed Gaussian quadrature rule combined with a recursive subdivision of the elements cut by the boundary of the physical domain.

In this process, every intersected element is subdivided into equal subcells until a pre-defined depth $k$ is reached. Quadrature points are then distributed on the domains of the leaf cells of this integration mesh. An example for the two-dimensional case is depicted in Figure 3.

In addition to spacetree-based schemes, numerous alternative approaches have been developed for the purpose of integrating through the discontinuous jump. These schemes reduce the integration error either by decomposing the integration domain into boundary-conforming integration cells [23, 24], or by modifying the quadrature weights, e.g. using moment-fitting equations [25].

The indicator function in Equation 3 has to be evaluted for every quadrature point. This requires the geometric model that represents $\Omega_{\text{phy}}$ to provide point-membership tests: given a quadrature point, does this point belong to $\Omega_{\text{phy}}$ or not?

Point membership tests can also be used to control the refinement process of the spacetree-based integration. To determine whether a cell is cut by the boundary of the physical domain, a set of test points is distributed inside the cell. Then, the inside-outside state of every test point is evaluated. The cell is definitely cut by the boundary of the physical domain if there is at least one pair of test points in the cell with changing inside-outside state. This idea is depicted in Figure 4.

The consequence of the above considerations is that, in the standard case, the only information that a geometric model needs to provide for the FCM is a robust point membership classification. Many geometric representations are able to answer such queries and have been successfully applied in combination with the FCM. Examples include voxel models from CT-scans [26], constructive solid geometries [17], boundary representations [24] and STL descriptions [27].

2.2 Inside-outside testing on point clouds

In point-cloud-based simulations, the domain $\Omega_{\text{phy}}$ is represented by a set of sample points $\mathbf{p}_{i}$ and their associated normal vectors $\mathbf{n}_{i}$. If no outliers are present, the set of pairs $S=\left\{\mathbf{p}_{i},\mathbf{n}_{i}\right\}$ represent a discrete sampling of the boundary $\partial\Omega_{\text{phy}}$ of the domain.

Each element in $S$ defines a hyperplane that separates the space in two half spaces: the open half-space $\Omega_{i}^{-}$ on the side of the hyperplane where the normal vector $\mathbf{n}_{i}$ points, and the closed half-space $\Omega_{i}^{+}$ on the other side. This concept is depicted in Figure 5.

For every $\mathbf{x}\in\Omega_{i}^{+}$, the following holds:

[TABLE]

Therefore, a simple way to estimate whether a quadrature point $\mathbf{q}$ lies inside or outside the domain is to find the $\mathbf{p}_{i}$ and the associated $\mathbf{n}_{i}$ in $S$ that lies closest to $\mathbf{q}$, and to evaluate the scalar product of Equation 7. While the heuristic nature of this approach may give the impression that it only works for simple shapes, our practical examples in Section 3.2 on more complex geometries show that the recovered indicator function is suitable to perform a finite cell analysis in these cases as well. The algorithm requires an efficient nearest neighbor query. In our examples, we use the k-d tree implementation of the C++ library nanoflann [28], while the clouds themselves are represented by the data structures of the Point Cloud Library [29]. The point membership classification method is summarized in Algorithm 1.

2.3 Point-based surfaces and Neumann boundary conditions

To apply boundary conditions in the weak sense, the contour integral in Equation 1 needs to be evaluated. For surface models, this is a relatively easy procedure, as they usually possess (or can be converted into) tessellations. Then, the integral over $\Gamma_{\text{N}}$ is computed as the sum of the integrals over the individual simplices in the tessellation.

However, for point-based geometries, no such tessellations exist. Although there are methods that are able to recover a triangulation from point cloud descriptions, their application would require to perform the same steps as the standard steps of the measurement-to-analysis pipeline in Section 1.

What is needed is an alternative formulation that allows for applying boundary conditions directly on point cloud-based surface representations. One possible solution to this challenge is to convert the contour integral into a domain integral by using the sifting property of the Dirac delta distribution [30]:

[TABLE]

with

[TABLE]

In numerical applications, the regularized variant of the Dirac delta distribution is employed. There are different choices available for the regularization, see e.g. [31, 32]. In our examples, we employ the following 1D formulation:

[TABLE]

where $\epsilon$ is a length scale parameter that controls the width of the regularization. Figure 6 depicts $\delta_{\epsilon}$ for different choices of $\epsilon$.

To extend the delta function to more dimensions, the distance function $d_{\Gamma}(\boldsymbol{x}):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is needed, which, for a given $\boldsymbol{x}\in\mathbb{R}^{n}$, returns the distance of $\boldsymbol{x}$ to the contour $\Gamma$. Then, the multi-dimensional regularized delta function, associated to the contour $\Gamma$, can be written as:

[TABLE]

In the point-cloud setting, $\Gamma$ is represented by the point set $S_{\Gamma}\subset S$. Instead of computing the distance to closest point $\boldsymbol{p}_{i}$, we compute the planar approximation of the n-neighborhood of $\boldsymbol{p}_{i}$ based on principal component analysis. Then, we evaluate the distance toward this approximant, as depicted in Figure 7.

Finally, the contour integral for the Neumann boundary condition in Equation 1 is formulated as:

[TABLE]

It is noted here that the transformation of boundary condition integrals into a domain integral using the delta function is an already existing concept in the context of FCM. Similar to point-cloud descriptions, phase-field models do not possess tessellations either. In this case, the approach outlined in this section can be applied to compute Neumann boundary conditions as well as Dirichlet boundary conditions formulated in the weak sense. For more details, refer to [33, 34].

2.4 Outliers

Sometimes, the cloud acquired by the shape measurement process carries outliers. As defined in [35], outliers are ,,observations that deviate so much from other observations as to arouse suspicion that it was generated by a different mechanism.” The detection and treatment of outliers – an aspect that lies beyond the scope of this article – is a well-studied problem in the literature, see e.g. the methods in [35, 36, 37, 38, 39, 40, 41]. These algorithms can be applied in a pre-processing step in order to recover a clean cloud where the majority of the outliers have been removed. The cleaning procedure is an essential step also in the traditional measurement-to-analysis pipeline, as most of the geometry recovery algorithms rely on a clean input cloud.

While the cleaning process is able to remove the bulk of the outliers, if isolated outliers remain in the cloud, they may have an effect on the point-membership classification process of Algorithm 1. To demonstrate this, Figure 8a shows the inside-outside state recovered from a point cloud consisting of 90 points and an outlier in the upper right quadrant of the domain. As can be seen in the figure, the presence of the outlier causes a region in the domain to be falsely detected as ,,inside”. To deal with isolated outliers, the point membership test of Algorithm 1 can be modified by testing in the n-neighborhood of the query point. In this process, instead of checking against a single closest point, the $n$ nearest points of $\boldsymbol{q}$ are found and the point membership with respect to each of them is computed. If $\boldsymbol{q}$ lies inside with respect to the majority of the points in the $n$ neighborhood, its membership is determined as inside, otherwise as outside. Following this idea, Figure 8b depicts the inside-outside state recovered using 3 nearest neighbor queries. The choice of $n$ for the number of nearest neighbor queries is a parameter that needs to be determined prior to the simulation. In the examples of Section 3.2, the parameter is chosen in the range $n=5...50$.

3 Numerical examples

This section demonstrates the proposed point-cloud-based FCM approach with numerical examples in two and three dimensions. First, we study an example with a known reference solution, where all the boundary conditions are aligned with the finite cell boundaries. This way, the modeling errors due to the approximate application of Neumann BC-s (Equation 12) can be ruled out. This is followed by an example where the performance of applying Neumann boundary conditions on non-conforming interfaces (as explained in Section 2.3) is tested. Finally, the point cloud-based FCM method is demonstrated on actual three-dimensional examples of historical structures.

3.1 2D studies

3.1.1 Perforated plate with circular hole

Figure 9a depicts a rectangular plate with a circular hole in the center. The plate is subjected to a constant traction along $\Gamma_{4}$, while symmetry boundary conditions are applied on $\Gamma_{1}$ and $\Gamma_{2}$. The Young’s modulus and the Poisson’s ratio of the material are $E=2.069\cdot 10^{5}\text{[MPa]}$ and $\nu=0.29$, respectively. Considering plane stress physics, the reference strain energy obtained by an overkill FEM analysis is $U_{\text{ref}}=0.7021812127$ [15]. For the FCM, the embedding domain is discretized into $2\times 2$ elements, and a value of $\alpha=10^{-12}$ is applied to scale the material parameters in the fictitious domain. The polynomial order of the shape functions is $p=12$. The continuous curve representing the circular hole is replaced by an oriented point cloud consisting of $n$ points, where $n$ is gradually increased in the range of $n=4...4096$. Refer to Figure 9b for an example of such a cloud.

The accuracy of the analysis is measured using the following error norm:

[TABLE]

where $U_{\text{num}}$ is the strain energy computed on the discrete point cloud representation of the circular hole.

The discretization of the circular interface into an oriented point cloud can be thought of as a replacement of the boundary by an $n$-sided polygon. This introduces a modeling error when integrating the term over $\Omega_{\text{phy}}$ in Equation 1, as the integration is not performed over an exact circle but rather over its polygonal approximation. Therefore, the error is expected to converge at the same rate as the area of an $n$-sided polygon converges towards the area of its inscribed circle, i.e. quadratically.

This expectation is confirmed by the error plots on Figure 10. It depicts the evolution of the error for increasing cloud densities, for different maximum levels of quadtree subdivision ${k=\left\{4,6,8\right\}}$, see Section 2.1. Initially, the polygonal approximation dominates the error, and the curves show quadratic convergence. Eventually, depending on the value of $k$, the convergence curves level off into a plateau. In the plateau region, the error due to the quadtree-based integration overtakes the polygonal discretization error: even though the interface is modeled by higher resolutions, the integration scheme is not able to resolve this increase in geometric accuracy.

3.1.2 Perforated plate with elliptical hole under internal pressure

In this example, the circular interface from the previous section is replaced by an elliptical curve (Figure 11). To investigate the effects of applying Neumann boundary conditions as described in Section 2.3, the elliptical hole is discretized into an oriented point cloud along which a constant internal pressure of $1\text{[MPa]}$ is applied. The reference value for the strain energy computed by p-FEM is $U=44.28375067893$. The domain is discretized into a regular grid of $6\times 6$ finite cells, where the polynomial order of the shape functions varies in the range $p=1...10$. The regularization parameter of the Dirac delta function in Equation 10 is chosen to be $\epsilon=0.0625$. The diffuse region obtained in this way (see Figure 7) is integrated using a composed Gaussian quadrature of order $10$ combined with a quadtree-based subdivision of maximum level $k=8$. While this integration depth seems to be excessively large for practical applications, it should be noted that the size of the surfaces for which weak boundary conditions need to be applied is usually significantly smaller than the overall sizes of the geometries of interest. Therefore, applying the boundary conditions in this manner does not lead to a significant performance penalty.

To rule out the errors associated to the discontinuity in the indicator function $\alpha(\boldsymbol{x})$, the stiffness matrix is integrated using the highly precise integration technique described in [23]. This way, the error due to the approximate application of the Neumann boundary condition is examined exclusively.

Figure 12 depicts the convergence of the error in the energy norm for different point cloud densities. As the figure shows, the expected exponential rate of convergence can be attained in the pre-asymptotic range. However, similar to the study conducted in Section 3.1.1, the convergence curves level off into a plateau, depending on the density of the point cloud. As higher cloud densities are able to represent the underlying elliptical contour with higher accuracy, the level-off location shifts towards lower errors for an increasing number of points in the cloud. It is noted here, however, that even a relatively low density (1000 points) is able to produce an error in the range of $1\%$, which is sufficient for most engineering applications.

3.2 3D examples

In the following, the proposed point-cloud-based FCM approach is demonstrated on three-dimensional structures represented by oriented point clouds.

3.2.1 Athlete

Figure 13a shows a statue from the museum ,,Glyptothek” located in Munich, Germany. Images of the object were taken from 36 different angles, using a cell phone camera. These input images were processed using the popular structure-from-motion toolbox VisualSFM [42, 43] and the multi-view reconstruction algorithm of [44], as demonstrated in Figure 13b. The resulting point cloud is depicted in Figure 13c.

The point cloud was embedded in a regular mesh of $325$ finite cells with polynomial order $p=5$, as shown in Figure 14a. Linear elastic material behavior is assumed, and the structure is loaded under its self-weight. The scaling factor $\alpha$ for the FCM was chosen as $10^{-6}$. Homogenous Dirichlet boundary conditions were applied on the bottom faces of the finite cell mesh, in order to rigidly fix the statue to the ground.

The resulting stress field is depicted in Figure 14b, while Figure 14c shows a detailed view of a cross-section of the left foot. As seen in the Figure, the peak stress occurs at the ,,ankle” region, which possesses the smallest cross-section over the entire structure. This phenomenon is in good accordance with other observations from the study of the structural behavior of stone statues: numerical computations conducted on ,,David” from Michelangelo showed a similar stress distribution, with the peak occurring in the ankle area [5]. Interestingly, there are also other areas of increased stresses – such as the area around the neck as well as areas in the upper thigh. Other intuitive candidates for high stresses, such as the left arms and shoulders can be disregarded.

While this example illustrates the main steps of the proposed point-cloud-based FCM pipeline, without further knowledge about the material parameters, the computed results merely allow for a qualitative assessment of the stress distribution in the statue, under the assumption that no internal cavities are present and that the material is uniform. The next example addresses this question on a structure fow which the material properties are known.

3.2.2 The cistern of the Hagia Thekla Basilica in Turkey

The archaeological site at Hagia Thekla (Meryemlik, today Turkey) was a major pilgrimage site in late antiquity [45]. The site features numerous above-ground structures.

The cistern of the Thekla Basilica is part of the water storage and distribution system of the main church of the site and its sacred surrounding area, which is enclosed by walls. It has a rectangular plan measuring approximately $12\times 14.6$ meters in the interior. The interior space is divided into three aisles by two rows of columns (Figures 15b and 15a). The columns in each row are connected by arches. Three barrel vaults cover the interior running in the north-south direction.

The columns supporting the upper structure are made of a pink calcareous stone and originally had a diameter of approximately $45$ cm. The columns have double capitals made of limestone. It is not possible to make observations about the condition of the column bases and the floor, due to the thick layer of earth accumulated inside the cistern over centuries. The outer walls were built with a multi-leaf masonry construction system. The outer facing of the walls are made of big limestone blocks, while the inner faces consist of brick and mortar. As seen in Figure 15a, the cross-sections of the columns have decreased remarkably. The exterior surfaces are flaking due to physicochemical effects; the erosion continues. In addition to surface erosion with a non-uniform pattern, there are deep cavities on the surfaces of the columns. One of the columns (Column 3) has already collapsed and was replaced by a concrete column in the 1960’s.

It is difficult to record the shape of the decayed column surfaces and cavities relying on manual measurement procedures. Therefore, a high definition surveying scanner was employed to document these elements. During the field campaign, the instrument was set up at a number of positions around each column at a distance of a few meters. Thus, a maximum point density of approx. $5$ mm was ensured to represent the highly decayed columns. More details on the measurement process can be found in [7]. Figure 16a depicts the measured point cloud, consisting of $10^{7}$ points.

The most vulnerable elements of the structure are the columns. As stress concentrations are expected at the cavities on the surfaces of the columns, a reduction of the discretization error by a refinement of the computational grid is needed. For reasons of efficiency, it is important to refine the grid only around the columns, where the stress field is expected to change rapidly. For the FEM and the FCM, such local refinement techniques have been well-studied recently. In our applications, we employ the multi-level $hp$-adaptivity technique of [46]. In the refinement procedure, those cells that are intersected by the points representing column 2 and 4 are recursively subdivided into eight equal subcells, until a subdivision depth of 4 is reached. A cross-sectional view of the refined mesh is depicted in Figure 17.

The material properties were defined to be linear elastic and isotropic, with an elastic modulus and a Poisson’s ratio of $E=2\cdot 10^{4}\text{MPa}$ and $\nu=0.2$, respectively. The specific gravity of the material was set to $27\text{ kN}/m^{3}$. In the fictitious domain, the material was given a stiffness of $2\cdot 10^{-4}\text{ MPa}$. The foundation of the structure was rigidly fixed to the ground. The maximum principal stress distribution computed by the FCM is depicted in Figure 18.

As expected, the highest compressive stresses occur in the columns. The stress values are in the range of $2...6\text{ MPa}$, while the peak value occurs at the connection between the columns and their capitals. This is in good agreement with the values computed in [7], following the traditional measurement-to-analysis procedure, i.e. reconstructing the surface, meshing the volume and performing a finite element analysis. A comparison of principal stresses along column 2 is given in Figure 19. As seen in the Figure, the stress patterns show high resemblance. However, the results computed by our technique neither required the recovery of a geometric surface model, nor the generation of a boundary-conforming finite element mesh, allowing for performing the analysis directly on the point cloud.

As the only requirement toward the geometric model is that it needs to be able to provide point-membership information, the finite cell method opens up a convenient way to investigate the effects of geometric changes, such as the removal of a specific column. Within the FCM, it is easily possible to integrate different geometric models following the idea of Constructive Solid Geometry (CSG) modeling [17]. In CSG, 3D objects can be created by combining geometric primitives into a tree structure, using boolean operations. To determine if a point lies within the model or not, the tree structure is traversed from the root towards the leaves, combining the inside-outside state according to the boolean operations at each level. Following this idea, we investigate the effect of removing Column 3 from $\Omega_{\text{phy}}$: using boolean difference, the bounding box of the points representing the column is subtracted from the point cloud, as shown in Figure 20.

The principal stress distribution in this scenario is depicted in Figure 21. As shown on the figure, the removal of the column causes a redistribution of the loads onto the neighboring supports, leading to an increase in the principal stresses in columns 2 and 4. The redistribution phenomenon can be demonstrated by comparing the principal stress trajectories in the configurations with and without the column, as depicted in Figure 22. Due to the removal of the column, a new ,,arc” forms between the two neighboring columns. This arc is in compression and carries the redistributed load. The newly formed stress state tells us why no structural failure occurred when column 3 collapsed: the structure was able to hold becuase stone is able to carry substantially higher loads in compression than in tension, and because the arc of principal stresses is predominantly in compression.

4 Conclusions and outlook

This contribution presented a method aiming at the numerical analysis of objects represented by oriented point clouds. Instead of relying on boundary-conforming computational meshes, the approach uses the finite cell method, which computes on a rectangular background grid and only requires inside-outside information from the geometric model. It was shown that oriented point clouds provide sufficient information to perform these point membership tests. This allows to avoid the difficult steps of geometry recovery and mesh generation in the usual measurement-to-analysis pipeline. The capabilities of the proposed method were presented on structures of engineering relevance, recorded by means of two shape measurement approaches: image-based reconstructions using photos taken with a hand-held device, and laser scanning.

It is noted here that image-based algorithms do not exclusively work with pictures stemming from hand-held devices. Recent developments in the technology of UAV-s have made it possible to take inexpensive but high-quality pictures of objects of virtually any size. In the field of cultural heritage preservation, this allows for an especially convenient way to assess the structural health of historical buildings. As an example, Figure 23 demonstrates a stress analysis conducted on the ruin of a medieval tower located at the ,,Hocheppan Castle” in South Tyrol. The results were computed following the point cloud-based FCM analysis pipeline, starting from 120 pictures taken of the ruin using a DJI Phantom 4 drone.

While the presented approach is a first step towards establishing seamless connections between shape measurement techniques and numerical analyses, future research should address further questions concerning at least the following three aspects:

The treatment of Dirichlet boundary conditions on non-conforming surfaces. For continuous surface descriptions, this can be resolved by applying these boundary conditions in the weak sense, using e.g. the penalty method or Nitsche’s method – allowing them to be formulated in terms of contour integrals that can be transformed into domain integrals, following the ideas presented in Section 2.3. 2. 2.

Especially for historical masonries, the anisotropic behavior introduced by the interaction between the mortar and the building blocks needs to be modeled. 3. 3.

Cracks, if present and visible in the original structure, should be incorporated in the corresponding numerical model.