# Spectral architectural geometry

Romain Mesnil, Kazuki Hayashi

► To cite this version:

Romain Mesnil, Kazuki Hayashi. Spectral architectural geometry. Computer-Aided Design, 2025, 188, pp.103927. 10.1016/j.cad.2025.103927. hal-05356066

HAL Id: hal-05356066

https://hal.science/hal-05356066v1

Submitted on 10 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Research Paper

# Spectral architectural geometry

Romain Mesnil<sup>a ID,*</sup>, Kazuki Hayashi<sup>a,b</sup>

<figure><img src="image_2.png" /></figure>

ARTICLE INFO

Keywords:

Architectural geometry

Spectral geometry

Discrete differential geometry

Shell structures

ABSTRACT

Spectral geometry is a mathematical field that links geometrical properties to eigenvalues of differential operators on surfaces. Although it is a well-established tool in geometry processing and has been used in many contexts, the structural engineering and architectural geometry communities have not yet adopted this framework for shape modeling. This paper aims to explore spectral methods for applications in architectural geometries. A novel methodology for generating anisotropic Laplacian operators based on regions of interest defined by the user is proposed. The potential of spectral methods in structural design is illustrated through design problems expressed on meshes and graphs.

1. Introduction

1.1. On architectural geometry

Free-form is an important trend in contemporary architecture and has fostered a significant amount of research in the last decades. It became apparent to the geometry processing community that the rationalization of the construction of free-form structures and envelopes featured many open research questions and was an interesting field of research [1]. The two main ingredients for architectural geometry are *discrete differential geometry* [2] and *numerical optimization*.

Discrete differential Geometry has led to a better understanding of Planar Quad (PQ) Meshes, and offset properties of meshes. Pottmann et al. have shown that the design space for discrete offsets yielding torsion-free nodes is a vector space [3]. The results have then been extended to non-standard meshes, for example, those based on semi-regular tessellations [4].

Concurrently with fundamental results on discrete differential geometry, several dedicated optimization schemes have been proposed to solve heavily constrained form-finding problems [5,6]. The main approach adopted in publications is to model a shape by solving an optimization problem, where each nodal coordinate of a mesh is a design variable [7]. In that sense, it can be compared to the “parameter-free” structural optimization paradigm [8,9]. The optimization approach is well-documented in practice but rarely applied by structural engineers or architects, who are not optimization experts.

Other researchers have promoted the use of restricted families of shapes fulfilling some geometrical constraints. Examples of this design approach are numerous, including the ruled surfaces used by Felix

Candela to build thin shell structures, translation surfaces or scale-
trans-surfaces [10], or shapes with planar lines of curvature for glazed
envelopes [11,12]. The main benefit of this bottom-up approach is the
mastery of the design constraints, with the disadvantage of a loss of
generality, making highly constrained problems difficult to solve.

1.2. *Current limitations*

Geometric exploration is made possible by numerical optimization, but this approach is rarely adopted in practice. Although several factors may explain this fact, the complexity of the architectural or struc- tural design task should be considered. Indeed, architectural design is a complex problem that involves the simultaneous consideration of multiple factors. Some factors are measurable, like cost, environ- mental footprint, lighting and thermal comfort. For example, Brown and Mueller have explored multi-criteria optimization of large span roofs, and parametrized their design space with six meaningful pa- rameters [13]. Other factors are intangible or boil down to designer preference. Mueller and Ochsendorf proposed a genetic algorithm to respect designer preference in a guided exploration [14].

In all those design scenarios, treating every nodal coordinate as a de-
sign variable is inappropriate due to its vast design space. Multi-criteria
optimization techniques scale poorly with the number of design pa-
rameters and cannot handle hundreds, if not thousands, of parameters
in conceptual design stages. Existing numerical methods considering
designer preference are based on metaheuristics, which also handle a
limited number of parameters [15,16]. Therefore, a current trend in
architectural design is dimensionality reduction to limit the number of

<footnote><sup>a</sup> Laboratoire Navier, UMR 8205, ENPC, UGE, CNRS, 6 et 8 avenue Blaise-Pascal Cité Descartes, Champs-sur-Marne, 77455, France</footnote>

<footnote><sup>b</sup> Department of Architecture and Architectural Engineering, Kyoto University, Kyotoodaigaku-katsura, Nishikyo, Kyoto, 6158540, Japan</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: romain.mesnil@enpc.fr (R. Mesnil), hayashi.kazuki@archi.kyoto-u.ac.jp (K. Hayashi).</footnote>
parameters and performance metrics to assess, which is no easy task, as observed in [17].

**Dimensionality reduction.** Dimensionality reduction is a field of machine learning; it is thus not surprising to see the adoption of state-of-the-art techniques from this field emerge in the last decade. Some authors have proposed Principal Component Analysis [18] to identify combinations of parameters that impact performance the most. Others have proposed techniques such as Self-Organizing Maps [19], or Bayesian inference [20]. More recently, some authors have focused on grammar rules to explore topological design spaces for complex structural patterns [21]. This exploration of meta-parametric design spaces prompts the question of efficient design space structuration to engage the designer in selecting appropriate topologies.

**Design robustness.** Sensitivity analysis can be included in a mathematical programming method to evaluate optimal criteria and obtain gradients analytically [22]. Analytical sensitivity analysis has been adopted for various structural optimization problems to boost optimization performance, including simultaneous optimization of geometry and topology of trusses [23,24]. However, the quality of the solution through gradient-based optimization is highly dependent on the initial solution, and optimization often converges at inferior local minima [25]. In addition, it is a known fact that treating independent nodal coordinates as design variables generates highly distorted shapes [26,27]. Therefore, reducing the number of design variables is key to successful optimization when using mathematical programming approaches.

### 1.3. Opportunity of spectral geometry processing

**Success of low-dimensional design space in practice.** There is, to our knowledge, no synthetic data on how engineers and architects explore parametric design spaces in real-life projects. The most documented examples favor low-dimensional design space, like gridshells by Jörg Schlaich and Hans Schober [10], or the famous courtyard of the British Museum. For the latter example, Chris Williams proposed to use only three functions to define the geometry of the roof [28], with each function controlling a key feature, like the curvature at the corners, or the rise of the gridshell.

**A signal processing point of view.** Generating a parametric design space for shape exploration is a long-established idea in the Computer-Aided Design community. In a seminal paper, Taubin proposed to adopt a signal-processing approach to geometry processing [29]. Akin to the Fourier analysis of a signal, Taubin showed that it is possible to construct a discrete version of the Laplacian that would serve as a basis for shape exploration analysis.

A multitude of operators can be used for spectral shape analysis, as reviewed in [30]. Particularly, the Laplacian operator is ubiquitous in the literature, as it plays a central role in geometry processing, due to its relation with the Dirichlet problem. The Dirichlet energy $\mathcal{E}$ associated with a scalar field $f$ on a surface $S$ is defined in Eq. (1).

$$ \mathcal{E}(f) = \int_S \|\nabla f\|^2 dS \qquad (1) $$

Loosely speaking, minimizing the Dirichlet energy results in the smoothest possible scalar field on a surface. This is equivalent to solving Laplace's equation.

$$ \mathcal{E} \text{ is minimal} \Leftrightarrow \Delta f = \nabla^2 f = 0 \qquad (2) $$

Previous studies have pointed out the limitations of spectral analysis in geometry processing due to the computational cost of computing eigenvectors. Some authors mention that thousands of modes on meshes, with several hundreds of thousands of vertices, are required to encode complex geometries correctly [31]. In addition, most mesh operations can be expressed as the solution of a linear system of equations, rather than an eigenvalue problem [32].

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Identification between structure and mesh data structure. Source: Adapted from a picture by Ivan Badura, distributed under CC BY 2.0 License.</figcaption></figure>

However, mesh processing in architectural geometry is distinct from conventional computer graphics applications. Indeed, the common practice in architectural geometry is to identify structural and envelope components with mesh attributes, as illustrated in Fig. 1. Common applications involve meshes of a few thousand vertices at most. For example, the roof over the Great Court of the British Museum spans over a rectangle of 73 m by 97 m and features only 1566 nodes [33]. The associated mesh is thus orders of magnitude smaller than meshes treated in Computer Graphics and does not pose complexity problems for spectral analysis.

### 1.4. Problem statement

This article discusses the potential offered by spectral mesh processing as a generic framework for modeling in architectural geometry. Although structural engineers are used to employ CAD-based parametrization of design spaces, very few of these approaches are compatible with fabrication constraints encountered in practice. The only attempts to do so [34,35] rely on bespoke descriptions of the parameter space and cannot be generalized, with the notable exception of [36].

Alternative operators derived from the Laplacian, such as the anisotropic Laplacian proposed by Andreux et al. [37] could open new design possibilities. However, most papers focus on shape analysis and do not consider user-defined features in the construction of operators. Nevertheless, we argue that designer preference should be considered in early design stages.

In order to demonstrate the relevance of spectral methods in architectural geometry, the following questions must be addressed:

* How to handle linear and affine constraints?

* How to modify the Laplacian shape operator to consider user-specified features?

* What design problems can be tackled?

The second section presents the modeling framework with spectral methods, including a novel methodology for defining ad hoc shape operators considering user-prescribed features. The third section focuses on the integration of linear and affine constraints. The fourth section discusses applications of spectral architectural geometry from structural optimization to surface parametrization.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Methodology for spectral exploration in architectural geometry.</figcaption></figure>

## 2. Methodology

### 2.1. Linear subspace exploration

In the following of this article, each structure is either modeled as a mesh, or a graph. The accepted modeling framework in architectural geometry is that edges correspond to members, faces correspond to panels, and vertices correspond to connectors.

Achieving a sound and efficient structural design requires tackling various design tasks:

* Definition of the nodal properties, like the position of each vertex.

* Definition of member properties, for example, the beam cross-section.

* Definition of the face properties, like panel thickness or opacity.

The proposed methodology in this article is depicted in Fig. 2. It expresses any design property from a linear combination of pre-computed scalar fields.

Formally, each occurrence of the design space is a linear combination of a restricted number of modal shapes ($v_i$), which are computed as eigenvectors of a matrix L, assumed to be diagonalizable. The vectors ($v_i$) should be understood here as a discretized version of scalar fields on a reference surface. Eq. (3) shows the formal expression of the combination, where $\alpha_i$ are scalar values chosen by the designer. A sufficient condition is, for example, that L is real symmetric.

$$ \mathbf{x} = \sum_{i=1}^{n} \alpha_i \mathbf{v}_i \qquad (3) $$

The scalars $\alpha_i$ are the parameters specified by the designer. The choice of the number of modes n thus defines the dimension of the design space. By definition, the eigenvectors are L-orthogonal:

$$ \forall i \neq j, \mathbf{v}_i \cdot \mathbf{L} \cdot \mathbf{v}_j = 0 \qquad (4) $$

The spectral theorem guarantees that the eigenvectors provide a basis for the space on which the operator acts. This design approach is thus exhaustive. The Laplacian, or discrete Laplace-Beltrami operators, have been chosen in numerous studies [30]. The proposed methodology can be applied to many design tasks, as scalar fields can be used to modify cartesian coordinates of vertices of a graph or a mesh.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Representation of scalar field u on a triangle.</figcaption></figure>

### 2.2. Scalar-field modeling

This article uses a piecewise linear description of scalar fields on triangulations, which is the most commonly used discretization. Fig. 3 displays a representation of a scalar field u on a triangle ABC. The scalar field of a point with barycentric coordinates $(\alpha, \beta, \gamma)$ is simply a linear combination of the nodal values, as written in Eq. (5).

$$ u(\alpha, \beta, \gamma) = \alpha \cdot u_A + \beta \cdot u_B + \gamma \cdot u_C \qquad (5) $$

Therefore, this discretization only requires nodal values of the scalar field. The same approach is adopted for graph modeling, where only nodal values of the scalar field are used. We denote the vector containing these nodal values as u.

**Gradient.** By definition, the gradient of piecewise-linear scalar fields is constant for each face. For a triangular mesh, the gradient in a face has a simple expression, depending on scalar field values u, edges vectors e, face normal vector N and face area A$_f$.

$$ \nabla u = \frac{1}{2A_f} \sum_i u_i (\mathbf{e}_i \times \mathbf{N}) $$

Notations are introduced in Fig. 4. When assembled for the whole mesh, the gradient of the scalar field is found by computing a matrix product Gu.

**Divergence.** The integrated divergence of a vector field X defined face-wise associated with a vertex i is given by:

$$ \nabla \cdot X = \frac{1}{2} \sum_j \cot \theta_1 (\mathbf{e}_1 \cdot X_j) + \cot \theta_2 (\mathbf{e}_2 \cdot X_j) $$

Notations are shown in Fig. 5. Similar to the gradient, the divergence is also a linear operator. It can be assembled with a matrix product DX, based on the previous equation.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Notations for computing gradient of scalar field **u**, after [38].</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Notations for computing divergence of vector field **X**, after [38].</figcaption></figure>

## 2.3. Modified Laplace-Beltrami operators

So far, we have mainly discussed applications of the Laplace-Beltrami operator in spectral mesh processing. However, other operators may be chosen, according to designers' need for editing specific regions of interest. To address this need, we extend the formulation of anisotropic Laplace-Beltrami operator $L_D$ initially proposed by Andreux et al. [37]. The Laplace-Beltrami operator is defined from the divergence (div) and gradient (grad) operators; the anisotropic Laplace-Beltrami operator is defined by introducing a 2 × 2 matrix A acting on tangent vectors.

$$L(u) = \operatorname{div}(\operatorname{grad}(u))$$

$$L_D(u) = \operatorname{div}(A \operatorname{grad}(u)) \qquad (6)$$

Andreux et al. proposed to construct A from principal curvature directions and showed that the obtained operator keeps key properties of the conventional Laplace-Beltrami operator: semi-definite positive-ness, symmetry, and the existence of a discrete spectrum. However, unlike the Laplace-Beltrami operator, their anisotropic operator is not isometry-invariant. While this extrinsic feature proves to be useful for shape analysis, this may be problematic for shape modeling. Meanwhile, there are various ways to construct alternative anisotropic Laplace-Beltrami operators using Eq. (6).

## 2.4. Construction of the anisotropic operator

### 2.4.1. Principle

As pointed out in [39,40], lowering the weight of specific vertices will increase their likelihood of being an extremum of the eigenvector. However, there seems to be little discussion on the appropriate way to assign such weights, especially a way to leave the intrinsic nature of the operator unchanged.

The simplest choice for the matrix A in Eq. (6) is the identity matrix $A = \alpha I_{2 \times 2}$, where $\alpha$ is a scalar value that may vary across the surface. It is then not difficult to see that the resulting operator can be expressed in its discrete form as the product between the true Laplacian and a diagonal matrix A, as shown in Eq. (7).

$$L_D = A \cdot L \qquad (7)$$

<figure><img src="image_8.png" /><figcaption>Fig. 6. Example of distance field $\phi$ computed with the heat method. The black dot is the source ($\phi = 0$) and the white curves are iso-values of $\phi$.</figcaption></figure>

When $\alpha > 0$, A is a diagonal matrix with strictly positive terms. It follows immediately that the matrix $L_D$ is symmetric positive semi-definite, as the product of two symmetric positive semi-definite matrices.

### 2.4.2. Construction of scalar field $\alpha$

The scaling factor $\alpha$ of the operator A may be computed by evaluating the geodesic distance $\phi$ between each vertex and the handle vertices H of the region of interest. Fig. 6 shows the geodesic distance field $\phi$ from a source point on a curved surface. This field has lower values next to the user-defined feature, which is one of the prerequisites for controlling modal shapes but is not directly usable. In this article, the scalar field $\phi$ is computed with the heat method, which conveniently uses the discrete Laplace-Beltrami operator [38].

We propose to define $\alpha$ by filtering the geodesic distance field $\phi$. The core idea is to lower the values of $\alpha$ near regions of interest while keeping comparatively larger values far from those. Let us write $\phi_{max} = \max(\phi)$, we consider a monotonic function $f$ so that Eq. (8) is satisfied.

$$\begin{cases} f(\phi_{max}) = 1 \\ f(0) = \epsilon \end{cases} \qquad (8)$$

The scalar $\epsilon < 1$ is the user-defined weight assigned to vertices in the region of interest. An example of such functions is given in Eq. (9), where R is the typical length of a region of interest (ROI).

$$\alpha = f(\phi) = \begin{cases} \frac{\epsilon}{2} & \text{if } \phi < 3R/4 \\ 1 - \frac{\epsilon}{2} \cos \frac{2(\phi - 3R/4)\pi}{R} & \text{if } \phi \in [3R/4, 5R/4] \\ 1 & \text{if } \phi > 5R/4 \end{cases} \qquad (9)$$

The smaller the value of $\epsilon$, the more likely local extrema occur in the region of interest. Another alternative is to use the following function defined in Eq. (10).

$$f(\phi) = \epsilon + (1 - \epsilon) \sin^k \left( \frac{\pi \phi}{2 \phi_{max}} \right) \qquad (10)$$

Fig. 7 shows the first modal shape for Laplacians defined in (7) with the weights defined in Eq. (10), with $k = 3$ and different values of $\epsilon$. The scalar field $\phi$ was obtained by specifying the outer and inner boundaries as source points, associated with the value $\phi = 0$. The initial geometry is the plane view of the courtyard of the British Museum, a rectangle of 73 × 97 m, with an inner circle of 44-meter diameter, shifted 3 m up north. The asymmetry of the initial shape yields a localization of high values in the south section of the courtyard for the first eigenvector of the true Laplace-Beltrami operator. The non-uniform operators tend to localize at the boundaries (low values of $\phi$): for $\epsilon = 10^{-3}$, high contrast of the eigenvector (visualized as high curvature of the height field) appears near the used-defined features. The first eigenmode for $\epsilon = 10^{-2}$ has non-zero values everywhere, and is closer to the shapes proposed
<figure><img src="image_9.png" /><figcaption>Fig. 7. First modal shape for several modified isotropic operators: top, front and perspective. The values of the eigenvectors are visualized both with color gradient and height field. The classical Laplace–Beltrami operator is on the left-hand side.</figcaption></figure>

by Chris Williams. This simple example illustrates that modifying the parameter $\epsilon$ allows for selecting relevant shape operators.

It shall be noted that the geodesic distance field $\phi$ is an intrinsic property of the shape, as it remains unchanged by isometries by definition. Therefore, the operators derived from the Laplace–Beltrami operator and $\phi$ are also intrinsic. This opens possibilities for several modeling approaches, discussed in the last section of this article.

## 2.5. Feature orientation

### 2.5.1. Principle

Feature alignment is an important component of shape and structural modeling. In structural applications, Li et al. have proposed a framework combining topology optimization and conformal maps to design efficient lattice structures [41]. In architectural geometry, Oval and coauthors have explored the design with quadrilateral patterns with boundary alignment properties [21,42]. The case of free edges in shell structures is particularly relevant since global stability is highly sensitive to the free-edge stiffness, as shown, for example, in Venuti et al. [43]. In practice, using a larger cross-section at the free edge or introducing localized curvature increases the geometric stiffness.

Tuning the operator’s anisotropy, for example, at free edges, seems thus like an interesting path. Formally, we assume that the designer can provide an input vector field aligned with anisotropic directions. For the sake of simplicity, we focus on vector fields that are gradients of geodesic distance fields $\phi$. We restrict the possible tensors $A$ to the form of Eq. (11), which is expressed in the orthonormal frame aligned with $(\nabla\phi, \nabla\phi^{\perp})$.

$$A = \alpha \begin{pmatrix} 1 & 0 \\ 0 & \beta \end{pmatrix} \qquad (11)$$

The parameter $\alpha$ plays the same role as in the previous paragraphs: it is considered a scaling parameter favoring local extrema of the eigenvectors in preferred regions. The parameter $\beta$ dictates the level of anisotropy; its effect is illustrated in the next section. The value of $\beta$ is supposed to be constant across the mesh in the following paragraphs, but it could be considered a variable.

### 2.5.2. Implementation

The anisotropic Laplace–Beltrami Operator is computed according to Eq. (6). The anisotropic matrix can be computed for each face: it requires a change-of-basis matrix $P_f$ from $(x, y, z)$ to the local frame of the face aligned with the guiding vector field.

$$A_f = P_f^\top \begin{pmatrix} 1 & 0 & 0 \\ 0 & \beta & 0 \\ 0 & 0 & 0 \end{pmatrix} P_f \qquad (12)$$

The anisotropic matrix is then obtained by the sequential matrix multiplication of a diagonal weight matrix $W$ whose diagonal entries correspond to nodal values of $\alpha$, the divergence matrix $D$, anisotropy matrix $A$ and gradient matrix $G$, as shown in Eq. (13).

$$L_D = W D A G \qquad (13)$$

Fig. 8 illustrates the influence of parameter $\beta$ on the first eigenvector, when $W$ is the identity matrix. The field $\phi$ is the geodesic distance field with $H$ on the boundaries. The value for $\beta$ decreases gradually from left to right, from $1$ to $10^{-5}$. The first eigenmode converges to a constant slope surface, whose gradient is aligned with $\phi$. Because $\phi$ is singular along the medial axis of the contour, the eigenvector exhibits gradient discontinuity at this location when $\beta$ becomes small. This simple illustration shows that aligning to design features is possible, but also that relative weight should be selected with great caution.
<figure><img src="image_10.png" /><figcaption>Fig. 8. Comparison of first eigenvector, visualized both as height field and color gradient for different values of β.</figcaption></figure>

There is a wide range of methods for computing eigenvectors of a matrix, some of which are reviewed in [44]. There are two families of methods: the ones aiming to retrieve all eigenvalues, and those which focus on determining only a handful of the smallest (or largest) ones. In the first case, the SVD of the matrix must be computed. The second case, which is of interest in our proposed approach, lends itself to several iterative methods. In this article, the eigenvalue problem is solved through the inverse power method, which results in the determination of the smallest eigenvalue. The implementation leverages the sparsity of $L_D$.

## 3. Modeling with spectral methods

### 3.1. Shape modeling

Now that the modeling background has been introduced, practical applications of spectral methods in architectural geometry shall be discussed. This section exemplifies shape modeling, including constraints that designers may often encounter.

The most evident application of spectral methods in the modeling of free-form structures is to interpret nodal values as nodal displacements. With this interpretation, the $(x, y, z)$ coordinates can be expressed as linear combinations of three families of modal shapes $(u, v, w)$:

$$ \left\{ \begin{aligned} x_i &= \sum_{j=1}^{N_x} \alpha_j u_j^j \\ y_i &= \sum_{j=1}^{N_y} \beta_j u_j^j \\ z_i &= \sum_{j=1}^{N_z} \gamma_j w_j^j \end{aligned} \right. \qquad (14) $$

In the following, the eigenvectors are represented as height field. Several form-finding processes make a distinction between z coordinates and the two others. Indeed, many architectural surfaces can be interpreted as a height field. Dedicated techniques have thus been adopted to design height-field surfaces. In architectural geometry, the *marionette technique* solves a linearly constrained set of equations to generate PQ-meshes [45]. In structural mechanics, the *thrust network analysis (TNA)* makes a distinction between horizontal and vertical equilibrium, also resulting in the solution of a linear system where z-coordinates are the only parameters [46].

The rest of this article looks at the handling of constraints. Without delving too much into the topic of optimization, we may recall that there are two ways to handle constraints, either as *soft constraints* which are enforced by penalization, but can be violated, or as *hard constraints* which cannot be violated. In the following, we limit our discussion to hard equality constraints, as the methods used with soft constraints are identical to the ones used in the unconstrained case.

### 3.2. Linear constraints

In practice, a designer has to impose certain constraints for design exploration. Eq. (15) is the formal definition of the generalized eigenvalue problem with linear equality constraints. Solutions of this problem satisfy $Lx = \lambda Bx$. With the classical cotan-operator, the true Laplacian $L$ is written as the product of the conformal laplacian and the inverse of the area matrix $A: L = A^{-1}L_C$. The conformal Laplacian is symmetric whereas the true Laplacian is not, the solution of the generalized eigenvalue problem with $B = A$ improves computational performance.

$$ \begin{array}{ll} \min & x^\top L x \\ \text{subject to} & x^\top B x = 1 \\ & C x = 0 \end{array} \qquad (15) $$

Poranne [36] proposed to solve the constrained eigenvalue problem by constructing a basis $N$ of the nullspace of $C$. Any vector satisfying $C \cdot X = 0$ can be written as $x = Nw$. Most numerical methods that compute the null-space of a matrix, like Francis' QR method, result in an orthonormal basis, meaning that $W^\top N^\top N W = W^\top W$. The constrained eigenvalue problem then becomes a classical eigenvalue problem:

$$ \begin{array}{ll} \min & w^\top N^\top L N w \\ \text{subject to} & w^\top N^\top B N w = 1 \end{array} \qquad (16) $$

The main drawback of this approach is that the matrix $N$ is usually not sparse, which makes the eigenvalue problem much more computationally expensive. Instead, augmented matrices $L^+$ and $B^+$ defined in (17) can be used for the generalized eigenvalue problem.

$$ L^+ = \begin{pmatrix} L & C^T \\ C & 0 \end{pmatrix}, \quad B^+ = \begin{pmatrix} B & 0 \\ 0 & 0 \end{pmatrix} \qquad (17) $$

The generalized eigenproblem is written in Eq. (18). It can be proven that this problem is equivalent to (15) by simply writing the Lagrangian of (15) and by writing the optimality conditions. $\lambda$ is the Lagrange multiplier associated with the constraint $x^\top B x = 1$, while $y$ is the Lagrange multiplier associated with the linear constraint.

$$ L^+ \begin{pmatrix} x \\ y \end{pmatrix} = \lambda \cdot B^+ \begin{pmatrix} x \\ y \end{pmatrix} \qquad (18) $$

The $L^+$ matrix is invertible when $C$ is full-rank, but not necessarily positive definite. The Inverse Power Method is thus still a valid approach to compute the first eigenmodes. In our implementation, a QR factorization is used instead of a Cholesky factorization, which results in a computational overhead that remains acceptable.

### 3.3. Assembly of the constraint matrix

The eigenvalue problem is solved far more efficiently when the matrix $C$ is full-rank. Rather than assembling a constraint matrix and determining its column space a posteriori, we directly assemble a full-rank matrix with the following method:
<figure><img src="image_11.png" /><figcaption>Fig. 9. Creation of affine subspace solution by addition of homogeneous eigenvector and the minimal norm particular solution. Here, the amplitude of the scalar field is visualized by moving the nodes in the vertical direction.</figcaption></figure>

1. The user defines constraints as a set of vectors ($\mathbf{v}_k$), so that $\mathbf{v}_k^\top \cdot \mathbf{x} = 0$

2. Normalize the vectors ($\mathbf{v}_k$)

3. Select a tolerance $\epsilon$. For an implementation with double precision, a value of $10^{-7}$ is satisfying.

4. At step $k$, compute the projection on the $h$ column vectors of $\mathbf{C}$, written as ($\mathbf{c}_i$):

$$ \tilde{\mathbf{v}}_k = \mathbf{v}_k - \sum_{i=1}^{h} (\mathbf{v}_k^\top \cdot \mathbf{c}_i) \mathbf{c}_i $$

5. If $||\tilde{\mathbf{v}}_k|| > \epsilon$, then write $\mathbf{c}_{h+1} = \mathbf{v}_k$.

6. Iterate on $k$

The algorithm recalls the Gram-Schmidt algorithm, however, we do not use the vector $\tilde{\mathbf{v}}_k$, so that the column vectors of $\mathbf{C}$ are not orthogonal. The main reason for this approach is to preserve the sparsity of the matrix $\mathbf{L}^+$.

This approach may seem computationally expensive at first glance. However, we shall point out that constraints usually consist of sparse vectors, and the projection cost is proportional to $h \cdot m$, where $m$ is the average number of non-null components in the constraint vectors. The complexity of the whole procedure is thus $m \cdot r^2$, where $r$ is the rank of the matrix $\mathbf{C}$.

### 3.4. Linear constraints with non-null vector

Consider now a linear constraint with constant terms, in the form $\mathbf{C} \cdot \mathbf{X} = \mathbf{c}$. The minimization problem of Eq. (15) can be rewritten and re-arranged owing to Lagrange multipliers

$$ \mathbf{L}^{+} \cdot \begin{pmatrix} \mathbf{x} \\ \mathbf{y} \end{pmatrix} = \lambda \cdot \mathbf{B}^{+} \cdot \begin{pmatrix} \mathbf{x} \\ \mathbf{y} \end{pmatrix} + \begin{pmatrix} \mathbf{0} \\ \mathbf{c} \end{pmatrix} \quad (19) $$

This problem is no longer an eigenvalue problem, and the solutions' space is no longer a vector space. Indeed, it can easily be seen that multiplying a solution of Eq. (19) by a scalar $\mu$ does not result in another solution. We suggest not solving Eq. (19), but instead computing the smallest solution of $\mathbf{C}\mathbf{x}^* = \mathbf{c}$. Then, we propose to solve the homogeneous eigenvalue problem of Eq. (18) and shift those eigenvalues by $\mathbf{x}^*$. Fig. 9 shows a particular example inspired by the British Museum, with Dirichlet boundary conditions (different altitudes for the outer and inner supports). Although this is not equivalent to the problem of Eq. (19), this approach maintains the vector space structure of the solution space, which seems desirable.

### 3.5. Relevant constraints in architectural geometry

**Dirichlet boundary conditions.** The Dirichlet boundary condition corresponds to the prescription of a value for a given node of index $i$.

$$ X_i = 0 \quad (20) $$

If the problem features only positional constraints, it is better to use the reduced Laplacian matrix [32], which consists of removing the $i^{th}$ rows and columns. Fig. 10 shows the modal shapes with the Dirichlet boundary conditions. The boundary is the one of the British Museum, described in [28].

**Neumann boundary conditions.** The Neumann boundary condition is a constraint on the gradient of the solution. In plate theory, this corresponds to constraints on rotations.

$$ \nabla \mathbf{X}_k \cdot \mathbf{n}_k = 0 \quad (21) $$

The gradient of a scalar field on a triangular mesh is written as follows, where $N$ is the face normal.

$$ \nabla X = \frac{1}{2A_f} \sum_i X_i (\mathbf{N} \times \mathbf{e}_i) \quad (22) $$

We may write $\mathbf{n} = \mathbf{N} \times \mathbf{e}$, where $\mathbf{e}$ is the tangent vector of the clamped edge. The condition (21) simplifies when $N$ is a unit vector.

$$ \sum_i X_i (\mathbf{e}_i \cdot \mathbf{e}) = 0 \quad (23) $$

**Periodicity and symmetry.** Symmetry may be an important feature to preserve during shape exploration. A symmetrical shape is inherently easier to construct due to the repetition of components. In a design context, we should expect symmetry to be enforced by the designer. In fact, this is the case in many studies on shape optimization of free-form structures [34,47,48]. In addition, optimization techniques seem to result in structures that maintain the symmetry of the geometry and loads [49]. We may provide a simple explanation for this fact. Most free-form structures are designed for buckling and are sensitive to geometrical imperfections. One of the most common heuristics is to modify the shape by adding the first buckling modal shape $x'$ which is generally non-symmetrical, as illustrated in [50].

Symmetry of eigenmodes can be enforced by restricting the problem to one half and imposing Neumann boundary conditions on the plane of symmetry, as shown in Fig. 11. Notice that only Neumann boundary conditions are applied and that the nodal position is not directly constrained. The remaining mesh can then be generated by applying symmetries.
<figure><img src="image_12.png" /><figcaption>Fig. 10. First two eigenmodes for the unconstrained (dashed) and constrained Laplacian with Dirichlet boundary conditions. Colors represent the magnitude of the eigenvector.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Boundary conditions for a symmetric mesh. The nodes on the axis of symmetry (purple curve) are clamped.</figcaption></figure>

Fig. 12 shows an example where symmetry is enforced. The operator is a Laplacian with the weighting proposed in Eq. (10) with $k = 3$ and $\epsilon = 10^{-2}$. Nodes on the inner ring and outer rectangle are fixed (Dirichlet homogeneous boundary conditions). The first eigenmodes are not symmetric. The enforcement of Neumann boundary conditions on the axis of symmetry before the generation of the symmetric mesh guarantees the smoothness of the result, illustrated by the continuity of contour lines and their perpendicularity to the axis of symmetry.

**Boundary value.** Another common constraint is setting a scalar field's value to be constant on a boundary $\partial\Gamma: X_i = c, c \in \partial\Gamma$. This constraint can be expressed as a linear constraint. Consider a boundary edge $e_{ij}$ from vertex $i$ to vertex $j$. The condition

$$x_i - x_j = 0 \qquad (24)$$

**Discrete offsets.** The other rationalization objective often treated in architectural geometry literature concerns discrete offsets of meshes. Indeed, it is preferable to construct torsion-free beam layouts, where beams are straight and without twist, meeting cleanly at nodes. Pottmann et al. have shown that this problem is equivalent to the construction of discrete parallel transformations, which is formulated as a linear problem [3]. The potential of this design space for support structure rationalization has been illustrated recently in [51]. Montagne et al. have already proposed to perform eigenmode decomposition to simplify the exploration of the linear space of parallel meshes [35].

<figure><img src="image_14.png" /><figcaption>Fig. 12. First three eigenmodes for a weighted Laplacian with (bottom) an without (top) enforcement of symmetry. The axis of symmetry is oriented along the vertical direction. White curves are contour lines of the eigenvectors.</figcaption></figure>

**PQ meshes.** One of the most common rationalization objectives in architectural geometry is to generate meshes with planar quadrilateral (PQ-meshes). Vaxman first proposed to study PQ-mesh generation through affine transformations of faces, with a linear compatibility constraint [52]. Poranne et al. [36] extended those results to other face transformations, followed by Mesnil et al. [53].

<figure><img src="image_15.png" /></figure>

There are several ways to enforce those constraints. A first direction could be to use quad meshes and adopt the polygonal Laplacian, [54]. However, the anisotropic Laplacian proposed in this article should be
<figure><img src="image_16.png" /><figcaption>Fig. 13. The first four eigenmodes of a mesh with a triangular plan with corner nodes fixed and highlighted in red. The planarity error is computed by the method of [52]. Enforcing planarity across all faces results in less smooth eigenmodes. In contrast, disregarding the planarity constraints for faces connected to the fixed nodes preserves smoothness comparable to the unconstrained case.</figcaption></figure>

adapted to quadrilateral meshes. Another option could be to triangulate the quadrilateral mesh with “shadow” edges represented as dashed lines in the Figure and used only for the computation of the operators, which is the approach used in the paper.

A practical example with the Marionette technique was implemented. With this method, the only variables are the z-coordinates of the vertices. It is thus sufficient to compute a single scalar field for shape editing. The planarity constraint for a quad ABCD reads as follows:

$$ ((d_{DB} + d_{BC} + d_{DC}) - d_{BC} - d_{DB} - d_{DC}) \cdot \begin{pmatrix} u_A \\ u_B \\ u_C \\ u_D \end{pmatrix} = 0 \quad (25) $$

where $d_{BC} = \det_{2D}(AB, AC)$, $d_{DB} = \det_{2D}(AD, AB)$, and $d_{DC} = \det_{2D}(AD, AC)$, as shown in [45]. Clearly, there is one planarity constraint for each face, which involves potentially ill-behaved behavior when overconstraining the eigenvalue problem.

Fig. 13 compares eigenmodes with and without marionette constraints, with Dirichlet boundary conditions. The eigenshapes without the marionette constraint of Eq. (25) are non planar and are not fit for applications to architectural envelopes. Adhering to PQ constraints on all quadrilateral faces leads to eigenshapes that differ significantly from the unconstrained case. Meanwhile, disregarding planarity of faces adjacent to fixed nodes yields smooth eigenmodes similar to the unconstrained solution. The planarity constraint allows real-time explorations of fabrication-aware design spaces compatible with manual modification, or optimization similar to preliminary work in [55].

# 4. Applications

## 4.1. Implementation

The implementation relies on the development of a C# library and of a dedicated plug-in for Grasshopper3D, a visual programming language that runs within Rhinoceros 3D CAD software, which is popular across the architectural engineering community. This implementation allows seamless integration with other engineering software, such as

finite element analysis, Computer Fluid Dynamics, daylight simulation, or robotic control. The design workflow consists in pre-computing eigenmodes, and to combine them. Combination can be done interactively, through slider manipulation, or automatically, through numerical optimization. The first scenario corresponds to the use-case of the British Museum, already mentioned in the literature review.

Grasshopper3D is a visual programming language which allows the designer to connect components in a canvas. Results can be visualized simultaneously in the CAD environment Rhinoceros3D. Fig. 14 displays a minimal Grasshopper3D Canvas which features some components allowing to perform spectral analysis. In this image, an input mesh is provided by the user. The components in the white groups allow to gather the input for the construction of the anisotropic operator. A component computes the geodesic distance field from a source point, here identified by its index value in the mesh. A second component computes its gradient, as a collection of 3D-vectors. The component on the far right of the image, in the blue group, solves the generalized eigenvalue problem. It outputs a collection of eigenvectors as a data tree. These eigenvectors can later be combined with basic operations.

Fig. 15 displays a canvas including the post-processing of the output. The eigenvalue problem results in eigenvector with a unit B norm. A re-normalization with infinity norm is performed (label Normalization (optional) in the image). Modes are then combined with number sliders (label scaling factors), which correspond to the numbers ($\alpha_i$) in Eq. (3). On the far right of the image, the resulting design vector x is used to modify nodal positions and create a new mesh from them.

## 4.2. Thickness and sizing

Notice that because in the absence of constraints, (1...11) is the first eigenvector of the true Laplace-Beltrami operator, it is always possible to enforce the positivity constraint. Consider controlling the thickness of a shell using the Laplacian. In this example, the scalar field is considered to be a parameter to define face-wise thicknesses.

Once the scalar values are obtained vertex-wise, the attribute of face j, written as $y_j$, is determined as the mean of the corner attributes.
<figure><img src="image_17.png" /><figcaption>Fig. 14. Grasshopper3D Canvas with components allowing the computation of eigenvectors.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15. Full Grasshopper3D Canvas including combination of eigenvectors and their visualization as a height field.</figcaption></figure>

Assuming that the specified initial volume is $\bar{V}$ and area of element j is $\bar{A}_j$, the thickness $t_j$ is determined as

$$t_j = y_j + \frac{(\bar{V} - \sum_j \bar{A}_j y_j)}{\sum_j \bar{A}_j} \quad (26)$$

The thickness computed by equation (26) ensures that the total structural volume remains $\bar{V}$.

To enhance the overall structural performance, the total linear strain energy is typically chosen as the objective function to be minimized. The total strain energy is calculated through structural analysis using discrete Kirchhoff triangle (DKT) elements. For shell structures, the thickness optimization problem is inherently nonlinear due to the coupling of membrane forces and bending moments, thus requiring nonlinear solvers. In the subsequent non-linear optimization examples, we employ the Sequential Least Squares Programming (SLSQP) algorithm in the SciPy library. SLSQP is an iterative method that addresses

nonlinearity by solving a sequence of sub-problems, in which the objective function is quadratically approximated and the constraint functions are linearly approximated with respect to the current variables.

The thicknesses of finite element shells can be optimized using Eq. (26), as depicted in Fig. 16. Here, the weights of the first $n_{mode} \in \{1, 5, 10\}$ eigenmodes are used as design variables. The optimization includes constraints to ensure the thickness values are positive. Smooth thickness distributions that minimize the strain energy can be obtained for the three cases.

### 4.3. Shape optimization

Another natural extension of our spectral framework is shape optimization of shell structures, regarding the per-vertex scalar field as nodal movements. Previous studies demonstrated that varying nodal coordinates individually often falls into local optima with corrugated
<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Thickness optimization results of a finite element shell with Young's modulus of 1.0 [N/m²]. The strain energy is minimized by adjusting the weights of $n_{\text{mode}}$ modes. Although the strain energy decreases with increasing $n_{\text{mode}}$, the contribution of higher eigenmodes is small, which is evident in the small difference between the cases $n_{\text{mode}} = 5$ and $n_{\text{mode}} = 10$.</figcaption></figure>

shapes [26]. This subsection demonstrates that spectral shape optimization is effective in obtaining smooth optimal shapes at a low computational cost.

Suppose a planar triangulated mesh shown in Fig. 17. This mesh contains 684 nodes, including 8 pin-supports; therefore, the degree of freedom is 676. The nodes are subjected to a 1 kN downward load. Assuming Young's modulus of 25000 N/mm², Poisson's ratio of 0.15, and a constant thickness of 20 cm across the shell, we aim to minimize its total linear strain energy under constraints on the total structural area and the non-negativity constraint for nodal heights. Optimization starts from a smooth height field obtained from the 1st eigenmode with its magnitude randomized with a random seed, and the maximum allowable mid-surface area is set to 1.2 times the projected area.

Fig. 18 highlights a key difference in optimization outcomes: independently varying nodal heights results in a highly corrugated shell, while our spectral approach, which expresses nodal heights as a weighted sum of low-frequency eigenmodes, yields a smooth optimal shape. Interestingly, the optimal shape, when utilizing all eigenmodes, showed little corrugation, even with the same number of design variables as individually varying the heights of free nodes.

Beyond shape smoothness, the spectral optimization approach also boasts considerable computational efficiency, especially when utilizing a limited set of low-frequency eigenmodes. Fig. 18 effectively depicts the critical trade-off between solution quality and computational efficiency as the number of eigenmodes used for optimization is varied. It is important to observe that solution quality plateaus after approximately

100 eigenmodes, thereby confirming the efficiency of focusing on low-frequency eigenmodes for structural optimization.

### 4.4. Explorations of graphs

Scalar fields can be utilized to design discrete structures represented as graphs, such as frames and trusses. Manipulating the geometry of the discrete structures based on smooth scalar fields has significant advantages in avoiding non-smooth shapes. Fig. 19 shows the first symmetrical eigenmodes for a graph representing a Pratt truss. The vertices of the bottom chord are treated as fixed points (Dirichlet boundary condition). Here, eigenmodes are treated as vertical displacement of the graph vertices. Notice that the first eigenmode corresponds to a uniform motion of the top chord, while further modes correspond to sinusoidal undulations.

Fig. 20 shows the assembly of a curved space frame using two scalar fields. One scalar field represents the nodal positions on the bottom layer, and the other scalar field represents the offsets of the nodes on the upper layer from the bottom layer.

The shape optimization of a space frame can be implemented by combining these two types of scalar fields and utilizing their few eigenmodes. In Fig. 21, the strain energy of a space frame with bar members of uniform cross-section is minimized by varying the weights of a pair of the first eigenmodes. This approach enables the attainment of a smooth optimal shape, avoiding inferior local optima.

### 4.5. Results on constrained anisotropic shape modeling

Fig. 22 shows an example of shape modeling using both the Laplace-Beltrami operator and an anisotropic Laplacian. The motivation was to replicate the free-edge feature observed in some designs by Heinz Isler. To do so, the anisotropic Laplacian was constructed from the scalar field $\phi$ with $\phi = 0$ on free edges. Since high localization was wanted, the filter proposed in Eq. (9) has been used, with $\epsilon = 10^{-3}$. Notice that the first eigenmodes of the anisotropic operators have very sharp features localized at the free-edge, akin to Isler's design.

The relevance of those features for mechanical performance is well documented in the literature. For example, parametric studies show that introducing high localized curvature increases the buckling load capacity of thin shell structures [34]. A simple sensitivity study was performed on the shells shown in 22, modeled with the finite element method as concrete shells (C40/50) of 4 cm thickness. The linear buckling load factor was increased by 68% and 71% respectively for the combination with first and third anisotropic eigenmodes while keeping total mass constant. This illustrates that meaningful design input may drastically improve overall performance and outperform conventional methods.

## 5. Discussion

### 5.1. Computational performance

Computational efficiency is one of the objectives of spectral methods in architectural modeling. We argue that a few eigenmodes (typically less than ten) are sufficient to describe complex shapes in many situations, which makes design space exploration possible with many computational techniques. In particular, subsequent optimization problems solved with gradient descents are expected to be solved in real-time.

The computational overhead for determining eigenmodes is also limited, due to the relatively small size of meshes used in architectural geometry. For example, the computation of 10 eigenmodes for a mesh with 3074 vertices and 366 linear constraints (see Fig. 12) takes on average 850 ms with Intel(R) Core(TM) i7-8650 CPU processor with 32 GB of RAM. The computation relies on our own C# implementation of the Inverse Power Method, which could be improved by using dedicated linear algebra packages using Arnoldi methods.
<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Initial mesh whose nodal heights are to be optimized to minimize the total strain energy (Left). Varying nodal heights individually often results in a highly corrugated local optimum, including bumps and dents (right).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Trade-off between the objective function and required optimization time when varying the number of eigenmodes used for optimization $n_{mode}$. Optimization was run 10 times with different random seeds, and the average results are plotted. For $n_{mode} = 1, 10, 100$, and 676, the optimal shapes with the least total linear strain energy are provided. Note that we incorporated sensitivity analysis, which analytically computes the derivatives of the objective function with respect to the design variables, into the SLSQP solver to improve its efficiency.</figcaption></figure>

<figure><img src="image_22.png" /></figure>

Reference graph

1st eigenmode

2nd eigenmode

3rd eigenmode

4th eigenmode

<caption>**Fig. 19.** First symmetrical eigenmodes for a graph representing a Pratt truss.</caption>

## 5.2. Interest of the method

The proposed workflow allows the definition of geometry-aware dimensionality reduction technique for parametrizing design spaces in architectural engineering. As demonstrated, many design criteria related to architectural geometry may be expressed with linear constraints and preserve the vector space structure of the design space. The extension of the anisotropic Laplace-Beltrami operator with predefined regions of interest considerably enhances control of the design space, with a low computational cost. The technique may be combined with other dimensionality reduction techniques mentioned in the literature for more efficiency.

<figure><img src="image_23.png" /><figcaption>**Fig. 20.** The shape of the space frame can be constructed by combining two scalar fields: one represents the rise of the bottom layer, and the other represents the offsets from the bottom layer.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. 21.** Initial shape (top) and the shape optimized using a pair of the first eigenmodes (bottom). A length constraint is applied so that the total member length is 1.2 times that of the initial flat state.</figcaption></figure>

### 5.3. Extensions

The proposed examples were computed with the *extrinsic* Laplacian which is constructed from nodal coordinates. It is well-known that it is possible to use the *intrinsic* Laplacian from an intrinsic Delaunay triangulations [56]. Recent advances in discrete differential geometry show that intrinsic mesh simplification could be performed to further improve computational efficiency when computing discrete Laplace-Beltrami operators [57]. Our explorations are restricted to graphs, or manifold meshes. It was shown in [58] that Laplacian can be extended to nonmanifold meshes, which would extend the applicability of the method proposed in this article.

### 5.4. Other applications

Applications of spectral methods are not necessarily limited to scalar field design. Vector fields and n-RoSy (Rotationally Symmetric) fields are useful means to construct discrete surface parameterizations. 2-RoSy fields are the natural framework to describe quadrangulations, or properties like principal curvatures or principal stress in a shell structure. Spectral methods for quadrangulations is a mature topic, as they have been proposed in early work by Dong et al. [40]. The exploration of quadrangulations in architectural and structural design has been proposed by Oval et al. [21], with discrete rules.

In architectural geometry applications, Boutillier et al. have proposed a framework for vector field generation as a gradient of scalar fields in [59]. Their approach consists of specifying Dirichlet boundary conditions at boundaries of a design domain, computing a minimal surface and creating a quadrangulation from the gradient and contour lines of this surface. Our spectral method handles Dirichlet boundary conditions and could be used to generate linear combinations eigenmodes. The main benefit of our approach is the reliance on smooth rather than discrete parametrization, which would open the door to many applications (gradient-based descent or sensitivity analysis).

However, restricting the search space to gradients of scalar fields does not capture the full extent of design possibilities. Using the Hodge decomposition theorem, we know that any vector field can be expressed as the sum of the gradient of a scalar field $\nabla\phi$, a purely rotational field $\text{rot}(u)$ and a harmonic field. Some methods focus on harmonic fields only [40]. However, the developments required to combine gradients and harmonic fields go far beyond the scope of this article.

### 5.5. Alternative shape operators

All the examples of this article were generated with the discrete Laplace-Beltrami operator or its weighted counterpart. The Laplace-Beltrami operator is an intrinsic operator, with eigenvalues independent of isometric transformations. Our modified operator has been constructed to satisfy this property. This seems desirable in many situations, especially if one considers robustness to mesh simplification, but extrinsic operators should also be explored.

The choice of the weighing function satisfying Eq. (8) was arbitrary. Finding the appropriate filters requires some insight, although the low computational cost for determining the first eigenmodes allows for manual tweaking of the parameters. Finally, other discretizations could be used to determine the Laplacian, especially ones using higher order basis functions that are known to offer better accuracy.

### 5.6. Limitations and open questions

Spectral methods are particularly well-suited for handling linear and affine constraints, as shown in the previous sections. However, non-linear constraints seem to be inherently incompatible with the creation of a vector space of solutions, which is a shortcoming of the proposed method. This limitation could be addressed by linearizing the constraints locally and exploring the neighborhood of an initial guess. A non-linear exploration could be achieved by recomputing and tracking eigenvectors from the new position and repeating the local exploration process ad libitum, recreating the manifold of constraints by linear approximations. For specific constraints, the theory of quadratic eigenvalue problems could also be leveraged. Extending spectral geometrical modeling to non-linear constraints is thus a vastly open research problem.

## 6. Conclusion

This article aimed to introduce spectral architectural geometry as a generic framework for modeling free-form structures. We identified three questions.

**How to handle linear and affine constraints?** It was shown that augmented matrices could incorporate linear and affine constraints. Furthermore, we have proposed combining the least-norm solution of the affine problem and eigenvectors of the homogeneous linearly constrained problem.

**How to modify the Laplacian shape operator to consider user-specified features?** A novel procedure to generate anisotropic operators considering user-defined geometrical features has been introduced. It is based on the computation of the geodesic distance field from feature curves or points, and its convolution with low-band filters.

**What design problems can be tackled?** Applications to shape modeling, thickness optimization and structural optimization problems have shown that spectral methods have a promising potential for a vast range of applications in architectural geometry and structural engineering.

The first explorations of spectral architectural geometry seem to point to promising applications by generating low-dimensional linear subspaces suitable for various design constraints. Yet, much remains to be explored, especially in mechanically optimal or funicular structures. Indeed, applying spectral methods to Airy stress functions could open novel design possibilities for the design of self-supporting shell structures.
<figure><img src="image_25.png" /><figcaption>Fig. 22. First eigenmodes for the conventional and an anisotropic Laplacian and linear combination. Color gradient corresponds to the amplitude of displacement.</figcaption></figure>

### CRediT authorship contribution statement

**Romain Mesnil:** Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Investigation, Conceptualization. **Kazuki Hayashi:** Writing – original draft, Validation, Software, Methodology.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Kazuki Hayashi reports financial support was provided by The Kajima Foundation. Kazuki Hayashi reports financial support was provided by Government of Japan Ministry of Education Culture Sports Science and

Technology. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

Kazuki Hayashi acknowledges partial financial support from the Kajima Foundation, Japan, for a Long-term Research Program, and from the Japan Ministry of Education, Culture, Sports, Science and Technology (MEXT) for the L-INSIGHT Project.

Romain Mesnil thanks Dr. Nicolas Montagne for valuable discussions about numerical implementation.
<table><caption>**Table A.1**
Comparison between theoretical and numerical results for different mesh resolutions. The number *n* corresponds to the number of vertices.</caption><thead><tr><th>n²+m²</th><th>λ<sub>theor</sub></th><th>λ, n = 121</th><th>λ, n = 441</th><th>λ, n = 2601</th></tr></thead><tbody><tr><td>2</td><td>19.74</td><td>19.25 (-2.6%)</td><td>19.26 (-2.5%)</td><td>19.56 (-0.9%)</td></tr><tr><td>5</td><td>49.35</td><td>46.83 (-5.4%)</td><td>47.57 (-3.7%)</td><td>48.65 (-1.4%)</td></tr><tr><td>5</td><td>49.35</td><td>49.04 (-0.6%)</td><td>48.15 (-2.5%)</td><td>48.99 (-0.7%)</td></tr><tr><td>8</td><td>78.96</td><td>75.79 (-4.2%)</td><td>76.46 (-3.3%)</td><td>78.13 (-1.1%)</td></tr><tr><td>10</td><td>98.70</td><td>92.72 (-6.4%)</td><td>94.82 (-4.1%)</td><td>98.35 (-0.4%)</td></tr><tr><td>10</td><td>98.70</td><td>90.94 (-8.5%)</td><td>95.76 (-3.1%)</td><td>96.97 (-1.8%)</td></tr><tr><td>13</td><td>128.30</td><td>118.02 (-8.7%)</td><td>124.25 (-3.3%)</td><td>126.17 (-1.7%)</td></tr><tr><td>13</td><td>128.30</td><td>118.72 (-8.1%)</td><td>123.18 (-4.2%)</td><td>127.92 (-0.3%)</td></tr><tr><td>17</td><td>167.78</td><td>145.98 (-14.9%)</td><td>159.64 (-5.1%)</td><td>166.75 (-6.6%)</td></tr><tr><td>17</td><td>167.78</td><td>147.77 (-13.5%)</td><td>162.35 (-3.3%)</td><td>166.74 (-0.6%)</td></tr></tbody></table>

<figure><img src="image_26.png" /><figcaption>**Fig. A.23.** Meshes used in the numerical validation. A Dirichlet boundary condition ($u = 0$) is applied to the contour (in purple). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

# Appendix

## A.1. Validation of numerical method

To validate our implementation of the eigenvector computation, we consider the classical problem of a membrane supported on a square of side length $a$. The solution for the eigenvalue problem $\Delta u = \lambda u$ with ($u = 0$) on the boundary is given by Eq. (A.1)

$$ \lambda = \frac{\pi^2}{a^2} (m^2 + n^2), \quad (m, n) \in \mathbb{N}^2 \qquad (A.1) $$

Several meshes of a square region with side length $a = 1$, shown in Fig. A.23, have been created, and the eigenvalue problem with Dirichlet constraints was solved. The first ten eigenvalues found with our implementations are listed in Table A.1. As expected, increasing mesh resolution results in a better approximation of the true eigenvalue. The numerical error is inferior to 5% with a few hundreds nodes, and below 1% for more than 2000 nodes.

## Data availability

Data will be made available on request.

## References

[1] Pottmann H, Eigensatz M, Vaxman A, Wallner J. Architectural geometry. Comput Graph 2015;47:145–64. http://dx.doi.org/10.1016/j.cag.2014.11.002.

[2] Bobenko AI, Suris YB. In: Discrete differential geometry: integrable structure, vol. 98, American Mathematical Soc.; 2008.

[3] Pottmann H, Liu Y, Wallner J, Bobenko A, Wang W. Geometry of multi-layer freeform structures for architecture. In: ACM SIGGRAPH 2007 papers. 2007, p. 65–es.

[4] Mesnil R, Douthe C, Baverel O. Non-standard patterns for gridshell structures: fabrication and structural optimization. J Int Assoc Shell Spat Struct 2017;58(4):277–86.

[5] Bouaziz S, Martin S, Liu T, Kavan L, Pauly M. Projective dynamics: fusing constraint projections for fast simulation. ACM Trans Graph 2014;33(4):1–11.

[6] Tang C, Sun X, Gomes A, Wallner J, Pottmann H. Form-finding with polyhedral meshes made simple. ACM Trans Graph 2014;33(4):1–9.

[7] Deng B, Bouaziz S, Deuss M, Kaspar A, Schwartzburg Y, Pauly M. Interactive design exploration for constrained meshes. Compu Aided Des 2015;61:13–23.

[8] Le C, Bruns T, Tortorelli D. A gradient-based, parameter-free approach to shape optimization. Comput Methods Appl Mech Engrg 2011;200(9):985–96. http://dx.doi.org/10.1016/j.cma.2010.10.004, URL: https://www.sciencedirect.com/science/article/pii/S0045782510002823.

[9] Hojat M, Stavropoulou E, Bletzinger K-U. The vertex morphing method for node-based shape optimization. Comput Methods Appl Mech Engrg 2014;268:494–513.

[10] Glymph J, Shelden D, Ceccato C, Mussel J, Schober H. A parametric strategy for free-form glass structures using quadrilateral planar facets. Autom Constr 2004;13(2):187–202.

[11] Mesnil R, Douthe C, Baverel O, Léger B, Caron J-F. Isogonal moulding surfaces: a family of shapes for high node congruence in free-form structures. Autom Constr 2015;59:38–47.

[12] Mesnil R, Douthe C, Baverel O, Léger B. Morphogenesis of surfaces with planar lines of curvature and application to architectural design. Autom Constr 2018;95:129–41.

[13] Brown NC, Mueller CT. Design for structural and energy performance of long span buildings using geometric multi-objective optimization. Energy Build 2016;127:748–61.

[14] Mueller CT, Ochsendorf JA. Combining structural performance and designer preferences in evolutionary design space exploration. Autom Constr 2015;52:70–82.

[15] Strobbe T, Pauwels P, Verstraeten R, Meyer R. Metaheuristics in architecture. Int J Sustain Constr design 2011;190–6.

[16] Wang L, Janssen P, Ji G. SSIEA: a hybrid evolutionary algorithm for supporting conceptual architectural design. Artif Intell Eng Des Anal Manuf 2020;34(4):458–76. http://dx.doi.org/10.1017/S0890060420000281.

[17] Brown NC, Mueller CT. Design variable analysis and generation for performance-based parametric modeling in architecture. Int J Archit Comput 2019;17(1):36–52.

[18] Yonekura K, Watanabe O. A shape parameterization method using principal component analysis in applications to parametric shape optimization. J Mech Des 2014;136(12):121401.

[19] Harding J. Dimensionality reduction for parametric design exploration. Adv Archit Geom 2016;2016:274–87.

[20] Conti ZX, Kaijima S. Enabling inference in performance-driven design exploration. In: De Rycke K, Gengnagel C, Baverel O, Burry J, Mueller C, Nguyen MM, Rahm P, Thomsen MR, editors. Humanizing digital reality: design modelling symposium Paris 2017. Singapore: Springer Singapore; 2018, p. 177–88.

[21] Oval R, Rippmann M, Mesnil R, Van Mele T, Baverel O, Block P. Feature-based topology finding of patterns for shell structures. Autom Constr 2019;103:185–201.

[22] Arora J, Haug E. Methods of design sensitivity analysis in structural optimization. AIAA J 1979;17:970–4. http://dx.doi.org/10.2514/3.61260.

[23] Ohsaki M, Hayashi K. Force density method for simultaneous optimization of geometry and topology of trusses. Struct Multidiscip Optim 2017;56(5):1157–68. http://dx.doi.org/10.1007/s00158-017-1710-8.

[24] Hayashi K, Ohsaki M. FDMopt: Force density method for optimal geometry and topology of trusses. Adv Eng Softw 2019;133:12–9. http://dx.doi.org/10.1016/j.advensoft.2019.04.002, URL: https://www.sciencedirect.com/science/article/pii/S0965997818310676.

[25] Allaire G, Jouve F, Toader A-M. Structural optimization using sensitivity analysis and a level-set method. J Comput Phys 2004;194(1):363–93. http://dx.doi.org/10.1016/j.jcp.2003.09.032, URL: https://www.sciencedirect.com/science/article/pii/S002199910300487X.

[26] Braibant V, Fleury C. Shape optimal design using B-splines. Comput Methods Appl Mech Engrg 1984;44(3):247–67. http://dx.doi.org/10.1016/0045-7825(84)90132-4, URL: https://www.sciencedirect.com/science/article/pii/0045782584901324.

[27] Haftka RT, Grandhi RV. Structural shape optimization—A survey. Comput Methods Appl Mech Engrg 1986;57(1):91–106. http://dx.doi.org/10.1016/0045-7825(86)90072-1, URL: https://www.sciencedirect.com/science/article/pii/0045782586900721.
[28] Williams CJ. The analytic and numerical definition of the geometry of the british museum great court roof. In: Mathematics & design 2001. Deakin University; 2001. p. 434–40.

[29] Taubin G. A signal processing approach to fair surface design. In: Proceedings of the 22nd annual conference on computer graphics and interactive techniques. 1995. p. 351–8.

[30] Wang Y, Solomon J. Intrinsic and extrinsic operators for shape analysis. In: Handbook of numerical analysis, vol. 20, Elsevier; 2019, p. 41–115.

[31] Vallet B, Lévy B. Spectral geometry processing with manifold harmonics. Comput Graph Forum 2008;27(2):251–60. http://dx.doi.org/10.1111/j.1467-8659.2008.01122.x.

[32] Sorkine O. Laplacian Mesh Processing. In: Chrysanthou Y, Magnor M, editors. Eurographics 2005 - state of the art reports. The Eurographics Association; 2005, http://dx.doi.org/10.2312/egst.20051044.

[33] Sischka J. Engineering the construction of the great court roof for the British Museum. Widespan Roof Struct 2000;199.

[34] Mesnil R, Douthe C, Richter C, Baverel O. Fabrication-aware shape parametrisation for the structural optimisation of shell structures. Eng Struct 2018;176:569–84. http://dx.doi.org/10.1016/j.engstruct.2018.09.026.

[35] Montagne N, Douthe C, Tellier X, Fivet C, Baverel O. Discrete voss surfaces: Designing geodesic gridshells with planar cladding panels. Autom Constr 2022;140:104200. http://dx.doi.org/10.1016/j.autcon.2022.104200, URL: https://www.sciencedirect.com/science/article/pii/S0926580522000735.

[36] Poranne R, Chen R, Gotsman C. On linear spaces of polyhedral meshes. IEEE Trans Vis Comput Graphics 2015;21(5):652–62.

[37] Andreux M, Rodola E, Aubry M, Cremers D. Anisotropic Laplace-beltrami operators for shape analysis. In: Computer vision-ECCV 2014 workshops: zurich, Switzerland, September 6-7 and 12, 2014, proceedings, part IV 13. Springer; 2015, p. 299–312.

[38] Crane K, Weischedel C, Wardetzky M. Geodesics in heat: A new approach to computing distance based on heat flow. ACM Trans Graph 2013;32(5):1–11.

[39] Koren Y, Carmel L, Harel D. Ace: A fast multiscale eigenvectors computation for drawing huge graphs. In: IEEE symposium on information visualization, 2002. IEEE, 2002, p. 137–44.

[40] Dong S, Bremer P-T, Garland M, Pascucci V, Hart JC. Spectral surface quadrangulation. In: Acm siggraph 2006 papers. 2006, p. 1057–66.

[41] Li D, Liao W, Dai N, Xie YM. Anisotropic design and optimization of conformal gradient lattice structures. Comput Aided Des 2020;119:102787. http://dx.doi.org/10.1016/j.cad.2019.102787, URL: https://www.sciencedirect.com/science/article/pii/S0010448519302386.

[42] Oval R, Mesnil R, Van Mele T, Block P, Baverel O. Two-colour topology finding of quad-mesh patterns. Comput Aided Des 2021;137:103030.

[43] Venuti F, Bruno L. Influence of in-plane and out-of-plane stiffness on the stability of free-edge gridshells: A parametric analysis. Thin-Walled Struct 2018;131:755–68.

[44] Kressner D. Numerical methods for general and structured eigenvalue problems. Berlin: Springer; 2005, http://dx.doi.org/10.1007/3-540-28502-4.

[45] Mesnil R, Douthe C, Baverel O, Leger B. Marionette meshes: Modelling free-form architecture with planar facets. Int J Space Struct 2017;32(3–4):184–98.

[46] Block P, Ochsendorf J. Thrust network analysis: a new methodology for three-dimensional equilibrium. J Int Assoc Shell Spat Struct 2007;48(3):167–73.

[47] Bletzinger KU, Ramm E. Form finding of shells by structural optimization. Eng Comput 1993;9:27–35.

[48] Tomás A, Martí P. Shape and size optimisation of concrete shells. Eng Struct 2010;32(6):1650–8. http://dx.doi.org/10.1016/j.engstruct.2010.02.013.

[49] Kiendl J, Schmidt R, Wüchner R, Bletzinger K-U. Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting. Comput Methods Appl Mech Engrg 2014;274:148–67. http://dx.doi.org/10.1016/j.cma.2014.02.001.

[50] Malek S, Wierzbicki T, Ochsendorf J. Buckling of spherical cap gridshells: A numerical and analytical study revisiting the concept of the equivalent continuum. Eng Struct 2014;75:288–98. http://dx.doi.org/10.1016/j.engstruct.2014.05.049.

[51] Tellier X, Mesnil R, Douthe C, Baverel O. On the design of torsion-free gridshells. J Int Assoc Shell Spat Struct 2024.

[52] Vaxman A. Modeling polyhedral meshes with affine maps. Comput Graph Forum 2012;31(5):1647–56. http://dx.doi.org/10.1111/j.1467-8659.2012.03170.x, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1467-8659.2012.03170.x.

[53] Mesnil R, Baverel O, Douthe C, Caron J-F, Léger B. Structural morphology and performance of plated structures with planar quadrilateral facets. J Int Assoc Shell Spat Struct 2017;58(1):7–22.

[54] Alexa M, Wardetzky M. Discrete Laplacians on general polygonal meshes. In: ACM SIGGRAPH 2011 papers. New York, NY, USA: Association for Computing Machinery; 2011, http://dx.doi.org/10.1145/1964921.1964997.

[55] Mesnil R, Douthe C, Richter C, Baverel O. Fabrication-aware shape parametrisation for the structural optimisation of shell structures. Eng Struct 2018;176:569–84.

[56] Bobenko AI, Springborn BA. A discrete Laplace-Beltrami operator for simplicial surfaces. Discrete Comput Geom 2007;38(4):740–56.

[57] Liu H-TD, Gillespie M, Chislett B, Sharp N, Jacobson A, Crane K. Surface simplification using intrinsic error metrics. ACM Trans Graph (ToG) 2023;42(4):1–17.

[58] Sharp N, Crane K. A Laplacian for nonmanifold triangle meshes. Comput Graph Forum 2020;39(5):69–80. http://dx.doi.org/10.1111/cgf.14069, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.14069.

[59] Boutillier R, Douthe C, Hauswirth L, Baverel O. Coupling grid topology generation and form-finding for the design of architectural meshes. In: Structures, vol. 68, Elsevier, 2024, 107119.