Rectifying Strip Patterns

Straight flat strips of inextensible material can be bent into curved strips aligned with arbitrary space curves. The large shape variety of these so-called rectifying strips makes them candidates for shape modeling, especially in applications such as architecture where simple elements are preferred for the fabrication of complex shapes. In this paper, we provide computational tools for the design of shapes from rectifying strips. They can form various patterns and fulfill constraints which are required for specific applications such as gridshells or shading systems. The methodology is based on discrete models of rectifying strips, a discrete level-set formulation and optimization-based constrained mesh design and editing. We also analyse the geometry at nodes and present remarkable quadrilateral arrangements of rectifying strips with torsion-free nodes.


INTRODUCTION
Classical and computational differential geometry offer methods to design and fabricate complex geometric structures from simple elements.This is of great significance, especially in architecture, where individual, large-scale designs demand bespoke solutions that cannot be mass-produced but need to stay within a given budget.A large amount of contributions to architectural geometry thus deals with the simplification of construction elements, such as panels, joints or beams.
The present paper is motivated by the recent emergence of various types of architectural structures that are fabricated by bending originally flat straight pieces of material, typically wood or metal.Fig. 2 shows three built examples of such gridshells.They differ in how the individual elements are placed with respect to an underlying reference surface and in the pattern that is formed by them.If a bent element lies tangential (Fig. 2 a) to the reference surface , it follows a shortest path (geodesic) on .If it is placed orthogonal (Fig. 2 b) to , it is aligned with an asymptotic curve of , i.e., a curve of vanishing normal curvature.Very recently, even the case of a constant angle (Fig. 2 c) (different from 0 and /2) of the bent elements against  has been considered.They follow so-called pseudo-geodesics on  and thus one speaks of pseudo-geodesic gridshells [Mesnil and Baverel 2023].(a) The Polydôme by Julius Natterer using stacked geodesic planks, and a diagonal geodesic cover [Natterer et al. 2000]; (b) Timber gridshell by Eike Schling et al. combining asymptotic and geodesic lamellas in a trihex grid (https://eikeschling.com).(c) Pseudo-geodesic gridshell by [Mesnil and Baverel 2023], using inclined strips.
The pattern of elements on the final structure may just be a single family of elements.Commonly, gridshells are constructed from two families of beams and a third set of elements that triangulate and thus brace the structure.This leads to three families of curved elements arranged in a so-called 3-web or Kagome pattern.In the present paper we are interested in all these potential patterns, thus generalizing the geodesic patterns of [Pottmann et al. 2010].
For some of the gridshells discussed above, there are efficient computational design methods.This includes asymptotic gridshells and geodesic gridshells.However, turning to pseudo-geodesic gridshells or hybrid asymptotic-geodesic shells, this is no longer true.There are computational models based on optimization, which however lack especially in the difficult initialization step.In other words, it is within some limits possible to modify a given design while keeping certain constraints, but it is hard to access the design space.The present paper aims to fill this gap.
While being motivated by specific gridshells and their underlying patterns, we present a more complete study of the problem area.We turn practical constraints on fabrication (such as torsion-free nodes) and function (e.g.shading) into geometric ones and show how to use them for computational design.

Contributions and Overview
Bent strips of material that are straight and flat prior to bending are already determined by their central curve or one of the boundary curves, as their binormals are strip surface normals (Sec.2.1).This reduces a major part of our work to a computational design problem of families of surface curves with constraints on their binormals.For its solution, we propose a novel level set formulation in the spirit of discrete differential geometry, by discretizing the geometry rather than equations.This methodology part of Section 3 is surveyed in Fig. 3.The second major contribution of our paper are individual applications, illustrated in Fig. 4. We conclude in Sec.7 with statistics Fig. 3. Overview of methodology.Given a surface that guides the shape of a strip structure, we compute base curve families of strips with a level set algorithm (Sec.3.1,3.2and 3.3), attach strips to each family (Sec.3.4) and optimize meshes which represent base curve families of strips (Sec.3.5).Fig. 4. Overview of applications.Gridshells formed by two or three families of rectifying strips are presented in Sec. 4 along with physical models.The difficult topic of gridshells with torsion-free nodes, where transversal strips meet at straight line segments, is addressed in Sec. 5.In Sec. 6 we design shading systems via a single family of rectifying strips.

Related Work
Differential Geometry.Geodesics and asymptotic curves on surfaces are a well-studied topic in classical differential geometry, and the same is true for their discrete counterparts [Bobenko and Suris 2008;do Carmo 1976].Wunderlich [1950] introduced and studied pseudo-geodesic curves on surfaces, especially cones and cylinders.The quadrilateral nets of asymptotic curves on negatively curved surfaces (A-nets) are very well understood, including their discrete versions [Bobenko and Suris 2008].A-nets with a constant intersection angle  of asymptotic curves represent surfaces with a constant ratio of principal curvatures [Wang and Pottmann 2022], which include minimal surfaces for  = /2.
Quad nets of geodesics (G-nets) provide more freedom than Anets.Remarkable are surfaces which possess conjugate G-nets, also known as Voss nets [Sauer 1970;Voss 1888].Discrete Voss nets are flexible quad meshes with planar faces, i.e., they become mechanisms when one considers faces rigid and places hinges at the edges [Bobenko and Suris 2008;Sauer 1970;Wunderlich 1951].Orthogonal G-nets exist exactly on developable surfaces.Their discrete models have been developed by Rabinovich et al. [2018] for geometric design with developable surfaces.Orthogonal nets where only one family of curves are geodesics exist on all surfaces.Discrete versions for design and fabrication are found in [Wang et al. 2019].

Geometry of webs.
In the present paper, we are mostly interested in triangular webs (3-webs), which consist of three families of curves on a surface, such that through each surface point there is a curve from each family.Such a 3-web can be represented by a parameterization x(, ) where the three families of web curves possess constant , constant  and constant  + , respectively.For the geometry of webs in general, we refer to the monograph by Blaschke and Bol [Blaschke and Bol 1938].In web geometry, these curve families are mostly considered as dense, but also discrete webs have been studied.There, one picks just iso-lines of integer values of , ,  + .
Our focus is on webs of asymptotic, geodesic and pseudo-geodesic curves.The best-studied ones are the 3-webs of geodesics.In the plane, these are webs of straight lines, which are formed by the tangents of an algebraic curve of class 3 [Graf and Sauer 1924].R. Sauer [1926] studied geodesic 3-webs on surfaces at hand of a discrete model.The characterizing partial differential equation for geodesic webs is found in [Mayrhofer 1931], but its explicit solution appears to be difficult.Geodesic webs exist on all surfaces of constant Gaussian curvature, since those can be mapped to the plane so that all geodesics get mapped to straight lines [Graf and Sauer 1924;Volk 1929].
Geometry Processing and Computational Design.There is a wealth of contributions on the efficient computation of geodesic curves (see [Crane et al. 2017] and references therein).Computational design of webs of geodesics and/or planar curves has been addressed in [Deng et al. 2011;Pottmann et al. 2010].Geodesic webs also appear as weaving patterns generalizing the traditional craft of basket weaving to more complex freeform shapes [Ayres et al. 2020[Ayres et al. , 2018;;Vekhter et al. 2019].The use of initially straight ribbons limits the space of potential shapes and thus this weaving technique has been extended to curved ribbons [Baek et al. 2021;Ren et al. 2021].[Pillwein et al. 2020;Pillwein and Musialski 2021] address the computational design of elastic geodesic gridshells, that are formed by an originally flat arrangement of straight planks, including solutions of the inverse problem of approximating a given shape with such a structure.Their shells are special X-shells [Panetta et al. 2019], for which a solution to the inverse problem is not known.[Liu et al. 2023] present structures that deploy from a stack of initially straight or circular elastic strips.They exhibit close relations to Chebyshev nets, which previously occurred in computational fabrication [Garg et al. 2014;Masson 2017;Sageman-Furnas et al. 2019].
Architectural Geometry and Gridshells.Gridshells have been of great interest in the architectural and engineering community for over a century.We point to early work by Vladimir Suchov (see [Schling and Barthel 2021]) and Frei Otto [Happold and Liddell 1975;Hennicke 1974].Natterer stacked flat planks along geodesic paths to create ribbed gridshells [Natterer et al. 2000].Geodesic gridshells recently gained popularity [Martín-Pastor and González-Quintial 2021;Montagne et al. 2020], also as deployable structures [Soriano 2017].Asymptotic gridshells are still a novelty in architecture.The tall orientation of planks creates a high stiffness, while offering an elastic erection process from flat to curved grid.Asymptotic gridshells have been built from steel, with constant node angle [Schling et al. 2018], and as web bisecting the principal curvature lines to allow for a simplified facade from developable or planar quad panels [Schling and Wan 2022].Asymptotic grids can be actuated to become highly controlled transformable structures [Schikore et al. 2021].Recently, hybrid asymptotic and geodesic 3-webs have been described [Schling et al. 2022], and implemented to build a robust trihex gridshell from timber planks (see Fig. 2 b).Mesnil and Baverel [2023] introduced pseudo-geodesic curves on rotational surfaces as a possible construction strategy for gridshells and recently built a large-scale prototype using timber planks and repetitive steel joints (Fig. 2 c).
Apart from practical applications as load-bearing gridshells, slender louvers are used in architecture for shading and ventilation, see Fig. 5.The controlled incline of louvers offers the possibility to redirect light and control views.

PRELIMINARIES 2.1 Bending straight strips: rectifying developables
Bending a rectangular strip  0 of inextensible flat material yields a developable surface strip  (Fig. 6).The straight central line c 0 of the strip  0 is bent to a geodesic curve c ⊂ .A geodesic on a surface is characterized by vanishing geodesic curvature.Equivalently, its principal normals are surface normals.This shows how to construct the strip when c is given: We consider the orthonormal Frenet frame (t, n  , b) along c, formed by the unit tangent vector t, principal normal vector n  and binormal vector b.Vectors t and n  span the osculating plane, and vectors t and b span the rectifying plane.For a geodesic on a surface, the rectifying planes are tangent to the surface.Therefore, the developable surface carrying the strip  has the rectifying planes of c as tangent planes.Since a developable surface contains only a one-parameter family of tangent planes, the surface  lies in the envelope  of rectifying planes, which is called the rectifying developable of c.As a developable surface,  contains a family of straight lines, called rulings (see Fig. 6). is also the rectifying developable of the strip boundaries.In the following, we call any developable strip with a straight development a rectifying strip.Sec.5.1 presents a closer local analysis of rectifying strips to better understand tolerances at nodes of a gridshell.
A discrete curve is a polyline  with vertices c  .Its edges can be seen as discrete tangents.The connecting planes   of 3 consecutive vertices c  −1 c  c +1 are discrete osculating planes at c  .The circle   through c  −1 c  c +1 is a discrete osculating circle.If edges are of constant length (discrete arc length parameterization), the tangent   ∈   of the circle   at c  is an outer bisector of the polyline, and the other bisector, containing the center of   , is a discrete principal normal.This yields a discrete Frenet frame (t  , n   , b  ).We use the definition via bisectors also for a general discrete parameterization.With unit edge vectors e  = (c +1 − c  )/∥c +1 − c  ∥, we have (1) The discretization of rectifying strips is addressed in Section 3.4.

Attaching strips at constant angle: pseudo-geodesics
A rectifying strip can be attached tangentially to a surface  exactly along a geodesic of .We also require strips  that are inclined under a constant angle  ≠ 0 against .Equivalently, we look at curves c ⊂  whose rectifying planes form a constant angle  ≠ 0 with the tangent planes of .In other words, the binormals of c form the constant angle  =  2 −  with the normals of .According to W. Wunderlich[1950], such curves are called pseudo-geodesics or shortly P-curves in the following.The case  = /2 characterizes asymptotic curves.
We complete the tangent vector t of the surface curve c and surface unit normal vector n via the side vector u = n × t to the Darboux frame.The curvature vector n  of c is decomposed into a tangential and normal component via n  =   u +   n, with   ,   as geodesic and normal curvature, respectively and n  as principal normal vector.For a P-curve c, the angle between n  and n equals , i.e. n  = sin  u + cos  n.Hence, a P-curve is characterized by a constant ratio of geodesic and normal curvature, tan  =   /  .Not only   , but also  and   have to be considered with a sign.The signs determine whether c turns to the left or right side of t.
Later we also consider surface curves along which the angle  is not constant, but varies according to the desired application.
For the main part of this paper, the surface  is given as a triangle mesh, and the curve is a polyline on it, namely a level set of a function  defined on the mesh.In the final steps of the individual applications, the curve may also be a parameter line in a constrained mesh (quad mesh for grid-like patterns and triangle mesh for 3webs).In all cases, the pseudo-geodesic property is expressed by a constant angle  between discrete binormal b and surface normal n, b • n − cos  = 0.
(2) [Jiang et al. 2019] proposed a tracing algorithm to find P-curves on triangle meshes, but their shape control is cumbersome, depending highly on the choices of  and tangent direction at the starting point.Moreover, it is hard to avoid intersections of neighboring P-curves or self-intersections.Therefore, we turn to a level-set formulation.

METHOD
In this section, we introduce our level-set based formulation and the optimization framework that forms the core of our computational design tools.The major part is devoted to the level-set algorithm.In Sec.3.4 we describe the computation of rectifying strips and Sec.3.5 deals with constrained mesh optimization for strip structures and user-guided shape modification.We found it easier to describe initially just the design of patterns formed by P-curves.In Sections 5 and 6 we show how to easily adapt our algorithms to similar constraints that occur in practical applications.

Curve families on surfaces via level sets
Given a surface  with disk topology and a function  defined on it, its level sets  = const.form a family of curves that cover the surface without intersections of different curves.If ∇ ≠ 0 everywhere, there are no self-intersections of individual curves and through each surface point there is exactly one level curve, orthogonal to ∇ .
It is well known (see e.g.[Pottmann et al. 2010]) that geodesic and normal curvature of level sets can be computed via where II is the second fundamental form of  and  is rotation by 90 degrees.Pseudo-geodesic level sets can then be expressed via the constant ratio   /  .Since we have to discretize later anyway and we have to include other constraints on binormals of level sets, we do not discretize equations, but take a geometric discrete approach: We express the binormal vectors of the discrete curves explicitly and control the behavior of the curves by applying constraints to the binormal vectors.
Controlling binormal vectors of the level curves.Let the surface  be represented by a triangle mesh (V, E, F ) that is topologically equivalent to a disk.We want to compute a scalar function  on , with non-vanishing gradient ∇ , whose level sets have binormal vectors that are constrained by the application.For pseudo-geodesic level sets, this is the constant angle constraint Eq. (2).
q q q q q q q q q q q q q q v 0 v Fig. 7. v 0 is an interior vertex of the triangle mesh .The level curve (red) passing through v 0 has value  0 , and has intersection points p, q with two edges (yellow) of the vertex star around v 0 .
The function  is defined by its values on the mesh vertices and linearly interpolated in the triangles.It is sufficient to express the constraint on binormals locally at each vertex v  as well as the level polyline passing through v  .This is done as follows.Consider the level polyline c passing through a vertex v 0 of  (Fig. 7) and let  0 be the function value at v 0 .If v 0 is not on the boundary of  there is a unique level curve c passing through v 0 , that intersects the boundary of the vertex star in two further points.These points p and q define two edges e p = v 0 −p and e q = v 0 −q .Assuming that e p and e q are linearly independent, they span the discrete osculating plane and the discrete binormal vector b of the curve c can be computed as b = ±e p × e q /∥e p × e q ∥.
Fig. 8.For a given angle  , Eq. ( 2) holds with appropriately oriented binormals for two osculating planes  1 and  2 .Equation (3) helps to select the desired osculating plane without having to care about the orientation of the binormal vector.
Handling the angle constraint between surface normals and binormals of level sets.Here we have to correctly select one of the two possible osculating planes  1 ,  2 that contain a given curve tangent and form an angle  with the tangent plane; see Fig. 8.In order to simplify the subsequent optimization, we want a criterion that is agnostic to the change b → −b.We use the side vector u and normal vector n.For one choice of osculating plane, the scalar products b • n and b • u have the same sign ( 1 in Fig. 8) and for the other choice they have a different sign ( 2 in Fig. 8), leading to the angle constraints (3) Clearly this depends on the sign of  and the side vector u.In Section 3.2.2,we show how to achieve a consistent side vector field.Fig. 9 illustrates the necessity of caring about the sign of  .
Validation of the discretization.We validate our discretization on a cylinder of revolution, whose pseudo-geodesics have known explicit parametric representations [Wunderlich 1951]; see Fig. 10.In Fig. 11, we also compare level sets on two models with the same shape, but different triangulations.The differences between the corresponding curves are minor, indicating consistency of our method under different triangulations.

Optimization framework
We apply Eq. (3) as our constraint on level sets to represent P-curves.
For that, we propose the following 3 steps to construct a function whose level curves are P-curves of an input angle  (see Fig. 12).The first step is to create a few guiding pseudo-geodesic curves on the reference triangle mesh using a tracing algorithm (Fig. 12 a).The second step initializes the scalar function  and optimizes it so that all level sets become pseudo-geodesics (Fig. 12 b-d).In the third step, we extract the level curves from the scalar function and  9 × 10 −6 ; see Eq. 12).The parametric P-curve (red) computed according to [Wunderlich 1951] matches the level set (black) well.The Hausdorff distance of the two curves is 0.47% of the curve length.11.Stability under re-meshing.The triangle mesh in (a) has 8960 vertices, and the level sets are P-curves of  = 60 • .The mesh is simplified to 657 vertices using fTetWild [Hu et al. 2020], but the shape is kept.We initialize  of (b) by linear interpolating the function values of (a), then optimize for P-curves (Section 3.2.2).We extract 10 level sets from each of (a) and (b) (Section 3.4) and compare them (c).The red and blue curves are for (a) and (b), respectively.The Hausdorff distance of the curves of (a) and (b) is 0.76% of the diagonal length of the bounding box of the mesh.
optimize the corresponding rectifying strips in a post-processing step (Fig. 12 e).

Tracing guide curves.
We first create a few guide curves on the reference triangle mesh with a robust version of the tracing algorithm in [Jiang et al. 2019].
Starting from a segment p  −1 p  on the triangle mesh (see Fig. 13), we look for the next point p +1 on an appropriate edge v  v  , such that the discrete osculating plane   spanned by p  −1 , p  , p +1 forms the desired angle  with the normal vector n  at p  .Equivalently, the binormal vector b  has to form the angle  = /2 − with n  .All planes through p  that form a constant angle  with n  are tangent to a rotational cone Γ  with axis vector n  and thus there are two candidate planes for   which also contain p  −1 .We compute the intersections of these two planes and the edges through the two vertices v 0 or v 2 in the following way.Parametrizing an edge where It is a quadratic equation in ; only one of the solutions yields the correct sign of the angle against a local Darboux frame (as explained in Sec.3.1).We select as next point p +1 the nearest intersection with an edge in the chosen neighborhood so that ∥p +1 − p  ∥ >  0 .The threshold  0 avoids too close vertices and unstable discrete osculating planes.The traced curves are used to guide level sets to desired directions, and quickly show if P-curves of angle  behave as desired for design.

Optimization of level sets towards P-curves.
To compute a function  on the given mesh whose level sets are P-curves, we first compute an initial guess using the traced polylines.Then, we optimize the level sets so that Eq. ( 3) is satisfied.
Initialization.Starting from a set C  of  ordered intersection-free traced curves c  ,  = 1, . . ., (Fig. 12a), we set the function value of each curve c  ∈ C  to a constant value   , e.g.  = .and optimizes the scalar function so that its level sets interpolate the guide curves (c).Then, we optimize the scalar function so that the level sets are pseudo-geodesic curves (d).The quality of the extracted curves and rectifying strips is improved through post-processing procedures (e).The detailed evaluation for (d) is shown in Table 2 and Fig. 40.The evaluation for (e) is shown in Table 4.  4), the next point p +1 of p  −1 and p  can be found.To improve stability, the edge containing p +1 (here v  v  ) is selected among the edges connecting v 0 or v 2 .
Having these function values   at the vertices p  of the curves c  , we initialize the scalar function  on a parameter domain.To map the mesh plus guide curves into the plane, we use a conformal mapping algorithm implemented in [Jacobson and Panozzo 2017].A point p ∈  is mapped onto (p) ∈ R 2 .To obtain non-intersecting level sets, we first initialize  as a linear function in 2D, as illustrated in Fig. 12 (b): , where c 1 0 and c  0 are the starting points of c 1 and c  , respectively.Then we have to optimize the values of  on all vertices v  of the mesh  so that the level sets are smooth and interpolate the guide curves.Using the notation in Fig. 13 for the local situation near a vertex p of a polyline c  ∈ C  , we define an objective function Here v 0 v 1 is the edge of  that contains p, and  0 and  1 are the function values at v 0 and v 1 , respectively. trace relates the value   at p by linear interpolation to the values  0 ,  1 at the end points of the edge v 0 v 1 .
To obtain smooth level sets, we adopt the squared Hessian energy  fair proposed by [Stein et al. 2018] and implemented in [Jacobson and Panozzo 2017]: A (v) is the area of the Voronoi cell of vertex v, and ∥ (v)∥ 2 is the squared Frobenious norm of the Hessian [Stein et al. 2018].This energy yields smoother level sets of  on the mesh boundaries than the Laplacian energy.Our final objective function for initialization is where  0 and  1 are properly chosen weights.Typically we choose  0 = 3,  1 = 1.A result of this step is shown in Fig. 12 (c).
Optimization for pseudo-geodesics.We optimize the scalar function  such that the level set through each vertex is locally a P-curve of angle  .Firstly, for each inner vertex v 0 of , we search among the boundary edges of its vertex star to find the edges that have intersections with the level curve passing through v 0 , as shown in Fig. 7.If v 0 is not a singularity of  , there are two edges v 1 v 2 and v 3 v 4 that intersect the level set at p and q, respectively.Let v 1 , v 2 , v 3 , v 4 be counterclockwise arranged w.r.t. the normal vector n at v 0 ,  1 >  2 and  3 >  4 , the direction of the tangent vector at v 0 of polyline {p, v 0 , q} can be prescribed as from p to q.This orientation is crucial to our method, since we need to construct a continuous Darboux frame along the curves to obtain continuous osculating planes.In view of Eq. ( 3), the angle constraint is written as where b is the unit binormal vector satisfying b and u is the unit side vector of the Darboux frame, constrained by The vectors b and u are used as auxiliary variables along with the values of  on mesh vertices.We note that the constraints in Eq. ( 10) cannot prevent u from inverting its orientation.In our implementation, in each iteration, we check if det(q − p, u, n) > 0; otherwise, we re-initialize u as the unit vector of n × (q − p).
Since we want to avoid self-intersecting level curves and singularities in the curve family, we add a term  grad that prevents a vanishing gradient of  , where ∇ (f) is the gradient of  on triangle f and  is a positive constant value.As mentioned in [Pottmann et al. 2010], this term also acts as a regularizer to control the spacing of level sets: ∥∇ ∥ −  = 0 means that the distance between two level sets  =  and  =  + 1 is roughly 1/ .Pseudo-geodesic level curves are obtained by minimizing The variables are the values of  on the vertices of the mesh , binormal vectors b, and side vectors u.We employ a Levenberg-Marquardt method as in [Madsen et al. 2004] to minimize  pg .The parameters  for the selected models are presented in Table 2.
To make sure that  angle better reflects the numerical errors in the angles regardless of the Voronoi cell areas of the vertices, we take such that the average of A () is 1.Signed angles are not necessary for geodesic and asymptotic level sets.These are optimized using energy terms  geo and  asy that slightly differ from  angle in Eq. ( 8) to avoid auxiliary variables, geo = 0 means that the osculating plane spanned by e p = v 0 − p and e q = v 0 − q contains the normal n at each vertex star, while  asy = 0 says that it is orthogonal to n and thus tangent to .
Interactive steering of pseudo-geodesic curves.The tracing algorithm produces P-curves which strongly depend on the starting point, starting direction, and the angle  .It requires the user to trace the curves one by one while avoiding intersections between them.To lift the burden of this tedious process, we also provide an interactive design method that allows the user to sketch on the surface and then solve for the best-fitting P-curves automatically.The minor adaptation required to extend this beyond P-curves, is also useful for the applications in Sections 5 and 6.
As shown in Fig. 14, non-intersecting curves represented as sequences of points are sketched by the user.Before optimization, we filter the curves such that on each triangle of  there is at most one sketched curve point.Then a smooth scalar function whose level sets follow the input curves is initialized via the method in Sec.3.2.2,but the tracing energy  trace is replaced by Here  0 ,  1 ,  2 are the function values at the vertices of the triangle  in which the sketched point p lies, (, , ) are the barycentric coordinates of p w.r.t. , and  c is a value which is constant along each sketched curve.In each iteration, we take  c as the average of the function values of the sketched points along the curve c in the previous iteration.
After obtaining a smooth scalar field whose level sets approximate the sketched curves, we compute the discrete binormal vectors on the vertices of the triangles where the sketched point lies and take the average  of the angles between the binormal vectors and the surface normal vectors as the target pseudo-geodesic angle.Then we solve for pseudo-geodesic level sets with angle  according to Sec. 3.2.2.
In our implementation, the interactive design of P-curves is fully automatic with prescribed parameters (see Table 2).To enhance the applicability of the parameters for triangle meshes with different scales or densities, we uniformly scale the meshes prior to optimization, so that their average edge length is 1.Even then we cannot   2 and Fig. 41).guarantee that the same parameters work for all models.The reasons are not only the unavoidable limitations of geometry, but also the avoidance of singularities (Section 7).Fig. 15 shows a failure example using optimization with default parameters.The problem can be resolved by changing parameters during optimization.
P-curve families with changing angles.For applications to gridshells, a small angle  between rectifying strips and the normals of the reference surface  can be an advantage.On a negatively curved surface , this is easy and one can even achieve  = 0 by placing the strips along asymptotic curves.However, for positively curved , P-curves to a constant angle  may be too curved and not form a family of non-intersecting curves, resulting in an optimization failure.Thus, we propose the following method which generates well-arranged P-curves with small angles  and allows for a variation of the angle  between different P-curves.
We aim at P-curves to sufficiently small angles  which arise as level sets of  .The computation of the angle assignment function  ( ) proceeds as follows.Given a smooth function  on , we compute the angle  between surface normal n and level set binormal b at each vertex.Then, the data (,  ) is fitted by a polynomial curve   , and the target angles on the vertices for the next iteration are obtained from   .Supervised by the fairness terms  fair and  grad , and iteratively updating   , a family of nearly parallel P-curves is generated.In our implementation, we use a polynomial of degree 3 to fit the data (,  ).To yield a smaller  , we decrease the target angle  by 5% each time we update the polynomial   .Fig. 16 shows   an example where we obtain a family of P-curves whose angles  vary from 36.2 • to 67.1 • .

Co-optimizing reference surface and level sets
In view of the targeted applications, we use a level-set formulation to obtain families of pseudo-geodesics where exactly one curve passes through each surface point.For a given surface and angle  , such families may not exist.As an example, take a spherical patch and a small angle  .P-curves then are circles of a small radius.They will have an envelope and the area beyond the envelope will not be covered at all.Thus, we implemented algorithms that can change the surface and the angle  .
Given a scalar function  on a triangle mesh , we can optimize  such that the level sets become the desired P-curves.Fig. 17 shows examples, where we changed the angles of P-curves by optimizing  while keeping the values of  at the mesh vertices unchanged.
The energy  pg is written as Here  is to provide smooth level sets on the surface.Taking the notations in Fig. 7, close is the approximation energy: where v * is the closest point of v on , n * is the normal vector of the initial surface on v * , and v ori is the original position of v.This term allows the vertices to move within the tangent planes of the reference mesh, preventing large-scale deviations from the initial surface.  edge restricts the changes in edge lengths, where  * e is the edge length of the original corresponding edge of e.    plays a role as a regularizer to prevent mesh shrinkage.4 show the evaluation for scalar field optimization and quality of the developable strips, respectively.
In this step we initially extract the level sets and improve the accuracy of angles between binormals and surface normals.The major part is devoted to the computation of the rectifying strips.
Having the scalar function  at mesh vertices, we linearly interpolate in the triangles and extract those level sets on which we build the strips.The  level sets correspond to the  scalar values sampled uniformly from [min ( ), max ( )].Energy  grad guides towards level sets uniformly distributed on the surface.However, the level polylines are initially uneven and binormals and the angle  are not highly accurate, see Fig. 18 (left).We re-sample the polylines and improve the angle  between binormals and surface normals by optimization.We still use Eq. 8 -10, but on vertices c  of the polylines instead of vertex stars of the triangle mesh.Next, we explain the subsequent computation of discrete rectifying strips; see also Fig. 18 (right).Preferred version of rectifying strips.We now turn to another discrete rectifying developable, which we did not see in prior work.It is a strip of planar quads that contains the given polyline  as a medial line and unfolds to a straight rectangular strip in the plane (see Fig. 20).The straight unfolding requires us to find discrete ruling lines d  through c  , so that the sum of angles  −  and  +  between d  and the two edges through   equals  (see Fig. Thus, one may choose one of these planes and get the others by successive reflection, resulting in a discrete rectifying developable.This degree of freedom, only present in the discrete model, is removed by optimization, aiming at proximity to the previously discussed model  (Fig. 19).We compute initial guesses for d  ∈   by averaging the unit direction vectors of the two discrete rulings   and  +1 in   (see Fig. 20).This initial choice does not guarantee coplanar lines d  −1 and d  .Thus, we finally optimize for planarity while keeping d  ∈   .To express that rulings d  , d +1 are coplanar with the edge vector e  , one introduces a unit normal vector n  17) Constraints ( 16) and ( 17) are quadratic in the variables c  , e  , d  , n   , n  ,+1 if the values of the denominators are taken from the previous iteration.The constraint system is formulated as a nonlinear least squares problem and solved by a Levenberg-Marquardt method, adding the standard fairness energy ∥c  −1 − 2c  + c +1 ∥ 2 as a regularizer with a small weight that is set to 0 in later iterations.For details on this type of optimization, see [Tang et al. 2014].The planarity of quads, straightness of the unfolded strips and approximation errors for selected models are presented in Table 4.

Optimization of base meshes for strip structures
If we have more than one family of strips in a structure, it is an advantage if the strips are attached to the two or three families of major mesh polylines of a quad mesh or a triangle mesh (in case of a 3-web) with regular combinatorics.Initial versions of such a base mesh are extracted from level sets, but we need to strictly enforce above constraints per mesh polyline that defines a strip.Using the notation of Fig. 21, we discuss the case of a quad mesh and consider a vertex star with center v 0 and neighbors v 1 , . . ., v 4 .As in Eq. ( 1), we add tangents t 1 , t 2 for mesh polylines from edge vectors via t 1 = (e 1 − e 3 )/∥e 1 − e 3 ∥, t 2 = (e 2 − e 4 )/∥e 2 − e 4 ∥.
The unit normal n at v 0 satisfies constraints, To express a discrete P-curve with angle  , e.g. in tangent direction t 1 , we use a binormal b 1 and the constraints as well as the angle constraints Eq. ( 8). S Again a Levenberg-Marquart algorithm is used.The parameters and statistics for mesh optimization are presented in Table .3. Similar to [Schling et al. 2022], this mesh optimization can serve for interactive editing of structures while keeping the constraints.Examples for designing rectifying strip structures are provided in Section 4, and Section 5 extends the presented mesh optimization to torsion-free structures.

APPLICATION TO GRIDSHELLS
The construction of gridshells by bending straight flat lamellas as in [Schling et al. 2022] and [Mesnil and Baverel 2023] has been one of the motivations for the present research.These highly constrained structures are not easily designed.While optimization algorithms for design modification are known in certain cases [Schling et al. 2022], it is especially difficult to access the design space through nontrivial shapes or to approximate a given surface by such a structure.It is not surprising that the pseudo-geodesic gridshells in [Mesnil and Baverel 2023] are just rotational shapes; a main contribution in that work is the development of construction elements such as appropriate nodes.Since these are now available, one needs effective digital design tools.We show how to properly apply our optimization framework to the design of various types of gridshells.
To simplify writing, we use the following abbreviations: A=asymptotic, G=geodesic, P=pseudo-geodesic (to  ≠ 0, /2).Gridshells may consist of two or three families of elements.In the former case, discussed in Subsection 4.1, we have quad nets that get labeled as AG, PG, PP and so on.In the latter case (Subsection 4.2), the three families are arranged in a 3-web and labeled as AAG, AGG, PPG or similar.A Kagome pattern as in Fig. 1 arises from a 3-web by shifting the evaluation of the level-set function for one family.Table 1 provides an overview of the net/web types and their constraints discussed in our paper.
Table 1.Different types of rectifying strip structures discussed in our paper.The , , and  constraints are the constraints applied to ,  parameter lines, and the diagonals of the structures; surf condition refers to conditions that the reference surfaces should satisfy.TFN denotes gridshells with torsionfree nodes (Sec.5.2.3),PN is for gridshells with parallel node axes (Sec.5.2.2), and ACCS is for gridshells with a family of asymptotic curves of constant slope (Sec.5.2.4).

Gridshells from nets of pseudo-geodesics
The methods of Section 3 are well suited for the design of PP-nets.One needs to compute two families of P-curves which intersect at angles that do not become too small.Hence, the design of the two families in two separate steps is not ideal.Thus, we simultaneously compute both families from two level-set functions , .In order to constrain the angle  between the level sets of  and , we have to constrain the scalar product of normalized gradients g  = ∇ /∥∇ ∥ and g  = ∇/∥∇ ∥ on each triangle, where ⊥ is rotation by 90 degrees.Moreover, using our methods, we can turn asymptotic or geodesic curves into P-curves.Results are shown in Figures 22,23,24,25 and 26.(a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 )

Gridshells from pseudo-geodesic webs
It is well-known (see e.g.[Blaschke and Bol 1938]) that a 3-web on a surface can be represented by the level sets of 3 functions , ,  which satisfy This simple constraint makes the level-set formulation very attractive for the design of 3-webs on surfaces.However, the webs we are dealing with, do not have enough degrees of freedom to represent any shape.Hence, it is important that our framework can also optimize for the underlying surface; this is most easily done through post-processing of a web initialized via the level-set approach.Geometrically, this means discretely sampling the two families of asymptotic curves to obtain an AA-net in which one family of diagonal curves are geodesics, which slide nicely in between the other two families of lamellas and create an efficient triangulation of the web.This will in general require a shape change of , which we want to keep minimal.Hence, we initialize  = − − and then simultaneously optimize all three functions , ,  so that ( +  +  ) 2 is minimized, and Aand G-properties are fulfilled as well as possible; see Fig. 28.This optimization brings level sets of  closer to geodesics, but results in some deviation of the level sets of  and  from asymptotic curves.Hence, in the final step we use the mesh optimization of Sec.3.4 to obtain a high-quality AAG-web (Figs. 28 and 29).

AGG-webs.
An AGG-web has more degrees of freedom than an AAG web, since only the asymptotic curves are already determined by the input surface and even there one can select one of the two families.Here and for AAG-webs we assume that the given surface is free of flat points so that there is no singularity in the net of asymptotic curves.
We initialize the design of an AGG-web by selecting a family of asymptotic curves and first optimize to an AG-net with a weak angle constraint between the A-and G-family.Later a third level-set function is added simultaneously we optimize to the AGG property, allowing a change of the reference surface as well; see Fig. 30.

Webs involving pseudo-geodesics.
In an analogous way we can compute webs which contain P-curves.For construction reasons, there should be at least one G-family in the web; see Figs. 1 and 27.

TORSION-FREE NODES IN RECTIFYING STRIP STRUCTURES
So far, we did not address the fabrication of structures from rectifying strips.If one has two or more families of strips involved, a key issue is how to connect them at the nodes where they meet.
We now turn to this problem and present special structures which simplify the construction of nodes.These are torsion-free support structures that are well-known in architectural geometry [Pottmann et al. 2015].However, the restriction to rectifying strips is essential and makes the problem more interesting and challenging.Then, a web is extracted from the level sets (c); it is not yet an AAG-web with sufficient precision.Hence, mesh optimization is applied as a post-processing procedure, resulting in an accurate AAG-web (d) that exhibits some deviation from the reference surface.(e) shows the rectifying strip structure defined by the web.The evaluation of the AAG web in (e) is shown in Table 3.

Nodes and the local geometry of rectifying strips
The rulings of a rectifying strip with central curve c are given by the Darboux vectors d = t+b.They are equal to the binormals b only at points with vanishing torsion .However, in all practical realizations of asymptotic gridshells, one assumes a small deviation from binormals and builds nodes with binormals as axes (see inset).In view of tolerances at nodes, this is a realistic assumption, drawing also on some non-isometric deformation behavior of the material being used for the strips (timber, metal, certain plastics).In other words, one models the strip by the binormal surface rather than the rectifying developable.Proposition 1, whose proof is given in the Appendix, quantifies the deviation of the binormal surface from developability.
Proposition 1.The ruled surface (, ) = c() + b() formed by the binormals of a space curve c has Gaussian curvature In particular, on each binormal ( = const.),the point with maximal | | belongs to  = 0, i.e., it is the curve point c().Thus, the maximum deviation from developability of the binormal surface is along the curve c itself.Just for a planar curve c ( ≡ 0), the binormal surface agrees with the rectifying developable and is a cylinder with rulings orthogonal to c's plane.
It is also good to understand how much the rectifying developable  deviates from the binormal.This requires the computation of the normal curvature of  in direction b, which follows easily from Euler's formula, as shown in the Appendix.
Proposition 2. The rectifying developable  has vanishing normal curvature  1 = 0 in direction of the Darboux vector d.At c, its second normal curvature  2 and normal curvature   (b) in direction of the binormal b are where  := / denotes the conical curvature.
The intersection curve c  between  and the normal plane of c has curvature  2 at the curve point c().The 2nd order Taylor expansion of c  in the normal plane reads b+ 1 2  2  2 n  , showing that a point on the binormal, at distance  from c, has an approximate distance | | 2  2 /2 to the rectifying developable .
By Euler's formula, the normal curvature   (a) in any direction a of the rectifying plane, that forms the angle  with d equals   (a) =  2 sin 2  .Therefore, Eq. ( 20) allows us to compute the deviation of the line c + a from  near c.This is relevant if one wants to place the node axis in direction a.
The singular point r() = c() +   d() on a ruling of  is found with   = − 1  ′ (see e.g.[Strubecker 1964], p. 216), where  ′ is the derivative of  w.r.t.arc length.For our application, the singular points have to be outside the strip.Even being too close to the strip boundary is an undesirable situation for fabrication and aesthetics.
Summary and implications on fabrication.The binormals are close to rectifying strips for small torsion  only.If binormals are used as node axes on a strip along an asymptotic curve c of a reference surface , the value  2 =  2 / should be small as well.By the way, the torsion  of an asymptotic curve is related to the Gaussian curvature  of  via  2 = −.Large local changes  ′ of the conical curvature  imply that singularities come close to the strip.Conical curvature  is the curvature of the spherical curve t() traced out by the unit tangent vectors.Hence,  can easily be estimated from a discrete model.These practical considerations on the strips can be incorporated into optimization.

Torsion-free rectifying strip structures
A node where two developable strips intersect along a straight line segment is called torsion-free (see inset).Of course, that line segment must lie on a ruling of both strips.Torsion-free nodes provide advantages for fabrication.Thus we now focus on such nets of rectifying strips which possess only torsion-free nodes and call them torsion-free rectifying strip structures (see Figs. 32 and 31).We first address a way to compute them by optimization and then turn to specific cases which possess advantages in design and fabrication.Those also constitute good initial models for design modification through optimization.

Torsion-free structures through optimization.
The mesh optimization algorithm of Sec.3.5 is already equipped with everything needed for the computation of torsion-free structures.One only needs to make sure that the ruling vectors of rectifying strips agree at the nodes.To initialize the computation, one can start from a net of rectifying strips and intersect discrete rectifying planes at nodes.The resulting lines are then initial guesses for the ruling vectors at nodes.We express the common ruling vectors d  on the vertices v  ,  = 0, . . ., 4 of each mesh vertex star (see Fig. 21), such that where n 0, is a discrete principal normal attached to the edge e  .The first row of the equations expresses d 0 as the common ruling, and the second row ensures developability.The model in Fig. 31, left, has been computed in that way.This approach will work well and cause only small shape changes if the ruling vectors of the two strips at the nodes are nowhere too different from each other.Another way of modeling uses initial structures which already possess torsion-free nodes, as the ones discussed next.
5.2.2 Torsion-free rectifying strip structures with parallel node axes.
The mathematical problem underlying the construction of torsionfree rectifying strip structures is to find surface parameterizations x(, ) along which there is a 2-parameter family of straight lines (, ) = x(, ) + a(, ) with the property that the generated ruled surfaces along parameter lines  = const.and  = const.are the rectifying developables of the corresponding parameter lines on .Hence, a(, ) has to agree with the Darboux vector of both parameter lines of x(, ).This appears to be a difficult problem, and we are not aware of any study of it in the geometry literature.
In the following, we discuss the much more accessible special case of parallel node axes (a = d =  .), which already shows that the search for torsion-free rectifying strip structures attached to a surface  is not a shape restriction on , at least locally.
In the following we imagine the parallel node axes to be vertical and choose the constant node axis vector (Darboux vector) as d = (0, 0, 1)  .Then, each strip  must be cylindrical with vertical rulings and the curve c is a geodesic on .The rulings of a cylinder  appear as parallel lines in the planar unfolding and the geodesic as a straight line segment, leading to a constant intersection angle between geodesic and rulings.Thus, c intersects the vertical rulings of  under constant angle  and is thus a so-called curve of constant slope (see e.g.[Pottmann and Wallner 2001]).The angle  between curve tangent t and Darboux vector d satisfies cot  = / = .Thus, curves of constant slope are characterized by constant conical curvature .This is in agreement with above-mentioned properties: The curve t() is a circle and thus has constant curvature  and the singular points on  are at infinity due to  ′ = 0.The slope is defined via the inclination angle  against horizontal planes as tan  = tan(/2 − ) = .Slope 0 belongs to curves in horizontal planes.Principal normals of a curve of constant slope are orthogonal to the vertical cylinder  and thus horizontal; binormals thus have the constant inclination angle  against the vertical direction.

Computational design.
Even after fixing the node axis direction d, there are still infinitely many choices for the curve net on .However, a curve of constant slope  ≠ 0 cannot reach those areas of  whose tangent planes have a slope smaller than the chosen one.
The boundary of reachable points is a socalled isophote   on , whose tangent planes possess the fixed slope .Curves of the same constant slope  already form a net of curves, which ends with cusps on   , since in areas with steeper tangent planes there are two directions of given slope  (see the inset for an illustration).
There is no need to choose the same slope.Here our level-set method is very useful, as it can easily be adapted to curves of constant slope.One just has to replace the angle condition for P-curves by one which expresses a constant angle between unit tangents t of level sets (gradients in triangles rotated by 90 degrees) and d: Now tracing gets simpler and user-guided patterns via strokes is a good design tool as well; see Fig. 33 for a result.Likewise, mesh optimization is easily adapted to this case; one just has to express the ruling vectors as described in Eq. 21.

5.2.4
Structures with parallel node axes and one family of strips orthogonal to the reference surface.Cylindrical strips may form small  angles with the reference surface , which could be undesirable.Since both families of rectifying strips in a torsion-free structure cannot be orthogonal to  unless  is a plane, the question arises whether at least one family of them can be attached orthogonal to .The answer is affirmative.If cylindrical strips are orthogonal to  they must follow asymptotic curves on .Hence, we are looking for surfaces where one family of asymptotic curves are of constant slope with respect to the same reference direction, but slopes vary across asymptotic curves.All these surfaces have been determined in [Brauner 1968], resulting in an explicit parameterization that involves a solution  (, ) of the heat equation, We summarize our findings: Proposition 3. Torsion-free rectifying strip structures with parallel node axes are formed by the rectifying cylindrical strips of a net of curves of constant slope on the reference surface .If, in addition, the strips of one family are orthogonal to , they follow one family of asymptotic curves in a class of surfaces , which in an adapted coordinate system can be parameterized as in Eq. ( 23).
The case with  (, ) =  − (cos  +  sin  + 2 cos ) is the basis of the structure in Fig. 32.This periodic surface, studied in detail in [Brauner 1968], is also a translation surface.The generating curves are those in horizontal planes and the -lines (in vertical planes).This could be interesting for a watertight skin using flat panels.
We note the following constructive advantage of all surfaces in [Brauner 1968]: Since the principal normals of the asymptotic curves of constant slope are horizontal tangents of , the asymptotic curves are also curves of steepest descent, forming an orthogonal curve net with the horizontal curves.Using the latter for the second family of strips (as in Fig. 32) results in repetitive torsion-free nodes.At each node axis, the tangent planes of the two strips are orthogonal and one can build nodes, aligned with strip boundaries, that are congruent along each non-vertical strip.

SHADING SYSTEMS
A family of appropriately arranged strips, mounted on a facade or transparent roof, can function as a shading system.The topic of shading systems for freeform skins has previously been addressed by [Wang et al. 2013].There, a torsion-free support structure composed of individually shaped flat faces is used, which is harder to fabricate than the simple system we propose.The construction of our system amounts just to fixing the originally flat straight strips of bendable material appropriately on the underlying structure.We do not address here potential ways of fixation, but clearly it is helpful to have control over the rulings of the strips or their deviation from straight lines as discussed in Sec. 5. Our focus is on the challenging computational problem of designing rectifying developable strips along a freeform surface , so that they provide shading at chosen times of the day; here  can be a freeform architectural skin or an auxiliary design surface in front of a flat facade so that the shading system becomes a design element.We also consider the case where in some areas the system does not block the direct sunlight, while in other areas it does.In order to function as a shading system, the rectifying developable strips have to block the light for selected light directions.We assume these to be given by unit vectors l  , oriented from the earth towards the sun.Vectors l  describe a set of points on the unit sphere  2 .During a single day, the points l  lie on a circle ⊂  2 , whose rotational axis is parallel to the earth's axis.Its radius depends on geographic latitude and date.Independent of latitude, it is a great circle at the equinox.
The user can select the light directions for which shading should be provided via a region   ⊂  2 .To simplify its representation, we use a Cartesian (, , ) system whose -axis is parallel to the earth axis and whose -plane is parallel to the meridian plane on which the building is located, as shown in the inset.In this system, we use angles (, ) with l = (cos  cos , sin  cos , sin ).Within one year,  ∈ [−23.5 • , 23.5 • ], and  ∈ [−90 • , 90 • ] during a day at the equinox.In our implementation,   corresponds to a choice of intervals for  and  and thus the constraint l = (  ,   ,   ) ⊂   amounts to where  min = sin  min ,  max = sin  max ,  min = tan  min ,  max = tan  max .Inequality constraints Eq. ( 24) are incorporated into our optimization framework as equality constraints through dummy variables  .An inequality  ≥  is written as  −  −  2 = 0. We study two different effects w.r.t. the incoming light, namely blocking the light for shading or letting it through in certain parts.In both cases, we minimize an objective function of the form where   controls the position of rectifying planes of level sets w.r.t. the light.This depends on the purpose and is now discussed in each of the two cases, together with an appropriate initialization.
Shading.A rectifying developable strip along a curve c on the given surface  blocks light in direction l most efficiently, if c's rectifying planes are orthogonal to l.This is expressed with the energy   =  block , where t and b are unit tangent and unit binormal vectors, respectively.If l was a constant light direction, minimizing   to 0 would result in a straight line c orthogonal to l, which will in general not exist on .However, we apply the constraint with variable l = l  for each vertex v  ∈  so that l  ∈   .
Our pipeline for the design of shading systems is seen in Fig. 34.To initialize optimization, we take an even simpler approximation by representing the region   as an arc of a great circle.For buildings not too far away from the equator, this great circle can be a circular arc in our local -plane.For simplicity of description, we confine ourselves here to this case.A strip  orthogonal to this choice of light must have tangent planes parallel to the -axis and thus lie on a general cylinder with -parallel rulings.It follows from Sec. 5.2.2 that the strip curve c is of constant slope.Hence, initialization is done with a family of curves of constant slope.The simplest choice would be given by planar curves, which in the adapted coordinate system are the level sets of the function  restricted to .While those curves are easiest to get, they may not be appropriate, since we need their normals to be in the selected light directions to be blocked.It is therefore important to have the entire family of curves of constant slope as candidates for initialization, as shown in Fig. 35.
Assuming   to be an arc of a great circle, we compute in each triangle of the mesh the range of vectors that are orthogonal to one light ray from   .Then, we compute their slopes and select the most common slope value for initialization (see Figs. 34 and 35).
Not blocking light.One may want to let the light pass the system in certain areas or at specific times of the day.This lighting effect requires that the light direction is coplanar with the rectifying planes of strip curves c, i.e., l is orthogonal to their principal normal vectors.Hence, we use   =  light , where n  is the principal normal vector that satisfies As an initialization, we choose a central light direction l  ∈   , corresponding to  = ( min +  max )/2,  = ( min +  max )/2.Initial curves c shall have rectifying planes parallel to l  , which means that they are geodesics on cylinders whose rulings are parallel to l  .As before, they are curves of constant slope, but now the angle  is measured against l  .The present case is a bit easier than for shading, since the constant angle  is the only constraint to be imposed.
Fig. 36 shows examples for different lighting effects.The model in Fig. 37 has been split into various parts with shading systems of different function, and also includes a PP-net.
We can apply different shading/lighting constraints on the same model without splitting it into separate meshes, as shown in Fig. 38.We first selected the east side as the part where sunlight should pass through.Following that is a transition area where only smoothness is applied so that there is a fair transition into the remaining part that is optimized for shading.

DISCUSSION
Our algorithms have been implemented in C++ and tested on an Intel Xeon E5-2687W 3.0 GHz processor.Tables 2, 3 and 4 as well as Figs.40 and 41 provide details on the used parameters and performance for selected models of various levels of complexity.Our implementation is available on our project web page.28: A, A, G.Note that the shading energy (see Fig. 36) needs not to have a very small value, since in general there is no exact solution.| Limitations.Some limitations are of a fundamental geometric nature.
Heavily constrained structures such as AAG-or AGG-webs do not exist on general negatively curved surfaces.Thus our algorithms (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 ) (a 1 )   may or may not find a result, even when allowing small deviations from the reference surface.Here, the avoidance of combinatorial singularities plays a role as well.Driven by the aim to avoid selfintersections or optimization failures of level sets, we also ruled out combinatorial singularities (Fig. 42).They are however allowed for The north side of the building (bottom row) is a PP-net with P-curves of 60 • , so that we can see through.We set the location to Vienna and the time of the bottom row to 8:00, 12:00, and 14:00, respectively, on August 1st.The shadow on the ground shows that the desired effects are achieved.the mesh optimization in Sec.3.4.For practical applications, the composition of combinatorially regular parts towards a final design may be sufficient, as illustrated in Fig. 26.A restriction of the current implementation is of a topological nature.The reference surface must be a topological disk, possibly with holes (as in Fig. 14).We did not go beyond that for a good reason: On a closed surface of any topology, there may be no patterns of the required type that extend smoothly over the entire model.That

Fig. 5 .
Fig. 5. Twisted louvers are used in architectural facades to control lighting and ventilation.(a) Soma Architects worked together with the University of Stuttgart to develop a kinetic louver facade for the One Ocean Pavilion at the EXPO 2012 (https://soma-architecture.com).(b) The Wave Car-Park by Scott Brownrigg (https://scottbrownrigg.com).(c) Aluminium wrapping by Amano Design in Omotesando, Japan (https://amanod.com/jingumae).
Fig.6.Rectifying developable strip  along a space curve c and its straight unfolding  0 .The ruling d through a point of c is given by the Darboux vector d = t + b, where  and  denote curvature and torsion, respectively.

Fig. 9 .
Fig. 9.The discrete curve (red) in (a) satisfies Eq. (2), but the osculating planes (blue rectangles) are not continuous, meaning non-continuity of the discrete Frenet frames.The curve in (b) has continuous Frenet frames since it follows Eq. (3).The normal vectors of the surface are marked in yellow.

Fig. 12 .
Fig.12.Overview of our algorithm for generating a family of pseudogeodesic curves: (a) Starting with tracing guide curves, our method first initializes the scalar function on a parameter domain of the surface (b) and optimizes the scalar function so that its level sets interpolate the guide curves (c).Then, we optimize the scalar function so that the level sets are pseudo-geodesic curves (d).The quality of the extracted curves and rectifying strips is improved through post-processing procedures (e).The detailed evaluation for (d) is shown in Table2and Fig.40.The evaluation for (e) is shown in Table4.

Fig. 13 .
Fig. 13.Local operation of the tracing algorithm.By solving a quadratic Eq. (4), the next point p +1 of p  −1 and p  can be found.To improve stability, the edge containing p +1 (here v  v  ) is selected among the edges connecting v 0 or v 2 .
(a) (a) (a) (a) (b) (b Fig. 14.Examples of interactive design of P-curves.The red curves drawn on the surfaces (a) are used as guide curves to initialize the scalar function  (b).We then average the angles  of the level sets around the sketched curves to determine the target angles and further optimize for P-curves (c).The level sets in the top and bottom figures in (c) are P-curves of  = 60.3 • and  = 75.36• , respectively (for details on optimization weights and detailed evaluation, see Table Fig. 15.Reaction to failures.This triangle mesh is the same as in Fig. 11 b.(a) The scalar field is initialized with hand-drawn strokes.(b) Optimization towards P-curves with  = 66.5 • using the default parameters (as in Fig. 14) fails.(c) We start from  fair =  grad =  angle = 10 −4 , then gradually increase  angle to 3 and decrease  grad to 0. The lower  angle at the beginning of optimization helps to avoid singularities.The error  angle is successfully reduced to 6.69 × 10 −7 for  = 66.5 • .

Fig. 16 .
Fig. 16.Starting with an initial scalar field  (a), we compute the angles  between b and n on the vertices, and fit (,  ) (green dots) by a polynomial curve (red, top left).The angles  are far from the fitted curve and the level sets are far from P-curves, since for each value of  , the angles are not constant.The final optimized P-curves and the corresponding angles for each level set are shown in (b) and bottom left.The extracted rectifying developable strips are shown in (c).
angle has the same expression as  angle , but now values of  at mesh vertices v are fixed, while vertices v ∈ V, binormal vectors b and side vectors u are variables.  fair is the fairness energy, which contains 3 terms:   fair =    +   uni +    .The first term    is the squared Hessian energy of the ,  and  coordinates of the mesh vertices [Stein et al. 2018], providing smooth surfaces during optimization.  uni is the uniform Laplacian energy [Taubin 1995] which prevents vanishing triangle areas.The third term

Fig. 17 .
Fig. 17.The level sets in (a) are P-curves of  = 60 • .Using our triangle mesh optimization method described in Section.3.3, we can optimize the underlying triangle mesh as shown in (b) so that  = 45 • .Similarly, (c) shows P-curves of  = 75 • .They are optimized as geodesics in (d).Observe the shape differences in the reference surface.

First
version of rectifying strips.The discrete Frenet frame in (1) is attached to the vertices of a discrete curve .The discrete rectifying planes   : (x − c  ) • n   = 0 define a discrete developable surface  with intersection lines   :=   −1 ∩   as discrete rulings.Note that the edges of the polyline  do not lie on .Fig. 19 illustrates the geodesic property of  w.r.t. by the planar unfolding of , where the vertices of  are nearly aligned along a straight line.

Fig. 19 .
Fig. 19.The outer bisecting planes at vertices of a discrete curve (polyline) define a discrete rectifying developable; unfolding it together with the vertices of  results in nearly collinear points.

Fig. 20 .
Fig. 20.Another version of a discrete rectifying developable strip, formed by a sequence of planar quads, with a precisely straight unfolding.This requires the discrete rulings to lie in the outer bisecting planes (blue) of the central polyline .
20).The extension of one edge beyond c  forms the angle  −  −  =  +  with d  , showing that d  must lie in the outer bisecting plane   at c  , or d  must be orthogonal to the discrete principal normals n   ; With unit edge vectors e  this leads to the constraints e  − c +1 − c  ∥c +1 − c  ∥ = 0, n   − e  − e  −1 ∥e  − e  −1 ∥ = 0, d  • n   = 0. (16) Moreover rulings d  −1 and d  must lie in a plane   −1, through c  −1 , c  .Reflecting the plane   −1, at   yields the next plane  ,+1 .

Fig. 21 .
Fig. 21.Notation for optimization of a quadrilateral strip structure.To make the polyline through v 3 v 0 v 1 a P-curve, the angle  between binormal b and surface normal n has to attain the assigned value.

Fig. 22 .
Two rotational meshes (a 1 ) and (b 1 ) with the orthogonal PG-net of circles and profiles are isometrically deformed to meshes (a 2 ) and (b 2 ) using the method of[Wang et al. 2019].Orthogonality of the net and geodesics are preserved during an isometry.While the P-curve property of circles (to different inclination angles  ) gets lost during isometric deformation, it is restored by optimization, yielding orthogonal PG-nets (a) and (b); for clarity only the P-strips are shown.

Fig. 24 .
Fig. 24.PP-gridshell on a shape with varying signs of Gaussian curvature.Here, an AA-net would not be possible and GG-nets result in lower stiffness[Schikore 2023].
Fig. 25.A cone (a,b) is represented as a discrete orthogonal geodesic net [Rabinovich et al. 2018] with a combinatorial singularity of valence 3 at the vertex.Three straight diagonals (red) separate the net into three patches of regular combinatorics.Optimizing the geodesics to become P-curves of angle  = 50 • results in a structure (c) whose prototypical model is shown in Fig. 26.

Fig. 26 .
Fig. 26.Prototypical model built from timber veneer strips at an inclination angle of constant 50°, creating a triple symmetric cantilevering roof structure.The constant incline allows for repetition in the joints and created a beneficial orientation of beams to carry gravity loads.

Fig. 27 .
Fig. 27.Digital and physical model of a PPG-3-web.The lamellas are connected laterally through 3D-printed joints and steel bolts.The pseudogeodesic lamellas on top (orange) and bottom (red) have  = 60 • inclination.The geodesic lamellas (yellow) are fitted in between.
4.2.1 AAG-webs.The approximation of a negatively curved surface  by an AAG-web, illustrated in Fig.28, is based on the following thoughts.Since the web consists of two families of asymptotic curves, initial representations of these two families by level-set functions  and  are easily found.However, the level-set functions are not unique; any reparameterizations F =  ( ) and H =  ( ) represent the AA-net as well, but increase the degrees of freedom to find a third function  = − −  whose level sets are close to geodesics.

Fig. 28 .
Fig. 28.Approximation of a negatively curved surface by an AAG-web.We first optimize the scalar functions  and  for the two families of asymptotic curves separately, shown as first and second image in (a), and initialize the function  = − −  , third image in (a).'s level sets are not geodesics yet.Then, the 3 scalar functions together are optimized towards a lower overall deviation from the AAG-property (b).Then, a web is extracted from the level sets (c); it is not yet an AAG-web with sufficient precision.Hence, mesh optimization is applied as a post-processing procedure, resulting in an accurate AAG-web (d) that exhibits some deviation from the reference surface.(e) shows the rectifying strip structure defined by the web.The evaluation of the AAG web in (e) is shown in Table3.

Fig. 31 .
Fig. 31.Torsion-free rectifying strip structures, generated by optimization from a PP-net with torsion in the nodes (left), and by editing part of the model in Fig. 32 so that node axes are no longer parallel (right).

Fig. 32 .
Fig. 32.Torsion-free rectifying strip structure with vertical node axes.Strips of one family (bottom left) are orthogonal to the reference surface; the other strips follow curves in horizontal planes.Tangent planes (blue) at node axes are orthogonal (bottom right); nodes along the same nonhorizontal strip are congruent.

Fig. 33 .
Fig. 33.The level-set method is used to design two families of curves of constant slope (right) as the basis of a torsion-free rectifying support structure with parallel node axes (left).The intersection angle at nodes is 65 • ; the detail shows the node axes and the arrangement of strips in two layers.

Fig. 34 .
Fig.34.Given a reference surface and the considered sunlight directions (great circular arc) (a), we select a region for placing a shading system and compute the best slope against the earth axis (red arrows in (b) show their direction).This serves for initializing the scalar function  (c).After optimizing the levelsets (d) with consideration of the more specific sunlight region, the rectifying strips of the shading system (e) are obtained as in Sec.3.4.The time and location of the scene is August 1st, 15:00, Los Angeles.

Fig. 35 .
Fig. 35.Built on the Tropic of Cancer, the two families of strip curves (a) and (b) have a constant slope w.r.t. the earth axis.The curves (a) are planar ( = 90 • ), the curves (b) have  = 30 • .On the summer solstice, both (a 1 ) and (b 1 ) provide proper shading, but on the winter solstice, (b 2 ) is more effective than (a 2 ).Thus the slope helps to provide shading for a larger variety of sun positions.

Fig. 37 .
Fig. 37. Strip patterns with different functions on the same model.The roof (top right) and the south patch (top middle) are designed to block light, and the east patch (top left) is designed to let the light pass through the building.The north side of the building (bottom row) is a PP-net with P-curves of 60 • , so that we can see through.We set the location to Vienna and the time of the bottom row to 8:00, 12:00, and 14:00, respectively, on August 1st.The shadow on the ground shows that the desired effects are achieved.

Table 4 .
Rectifying strip quality evaluation to selected figures.|Strip| is the number of strips, max(diff  ) is the maximal deviation of  from the input binormal strips.max( straight ) is the ratio between the maximal deviation of an unfolded strip from a straight line and the strip length. planar is the total energy on quad planarity, expressed as  det(d  , d +1 , e  ) 2 , where the unit vectors d  , d +1 and e  are defined as in Eq. (17).To avoid optimization failures for strips defined by curves with large torsion and small curvature (see Sec. 5.1), our optimization for developable strips (Sec.3.4) can slightly change the base curve of a strip to obtain better rulings, leading to a possible loss of accuracy in  .This deviation is usually minor as shown by max(diff  ) values. is the ratio (in percentage) of ℎ to , with ℎ as Hausdorff distance between input surface and base curves of strips and  as bounding box diagonal.

Fig. 38 .
Fig.38.This shading system has been designed so that it allows light to pass through during the morning, but blocks it at noon and in the afternoon.

Fig. 39 .
Fig. 39.Outside view onto the shading system of Fig. 38.Less extremely shaped systems may not only serve for shading a flat facade, but also as a design element.The date and location of this scene are August 15th in London.

Table 2 .
Level set optimization statistics to selected figures.| | is the number of vertices, bbd is the diagonal length of the axis-aligned bounding box of the surface.The multiple  angle values of one model are for: Fig. 1:  1 . | | bbd  fair  edge  close  angle T/iter