Simulation of flow and view with applications in
computational design of settlement layouts
Christian Valdemar Hansena , Anders Loggb , Carl Lundholmb

arXiv:1610.02277v1 [] 7 Oct 2016

a

Department of Mathematics and Computer Science,
University of Southern Denmark, Odense.
b
Department of Mathematical Sciences,
Chalmers University of Technology and University of Gothenburg

Abstract
We present methodology, algorithms and software for evaluating flow and
view for architectural settlement layouts. For a given settlement layout consisting of a number of buildings arbitrarly positioned on a piece of land,
in the present study an island situated on the west coast of Sweden, the
methodology allows for evaluation of flow patterns and for evaluating the
view experienced from the buildings. The computation of flow is based on a
multimesh finite element methods, which allows each building to be embedded in a boundary-fitted mesh which can be moved around freely in a fixed
background mesh. The computation of view is based on a novel and objective
measure of the view which can be efficiently computed by rasterization.
Keywords: Settlement layout, computational architecture, flow, view,
Stokes, multimesh, FEniCS.
1. Introduction
When designing settlement layouts, architects need to take a large number
of variables into consideration, such as economic interests and connections
to infrastructure (roads, water and electricity), the experienced quality of
view from the buildings, wind conditions, and many more. In this study, we
examine how to efficiently compute wind patterns and how to evaluate view.
The aim is to provide architects with a computational tool that can be used
Email addresses: cvh@imada.sdu.dk (Christian Valdemar Hansen),
logg@chalmers.se (Anders Logg), carlun@chalmers.se (Carl Lundholm)
Preprint submitted to Computers and Mathematics with Applications September 30, 2016

as a guide as part of an iterative architectural design process. For a current
example of a challenging design problem in architectural settlement layout,
see Johansson et al. (2014). There are several examples of CFD simulations in urban environments in the literature; see, e.g., Baskaran and Kashef
(1996); Blocken et al. (2007, 2011); Heuveline et al. (2011). A recent study
by Ingelsten et al. (2016) investigates the simulation of flow on geometries
directly defined by point clouds.
A central issue when designing a computational tool for settlement layout
design is that the tool should be able to quickly evaluate a multitude of
suggested designs, either as part of an optimization loop or as part of a
manual (artistic) iterative design process. Standard numerical methods for
computation of flow require that a computational mesh or grid is generated
around both buildings, ground and other objects. Generating such a bodyfitted mesh is a costly procedure and even more so when a large number of
different meshes must be created, one for each configuration of the buildings.
Instead, we examine in this work the use of multimesh finite element
methods. Multimesh finite element methods allow a problem to be posed
not on a single body-fitted mesh but instead on a collection of meshes that
may overlap arbitrarily and which together define the computational domain.
Figure 1 gives an example of such a configuration.
By allowing the computational domain to be discretized by multiple overlapping meshes, one may freely move the overlapping meshes relative to one
another, which allows a multitude of configurations to be computed and evaluated without the need for costly mesh generation. However, this flexibility
comes at a price. First, one must ensure that the finite element discretization
remains stable and that convergence is retained independently of the relative
positions of the meshes. If the finite element method is not carefully designed and does not incorporate the correct stabilization terms, certain configurations may lead to very ill-conditioned systems, low accuracy and even
blow-up. Another concern is that the formulation of multimesh finite element
methods requires integration over cut cells and interfaces. This means that
an implementation faces challenging problems in computational geometry,
when intersections and quadrature points must be computed efficiently and
robustly.
For the computation of view, we develop in this work a novel and objective measure that can be used to quantify the view for any given design.
This measure may be computed efficiently using a technique from computer
graphics known as rasterization. The measure itself allows for easy incorpo2

ration of weights that may be tweaked to give different weights to view of
air, water, ground, buildings, herbage and other objects.
In the remainder of this paper, we first review in Section 2 multimesh
finite element methods and their relation to existing methodologies for discretization of multiphysics problems and complex geometries. We then discuss in Sections 3 and 4 details of the algorithms and the implementation.
Results are presented in Section 5 and finally we present our conclusions and
discuss current limitations and future work in Section 6.

Figure 1: A collection of arbitrarily overlapping meshes together definining the computational domain, here the surroundings of a 3D propeller immersed in water. Images courtesy
of Johansson et al. (2015).

2. Multimesh finite element methods
A number of methodologies have been proposed to circumvent the limitations of interface-fitted discretization. Notable examples are the fictitious
domain (FD) method by Glowinski et al. (2001) and the extended finite
element method (XFEM) by Belytschko et al. (Moes et al., 1999). Both
methods have been successful in extending the range of problems that can
be simulated, but both suffer from limitations in that the conditioning of the
3

discretization cannot be guaranteed, and a theoretical framework for convergence analysis and error estimation is lacking. In particular, time-dependent
multiphysics problems on evolving geometries are typically discretized using
ad hoc low order discretization methods, which cannot easily be analyzed,
nor extended to higher order.
Over the past decade, a theoretical foundation for the formulation of
stabilized cut FEM has been developed by extending the ideas of Nitsche
(1971) to a general weak formulation of the interface conditions, thereby
removing the need for interface-fitted meshes. The foundations of cut FEM
were presented in Hansbo and Hansbo (2002); Hansbo et al. (2003) and the
methodology has since been extended to a number of important multiphysics
problems; see Burman and Fernández (2007); Burman and Hansbo (2007);
Becker et al. (2009); Massing et al. (2013).
We refer to multimesh finite element methods as finite element methods
based on the stabilized finite element formulations of cut FEM in combination
with a flexible and general treatment of multiple and arbitrarily overlapping
meshes. In addition to the theoretical foundation in cut FEM, the implementation and application of multimesh finite element methods relies on efficient
and robust computation of mesh-mesh intersections in 3D. This puts high
demands on the implementation and requires much more sophisticated algorithms from computational geometry than what is normally the case for finite
element problems. Figure 2 illustrates a collection of body-fitted meshes that
move independently through a fixed background mesh to discretize the flow
around a collection of moving bodies immersed in a fluid.

Figure 2: Multimesh discretization of the flow around a collection of bodies moving through
a channel.

4

Support for the formulation of multimesh finite element methods has
recently been added to the popular open-source finite element package FEniCS (Logg and Wells, 2010; Logg et al., 2012). With version 2016.1 of FEniCS, users may formulate and automatically discretize basic multimesh finite
element formulations of systems of PDE such as the Stokes problem. The
implementation in FEniCS relies on a novel implementation of mesh-mesh intersections, based on efficient generation and traversal of axis-aligned bounding box trees (AABB trees), robust low-level operations for computing and
representing the intersections of triangles (2D) and tetrahedra (3D), and generation of quadrature points on cut cells. The implementation is integrated
with the automatic code generation of FEniCS which allows multimesh discretizations to be formulated in (close to) natural mathematical language,
as will be demonstrated in Section 4. The present study constitutes the first
application of the newly implemented multimesh framework of FEniCS.
For the flow problems studied in the present work (modeled by the Stokes
equations), the finite element formulation is based on the stabilized Stokes
discretization analyzed by Johansson et al. (2015) as an extension to higherorder function spaces of the discretization previously analyzed by Massing
et al. (2014).
3. Algorithms
We here present an overview of the methodology used to simulate the
flow of air and the view for the application under consideration: the design
of a settlement (placement of houses) on a small island on the west coast of
Sweden.
3.1. Computation of flow
To model the flow over the island and houses, a finite element method
for Stokes equations on overlapping meshes is used. The method is a slight
modification of the cut finite element method for Stokes equations, presented
by Johansson et al. (2015). For a bounded domain Ω ⊂ Rd with boundary
∂Ω, the strong problem formulation for Stokes equations reads: Find the
velocity u : Ω → Rd and the pressure p : Ω → R such that


 −∆u + ∇p = f in Ω,
∇ · u = 0 in Ω,
(1)


u = 0 on ∂Ω,
5

where f : Ω → Rd is a given right-hand side. To obtain a finite element
formulation for Stokes equations on two overlapping meshes, we start by
b 0 and Ω
b 1 , called predomains.
considering two separate bounded domains Ω
b 1 on top of Ω
b 0 . Let Ω0 := Ω
b0 \ Ω
b 1 and Ω1 := Ω
b 1 , with
We may place Ω
boundaries ∂Ω0 and ∂Ω1 , respectively. We define the solution domain by
Ω := Ω0 ∪ Ω1 and the joint boundary between Ω0 and Ω1 by Γ := ∂Ω0 ∩ ∂Ω1 .
Proceeding in the usual way, by introducing suitable function spaces, multiplying the equations with test functions, and integrating by parts, one arrives
at the variational formulation of the problem. To get the corresponding finite
b 0 and Ω
b 1 , are tessellated to create the
element formulation, the predomains, Ω
b
b
b
meshes Kh,0 and Kh,1 , respectively. We call Kh,0 the background mesh, and
b h,1 the overlapping mesh. This gives us a mesh hierarchy on which we define
K
a finite element space for the velocity, Vh , and a finite element space for the
pressure, Qh . Here, Taylor-Hood elements are used to ensurethe stability of
the solution, i.e., polynomials of degree two for the velocity and polynomials
of degree one for the pressure. The finite element formulation for Stokes
equations on two overlapping meshes reads: Find (uh , ph ) ∈ Vh × Qh such
that
(Duh , Dv)Ω0 + (Duh , Dv)Ω1
− (h(Duh ) · ni, [v])Γ − ([uh ], h(Dv) · ni)Γ
+ βh−1 ([uh ], [v])Γ + ([Duh ], [Dv])Ωh,0 ∩Ω1
− (∇ · uh , q)Ω0 − (∇ · uh , q)Ω1 + ([n · uh ], hqi)Γ
− (∇ · v, ph )Ω0 − (∇ · v, ph )Ω1 + ([n · v], hph i)Γ
+ γ([ph ], [q])Ωh,0 ∩Ω1 + h2 (∆uh − ∇ph , ∆v + ∇q)Ωh,0 \ωh,0
=

(2)

(f , v)Ω0 + (f , v)Ω1 − h2 (f , ∆v + ∇q)Ωh,0 \ωh,0 ,

for all (v, q) ∈ Vh × Qh . Here, n is the unit normal to Γ exterior to Ω1 ,
hvi = (v0 + v1 )/2 is the average of v on Γ (vi is the limit of v on Ωi as we
approach Γ, for i = 0, 1), [v] = v1 − v0 is the jump in v on Γ, Ωh,0 ∩ Ω1 is the
overlapped part of all the background cells that are cut by Γ, Ωh,0 \ ωh,0 is
the visible part of all the background cells that are cut by Γ, h is the mesh
size, and β and γ are stability parameters. For more details on the method,
see Johansson et al. (2015).
The stability term γ([ph ], [q])Ωh,0 ∩Ω1 has been added to the original formulation of the method presented in Johansson et al. (2015), and is the
6

slight modification mentioned earlier. Numerical tests have indicated that
the presence of this stability term increases the robustness of the method, in
particular for the pressure. This indicates that further analysis of the original
formulation may be of interest.
For the application studied in this work, the background mesh is also
referred to as the air mesh, since it discretizes the air above the island. The
overlapping meshes are also referred to as house meshes, since they contain
the houses. We are now ready to formulate the algorithm for obtaining the
flow over the island and houses by solving (2).
Algorithm 1 Flow model
1: Geometries for the island and houses are imported. See Figure 3.
2: Meshes are generated around the geometries. For the house meshes, cells
are also generated inside the houses. See Figures 4 and 5.
3: The house meshes are placed inside the air mesh. See Figures 6- 8.
4: Boundary conditions (inlet, outlet and no-slip) are set on the air mesh.
5: A no-slip boundary condition is applied to all facets belonging to house
mesh cells that are not entirely inside the air mesh and also to facets
belonging to house mesh cells that are inside the house. For more details
see Figure 9.
6: The linear system of equations, resulting from (2), is solved.

Figure 3: Left: Geometry of the island. Right: Geometry of a simplified house.

7

Figure 4: Air mesh viewed from the outside and from the inside.

Figure 5: House mesh viewed from two different angles.

Figure 6: House mesh inside air mesh viewed from two different angles.

Figure 7: House mesh placed on the surface of the island viewed from two different angles.

8

Figure 8: Closeup on the placed house mesh.

Figure 9: Cross sections of a house mesh with a cross section of the air mesh added
for reference. Red facets are marked for no-slip, blue facets are not. Top left: Initially
all facets are marked for no-slip. Top right: First facets intersected by the air mesh are
removed. Bottom left: Then facets cut by the boundary of the air mesh are added. Bottom
right: Finally facets of the house are added.

3.2. Computation of view
We here present a novel and objective measure for evaluating the view
from a location such as the window of a building. The measure assigns a
value V between zero and one, zero being the worst possible view and one
the best possible view.
9

It is interesting to first consider the best and the worst cases. In the
present case, the worst case would occur if the view is nothing else than
another house; see illustration in Figure 10. Thus in this case the view
should evaluate to V = 0. The best case V = 1 is a view consiting entirely
of sea and sky; see Figure 11.

Figure 10: Worst case: V ≈ 0.

Figure 11: Best case: V = 1.

Examining again the worst case, the distance to a neighbouring house
should influence the value of the view. In general, the negative impact of
objects on the view should decrease by the distance. Our proposed measure
of view is expressed as an integral over the integration domain ω = ωφ × ωθ
of size |ω|:
ZZ
1
σ(φ, θ) dφ dθ .
(3)
V =
|ω|
ω

Here, σ(φ, θ) is the weight of the object viewed
weight must take a value between zero and one.
have used the following weights:


 1,
1, 
σ(φ, θ) =


 2 · s w(φ, θ) l(φ,θ) − 1,
L

at the angle (φ, θ). This
In the present study, we
if water,
if sky,

(4)

otherwise,

where s(t) is the Sigmoid function
s(t) =

1
,
1 + e−t

(5)

and l(φ, θ) is the distance to the nearest object viewed at the angle (φ, θ).
The specific element weight w(φ, θ) is set to be 0.1 and 0.7 if there is a house
10

or ground viewed at (φ, θ), respectively. Based on the premise that σ(φ, θ) =
0.9 for a house viewed at (φ, θ) and placed at the horizon approximately 5
kilometers away, the constant L is found to be 0.17 km.
Equation (3) weights the view independent of the cardinal direction. However, in northern countries like Sweden, a southern view is often weighted
higher than northern view because of the sunlight. Thus, to obtain a view
formula which can depend on the cardinial direction, we let θ ∈ [0, 2π] be
the angle in the horizontal plane. Let south be at θ = 0, consequently θ = π
is the northern direction. Assuming that one would weight the view in the
south direction three times as high as the view in the north direction, the
cardinal direction weight function could be expressed by
D(θ) = 1 +

1
sin(θ − 3π/2) .
2

(6)

For a 360◦ horizontal view valuation, we now introduce D(θ) in (3). We may
thus make the following modification to the measure of view:
Z 2π
Z
1
V360 =
D(θ)
σ(φ, θ) dφ dθ .
(7)
|ω| 0
ωφ
Note that V360 ∈ [0, 1].
4. Implementation
4.1. Implementation of flow computation
The flow model has been implemented with the open-source finite element
software FEniCS. Simulations of the flow model have been run with FEniCS
v2016.1 on a MacBook Pro with operating system OS X Yosemite version
10.10.5, a 3.1 GHz Intel Core i7 processor and 16 GB RAM.
The island geometry has been imported as an STL-file and placed on the
bottom of a sufficiently large box. The air mesh has then been generated and
stored as an XML-file by the use of functionality provided by the FEniCSpackage mshr. The geometries for the simplified houses have been made
in Gmsh (Geuzaine and Remacle, 2009). The house meshes have also been
generated with Gmsh and then converted to XML-files. To demonstrate how
the FEniCS multimesh functionality is used to implement the flow model,
we present the following code-snippets.

11

Python code
from dolfin import *
# Read the meshes from file
mesh_air = Mesh ( " mesh_air . xml " )
mesh_house1 = Mesh ( " mesh_house1 . xml " )
mesh_house2 = Mesh ( " mesh_house2 . xml " )
# Initialize multimesh , add meshes and build
multimesh = Multimesh ()
multimesh . add ( mesh_air )
multimesh . add ( mesh_house1 )
multimesh . add ( mesh_house2 )
multimesh . build ()

The order in which the meshes are added will create a mesh hierarchy, where
the latest added mesh will be on the top. A multimesh may then be used
to create finite element function spaces that are used to define the trial and
test functions.
Python code
# Create function space
P2 = VectorElement ( " P " , tetrahedron , 2 )
P1 = FiniteElement ( " P " , tetrahedron , 1 )
TH = P2 * P1
W = M u l t i M e s h F u n c t i o n S p a c e ( multimesh , TH )
# Define trial and test functions
(u , p ) = TrialFunctions ( W )
(v , q ) = TestFunctions ( W )

Here, Taylor-Hood elements of degree 2/1 are created. The trial and test
functions are used to define the bilinear and linear forms, which are in turn
are used to assemble the system matrix and load vector, respectively, for the
linear system of equations.
Python code
# Define facet normal , mesh size and stability parameters
n = FacetNormal ( multimesh )
h = 2 . 0 * Circumradius ( multimesh )
beta = Constant ( 1e1 )
gamma = Constant ( 1e8 )
def tensor_jump (v , n ) :
return outer ( v ( ’+ ’) , n ( ’+ ’) ) + outer ( v ( ’ - ’) , n ( ’ - ’) )

12

def a_h (v ,
return
+

w):
inner ( grad ( v ) , grad ( w ) ) * dX \
inner ( avg ( grad ( v ) ) , tensor_jump (w , n ) ) * dI \
inner ( avg ( grad ( w ) ) , tensor_jump (v , n ) ) * dI \
beta / avg ( h ) * inner ( jump ( v ) , jump ( w ) ) * dI

def b_h (v , q ) :
return - div ( v ) * q * dX + jump (v , n ) * avg ( q ) * dI
def l_h (v , q , f ) :
return inner (f , v ) * dX
def s_O (v , w ) :
return inner ( jump ( grad ( v ) ) , jump ( grad ( w ) ) ) * dO
def s_C (v , q , w , r ) :
return h * h * inner ( - div ( grad ( v ) ) + grad ( q ) ,
- div ( grad ( w ) ) - grad ( r ) ) * dC
def l_C (v , q , f ) :
return h * h * inner (f , - div ( grad ( v ) ) - grad ( q ) ) * dC
def s_P (q , r ) :
return gamma * inner ( jump ( q ) , jump ( r ) ) * dO
# Define bilinear form
a = a_h (u , v ) + b_h (v , p ) + b_h (u , q ) + s_O (u , v ) \
+ s_C (u , p , v , q ) + s_P (p , q )
# Define linear form
L = l_h (v , q , f ) + l_C (v , q , f )
# Assemble linear system
A = as se mb le _m ul ti me sh ( a )
b = as se mb le _m ul ti me sh ( L )

Note the resemblence of the code to the mathematical notation used in the
finite element formulation (2). Then boundary conditions for the air mesh
are created and applied to the system matrix and load vector.
Python code
# Mark boundaries for air mesh
face t_mark ers_ai r = FacetFunction ( " size_t " , mesh_air )
face t_mark ers_ai r . set_all ( 0 )

13

noslip_boundary . mark ( facet_markers_air , 1 )
inflow_boundary . mark ( facet_markers_air , 2 )
outflow_boundary . mark ( facet_markers_air , 3 )
# Create boundary conditions
bc0 = M u l t i M e s h D i r i c h l e t B C (V , noslip_value ,
facet_markers_air , 1 , 0 )
bc1 = M u l t i M e s h D i r i c h l e t B C (V , inflow_value ,
facet_markers_air , 2 , 0 )
bc2 = M u l t i M e s h D i r i c h l e t B C (Q , outflow_value ,
facet_markers_air , 3 , 0 )
# Apply boundary conditions
bc0 . apply (A , b )
bc1 . apply (A , b )
bc2 . apply (A , b )

The no-slip condition for the concerned facets of the house meshes are created
and applied in a similar fashion. Finally, the linear system is solved and the
solution components on the different meshes are extracted.
Python code
# Compute solution
w = MultiM eshFun ction ( W )
solve (A , w . vector () , b )
# Extract solution components
u_air = w . part ( 0 ) . sub ( 0 )
p_air = w . part ( 0 ) . sub ( 1 )
u_house1 = w . part ( 1 ) . sub ( 0 )
p_house1 = w . part ( 1 ) . sub ( 1 )
u_house2 = w . part ( 2 ) . sub ( 0 )
p_house2 = w . part ( 2 ) . sub ( 1 )

4.2. Implementation of view computation
Given the mesh constructed for the flow simulation, it is now possible to
visualize the view from a given point. In practice we do this by rendering an
image from the 3D mesh with the use of the rasterization rendering technique.
This technique goes back to Catmull (1974); Pineda (1988); Heckbert (1989).
The rasterization algorithm projects the triangles from the 3D mesh onto a
2D image.
The projecting principle is sketched in Figure 12, where the vertices which
are mapped onto the 2D image plane are used to check which pixels the
14

Figure 12: Projection of triangles from the 3D mesh onto the 2D image.

triangle covers. The naive idea is to loop across all the pixels in the image
and check if they are inside the projected triangle or not. The efficiency
of this approach depends on the size of the triangles. To account for small
triangles, one may optimize the search by only searching the pixels which lie
inside the bounding box of the triangle. In Figure 13, the blue square around
the triangle illustrates the bounding box for which the corner coordinates are
rounded to the nearest pixel.
This can be optimized further by not checking all the pixels inside the
bounding box, but starting in the pixel which contains the top-point of the
triangle. Then go stepwise down and check the pixels to the left and right of
the reference point. If we already visited one or more pixels which are inside
the projected triangle and then come to a pixel which is not in the triangle,
the search stops in that direction. The search algorithm is illustrated with
the green lines in Figure 13. There are several ways to optimize this, see
Pineda (1988) for further reading.
Looping over all the triangles to do the projections one by one, can cause
two or more projected triangles to overlap. To decide which one that should
be shown in the image we have to look at the distance to the them. The
distances to the triangles which are already shown in the image are stored in
a two dimensional array with the same dimension as the image. Thus only
the element with the shortest distance is shown in the image.
The rasterization algorithm is implemented in C++ and by the use of the
15

Figure 13: To check which pixels that lies inside the triangle, we start from the top-point
of the triangle (the green dot) going down one step at a time to check if the pixels to the
left and right are contained in the triangle.

SWIG interface compiler it is possible to access the rasterization algorithm
from the C++ code in a Python script. The rasterization algorithm needs as
input a list of FEniCS (.xml) meshes, where the first mesh in the list should
be the main island mesh, and the rest house meshes. With the scene set, the
algorithm needs to know the size of the image, both the size in pixels and the
real size measured in the same units as the meshes. Also the position of the
camera and the distance between the camera and the image is needed. The
direction for the camera expressed by a vector should also be given. With
these inputs, the algorithm generates the image and a matrix S which has
the same size as the image and contains a value for σ for each pixel. The
view valuation V is then found from S and given as output.
The view V may thus be computed from the image of the view generated
by rasterization. If we want to compute the 360◦ view from (7), the domain ω
is assumed to be either a cylinder or a sphere. As the rasterization algorithm
generates flattened images, (7) cannot be used directly. However (7) could
be estimated by
N
−1
X
D (2πi/N )
Vi ,
(8)
V360 =
N
i=0
where N ≥ 3 is the number of images and Vi is the view valuation for image
16

number i. It is important that no images are overlapping and that they, when
joined together, form the boundary of a convex polygon with N edges, when
seen from above. Thus Vi is evaluated in the direction with angle 2πi/N with
respect to south and an image width of 2d tan(π/N ), where d is the distance
between the camera and the image.
5. Results
5.1. Flow
Results from two simulations using the multimesh finite element method
for Stokes equations on overlapping meshes (2) are presented. The first
simulation was run with one overlapping house mesh and the results can be
seen in Figures 14 and 15. The other simulation was run with two overlapping
house meshes and the results can be seen in Figures 16 and 18. The results are
in the form of tubed streamlines from the velocity field of the finite element
solution. Tubes from the solution on the house mesh are made thicker than
the ones from the air mesh for illustration purposes. The flow in all the
figures goes from left to right.

Figure 14: Tubed streamlines from the simulation with one overlapping house mesh. With
and without the house mesh visible. Notice the continuous velocity field going between
the two different meshes used for discretization.

17

Figure 15: Tubed streamlines from the simulation with one overlapping house mesh. In
front of and behind the house.

Figure 16: Tubed streamlines from the simulation with two overlapping house meshes.
Left: Overview with the house meshes visible. Right: Overview without the house meshes
visible but with some extra houses added for illustration.

Figure 17: Tubed streamlines from the simulation with two overlapping house meshes.
Closeup on the first house with and without the house meshes visible.

18

Figure 18: Tubed streamlines from the simulation with two overlapping house meshes.
Closeup on the second house with and without the house mesh visible.

5.2. View
To test the computation of view, houses were arbitrarily positioned on
the island as seen in Figure 19.

Figure 19: Overview of the island with ten houses. The yellow dot represents the position
of the camera for the view computation.

For this test case we would like to examine how attractive it is to place
a new house at the yellow dot, just next to the already existing house. The
views in Figures 20–23 are all views from the camera located at the yellow
dot in Figure 19.
Figures 20–23 show a variety of views from the selected location. Figure
20 shows a fair view with a lot of sky and sea view but also three houses
and a part of the island can be seen. Therefore the view valuation becomes
V = 0.74. To get a better view there should be as much sea or sky as
possible and less houses and ground. A better example can be seen in Figure
21, where V = 0.87. Examples of views from the other end of the view scale
can be found in Figure 22 and 23, where a neighbouring house blocks the
view. In Figure 23 one can almost only see the neighbouring house, which
also results in a very low value of the view.
19

Figure 20: Fair view: V = 0.74

Figure 21: Good view: V = 0.87

Figure 22: Poor view: V = 0.38

Figure 23: Bad view: V = 0.02

The 360◦ horizontal view measure (8) naturally depends on the cardinal
direction. We compute the value based on 32 computed rasterizations. If
the south direction is chosen to be in the direction of the view, we obtain
the value V360 = 0.68 (as a result of looking straight at a neighboring house),
whereas if we choose the south direction to be in the opposite direction of
the neighbouring house we obtain V360 = 0.73.
6. Conclusions and future work
We have presented a generic framework for evaluating flow and view for
settlement layouts. The framework allows multiple configurations to be computed and evaluated with relative ease. The current proof-of-concept implementation has several limitations that will be addressed in future work.
These limitations and extensions fall into three different categories: efficiency,
robustness and ease of use.
Regarding efficiency, the current implementation of multimesh methods
adds a significant overhead in computational time compared to standard fi20

nite element methods. In particular, the assembly of multimesh finite element
variational forms is significantly slower than standard assembly. Another limitation is the lack of a properly preconditioned iterative method (the current
implementation uses the direct sparse solver UMFPACK). Both these limitations are the focus of ongoing work. Future work will also consider the
extension of the current multimesh implementation to parallel architectures.
Regarding robustness, we noticed during this study that some particular
configurations of buildings resulted in numerical instabilities. This is likely
due to bugs or untreated corner cases in the computational geometric framework of the FEniCS multimesh implementation. This issue is also the focus
of ongoing work.
To be a useful tool in an iterative architectural process, the framework
must not only be efficient and robust. It must also be easy to use. Future
work will also consider the creation of a user-friendly interface. In particular,
it would be highly relevant to consider the creation of VR or AR interfaces
to our framework.
Acknowledgements
This work was supported by the profile area Virtual Cities as part of
Building Futures, an Area of Advance at Chalmers. We are also greatful towards CREAM Architects in Gothenburg for ideas, feedback and inspiration
during this project.
References
References
Baskaran, A., Kashef, A., 1996. Investigation of air flow around buildings
using computational fluid dynamics techniques. Engineering Structures
18 (11), 861–875.
Becker, R., Burman, E., Hansbo, P., 2009. A Nitsche extended finite element method for incompressible elasticity with discontinuous modulus of
elasticity. Comp. Methods Appl. Mech. Engrg. 198 (41), 3352–3360.
Blocken, B., Carmeliet, J., Stathopoulos, T., 2007. Cfd evaluation of wind
speed conditions in passages between parallel buildings—effect of wallfunction roughness modifications for the atmospheric boundary layer flow.
Journal of Wind Engineering and Industrial Aerodynamics 95 (9), 941–962.
21

Blocken, B., Stathopoulos, T., Carmeliet, J., Hensen, J. L., 2011. Application
of computational fluid dynamics in building performance simulation for
the outdoor environment: an overview. Journal of Building Performance
Simulation 4 (2), 157–184.
Burman, E., Fernández, M. A., 2007. Stabilized explicit coupling for fluidstructure interaction using Nitsche’s method. C. R. Math. Acad. Sci. Paris
345 (8), 467–472.
Burman, E., Hansbo, P., 2007. A unified stabilized method for Stokes’ and
Darcy’s equations. J. Comput. Appl. Math. 198 (1), 35–51.
Catmull, E. E., 1974. A subdivision algorithm for computer display of curved
surfaces. Ph.D. thesis, The University of Utah.
Geuzaine, C., Remacle, J.-F., 2009. Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities. International Journal
for Numerical Methods in Engineering 79 (11), 1309–1331.
Glowinski, R., Pan, T. W., Hesla, T. I., Joseph, D. D., Periaux, J., 2001.
A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate
flow. J. Comput. Phys. 169 (2), 363–426.
Hansbo, A., Hansbo, P., 2002. An unfitted finite element method, based on
Nitsche’s method, for elliptic interface problems. Comp. Methods Appl.
Mech. Engrg. 191 (47), 5537–5552.
Hansbo, A., Hansbo, P., Larson, M. G., 2003. A finite element method on
composite grids based on Nitsche’s method. ESAIM, Math. Model. Numer.
Anal. 37 (03), 495–514.
Heckbert, P. S., 1989. Fundamentals of texture mapping and image warping.
Tech. rep., University of California at Berkeley, Berkeley, CA, USA.
Heuveline, V., Ritterbusch, S., Ronnas, S., 2011. Augmented reality for urban
simulation visualization. Preprint Series of the Engineering Mathematics
and Computing Lab 0 (16).
Ingelsten, S., Mark, A., Edelvik, F., Logg, A., Österbring, M., 2016. Urban
cfd-simulation using point cloud data. In: Proceedings of NSCM-26: the
26th Nordic Seminar on Computational Mechanics.
22

Johansson, A., Larson, M. G., Logg, A., 2015. High order cut finite element
methods for the Stokes problem. Advanced Modeling and Simulation in
Engineering Sciences 2 (1), 1–23.
URL http://dx.doi.org/10.1186/s40323-015-0043-7
Johansson, G., Karlén, F., Stark, M., 2014. Lilla fjellsholmen, varsamt byggande i en unik skärgårdsmiljö. Master’s thesis, Chalmers University of
Technology, Gothenburg, Sweden.
Logg, A., Mardal, K.-A., Wells, G. N., et al., 2012. Automated Solution of
Differential Equations by the Finite Element Method. Springer.
Logg, A., Wells, G. N., 2010. DOLFIN: Automated finite element computing.
ACM Transactions on Mathematical Software 37 (2).
URL http://www.dspace.cam.ac.uk/handle/1810/221918/
Massing, A., Larson, M. G., Logg, A., Rognes, M. E., 2013. An overlapping
mesh finite element method for a fluid-structure interaction problem. submitted to Communications in Applied Mathematics and Computational
Science.
Massing, A., Larson, M. G., Logg, A., Rognes, M. E., 2014. A stabilized
nitsche overlapping mesh method for the stokes problem. Numerische
Mathematik, 1–29.
Moes, N., Dolbow, J., Belytschko, T., 1999. A finite element method for
crack growth without remeshing. Int. J. Numer. Meth. Engng 46, 131–150.
URL http://venus.usc.edu/PAPERS/MultiScaleMechanics/XFEM.pdf
Nitsche, J., 1971. Über ein Variationsprinzip zur Lösung von DirichletProblemen bei Verwendung von Teilräumen, die keinen Randbedingungen
unterworfen sind. In: Abhandlungen aus dem Mathematischen Seminar
der Universität Hamburg. Vol. 36. Springer, pp. 9–15.
Pineda, J., 1988. A parallel algorithm for polygon rasterization. In: In Proceedings of Siggraph ’88. pp. 17–20.

23

