Interrogation of spline surfaces with application to isogeometric design and analysis of lattice-skin structures

by   Xiao Xiao, et al.
University of Cambridge

A novel surface interrogation technique is proposed to compute the intersection of curves with spline surfaces in isogeometric analysis. The intersection points are determined in one-shot without resorting to a Newton-Raphson iteration or successive refinement. Surface-curve intersection requires usually the solution of a system of nonlinear equations. It is assumed that the surface is given in form of a spline, such as a NURBS, T-spline or Catmull-Clark subdivision surface, and is convertible into a collection of Bézier patches. First, a hierarchical bounding volume tree is used to efficiently identify the Bézier patches with a convex-hull intersecting the convex-hull of a given curve segment. For ease of implementation convex-hulls are approximated with k-dops (discrete orientation polytopes). Subsequently, the intersections of the identified Bézier patches with the curve segment are determined with a matrix-based implicit representation leading to the computation of a sequence of small singular value decompositions (SVDs). As an application of the developed interrogation technique the isogeometric design and analysis of lattice-skin structures is investigated. Current additive manufacturing technologies make it possible to produce up to metre size parts with designed geometric features reaching down to submillimetre scale. The skin is a spline surface that is usually created in a computer-aided design (CAD) system and the periodic lattice to be fitted consists of unit cells, each containing a small number of struts. The lattice-skin structure is generated by projecting selected lattice nodes onto the surface after determining the intersection of unit cell edges with the surface. For mechanical analysis, the skin is modelled as a Kirchhoff-Love thin-shell and the lattice as a pin-jointed truss. The two types of structures are coupled with a standard Lagrange multiplier approach.



There are no comments yet.


page 1

page 2

page 3

page 4


Concurrent infill topology and shape optimisation of lattice-skin structures

Lattice-skin structures composed of a thin-shell skin and a lattice infi...

Toric degenerations of Bezier patches

The control polygon of a Bezier curve is well-defined and has geometric ...

Interpolation of a spline developable surface between a curve and two rulings

In this paper we address the problem of interpolating a spline developab...


BSplines are one of the most promising curves in computer graphics. They...

Multi-Axis Support-Free Printing of Freeform Parts with Lattice Infill Structures

In additive manufacturing, infill structures are commonly used to reduce...

Evolution of natural patterns from random fields

In the article a transition from pattern evolution equation of reaction-...

A non-iterative method for robustly computing the intersections between a line and a curve or surface

The need to compute the intersections between a line and a high-order cu...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Almost all geometry models in isogeometric analysis are based on parametric rational spline curves and spline surfaces, such as NURBS, T-splines or subdivision surfaces Hughes et al. (2005); Bazilevs et al. (2010); Cirak et al. (2000); Wei et al. (2015); Zhang et al. (2018). The process of interacting with and exploring such geometries is commonly referred to as shape interrogation in computer-aided design (CAD) literature. Shape interrogation of curves and surfaces includes, for instance, the plotting of their differential geometric properties, like iso-curvature lines or geodesics, and the computing of their intersections with rays and planes as in ray tracing and contouring. Due to its importance in CAD there is copious literature and a number of ingenious techniques on shape interrogation, see e.g. the monographs Farin et al. (2002); Patrikalakis and Maekawa (2009) and the references therein. All shape interrogation involves usually, at some level, the computation of the intersection between curves, surfaces with curves, and surfaces with surfaces. The focus of this paper is on curve-surface intersection, but the presented techniques can be extended to curve-curve and surface-surface intersections. Evidently, interrogation of spline surfaces is critical in a number of isogeometric analysis applications as well. Robust intersection computations are essential, for instance, for enforcing non-penetration constraints between spline curves or surfaces Temizer et al. (2011); De Lorenzis et al. (2011), integrating over cut-cells Rüberg and Cirak (2012); Schillinger et al. (2012); Upreti and Subbarayan (2017), coupling trimmed shell patches Marussig and Hughes (2017); Breitenberger et al. (2015); Guo et al. (2018), or, as in the present paper, for coupling shells with lattices. The number of finite elements, i.e. spline patches, in analysis is usually significantly larger than the number of patches in usual CAD models so that the efficiency of the intersection computations becomes paramount in isogeometric analysis.

The intersection between a spline surface and a curve leads to a nonlinear root-finding problem, which is easy to formulate but extremely hard to solve robustly. In particular, multiple intersections and curves that intersect the surface tangentially cause difficulties to most intersection algorithms. The prevalent techniques to determine intersections include triangulation Piegl and Tiller (1998), Newton-Raphson iteration Toth (1985), subdivision Houghton et al. (1985); Sederberg and Parry (1986); Sabin (2002), implicitisation, to be used in this paper, or a combination thereof; see also the recent review Marussig and Hughes (2017) with a focus on trimming and isogeometric analysis. In triangulation based techniques the surface is approximated with a sufficiently fine triangle mesh and the curve with a polygon. The intersection between the triangle mesh and polygon is straightforward to compute and can serve as an approximation to the true intersection point. Subdivision based methods rely on the refinability and convex-hull properties of spline surfaces. The convex-hull property guarantees that a spline surface lies always within the convex-hull of its control points, and the refinability property enables to create a hierarchy of nested patches with decreasing convex-hull size. Hence, the intersection problem can be reformulated as the intersection between two bounding volume hierarchies for the convex-hulls of the spline patches and the curve. To simplify the intersection computation between convex-hulls it is expedient to approximate them with axis-aligned bounding boxes, spheres or, as in the paper, with k-dops (discrete orientation polytopesLin and Gottschalk (1998); Klosowski et al. (1998). In both triangulation and subdivision methods the smallest element or bounding volume size has to be chosen similar to the required accuracy for the intersection computation. Hence, the computing time and memory needed depend on the required accuracy. In contrast, in Newton-Raphson iteration the patches can be relatively large and convergence is quadratic sufficiently close to the intersection points. The principal difficulties encountered with Newton-Raphson iteration, mainly divergence for not carefully selected starting points, are the same ones as in path-following in nonlinear static analysis Ramm (1981); Wriggers (2008).

In the present paper the intersections between spline surfaces and curves are computed with an implicitisation technique proposed by Busé et al. Busé (2014); Shen et al. (2016). Prior to computing the intersections, the surface patches close to a curve are identified using a hierarchical k-dop bounding volume tree Klosowski et al. (1998)

. As our numerical experiments indicate a conventional axis-aligned bounding volume tree can be significantly slower than a k-dop bounding volume tree depending on the orientation of the surface and the curve. Subsequently, the intersections of the identified patches with the curve are determined by computing the singular values and null vectors of an implicitisation matrix, without resorting to a Newton-Raphson iteration or subdivision refinements. In implicit form a parametric surface,

with and is represented as the zero iso-contour of a scalar function, i.e. . Implicitisation naturally leads to algebraic geometry, which informally deals with systems of implicit polynomial equations of the form with and the geometric objects defined by them. For instance, a planar parametric curve, and , can be implicitised with resultants from algebraic geometry (see Hoffmann (1989); Sederberg and Zheng (2002); Cox et al. (2007); Sederberg (2012) for an introduction to resultants). A resultant is a polynomial expression for and , which is if and only if both polynomials satisfy and ; or in other terms only for points on the parametric curve. For a planar parametric curve the resultant can be computed as the determinant of the classical Sylvester or Bézout matrices. The moving lines introduced by Sederberg et al. Sederberg et al. (1994); Sederberg and Chen (1995) provide a related and more intuitive geometric approach to implicitisation. According to Sederberg et al. a planar parametric curve can be described as the intersection of two lines and that move with the parameter . The extension of the moving lines idea to surfaces is known to be challenging due to the presence of so-called base points and the impossibility of describing a surface with exactly three moving planes. The matrix-based implicitisation introduced by Busé Busé (2014) can be understood as the generalisation of the moving lines idea to surfaces. Instead with a closed-form implicit equation the surface is represented using the singular values of an implicitisation matrix. There is a drop in the rank of the implicitisation matrix when evaluated on the surface. The singular values of the rectangular implicitisation matrix are determined with a singular value decomposition (SVD), see e.g. Strang (2016)

. In intersection computations the implicitisation matrix is used to formulate a generalised eigenvalue problem for directly computing the spatial and parametric coordinates of the intersection point. The relatively small SVD and eigenvalue problems can be solved with robust and efficient solvers available in most modern numerical linear algebra packages. For completeness, there are a few other techniques for computing intersections 

Manocha and Canny (1991); Dokken (2001) which also use a SVD but are different from the implemented approach.

(a) Lattice structure with BCC unit cells
(b) Finite element computation
Figure 1: Stanford bunny with an internal lattice and its deformation under applied loading. The intersection points of the lattice with the surface are computed using the implicit matrix representation after identifying possible intersections with a k-dop bounding volume tree.

As an application of the developed shape interrogation technique the isogeometric design and analysis of lattice-skin structures is considered, see Figure 1. The advances in additive manufacturing, or 3d printing, make it possible to produce up to metre size macroscopic lattice structures with designed geometric features reaching down to submillimetre scale Gibson et al. (2010). Structures with lattice-like components can achieve significantly higher stiffness-to-weight and strength-to-weight ratios than ordinary structures with comparable weight Gibson and Ashby (1999); Fleck et al. (2010). To achieve this the struts within each unit cell must primarily be subjected to stretching (rather than bending). One of the most expedient application areas of lattice-like cellular solids is the lightweight design of structures where they are used to provide stiffness to plates and freeform shell skins Wang et al. (2013). It is usually not feasible to model the entire lattice-skin structure in a CAD system due to the required very large number of geometric entities and poor scalability of CAD systems Gibson et al. (2010), although the size of lattices that can be represented is continuously increasing. In practice, this limitation is overcome by designing only the skin in the CAD system. After the skin design is completed it is triangulated and the lattice is added outside the CAD system, see e.g. Weeger et al. (2017)

. The concurrent use of several representations leads to the usual proliferation of geometry models for representing the same part and the attendant interoperability and accuracy issues. The situation is further complicated when the generated design is to be analysed with finite elements. The interrogation technique presented in this paper allows the combination of any periodic lattice with a spline surface without the need for triangulation. After computing the intersection of the lattice with the spline surface the closest joints are projected onto the spline surface. The periodic lattice is modelled as a truss structure with pinned joints which do not transfer moments. This approximation is sufficient for lattices which are stretch dominated, see e.g. 

Deshpande and Fleck (2001). The shell skin is modelled as a Kirchhoff-Love thin-shell and is discretised with spline basis functions Cirak et al. (2000); Cirak and Long (2011); Zhang et al. (2018). The truss and Kirchhoff-Love thin-shell finite elements are combined with the standard Lagrange multiplier approach, see e.g. Guidault and Belytschko (2009). Although not pursued further here, after the design and analysis of a part is completed it can be triangulated for slicing for 3d printing.

The outline of this paper is as follows. In Section 2 the developed interrogation approach is introduced. After a cursory review of Bézier patches, first the proximity search between Bézier patches and curve segments using the hierarchical k-dop bounding volume tree is presented. Subsequently, the matrix-based implicitisation of Bézier patches and its use for computing their intersection points with curve segments are discussed. In Section 3 briefly the need and application of robust interrogation techniques in the isogeometric design and analysis of lattice-skin structures are illustrated. The efficiency and accuracy of the developed approach are studied in Section 4 with detailed timing studies for generating lattice-skin structure geometries and finite element analysis of simply-supported sandwich plates. In addition, the design of an optimally stiff sandwich cap with a lattice core is studied, which demonstrates that there is an optimal material distribution between the lattice and the skin. In Appendix, amongst others, an illustrative example for computing the intersection between two lines using the implicit matrix representation approach can be found.

2 Geometry representation and interrogation

2.1 Spline surfaces

A surface is assumed to be given in form of a spline, such as a NURBS, T-spline or Catmull-Clark subdivision surface. Each of these splines can be converted into a set of rational tensor-product Bézier patches, see e.g. 

Piegl and Tiller (1997); Borden et al. (2011); Stam (1998). In homogeneous coordinates each Bézier patch is given by


with the coordinates

the bivariate Bézier basis functions , the parametric coordinates , the prescribed control point coordinates and the scalar weights . The Bézier control points are obtained from the spline control points by means of a projection, referred to as a Bézier projection in case of tensor-product b-splines Borden et al. (2011). The basis functions and the control points are labelled with the multi-index . The tensor-product Bézier basis functions of bi-degree are defined by


with the Bernstein polynomials


where and . The Bernstein polynomials are nonnegative and sum up to one. In applications the same degree is usually used in both directions and in (1) the weights are except when quadric surfaces, like cylinders or spheres, are represented.

Although it would be feasible to re-express a spline surface using any polynomial basis, the Bernstein polynomials and the resulting Bézier patches have the advantage of numerical stability and convex-hull property. The convex-hull property guarantees that the Bézier surface lies always in the convex hull of its control points (when, as usual, Farin (2002). This is a direct consequence of the nonnegativity and partition of unity property of the Bézier basis functions. As will be discussed in Section 2.2, the intersection of convex-hulls can be conveniently used as a preliminary check to identify possible intersections.

It is worth noting that in case of Catmull-Clark subdivision the nested nature of spline rings around extraordinary vertices, with other than four connected faces, gives rise to rings of Bézier patches of decreasing size Stam (1998); Zhang et al. (2018). The size of the Bézier patches tends to zero when the extraordinary vertex is reached.

2.2 Surface-curve intersection

The computation of the intersection point between a Bézier patch and a given curve is an important high-level primitive task in the design of lattice-skin structures. For both design and analysis purposes the coordinates of the intersection point and its parametric coordinates on the surface and the curve are required. This leads to a system of nonlinear equations which are hard to solve robustly and efficiently. Although it is possible to check the intersection of a curve with each Bézier patch, this has evidently a complexity of , where is the number of patches, and becomes quickly inefficient for a large number of curves. Therefore, a hierarchical bounding volume tree and k-dops (discrete orientation polytopes) are used to first identify all the Bézier patches that may be intersected by the curve. This reduces the complexity, under some mild assumptions, to . Subsequently, the intersection points of the identified Bézier patches with the curve are computed one by one and valid intersections are determined. The intersections are computed with a matrix-based implicit representation without resorting to a Newton-Raphson iteration.

2.2.1 Bounding volume trees and k-dops for proximity search

Bézier patches that may be intersected by a given curve can be efficiently identified with a hierarchical bounding volume tree. Due to the convex-hull property of Bézier patches it is expedient to base such a tree on the convex hull of patch vertices. For ease of implementation and efficiency, however, convex-hulls can be approximated with k-dops Sabin (1992); Klosowski et al. (1998), see Figure 2. K-dops are the generalisation of axis-aligned bounding boxes and their bounding surfaces come from a set of fixed orientations ( in 3D). Axis-aligned bounding boxes with usually provide only a loose fit to convex-hulls and give many false positives. For a relatively high number of orientations k-dops can provide a much tighter fit and their intersections are easy to compute.

(a) Quadratic Bézier patch
(b) Convex hull of the vertices
(c) 8-dop of the control polygon
Figure 2: Bézier patches lie entirely within the convex hull and the k-dop formed by their vertices. Two-dimensional illustration with all control points .

In the implemented k-dop bounding volume tree each node represents the k-dop of either one single or a set of neighbouring patches, see Figure 3. The tree is constructed by recursively splitting the patches into eight subsets starting from the entire surface. In each step the splitting is based on the coordinates of the patch centroids (computed as the mean of vertex coordinates). To obtain a balanced tree all the nodes on the same refinement level have to contain approximately the same number of patches. After the patches have been assigned to tree nodes, it is straightforward to determine the k-dop of a set of patches either from their individual k-dops or from the union of their vertices.

(a) Root
(b) Level 1
(c) Level 2
Figure 3: An 8-dop bounding volume tree is constructed by successively splitting the Bézier patches into subsets. Two-dimensional illustration for a control polygon (in dark blue) consisting of four linear Bézier segments.

Indeed, for intersection checks there is no need to explicitly construct any k-dops. To determine whether two k-dops intersect it is sufficient to consider the projections of their control vertex coordinates along the prescribed directions with . As an example, in Figure 4 the intersection of two k-dops belonging to a linear and a cubic Bézier curve is illustrated. The support heights and of the cubic Bézier curve with the vertex coordinates in the direction are defined as


The support heights and for the linear Bézier curve are determined in a similar way. The two k-dops intersect only when


and they do not intersect when


Hence, in the k-dop bounding volume tree data structure it is sufficient to store only the support heights. The support heights for each tree node can be efficiently computed by starting from the leaf nodes and determining the support heights of parent nodes from children nodes according to

(a) Support heights in
(b) Support heights in
Figure 4: Intersection detection between a linear Bézier curve and a cubic Bézier curve. There is an intersection in direction and no intersection in direction.

2.2.2 Implicit matrix representation for intersection computation

Next, the matrix-based implicitisation of a single Bézier patch of bi-degree as defined in (1) is developed. To this end, consider the auxiliary vector of polynomials


with Bézier basis functions , coefficients to be yet determined and the multi-index . According to Busé (2014), the bi-degree of the Bézier basis can be chosen with or higher; alternatively or higher is also possible. The bi-degree of has to be sufficiently high in order to be able to discover later all the intersection points. The two vectors and , with each having four components, are required to be orthogonal


that is,


Evidently, the product of two Bézier basis functions and can again be expressed as a Bézier basis function, see Appendix 6.1. Hence, (10) can be written with a new Bézier basis of bi-degree


where are the coefficients of a matrix . Note that there are basis functions in the basis with bi-degree , basis functions in the basis with bi-degree and each has four coefficients. The matrix contains the control point coordinates, associated weights and the basis transformation coefficients resulting from the Bézier basis multiplication.

The orthogonality constraint (9) implies equations for determining unknown components of ’s so that the matrix must be rank deficient. The non-trivial ’s satisfying (11) must lie in the null space of , which is straightforward to compute with a singular value decomposition (SVD) Strang (2016). The singular value decomposition of reads


where is a matrix of left singular vectors, is the diagonal matrix of singular values and is a matrix of right singular vectors. The singular values in are usually sorted in descending order starting from the top. The right null vectors are the columns of which correspond to zero diagonal entries in . For a specific Bézier patch the number of right null vectors of depends on its control vertex coordinates  and the corresponding weights .

Denoting the right null vectors of with and introducing them into the auxiliary vector of polynomials (8) yields


To advance an intuitive geometric interpretation to implicitisation consider, for a fixed parametric coordinate , the set of planes defined by the vectors  with the implicit equations


Due to the orthogonality condition (9) the surface point with the homogeneous coordinate and the surface coordinate

satisfies all plane equations and, hence, is present on all planes, i.e. . For a point to be present on two or more planes it must lie on their intersection. Note that the intersection of two planes defines a line and the intersection of a plane with a line defines a point. These observations motivate the definition of the surface  as the intersection of a set of planes. Evidently, the planes must move, or change their inclination, while the parametric coordinate is varied. Hence, each for a fixed represents a family of planes and is traditionally referred to as a moving plane or a pencil of planes Sederberg et al. (1994). Although three planes are sufficient to define a point, according to Busé (2014) more than three planes must be used, in general, to describe an entire patch. As an example in Figure 5 the description of a Bézier patch with four moving planes is illustrated. For the sake of clarity in Figure (b)b only three of the planes are depicted.

(a) scale=1.0
(b) scale=1.0
Figure 5: A Bézier patch and its four moving planes at . For the sake of clarity, in (a) instead of the moving planes their pairwise intersection lines are shown and in (b) only three of the four moving planes are shown.

Instead of the geometrically motivated moving plane interpretation it is possible to advance an equivalent matrix-based linear algebra explanation. To this end, the plane equations (14) are, after introducing (8), rewritten as


where constitutes the -th column of a matrix , that is,


The number of columns of  is equal to the number of planes or null vectors , and its number of rows is . Again, due to the orthogonality condition (9) for a point with the homogeneous coordinate  and the surface coordinate  the set of equations (15) are satisfied for all . This implies that there has to be a change in the rank of  when evaluated at . Furthermore, the minimum degree suggested for the basis  ensures that  is full rank except on the surface Busé (2014). This motivates the use of the change of rank of to define the surface , which gives rise to the notion of matrix representation of a spline surface.

The matrix allows the formulation of the interrogation of the Bézier patch  as a linear algebra problem. As an example the intersection of a patch with a line is discussed next. For the intersection of a patch with a higher-order curve see Appendix 6.3. The line is assumed to be given in parametric form with


where are two prescribed vectors and is the parametric coordinate. It is assumed that the patch  intersects the line at the parameter value(s) , which can be several. To determine the line (17) is introduced into (16) yielding


where the matrices and contain the components of  independent and linear in , respectively. As discussed, at the intersection points  must be rank deficient or (15) must be satisfied. Hence, the generalised eigenvalue problem


is considered for computing all the intersection points 

of the line with the patch, where only the real eigenvalues are relevant. Here, the eigenvectors 

are of no significance for the intersection computations.

The matrix , and hence the matrices  and , are usually not square. Therefore, the eigenvalue problem (19) needs to be first transformed into a block column echelon form by applying a sequence of orthogonal transformations. Note that the eigenvalues are invariant with respect to orthogonal transformations and that in block column echelon form the eigenvalues depend only on the diagonal blocks. Following the Algorithm 4.1 given in Van Dooren (1979), first the SVD of the matrix is computed


This decomposition has the usual orthogonality properties  and , where

is the identity matrix, and the singular eigenvalues in 

are sorted in descending order. The matrices  considered are rank deficient so that  contains a number of zero diagonal entries. Hence, the right singular vectors  can be used to partition the matrices and such that


The number of columns of corresponds to the rank of and

denotes a zero matrix.

is chosen to have the same number of columns as . Next the SVD of the matrix is computed


The corresponding left singular vectors are used to further partition the matrices  and  such that


The number of rows of the upper blocks corresponds to the rank of . With these partitioned matrices the generalised eigenvalue problem (19) can be rewritten as


Permuting the rows yields the desired column echelon form


The eigenvalues  of this problem are the same as the ones for the original eigenvalue problem (19) because only orthogonal transformations have been applied. Furthermore, to find a set of eigenvalues  and eigenvectors that satisfy (25) it is sufficient to consider the eigenvalue problem stipulated by the top diagonal . The top diagonal may however not yet be a square matrix in which case the above sketched transformations need to be repeated. These transformations can also be realised with other methods, for example LU decomposition Thang et al. (2009)

and QR decomposition 

Beelen and Van Dooren (1988), which may potentially improve the computation speed.

Among the eigenvalues of (25) only the non-complex ones correspond to intersection points with the patch. The coordinates of the intersection points  are obtained by evaluating the line equation . The corresponding parametric coordinates  on the patch can be determined with a final singular value decomposition . According to (15) the vector of the basis function values  must be in the left null space of . Multiple left null vectors, i.e. multiple zero diagonal entries in , indicate a self intersection of the surface at . The parametric coordinate is determined from the set of nonlinear equations


where represents the coefficients of one of the left null vectors. Although these equations are nonlinear the two components of  are determined by simply considering their ratios, like


which yields two linear equations that can be trivially solved. All intersection points with the parametric coordinates in the range are valid.

Finally, a criterion is needed to determine the rank of a matrix in computing the singular value decompositions with finite precision arithmetic. The numerical rank  of a matrix can be defined considering the ratios of consecutive singular values according to


where is the -th singular value (sorted in descending order) and is a prescribed tolerance. In our computations the tolerance is chosen with .

3 Lattice-skin structures

3.1 Geometry generation

The closed surface and the specific lattice structure that is to be added to support the surface are assumed to be given, see Figure 6. In this paper only periodic lattices with cubic cells of the same size are considered. Each cubic cell contains an arrangement of a small number of struts that is periodically repeatable over all cells. Depending on the placement of the joints within a cell and their connections with struts different structures are possible, like the octet, pyramidal or tetrahedral truss lattice. The assumption of cubic cells is here not overly restrictive, considering that amongst polyhedra with a small number of faces only the cube and the dodecahedron are able to uniformly fill the space Fleck et al. (2010). Lattices containing several different types of polyhedral cells or with gradually changing cell size are possible but will be not considered for the sake of simplicity.

Figure 6: Surface and the lattice and the projection of unit cell corners onto the surface.

To create the lattice-skin structure first the intersections of the lattice lines with the surface are determined. This is achieved by looping over all lines of the lattice and first determining the Bézier patches which might intersect the line segment and subsequently computing the exact intersection points. In design and analysis both the spatial and parametric coordinates of the intersection points are required. In most applications the size of cells is significantly smaller than the characteristic size of the surface so that there are a very large number of intersection points.

After the intersection points have all been determined the shapes of the cells closest to the surface are modified by projecting their corner vertices onto the surface. The lattice joints inside the space enclosed by the surface can be identified by examining their relative positions on the lattice line. The joints between the first and second, third and fourth, and so on intersection points have to be inside the surface. As indicated in Figure (a)a, always the corner vertices which lie within the solid enclosed by the surface are projected. It is alternatively possible to choose the corner vertices which lie outside. After a lattice matching the given surface has been created it becomes straightforward to create the truss lattice structure by first creating the joints and then the struts which connect them, see Figure (b)b.

(a) Geometry fitted lattice (b) Truss lattice structure
Figure 7: Truss lattice generation from a geometry fitted lattice.

3.2 Structural model

3.2.1 Skin

We briefly review the governing equations for the linear Kirchhoff-Love model which is used for modelling the skin, for more details see, e.g., Ciarlet (2005); Cirak et al. (2000). The mid-surface of the shell is parameterised with the curvilinear coordinates . The position vector of material points on the surface is denoted with . The superscript s has been introduced to distinguish between the variables for the shell and the lattice.

The covariant basis vectors and the unit normal vector of the mid-surface are defined with


Here and in the following the Greek indices take the values and the summation convention is used. As usual, the contravariant base vectors are defined through the relation , where is the Kronecker delta. In case of small mid-surface displacements the linearised membrane strain tensor and the bending strain tensor are, according to Cirak et al. (2000), derived to be




with .

The total potential energy of the displaced thin-shell is given by


where and denote the membrane and bending energy densities, and is the surface load applied to the shell. For an elastic and isotropic material the membrane and bending energy densities are given with


where is the shell thickness, and and are the Young’s modulus and the Poisson’s ratio. The fourth-order auxiliary tensor is defined as


with the contravariant metric .

3.2.2 Lattice

The lattice structure is modelled as a pin-jointed truss consisting of struts connected by pins that do not transfer moments. Each strut deforms only by stretching and does not bend. We denote the coordinates and the displacements of the joints with and . In case of small displacements the axial strain in a strut with the index is given by


where and are two “picking matrices” for selecting the displacements of the two joints attached to the strut with the index , is the strut length and is the unit tangent to the strut.

The total potential energy of the lattice, with no externally applied loading at the joints, is composed of the potential energies of the individual struts


where is the cross-section area of the strut with the index . The external potential energy has been omitted as the loading is usually applied through the skin. For an elastic material with the Young’s modulus the energy density is given by


3.2.3 Lattice-skin coupling and discretisation

Some of the lattice joints are by design located on the shell mid-surface . The displacements of these joints with  have to be compatible with the shell displacements


where are the parametric coordinates of the joints on the surface. This condition can be imposed with Lagrange multipliers and considering the total potential energy for the lattice and the skin


The stationarity condition for this potential yields the equilibrium equations


These equations lead to a linear system of equations after discretising the shell. In isogeometric analysis the shell surface and displacements in (40b) and (40e) are both represented with the same type of basis functions. Hence, we use the splines introduced in Section 2 for geometry also for approximating the displacements. For displacements there is usually no need for rational splines as the exact representation of quadric sections is not relevant. Moreover, both the geometry and displacements can be represented with Bézier instead of the spline basis functions as discussed in Section 2.1. If Bézier basis functions are used the computed element stiffness matrices and force vectors have to be projected back to spline control points Borden et al. (2011). Informally, this projection ensures that the Bézier basis functions on neighbouring patches remain smoothly connected for the displaced structure.

4 Examples

4.1 Freeform surface with an internal lattice

In our first example we investigate the time needed for generating a lattice within a given surface. Figure 8 shows the steps involved in generating an internal lattice for the Stanford bunny. The bunny surface is a Catmull-Clark subdivision surface and consists of facets. The number of unit cells for the lattice is chosen with . This leads to intersection points between the lattice and the surface. Each lattice cell is tessellated with struts to give a body-centred cubic unit cell as shown in Figure (a)a. The deformation of the lattice-skin coupled structure under a uniform loading is depicted in Figure (b)b.

(a) Control mesh (b) Immersed lattice (c) Intersection points
Figure 8: Generation of an internal lattice for the Catmull-Clark subdivision surface of the Stanford bunny.

All the intersection points between the surface and lattice are computed using the intersection algorithm developed in Section 2. The SVD and eigenvalue computations involved are computed with the open source library Eigen Guennebaud et al. (2010). For intersection detection we use 14-dops, consisting of six directions orthogonal to the coordinate planes, two directions parallel to the average normal of all the Bézier patches (evaluated by sampling), and six directions parallel to the lattice. There is no need to ensure that these 14 directions are always unique. In Figure 9 the 14-dops and the axis-aligned bounding boxes (with k = 6) are compared for a lattice with unit cells in terms of efficiency of the intersection detection. Different lattice orientations with respect to the coordinate axis have been considered. As to be expected, 14-dops are far more efficient than the axis-aligned bounding boxes when the lattice is not aligned with the coordinate axis. Lattice lines not aligned with the coordinate axis have much larger axis-aligned bounding boxes.

(a) Lattice orientation 0 (b) Lattice orientation 45
(c) Computation times
Figure 9: Efficiency comparison between k-dops and axis-aligned bounding boxes in dependence of the lattice orientation.

For evaluating the accuracy and efficiency of the implicitisiation, we compare it with the widely used subdivision, also referred to as divide-and-conquer, method for intersection computation Houghton et al. (1985). In the subdivision method the Bézier patches are successively refined until they are sufficiently flat so that the intersection computation reduces to the trivial intersection between a triangle and a line segment. In the chosen flatness criterion the difference between the maximum and minimum support heights along the direction of the average normal of a Bézier patch are measured, cf. (4). This difference controls the accuracy of the intersection computation and has to be less than a prescribed tolerance. The tolerance becomes, however, ineffective when the line is (nearly) tangential to the surface, which is not a trivial problem for the subdivision method Hasle et al. (2007). The performance comparison between the implicitisation and the subdivision methods is shown in Figure 10. Both the computation time and memory consumed increase linearly with the decreasing tolerance for the subdivision method. This reflects the linear relationship between the depth of the bounding volume tree and the required accuracy. In contrast, the implicitisation exhibits a high accuracy (more than ) but takes less time and memory compared with the subdivision method. The computations were performed on a computer equipped with an Intel Xeon(R) CPU E5-2623 v3 @ 3.00GHz processor and 32GB memory.

Figure 10: Computation time and memory usage for the implicitisation and the subdivision methods for different flatness tolerances (lattice orientation is  and 14-dops are used).

One advantage of the subdivision method is that it is, without any modification, applicable around extraordinary vertices, where each spline consists of a sequence of nested patches Stam (1998); Peters and Reif (2008); Zhang et al. (2018). Therefore, in the present paper around extraordinary vertices the intersections are computed with the subdivision method. In the considered Catmull-Clark mesh for the Stanford bunny the subdivision method was applied in out of patches always.

4.2 Sandwich plate

Next we consider a sandwich plate with a pyramidal core to verify the accuracy and convergence of the finite element computations, see Figure (a)a. Although it is trivial to determine the intersection points between the flat skin plates and the lattice, the parametric coordinates of the intersection points still require the solution of nonlinear systems of algebraic equations. The size of the sandwich plate is chosen with  and the distance between the top and bottom skins is . The pyramidal lattice core consists of cubic unit cells of side length  and the struts are inclined by  with respect to the bottom skin, see Figure (b)b. The both skins and the struts are made of the same material with a Young’s modulus  and a Poisson’s ratio . The top skin is subjected to a uniform pressure loading of . The bottom skin is at its edges simply supported and each skin is discretised with  Kirchhoff-Love shell elements and cubic b-splines.

(a) Sandwich plate with pyramidal lattice core
(b) Unit cell with pyramidal tessellation
Figure 11: Geometry of the sandwich plate.

The homogenised material properties for the pyramidal lattice core can be found in Deshpande and Fleck (2001) . Its in-plane Young’s modulus is zero and its out-of-plane shear modulus is


The relative density  of the core is


where is the strut diameter and is the strut length. Note that the two out-of-plane shear moduli have the same value  due to the symmetry of the pyramidal core. Analytic expressions for the displacements of a simply supported isotropic sandwich plate under uniform loading can be found in Allen (1969).

(a) Comparison with analytic results
(b) Convergence of the maximum displacement
Figure 12: Maximum displacement of the sandwich plate under uniform pressure loading.

Three different skin thicknesses and eight different strut diameters are considered. The comparison of the numerical and the analytic maximum displacements is plotted in Figure (a)a. It can be seen that the numerical results agree well with the analytic ones for a range of different strut diameters and skin thicknesses. As depicted in Figure (b)b the maximum displacement exhibits a quadratic convergence rate (for a plate with strut diameter and skin thickness ). The convergence rates for different strut diameters and skin thicknesses may fluctuate around the quadratic rate since in the reference solution Allen (1969) the lattice core is approximated with an isotropic bulk material, and the skin-lattice coupling has also an effect on the convergence rate.

4.3 Spherical sandwich cap

As a further application of the developed analysis technique the optimal design of a spherical sandwich cap is considered. The sandwich cap has a BCC lattice core as shown in Figure 13 and the lattice edges are aligned with the principal axes of the cap. Each unit cell has a side length of  and all struts have a diameter depending on the prescribed lattice volume ratio. The size of the sandwich cap projected onto the horizontal plane is . The distance between the top and bottom skin is . Each skin is discretised with  Kirchhoff-Love shell elements and cubic b-splines. The Young’s modulus and the Poisson’s ratio of the thin-shell and the lattice material are and . A uniform loading with the magnitude is applied in the region around the centre of the upper skin. The four corners of the top and bottom skins are fixed.

Figure 13: Doubly-curved sandwich panel.

The total volume of the sandwich cap is prescribed and the optimal lattice to total volume ratio is sought. In Figure 14 the structural compliance for different lattice to total volume ratios are plotted. Three different total volumes are considered. For each total volume there exists an optimal lattice to total volume ratio as can be inferred from Figure 14. Structures consisting either only of a lattice or only one shell skin are non-optimal. For a relatively high lattice volume ratio the two shell skins become ineffective and the compliance becomes very large. In contrast, for a relatively low lattice volume ratio the lattice is not sufficiently stiff to carry the shear forces occurring between the top and bottom skins. As also recently discussed in Sigmund et al. Sigmund et al. (2016) and illustrated in Figure 14 a lattice structure is not per se the most effective structure.

Figure 14: Compliance for different lattice to total volume ratios.

5 Conclusions

A surface interrogation technique for computing the intersection points of curves with spline surfaces has been introduced. The two key components are the hierarchical k-dop bounding volume tree for fast proximity search and the implicit matrix representation for intersection computations. As demonstrated, in contrast to conventional axis-aligned bounding boxes the k-dop bounding volumes can significantly shorten the search times and are, most crucially, insensitive to the orientation of the considered geometry. The matrix-based implicitisation of surface patches reduces intersection computations to the solution of SVD and generalised eigenvalue problems. The obtained non-complex eigenvalues represent all the intersection points between a surface patch and a curve. Multiple real eigenvalues indicate multiple intersections between the curve and the patch. While establishing the implicitisation matrix and solving the generalised eigenvalue problem the null spaces of a range of matrices need to be computed. The largest matrix is of size , i.e. for the cubic patches used in this paper. All the required null spaces and eigenvalues are computed with a high-performance linear algebra software library. The main advantage of the proposed interrogation technique is its robustness in comparison to Newton-Raphson, or its variants, and its efficiency in comparison to subdivision based methods. In none of the introduced examples, with up to intersections, were any robustness issues observed. Subdivision based methods are not competitive due to strict accuracy requirements in finite element computations. This leads to very deep bounding volume trees, with levels, which may consume too much memory and take too long to construct and traverse.

The developed interrogation technique has been applied to the isogeometric design and analysis of lattice-skin structures to be manufactured with 3d printing. The isogeometric design and analysis of such structures is currently hampered by the poor scalability of prevalent CAD systems and their limitations in representing non-solid entities, like lattices. The implemented approach combines a lattice, consisting of a set of nodes and their connectivity, with a surface for the skin. The lattice is discretised as a pin-jointed truss structure and the skin as a Kirchhoff-Love shell. This makes it feasible to design and analyse structures with several millions of entities and the entire approach is amenable to parallelisation. Although not pursued in this paper, the presented interrogation technique is also suitable for slicing the geometry model to generate machine instructions for the 3d printer Ma et al. (2004); Starly et al. (2005).

In closing, we reiterate the importance of surface interrogation in a gamut of isogeometric analysis applications offering opportunities for research, including in enforcement of non-penetration constraints in contact, tessellating cut-cells in immersed or embedded boundary methods and coupling trimmed shell patches. Beyond smooth splines, the proposed interrogation technique is applicable to any polynomial tensor-product basis, including a Lagrange basis.

6 Appendix

6.1 Products of Bernstein polynomials

According to the definition of the univariate Bernstein polynomials (3), the product of two univariate Bernstein polynomials with and with can be expressed with


Similarly, the product of two bivariate Bernstein polynomials and can be expressed with