Acta Polytechnica CTU Proceedings


doi:10.14311/APP.2018.15.0025
Acta Polytechnica CTU Proceedings 15:25–30, 2018 © Czech Technical University in Prague, 2018

available online at http://ojs.cvut.cz/ojs/index.php/app

ISOGEOMETRIC BEAM ELEMENT EXTENDED FOR GRIDSHELL
ANALYSIS

Edita Dvořáková∗, Bořek Patzák

Department of Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7,
166 29 Prague 6, Czech Republic

∗ corresponding author: edita.dvorakova@fsv.cvut.cz

Abstract. The design of gridshell structures is a complicated process, as the resulting shape of
structure depends on initial undeformed grid geometry as well as on the history of boundary conditions
used to form the final structure. In practice, the physical model is often used to determine the shape of
the structure and the initial grid at the same time. Introducing Isogeometric analysis into a design of
gridshells simplifies the design process as the problem can be easily recalculated when initial grid or
boundary conditions change and the resulting shape can be immediately illustrated. The presented
paper discusses possibilities in isogeometric gridshell modeling and proposes possible solutions of
identified problems.

Keywords: Curved beam element, gridshell structure, isogeometric beam element, NURBS.

Figure 1. Realization of a multipurpose gridshell
pavilion in Třebešov [1].

1. Introduction
Gridshell is a light-weight structure gaining its
strength from a doubly-curved geometry. Such a struc-
ture is typically made of wood and it is constructed
from initially straight grid of laths. The laths are
mutually connected by joints allowing only for the
in-plane rotation. The final shape of the structure
is reached by subsequent introduction of a load or
a prescribed displacement at selected nodes. Once
the desired shape is reached, the joints are tightened
in order to fix the in-plane rotation and stiffen the
structure. See Fig. 1 for an example of a gridshell
structure.

It is important to design the structural shape prop-
erly to reach high low-bearing capacity. The design
process of a gridshell can be divided into three main
parts: form-finding, analysis, and construction. Form-
finding is typically performed by an architect who is
seeking for the optimal shape of the gridshell taking
into the account the aesthetics and structural behavior

of the structure. Either physical modeling or com-
puter simulation is used for form-finding, nevertheless
these methods (and available tools) do not provide
knowledge about structural behavior. When the de-
sired shape is defined the structure is assessed in a
Finite Element Method (FEM) solver, however the
analysis considers only the final shape of the struc-
ture. There are no available tools which would take
into account the history of the construction steps and
would allow to assess all the intermediate construc-
tion stages. Currently the construction itself relies on
trial-error methods supported by a human intuition
and experience.
Our goal is to develop a complex tool for analysis

assisted form-finding which would allow to find the
optimal shape of gridshell and which would help to
propose sequence of construction steps. The idea is to
use the isogeometric approach which eases incorpora-
tion of the analysis within available design tools and
which may bring higher accuracy into the analysis of
doubly-curved geometry thanks to the more suitable
geometry description. The presented paper discusses
the advantages of isogeometric gridshell modeling and
also points out some issues connected to its applica-
tion.

2. Isogeometric analysis of beams
Isogeometric analysis has been introduced by Thomas
Hughes in 2005 [2] with the aim to bridge the gap
between the Computer-Aided Design (CAD) and the
Finite Element Analysis (FEA). This goal has been
achieved by using the same CAD basis functions for
geometry description and unknowns approximation in
the analysis. Thus one model can be shared between
the design and the analysis.

25

http://dx.doi.org/10.14311/APP.2018.15.0025
http://ojs.cvut.cz/ojs/index.php/app


Edita Dvořáková, Bořek Patzák Acta Polytechnica CTU Proceedings

Figure 2. Description of B-spline (NURBS) finite element geometry. The geometry is modeled by one patch
consisting of two knotspans with quadratic NURBS approximation.

2.1. NURBS-based analysis
The most wide-spread technology in CAD indus-
try are NURBS (Non-Uniform Rational B-Splines).
The NURBS curve is obtained by linear combination
of Cartesian coordinates of the control points P using
NURBS functions Rpi

C(ξ) =
n∑
i=1

R
p
i (ξ)Pi. (1)

Each NURBS function is generated by weighting B-
spline functions Npi with a given weight wi associated
with control point Pi

R
p
i (ξ) =

Ni,p(ξ)wi∑n
j=1 Nj,p(ξ)wj

. (2)

B-spline functions form a special subset of NURBS
functions (corresponding to constant wi) and are de-
rived recursively starting with a piecewise constant
functions

Ni,0(ξ) =
{

1 if ξi ≤ ξ < ξi+1
0 otherwise , (3)

Ni,p(ξ) =
ξ − ξi

ξi+p − ξi
Ni,p−1(ξ) +

ξi+p+1 − ξ
ξi+p+1 − ξi+1

Ni+1,p−1(ξ), (4)

where p is the degree of the B-spline function, ξi is
the coordinate of the ith-knot and ξ ∈ 〈0, 1〉 is a
parametric coordinate.

The computational domain in Isogeometric analysis
(IGA) is firstly divided into patches, which are further
divided into knotspans (often referred to as elements

in IGA). Understanding of knotspans is similar to
elements in standard FEM, nevertheless higher conti-
nuity up to Cp−1 between knotspans can be achieved
naturally using NURBS, unlike C0 in standard FEM.
Moreover, the NURBS basis functions are generally
non-iterpolatory. In the knots (points which are di-
viding patch into knotspans), the continuity can be
locally decreased up to C0 by knot multiplication,
which is a standard technique in NURBS technology.
This procedure imposes the interpolatory behavior of
the functions at the particular point.
For better understanding of NURBS see Fig. 2,

where mapping between real and parametric geometry
is illustrated. For more information the reader could
refer to [3].

2.2. Beam element formulation
The formulation of the three-dimensional beam ele-
ment [4] is derived in local coordinate system (t,n,b),
where t,n,b are tangent, normal and binormal vec-
tors and s is the curvilinear coordinate measured along
the centerline of the beam. Presented formulation is
based on Timoshenko beam theory and therefore, in
three-dimensional space, there are six independent un-
knowns: tangential displacement ut, normal displace-
ment un, binormal displacement ub and rotations θt,
θn, θb (see Fig. 3).

A strain-displacement matrix B is defined as

ε = Br, (5)

where ε = {εm,γn,γb,χt,χn,χb}T and r =
{ut,un,ub,θt,θn,θb}T . B is derived from the follow-
ing relations for membrane strain εm, shear strains
γn and γb, torsional strain χt, and bending strains χn

26



vol. 15/2018 Isogeometric Beam Element Extended for Gridshell Analysis

b(s)

r(s)

t(s)n(s)

ut
θt

θb

ub
θn

un

Figure 3. Local coordinate system and degrees of
freedom of a three-dimensional beam element.

and χb

εm = u′t −κun,
χt = θ′t −κθn,
γn = κut + u′n − τub −θb,
χn = κθt + θ′n − τθb, (6)
γb = τun + u′b + θn,
χb = τθn + θ′b.

Curvature κ is given as

κ =
∣∣∣∣
∣∣∣∣d2r(s)ds2

∣∣∣∣
∣∣∣∣ (7)

and torsion τ

τ =
dn(s)
ds

b(s), (8)

where r(s) is the position vector. The material stiff-
ness matrix D is given by constitutive relations

N = EAεm,
Mt = GIkχt,
Qn = GAnγn,
Mn = EInχn, (9)
Qb = GAbγb,
Mb = EIbχb,

where E is the Young’s modulus, G is the shear modu-
lus, Ik is the torsional moment, In,Ib are the moments
of inertia, A is the cross-sectional area, Ab,An are the
reduced shear areas. Finally, the element stiffness
matrix K is defined as

K =
∫ L

0
BT DB ds, (10)

where the integral is evaluated using Gaussian quadra-
ture, which is not exact for NURBS functions, but it
has a sufficient accuracy.

2.3. Numerical locking
Well-known problem of numerical locking of Timo-
shenko beam elements in standard FEM occurs also

0 0.5 1

x

-0.5

0

0.5

Q
n

Shear force Q

0 0.5 1

x

0

0.1

0.2

0.3

M
b

Bending moment M

Figure 4. Simple gridshell structure modeled by four
patches with cubic approximation. Joint is modeled
by connection of patches. Resulting internal forces
overlap exact solution.

for the isogeometric formulation. The same approxi-
mation order for displacements and for rotations leads
to the shear strains of higher order than bending mo-
ments, while it should be vice-versa. Moreover the
constant strains cannot be approximated along the
entire patch due to the field inconsistency in the shear
terms.

The reduced integration, popular locking-treatment
in standard FEM for its computational efficiency, is
not convenient for IGA, as the optimal reduced inte-
gration rule depends on the choice of approximation
order and continuity and such a rule is hard to be de-
termined. On the contrary, Discrete shear gap (DSG)
method [5, 6] or B̄-method [5] can be used indepen-
dently of the choice of approximation, the latter being
used in the presented implementation of the element.

The main idea of the B̄-method is to project mem-
brane strain εm and shear strains γn,γb onto a basis
of lower order. For example, the contribution of mem-
brane strain εm to the B̄ matrix is

B̄m =
[
B̄

1
m B̄

2
m 0 0 0 0

]
, (11)

where for example

B̄
1
mij

=
∫ L

0
ÑiN

′
j ds (12)

and where Nj is the jth basis function of the original
basis and Ñi is the ith basis function of the lower
order basis. The membrane contribution to stiffness
matrix Km is

Km = EA B̄
T

mM̃
−1
B̄m (13)

where M̃ is the element “mass matrix” calculated in

27



Edita Dvořáková, Bořek Patzák Acta Polytechnica CTU Proceedings

0 0.5 1

x

-0.5

0

0.5

Q
n

Shear force Q

0 0.5 1

x

-0.2

0

0.2

0.4

M
b

Bending moment M

0 0.5 1

x

-1

0

1

Q
n

Shear force Q

0 0.5 1

x

-0.2

0

0.2

0.4

M
b

Bending moment M

a) C3-continuity b) C2-continuity

0 0.5 1

x

-1

0

1

Q
n

Shear force Q

0 0.5 1

x

-0.2

0

0.2

0.4

M
b

Bending moment M

0 0.5 1

x

-0.5

0

0.5

Q
n

Shear force Q

0 0.5 1

x

0

0.1

0.2

0.3

M
b

Bending moment M

c) C1-continuity d) C0-continuity

Figure 5. Resulting internal forces (blue continuous line) and exact solution (red dash-dot line) for different continuity
at the position of concentrated force load.

the lower order basis

M̃ij =
∫ L

0
ÑiÑj ds. (14)

The analogical procedure can be followed also for the
shear strains γn and γb.

3. Application to gridshells
Isogeometric formulation of beam element allows to
perform the analysis directly within the CAD system,
however the modeling of gridshells with presented iso-
geometric beam element formulation is not straightfor-
ward. There are two possible approaches to modeling
gridshells: a) multiple patches per one lath, b) single
continuous patch per lath.

3.1. Lath modeled by multiple patches
This approach is similar to the use of standard FEM,
where the joint is usually modeled by inserting a node.
In IGA, it means to assign at least one patch to
each lath segment between two joints leading to the
interpolatory degrees of freedom at particular points
(joints). This approach is simple and does not require
any additional implementation, nevertheless it leads
to loss of some advantages of isogeometric formulation,
namely loss of higher inter-element continuity along
the entire lath and loss of possibility to handle one lath
as a single entity, which might be highly convenient
for the development of a complex design tool.
A simple test has been performed to verify the

"multiple-patches" approach. Two laths mutually con-
nected in the middle, each modeled by two patches
with one knotspan and cubic B-spline approximation,
have been subjected to the transverse force placed at
the joint. The model and resulting shear force γn and
bending moment χb are shown in Fig. 4. The results

0 0.5 1

-0.1

0

0.1

N
*

0 0.5 1

0

0.05

0.1

N
v

*

Figure 6. Additional basis functions to support case
of concentrated force load on 2D straight beam with
NURBS approximation.

overlap the exact solution and prove the possibility to
use this modeling technique for gridshell analysis.

3.2. Lath modeled by a single patch
Isogeometric approach offers a possibility to model
each grid lath by a single patch. Such a model requires
to add additional constraints to enforce the compati-
bility of displacements and rotations in joints as the
compatibility cannot be achieved by sharing particular
degrees of freedom due to the non-interpolatory nature
of isogeometric basis functions. A special boundary
condition based on the Lagrange multipliers method
has been implemented to support this approach. The
potential energy is enhanced by an additional term
yielding the required continuity

Wint +
∑
k

λk(r̄ki − r̄
k
j ) = Wext, (15)

where λk is vector of Lagrange multipliers for kth joint,
r̄ki and r̄

k
j are the vectors of constrained unknowns at

the kth joints on the ith and jth beams, respectively.
Note, that due to the non-interpolatory nature of
NURBS basis functions, the vectors r̄ki and r̄

k
j have

to be expressed as a linear combination of the control
points values using NURBS basis functions evaluated

28



vol. 15/2018 Isogeometric Beam Element Extended for Gridshell Analysis

0 0.5 1

x

-1

-0.5

0

0.5

1

Q
n

Shear force Q

0 0.5 1

x

0

0.1

0.2

0.3

M
b

Bending moment M

0 0.5 1

x

-2

-1

0

1

2

Q
n

Shear force Q

0 0.5 1

x

0

0.1

0.2

0.3

0.4

M
b

Bending moment M

Figure 7. Resulting internal forces for simply supported beam subjected to a) a concentrated force, b) a concentrated
force and continuous force loading. NURBS approximation enriched by the functions tailored to support case of
concentrated force has been used. The results overlap the exact solution.

at the joint position. The local coordinate ξi (ξj)
of the joint on the beam i (j) has to be calculated
according to the parametrization.
The physical meaning of Lagrange multipliers λ

corresponds to the forces and moments required to
enforce the compatibility and thus such an approach
results in concentrated force and moment loadings
at the joints positions. The problem of such a for-
mulation is the inability of NURBS basis functions
to represent exact solution with discontinuities in
strains and internal forces corresponding to concen-
trated loadings. This inability is a consequence of the
higher continuity of NURBS basis along the entire
patch.
This problem is documented on the similar test

as the one used in Section 3.1 (see Fig. 4). Unlike
in previous section, each lath is now modeled by a
single patch using just a single span with cubic B-
spline approximation, which results in C3-continuity
at the joint position causing the oscillations in internal
forces, see Fig. 5a for the results. The continuity can
be lowered to C2 by placing a knot at the particu-
lar position, the continuity can be further reduced
up to C0 by increasing the knot multiplicity. The
resulting internal forces for all choices of continuity
are shown in Fig. 5. The oscillations in the numerical
solutions are sufficiently removed only in the case of
C0-continuity. Note, that the procedure of decreasing
the continuity to C0 results in the same model as
in Sec. 3.1 and thus the advantages of isogeometric
analysis are eliminated.

The main idea, how to preserve the higher continu-
ity while avoiding oscillations in numerical solutions,
is to enrich the NURBS basis by the special basis

functions tailored to match the exact solution. In this
paper, the derivation of the additional basis functions
is illustrated on the case of two-dimensional straight
beam subjected to concentrated force.

From the knowledge of the exact solutions of bend-
ing moment Mex and shear force Qex under the unit
force

Mex = EIβ′(x), (16)
Qex = kGA (v′(x) + β(x)) , (17)

where β is a rotation and v is a transverse displace-
ment, the additional basis functions N∗β, N

∗
v can be

derived

N∗β = β =
∫

1
EI

Mex dx, (18)

N∗v = v =
∫

1
kGA

Qex −β dx. (19)

See Fig. 6 for the resulting basis functions. The final
approximation for the enhanced solution is

v(x) =
∑

Nivi + N∗vv
∗, (20)

β(x) =
∑

Niβi + N∗ββ
∗. (21)

The integration constants were determined with re-
spect to the boundary conditions of a simply sup-
ported beam.
The performance of the enhanced approximation

has been tested by means of two simple test. Firstly,
simply supported beam subjected to the concentrated
force load was considered, in the second test the con-
centrated force has been combined with continuous
loading. The resulting internal forces for both load
cases overlap exact solutions, see Fig. 7.

29



Edita Dvořáková, Bořek Patzák Acta Polytechnica CTU Proceedings

4. Conclusions
Two possible approaches to gridshell modeling have
been presented. The simple FEM-like approach, where
the joints are modeled by connection of two patches,
provides good results, however the advantages of isoge-
ometric approach are lost. In order to keep the higher
continuity along the entire lath, it has to be modeled
by a single patch. In such a case, the enforcement of
continuity of displacements and rotations in joints by
Lagrange multipliers inevitably leads to application of
concentrated forces and moments. It has been shown,
that enhancement of NURBS basis functions by a
set of the special functions is sufficient to support
concentrated loadings.

The paper documents the first necessary steps of the
development of a complex design tool for gridshells.
The presented idea of enhancement of the NURBS
basis has to be combined with locking treatment to
provide robust element formulation, moreover due
to the large deformations, the formulation has to be
extended for geometrically non-linear analysis.

Acknowledgements
The financial support of this research by the Grant Agency
of the Czech Technical University in Prague (SGS project
No. SGS18/037/OHK1/1T/11) is gratefully acknowledged.

References
[1] A. Vaňek. Multipurpose gridshell pavilion in Třebešov,
Czech Republic, physical model as a groundwork for
realization, harvest super, 2016.

[2] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs.
Isogeometric analysis: CAD, finite elements, NURBS,
exact geometry and mesh refinement. Comput Methods
Appl Mech Engrg 194:4135–4195, 2005.

[3] L. Piegl, W. Tiller. The NURBS Book.
Springer-Verlag Berlin Heidelberg, New York, 1997.

[4] G. Zhang, R. Alberdi, K. Khandelwal. Analysis of
three-dimensional curved beams using isogeometric
approach. Engineering Structures 117:560–574, 2016.

[5] R. Bouclier, T. Elguedj. Locking free isogeometric
formulations of curved thick beams. Comput Methods
Appl Mech Engrg 245-246:144–162, 2012.

[6] R. Echter, M. Bischoff. Numerical efficiency, locking
and unlocking of NURBS finite elements. Comput
Methods Appl Mech Engrg 199:374–382, 2010.

30


	Acta Polytechnica CTU Proceedings 15:25–30, 2018
	1 Introduction
	2 Isogeometric analysis of beams
	2.1 NURBS-based analysis
	2.2 Beam element formulation
	2.3 Numerical locking

	3 Application to gridshells
	3.1 Lath modeled by multiple patches
	3.2 Lath modeled by a single patch

	4 Conclusions
	Acknowledgements
	References