1 Introduction
Almost all geometry models in isogeometric analysis are based on parametric rational spline curves and spline surfaces, such as NURBS, Tsplines 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 computeraided design (CAD) literature. Shape interrogation of curves and surfaces includes, for instance, the plotting of their differential geometric properties, like isocurvature 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 curvesurface intersection, but the presented techniques can be extended to curvecurve and surfacesurface 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 nonpenetration constraints between spline curves or surfaces Temizer et al. (2011); De Lorenzis et al. (2011), integrating over cutcells 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 rootfinding 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), NewtonRaphson 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 convexhull properties of spline surfaces. The convexhull property guarantees that a spline surface lies always within the convexhull of its control points, and the refinability property enables to create a hierarchy of nested patches with decreasing convexhull size. Hence, the intersection problem can be reformulated as the intersection between two bounding volume hierarchies for the convexhulls of the spline patches and the curve. To simplify the intersection computation between convexhulls it is expedient to approximate them with axisaligned bounding boxes, spheres or, as in the paper, with kdops (discrete orientation polytopes) Lin 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 NewtonRaphson iteration the patches can be relatively large and convergence is quadratic sufficiently close to the intersection points. The principal difficulties encountered with NewtonRaphson iteration, mainly divergence for not carefully selected starting points, are the same ones as in pathfollowing 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 kdop bounding volume tree Klosowski et al. (1998)
. As our numerical experiments indicate a conventional axisaligned bounding volume tree can be significantly slower than a kdop 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 NewtonRaphson iteration or subdivision refinements. In implicit form a parametric surface,
with and is represented as the zero isocontour 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 socalled base points and the impossibility of describing a surface with exactly three moving planes. The matrixbased implicitisation introduced by Busé Busé (2014) can be understood as the generalisation of the moving lines idea to surfaces. Instead with a closedform 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.As an application of the developed shape interrogation technique the isogeometric design and analysis of latticeskin 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 latticelike components can achieve significantly higher stiffnesstoweight and strengthtoweight 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 latticelike 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 latticeskin 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 KirchhoffLove thinshell and is discretised with spline basis functions Cirak et al. (2000); Cirak and Long (2011); Zhang et al. (2018). The truss and KirchhoffLove thinshell 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 kdop bounding volume tree is presented. Subsequently, the matrixbased 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 latticeskin structures are illustrated. The efficiency and accuracy of the developed approach are studied in Section 4 with detailed timing studies for generating latticeskin structure geometries and finite element analysis of simplysupported 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, Tspline or CatmullClark subdivision surface. Each of these splines can be converted into a set of rational tensorproduct 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(1) 
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 tensorproduct bsplines Borden et al. (2011). The basis functions and the control points are labelled with the multiindex . The tensorproduct Bézier basis functions of bidegree are defined by
(2) 
with the Bernstein polynomials
(3) 
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 reexpress a spline surface using any polynomial basis, the Bernstein polynomials and the resulting Bézier patches have the advantage of numerical stability and convexhull property. The convexhull 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 convexhulls can be conveniently used as a preliminary check to identify possible intersections.
It is worth noting that in case of CatmullClark 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 Surfacecurve intersection
The computation of the intersection point between a Bézier patch and a given curve is an important highlevel primitive task in the design of latticeskin 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 kdops (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 matrixbased implicit representation without resorting to a NewtonRaphson iteration.
2.2.1 Bounding volume trees and kdops 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 convexhull 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, convexhulls can be approximated with kdops Sabin (1992); Klosowski et al. (1998), see Figure 2. Kdops are the generalisation of axisaligned bounding boxes and their bounding surfaces come from a set of fixed orientations ( in 3D). Axisaligned bounding boxes with usually provide only a loose fit to convexhulls and give many false positives. For a relatively high number of orientations kdops can provide a much tighter fit and their intersections are easy to compute.
In the implemented kdop bounding volume tree each node represents the kdop 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 kdop of a set of patches either from their individual kdops or from the union of their vertices.
Indeed, for intersection checks there is no need to explicitly construct any kdops. To determine whether two kdops 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 kdops 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
(4) 
The support heights and for the linear Bézier curve are determined in a similar way. The two kdops intersect only when
(5) 
and they do not intersect when
(6) 
Hence, in the kdop 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
(7) 
2.2.2 Implicit matrix representation for intersection computation
Next, the matrixbased implicitisation of a single Bézier patch of bidegree as defined in (1) is developed. To this end, consider the auxiliary vector of polynomials
(8) 
with Bézier basis functions , coefficients to be yet determined and the multiindex . According to Busé (2014), the bidegree of the Bézier basis can be chosen with or higher; alternatively or higher is also possible. The bidegree 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
(9) 
that is,
(10) 
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 bidegree
(11) 
where are the coefficients of a matrix . Note that there are basis functions in the basis with bidegree , basis functions in the basis with bidegree 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 nontrivial ’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
(12) 
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
(13) 
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
(14) 
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.
Instead of the geometrically motivated moving plane interpretation it is possible to advance an equivalent matrixbased linear algebra explanation. To this end, the plane equations (14) are, after introducing (8), rewritten as
(15) 
where constitutes the th column of a matrix , that is,
(16) 
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 higherorder curve see Appendix 6.3. The line is assumed to be given in parametric form with
(17) 
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
(18) 
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
(19) 
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
(20) 
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(21) 
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(22) 
The corresponding left singular vectors are used to further partition the matrices and such that
(23) 
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
(24) 
Permuting the rows yields the desired column echelon form
(25) 
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 noncomplex 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
(26) 
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
(27) 
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
(28) 
where is the th singular value (sorted in descending order) and is a prescribed tolerance. In our computations the tolerance is chosen with .
3 Latticeskin 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.
To create the latticeskin 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.
3.2 Structural model
3.2.1 Skin
We briefly review the governing equations for the linear KirchhoffLove model which is used for modelling the skin, for more details see, e.g., Ciarlet (2005); Cirak et al. (2000). The midsurface 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 midsurface are defined with
(29) 
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 midsurface displacements the linearised membrane strain tensor and the bending strain tensor are, according to Cirak et al. (2000), derived to be
(30) 
and
(31) 
with .
The total potential energy of the displaced thinshell is given by
(32) 
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
(33) 
where is the shell thickness, and and are the Young’s modulus and the Poisson’s ratio. The fourthorder auxiliary tensor is defined as
(34) 
with the contravariant metric .
3.2.2 Lattice
The lattice structure is modelled as a pinjointed 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
(35) 
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
(36) 
where is the crosssection 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
(37) 
3.2.3 Latticeskin coupling and discretisation
Some of the lattice joints are by design located on the shell midsurface . The displacements of these joints with have to be compatible with the shell displacements
(38) 
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
(39) 
The stationarity condition for this potential yields the equilibrium equations
(40a)  
(40b)  
(40c)  
(40d)  
(40e) 
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 CatmullClark 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 bodycentred cubic unit cell as shown in Figure (a)a. The deformation of the latticeskin coupled structure under a uniform loading is depicted in Figure (b)b.
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 14dops, 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 14dops and the axisaligned 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, 14dops are far more efficient than the axisaligned bounding boxes when the lattice is not aligned with the coordinate axis. Lattice lines not aligned with the coordinate axis have much larger axisaligned bounding boxes.
For evaluating the accuracy and efficiency of the implicitisiation, we compare it with the widely used subdivision, also referred to as divideandconquer, 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 E52623 v3 @ 3.00GHz processor and 32GB memory.
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 CatmullClark 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 KirchhoffLove shell elements and cubic bsplines.


The homogenised material properties for the pyramidal lattice core can be found in Deshpande and Fleck (2001) . Its inplane Young’s modulus is zero and its outofplane shear modulus is
(41) 
The relative density of the core is
(42) 
where is the strut diameter and is the strut length. Note that the two outofplane 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).


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 skinlattice 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 KirchhoffLove shell elements and cubic bsplines. The Young’s modulus and the Poisson’s ratio of the thinshell 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.
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 nonoptimal. 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.
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 kdop bounding volume tree for fast proximity search and the implicit matrix representation for intersection computations. As demonstrated, in contrast to conventional axisaligned bounding boxes the kdop bounding volumes can significantly shorten the search times and are, most crucially, insensitive to the orientation of the considered geometry. The matrixbased implicitisation of surface patches reduces intersection computations to the solution of SVD and generalised eigenvalue problems. The obtained noncomplex 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 highperformance linear algebra software library. The main advantage of the proposed interrogation technique is its robustness in comparison to NewtonRaphson, 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 latticeskin 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 nonsolid 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 pinjointed truss structure and the skin as a KirchhoffLove 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 nonpenetration constraints in contact, tessellating cutcells in immersed or embedded boundary methods and coupling trimmed shell patches. Beyond smooth splines, the proposed interrogation technique is applicable to any polynomial tensorproduct 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
(43) 
Similarly, the product of two bivariate Bernstein polynomials and can be expressed with
Comments
There are no comments yet.