Overview of the methods for the modelling of historical masonry structures

An overview of methods for the modelling of historic masonry structures, starting from simple traditional to modern methods suitable for use in computer programs, is presented. Main characteristics of masonry structures, as well as principal features of methods for their analysis, such as dynamic analytical methods, limit state methods, finite element methods, and discrete element methods, are presented. A review of application of individual methods with regard to the required level of accuracy, quantity of input data, and calculation time, is presented.


Introduction
Construction of buildings with stone or clay bricks that are held together by mortar is one of the oldest building techniques that is still in use today. Some of the oldest manmade structures are masonry stone huts in form of a circle found near Lake Hullen in Israel that date back to the period between 9,000 BC and 8,000 BC [1]. These first masonry structures were amazing piles of natural stone [2]. As humans became more skilled and began using tools, masonry structures started to be more symmetrical. Later on humans learned to make bricks moulding them from mud or clay and leaving them to dry in the air and, later, baking them in ovens. Bricks were strong, uniform and easy to make. In addition to improvements in block-making techniques, various cultures started to use architectural features such as pillars to obtain the height they needed, or beams, arches, and domes to bridge various distances. Due to its simplicity, masonry has a long worldwide tradition of use in construction. The durability of masonry structures is evidenced by the number of structures that are still in use after hundreds or even several thousands of years. Some examples of structures that have become symbols of certain cultures are the Egyptian pyramids that originate from the period between 2800 BC and 2000 BC, the Parthenon in Greece from the fifth century BC, The Great Wall of China whose construction began in the fifth century before Christ, and the Colosseum in Rome dating back to the first century AD (Figure 1). Despite the simplicity that is manifested during construction of masonry structures, the understanding and describing mechanical behaviour of those structures, especially in conditions of seismic loading, represents a true challenge due to the nature of masonry structures which are characterized by a complex and particularly nonlinear behaviour, because of the presence of joints among blocks that can but need not be filled with mortar. Many masonry structures are located in seismically active zones where earthquakes have revealed on occasion their vulnerability. These earthquakes often inflict damage to buildings and monuments that are classified as cultural heritage, but also to contemporary masonry structures. In order to reduce the number of human victims and level of damage to such structures, a better look must be taken into Overview of the methods for the modelling of historical masonry structures behaviour of masonry structures under seismic load. The main motivation for this paper is to describe the methods that have been developed in relation to the seismic analysis of masonry structures. The goal of the paper is to provide scientists with the current state-of-the-art about the behaviour of masonry structures under seismic load, and with methods for modelling the effects of seismic activity on these structures.

Traditional and modern approaches to the analysis of historical structures
The calculation of masonry structure had not been mentioned for a long period of time during recorded history. The art of constructing masonry structures was something gained by experience and passed down from one generation to another. Roman architect Vitruvius in his work Ten Books on Architecture [3] compares the quality of stone and wood from various locations and speaks about proportions of various construction elements and structures, but he does not say anything about calculations.
A structure that has right proportions was considered structurally adequate and that way of thinking was kept during the Middle Ages. The characteristic of this time is strictly keeping knowledge of verified proportions that are passed on from one generation to another. Many impressive structures made in these times that are still here today show that the experience and knowledge about stability and distribution of forces within the masonry structure was not negligible. Many of the buildings constructed during the Renaissance in the fifteenth century became more slender. This required a proper theoretical basis for their construction. In the second half of seventeenth century, Robert Hooke noticed that parts of stone arch are shaped as a reverse catenary. The form of a catenary was mathematically described by David Gregory who independently came to the Hooke's theory and expanded it so that it can be applied to arches of final thickness. According to Gregory, arches are stable when a catenary can be placed within their thickness. The analogy with catenary was used throughout the eighteenth and nineteenth centuries for the design and analysis of stone bridges and domes. One of the most characteristic examples is the analysis of the Saint Peter's dome in Vatican made by Polenius [4] (Figure 2). A different view about the same problem was developed during the eighteenth century in France. La Haire, Couple and Coulomb described the arch as a series of rigid blocks that can have relative movements. According to Couple, a breakdown occurs when the number of joints sufficient to activate a mechanism appear inside of an arch [5]. The first general theory about the stability of arches was published by Coulomb in 1773. [6]. Coulomb developed a mathematical base to describe various types of arch collapse taking in consideration relative rotations and sliding among blocks. Coulomb considered that sliding among blocks happens rarely, and so he suggested that only the shapes of fracture caused by relative rotations of blocks should be considered. A further progress in the analysis of arches happened in the nineteenth century with the appearance of graphic statics and the theory of thrust line. The graphic statics was used to analyze many types of stone bridges and structures all the way to the beginning of the twentieth century. For instance Snell [7] used the graphic method to analyse stability of stone arches. The force polygon that was used by Snell to make a stone arch thrust line is shown in Figure 3. Similarly, Rubio analysed the Mallorca Cathedral using the theory of thrust line [8] ( Figure 4). In order to be able to talk about the calculation and modelling of masonry structures it is necessary to consider their main feature, i.e. their composite nature due to the fact that they are made of blocks divided by joints that can but do not need to be filled with mortar. The presence of joints, which in masonry structures represent the weakest link, is the cause of their nonlinear and complex behaviour that creates difficulties in numerical modelling. For that reason we now have a variety of methods and numerical models for the analysis of masonry structures based on the degree of complexity, volume of input data, and required accuracy. As every method has its own area of use, the best method depends on the structure we wish to analyze, input data that are available, and the results that we wish to obtain, which is all coupled with the experience and level of expertise of the researcher [9]. The best method is the one that provides required information with an acceptable level of error in the most efficient manner. Each of these methods will not be further described in this paper. The paper focuses on main features and the area of use of some methods. Two basic approaches are used in numerical modelling of masonry structures: idealization using the continuum and discontinuum. The adoption of the continuum hypothesis assumes that strains and deformations across the observed structure are described by continuous functions. Connections between strains and deformations are given by the constitutive law of material behaviour. By combining the constitutive law of material behaviour with the equation of equilibrium, it is possible to obtain differential equations. The solving of these equation by satisfying boundary conditions gives the solution to the problem in terms of dislocation and strains.
Since it is mainly in differential equations that the solution is not known in analytical form, these equations, together with boundary conditions, are frequently translated using the balance principles into a variation problem, which is a fragile formulation that can be solved with approximate numerical procedures, where the most frequently used method is the finite element method.
On the other hand, the idealization using discontinuum observes the structure as a group of discrete elements which can be separated during the analysis, move freely and meet again in a mutually dynamic interaction. In this approach, discrete elements are mostly taken to be absolutely rigid, and numerical integration of equation of motion of blocks in time is most often conducted in an explicit way. Numerical models that use this approach are included in discrete element methods.
Recently there are more and more numerical models that combine the advantages of idealization of structures using the continuum and discontinuum approach. Some of these methods have been developed from finite element methods, while others have emerged from discrete elements methods. Those that have developed from finite element methods realize discontinuum advantages with contact elements that can be implemented between the finite elements mesh. Contact elements are used to describe discontinuity in the field of movement when cracks appear. On the other hand, numerical models that have developed from the discrete element method frequently show advantages of idealization with continuum in such a way that every discrete element is discretized with its own finite element mesh, so that the deformability of elements can be taken in consideration.
No matter what approach we are talking about when taking into consideration the level of simplicity and accuracy, we can talk about micro-modelling, simplified micro-modelling, and macro-modelling [1] ( Figure 5).  Overview of the methods for the modelling of historical masonry structures In macro-modelling, all points of the structure have equal characteristics (i.e. the same constitutive elements) that are obtained from the sufficiently big representative volume according to the theory of averaging or homogenization. In this case, the characteristics of the mortar and characteristics of blocks are uniformly spread across the structure, which itself is treated as a homogeneous anisotropic continuum. This approach is appropriate when analyzing larger structures because the calculation is less demanding. If the representative volume is the size of a block or smaller then micro-modelling is used because the block and mortar are discretized and modelled with finite elements, while the connection between the block and mortar is presented with contact elements. In this approach, constitutive laws of mortar and blocks are observed separately. This approach is appropriate for modelling smaller structures but not real structures because of big calculation demands. In simplified micro-modelling, the expanded block is modelled with finite elements while the mortar to block connection is described with contact elements that also represent potential cracking points. This approach is characterized by lower accuracy due to the fact that the Poisson's ratio of mortar, which greatly affects the compressive strength, is not taken in consideration.
As for the ways of treating seismic load, methods designed for the seismic analysis of masonry structures can be classified into linear, one of which is the simplified equivalent static analysis and modal analysis, and nonlinear, one of which is the static method of gradual pushing ("pushover analysis") and methods for analyzing response of structure over time. In the simplified equivalent static analysis the seismic loading is approximated with an equivalent static charge in two ways. The first way is to expose the structure to constant horizontal acceleration. This approach does not take into account the fact that during the earthquake the constant acceleration of base lasts only for a short period of time. Vibration effects that appear in the structure during the earthquake due to elastic features of materials are also neglected. This approach is appropriate for analysing stability of some types of older stone structures, such as arches, where the role of elastic features of materials is not significant. The second way is to set a distributed horizontal load along the height of a structure that increases from bottom towards the top, taking into consideration the distribution of horizontal forces caused by dynamic response of the structure. The seismic calculation of a structure using a modal analysis is mostly made using a finite element method by which eigenvectors and eigenvalues can be calculated. Considering the shape of eigenvectors and their contribution to the dynamic response of the structure, a horizontally distributed load is applied on the height of the structure that later gets combined using various procedures. The pushover analysis is based on gradual amplification of the horizontal force amplitude, in parallel with the observation of the structure's response [10]. The method of response over time consists of the calculation of the structure's response in time and considers the shape of strain, deformations and dislocations for certain entry of acceleration at the base.
The main features and the possibility of using numerical models that are most frequently applied in seismic analysis of masonry structures will be described later on in the paper.

Dynamic Analytic Methods
The purpose of these methods is to predict the structure's response during dynamic impulse, or to predict the smallest rate of horizontal acceleration of the base that will cause the collapse of a structure, using an analytical approach.
Due to the fact that in these methods structural elements are assumed to be absolutely rigid, it is assumed that the structure's collapse will not occur because the strength of material is exceeded, but due to the loss of stability only.
Because of the complexity of analytic equations, these methods are restricted to the analysis of simpler structures such as blocks on horizontal base, doorway frames, and arches.
Housner was the first to use dynamic laws to analyze the overturning of the rigid block on horizontal surface exposed to constant horizontal, sine and seismographic acceleration of the base [11]. By adopting the assumption that there is no sliding between the block and the base, that the block does not leap from the base, that the block is thin enough, i.e. that the angle α ( Figure 6) is less than twenty degrees, and that the angle of rotation Ө during the oscillation is small, Housner analyzed the necessary duration of the rectangle and sine impulse that would cause the block to overturn. Housner took into consideration the loss of energy during collision of the block with the base, and introduced the hypothesis that during the alternating rotation of blocks around one and other edge the moment of the quantity of motion is preserved. Housner also showed how the stability of a high thin block during seismic charge is a lot greater from the one that it has during constant acceleration of the base. Housner's work represented the basis for many other studies about rocking motion of blocks. Continuing on the Housner's work, Yim [12] approached the problem of the rocking motion of block on the horizontal base without introducing the presumption of small rotation angles. By solving the equation of motion using a numerical procedure, Yim tackled the problem of overturning of the block caused by seismic impulse from a probabilistic point of view, because it was observed that dynamic responses of the block to various accelerations of the base vary a lot at identical amplitude: from small oscillations to the complete overturning. While Yim [12] focused on the dynamic block response during accidental acceleration of the base, Spanos and Koh [7] placed emphasis on the dynamic block response during harmonious acceleration of the base. Following Spanos and Koh, many scientists conducted research on dynamic response of the block during harmonious acceleration of the base [13][14][15][16][17]. All these studies are based on the assumption that dynamic motion can be described over a spectrum of responses.
In order to test a much complex block behaviour during acceleration of the base, may researchers speculated about the possibility of sliding and bounding the block away from the base [18][19][20][21][22][23].
Several researchers studied response of the system made of multiple blocks, and so Sinopoly and Sepe [24] studied response of a frame structure made of three blocks and exposed to horizontal acceleration of the base. In a similar way, Spanos et al. [25] analysed dynamic behaviour of the structure made of two blocks. However, equations used to analyse dynamic behaviour of these two blocks were so complex that the scientists came to the conclusion that the discrete element method is much more appropriate for dynamic analysis of multi blocks systems. The discrete element method will be discussed in greater detail in a separate section of the paper. Besides analyzing blocks, some scientists use analytic methods to analyze the stability of arches exposed to constant horizontal acceleration of the base [26,27]. Both authors use the equivalent static analysis to determine places where joints in arch are formed, and assume that these places will remain unchanged during the dynamic response of a structure. Simplification of this system that is made of multiple blocks in the scope of the system with one degree of freedom (figure 7) enables analytic solution. The complexity of analytic methods, namely in calculation of the dynamic seismic response of more complex masonry structures, has animated scientists to develop numerical methods. The most appropriate numerical method for calculating stability of masonry structures, capable of analysing dynamic behaviour of rigid bodies, is the discrete element method.

Limit Analysis Method
Limit analysis methods are based on assumptions made by Couple in 1730: (1) Masonry structures do not have a tensile strength, (2) masonry structures have an infinite strength in compression, and (3) sliding can not occur between joints. Heyman [28] was among the first to use these assumptions in the stone arch stability analysis. Adoption of these assumptions enables the use of cinematic and static theorems of plasticity which are used to determine, for the set vector of outside loading F, the factor of charge α which is necessary to increase external load until structural collapse. According to the static theorem or low boundary theorem, a structure is stable and will not collapse if a statically acceptable field of internal force can be found for outside loads. The highest charge value for which the structure is still in balance represents the lowest level of safety factor. According to cinematic theorem, the structure is unstable if the mechanism for which the work of external forces is greater than or equal to zero can be found. A bottom level of load factor in which virtual work of outside forces is zero represents the bottom upper level of safety factor for a structure. In conclusion, the safety coefficient of a structure, obtained using either static or cinematic approach, has to be the same. Traditional limit analysis methods made use of a static approach that was based on the use of graphic statics during graphic interpretation of thrust lines of stone arches [30,31]. If the thrust line was located inside the outline of an arch, the arch would be considered statically stabile.
With the advent of computer era, Harvey and Maunder [32] used tabulated calculations to obtain three dimensional shape of thrust line while Block et al. [33,34] developed an interactive computer analysis based on the combination GRAĐEVINAR 65 (2013) 7, 603-618 Overview of the methods for the modelling of historical masonry structures of static and cinematic theorems to obtain the thrust line when treating three dimensional issues.
In addition to analysing arches, the limit analysis method was also used for other types of structures. Based on the collapse mechanism observed for old masonry structures, Giuffré [35,36] and Carocci [29] used cinematic approach, after decomposition of a structure into rigid blocks, to estimate seismic resistance of structures ( Figure 8). This method is appropriate for the analysis of structures in which inter-terrace structures are not rigidly connected with walls. Giuffré also attracted considerable attention with his proposal to combine block representation of structure with carrying capacity methods [37,38] to estimate seismic resistance of masonry structures. Roca [39] also suggested a method based on the static theorem for the analysis of reinforced masonry structures. Ochsendorf [40] used the limit analysis method for the analysis of arches with deformed base, while De Luca [41] used a finite element method combined with the limit analysis method to analyze seismic resistance of stone arches. After he found the places where joints originate using a finite element method based on the linear elastic analysis, De Luca used the limit analysis method to determine an ultimate load that would cause collapse of a structure. In recent times, many computerbased limit analysis methods have been developed [42][43][44][45] using mostly the cinematic approach. Most of these methods are based on the following assumptions: 1. masonry structures do not have a tensile strength 2. masonry structures have an infinite strength in compression 3. shear behaviour in contact among blocks is perfectly plastic 4. limit load occurs with small displacements.
Similar to limit analysis methods based on graphic approach, these methods assume rigid perfectly plastic behaviour of materials and are used to estimate the carrying capacity and to provide a better insight into the mechanisms of structural collapse. To enable mathematical use of the model of rigid-perfectly plastic behaviour of material, the flow function φ defined in the unit of strain was appropriated and in this respect the following is valid: if φ < 0 material remains rigid, if φ = 0 material becomes plastic, if φ > 0 inadmissible state of strain is in course. The group of conditions for which φ > 0 forms the surface of flow. All conditions located inside or on the surface of flow gratify the flow criteria, while conditions that fall outside of the surface of flow are considered inadmissible. Material becomes plastic for the strain conditions that are located on the surface of flow, which means that the flow direction, determined with flow function, must be defined. If the flow direction towards the surface of flow is vertical, then the associative law of flow used in traditional limit analysis methods is applied. The associative law of flow implies that the angle of dilatation is equal to the angle of friction, which is inadmissible for most stone structures in which the angle of dilatation is approximate to zero. This problem can be solved with appropriation of the non-associative law of flow that leads to a non-standard limit analysis method in which theorems of limit conditions (cinematic and static) can not be strictly used. A model developed by Orduña and Lourenço [46][47][48], which takes into consideration limited compression firmness, and Gilbert's model [49], are only some of the models based on the non-standard limit analysis method. Figure 9 shows a brick pillar whose mechanism of collapse is analysed by Orduña and Lourenço [46] using the limit analysis method and finite element method, from which it can be observed that both methods give the same results.

Finite Element Method
Due to its long tradition, the finite element method is the most used method not only for calculation of masonry structures but also for the calculation of structures in general. So far a number of numerical models based on finite element methods have been developed, and they are distinguished according to the type of finite elements with which the structure is discretized, and according to the constitutive law of materials that can be linear and nonlinear.  When it comes to calculation, the simplest way of modelling masonry structures using the finite-element method is the discretization of structures using the skeletal system and linear finite elements. Molins and Roca [50] have developed a numerical model for the analysis of space structures that is made of linear space elements of variable cross section. The material and geometric nonlinearity and Mohr-Coulomb criterion of failure in shear, are included in the model. Several simplified models that approximate a structure with an equivalent frame system [51,54] have been developed for the analysis of masonry structures. In recent times, many studies have been aimed at modelling masonry structures using macro elements [55][56][57][58][59], which contributed to the reduction in the number of degrees of freedom and duration of calculation. Every macro element can represent an entire wall or any wall with openings. The element can be approximated by macro elements that are placed so as to ensure that the connection between two macro elements is established at the place from which the crack is expected to originate ( Figure 10). Although this approach is suitable for determining a collapse mechanism and carrying capacity, it can not be used to describe behaviour of a particular structure element.
Due to difficulties in discretization of old stone structures using structural elements, but also because of the need to make a more detailed analysis, two-dimensional and three-dimensional finite elements are used for modelling masonry structures on a macro level. Macro-modelling is the most frequently used approach for the analysis of masonry structures in practical work because it provides the best balance between the calculation price, and the level of accuracy. The simplest numerical models of this type, based on linear elastic behaviour of materials, have often been used to analyse big masonry structures. This is due to the scarcity of better quality models and also to calculation demands. Due to the fact that masonry structures, because of very small ultimate tensile strengths, demonstrate an explicit nonlinear behaviour already at very small loads, the use of linear analysis in modelling of masonry structures is considered unacceptable, as it can lead to wrong results and inappropriate conclusions. However, its use can be justified in some cases when the intention is to observe behaviour of a structure until appearance of the first cracks, or to evaluate the places where first cracks could appear, which will be further analyzed in detail. All effects that appear in masonry structure, starting from initiation and propagation of cracks, and all the way to final collapse, can be determined by nonlinear analysis only. Suitable nonlinear macro models that are used to analyse masonry structures take in consideration different ultimate tensile and compressive strength values and different elastic and inelastic characteristics along the material axis and, in this way, the structure is treated as a homogeneous orthotropic continuum. Elastic and inelastic parameters Overview of the methods for the modelling of historical masonry structures of that continuum are most frequently determined on the basis of experimental study of sufficiently big samples that are exposed to the homogeneous state of strain. As an alternative to complex experimental testing conditions, single components of masonry structure (mortar and blocks) can be determined and used as input data for numerical homogenization techniques. A comprehensive survey of numerical homogenization techniques is presented in [64].
The plasticity theory and damage mechanics are the most widely known theories for formulating nonlinear constitutive law of behaviour of materials. The plasticity theory tries to describe plastic behaviour of materials when permanent deformation occurs. Although it was at first used for modelling ductile materials, it is nowadays also intensively used for other materials such as soil, concrete and masonry structures [1,[65][66][67][68]. Due to the fact that the plasticity theory is suitable for monotonic load increase only, most of the named models can not take into account the cyclic behaviour. In order to eliminate that imperfection, some scientists have implemented in the traditional theory of plasticity the most significant features of materials that characterize cyclic behaviour, such as the energy dissipation due to hysteresis, and stiffness degradation [2,69].
The main feature of the damage mechanics is the concept of damage that can be defined as degradation of elastic features of materials due to initiation and evolution of micro cracks that result in gradual decease of the surfaces that transmits inside forces. As a result of the damage process the elastic features of the material decrease. Because of its mathematical complexity, this approach is basically used by adopting the hypothesis about isotropic material. There are also a lot of numerical models that take into consideration orthotropic behaviour of materials [71][72][73]. Macro-modelling has been intensively used to analyze seismic response of complex structures like arch bridges [70,74] (Figure 12), historical buildings [75], and cathedrals [72].
A drawback of most macro models is the fact that they are not able to simulate discontinuity that appears between blocks or parts of a masonry structure. These discontinuities that are already determined, like in case of old stone structures, or that can appear later on in form of cracks, can lead to various problems such as the sliding or rotation of certain parts of a structure, separation of blocks, etc. All these effects can not be taken into consideration using a traditional finite-element method that is conceived for presentation of a structure with the use of continuum. One of the ways of solving this problem is to insert contact elements between the finite-element mesh points [76][77][78][79], which represent the places of potential cracking. In this approach, the finite element mesh must be conceived in such a way that every block is discretized with its own mesh that leads towards modelling on the micro or simplified micro level. A material nonlinearity is concentrated in contact elements, while the behaviour in finite elements is most frequently linear elastic. Nonlinear behaviour of contact elements is based on the plasticity theory [80,81] or damage mechanics [82][83][84]. If modelling on a simplified micro level is considered, then contact elements describe a nonlinear behaviour of mortar and connections between blocks and mortar, while it is assumed that there is a linear elastic behaviour inside of a block, which is described with finite elements. In case of modelling on the real micro level, the block and mortar are modelled with different finite elements between which contact elements that describe a material nonlinearity are implemented. In that case, nonlinear behaviour of mortar is described with one contact element, block behaviour with others, and the connection between block and mortar with thirds ones. Due to its excessive calculation requirements, this calculation approach is mostly used to analyse smaller details of a structure that are exposed to heterogeneous conditions of strain, or in homogenization techniques where attempts are made to obtain basic mechanical characteristics of mortar and blocks based on the constitutive law of masonry structures. Almeida [85], Pegon and Pinto [76] (Figure 13) are some of the researchers that have used the finite-element method in combination with contact elements for the analysis of masonry structures.

Discrete Element Method
The discrete element method is a group of methods defined by Cundall and Hart [86] as a computing approach that: (1) enables finite displacements and rotations of discrete bodies including their complete separation, and (2) automatically recognizes new contacts among bodies as the calculation progresses. Cundall [87] has also developed a method known as the Distinct Element Method (DEM) whose original purpose was to simulate sliding and separation of connected rock mass along already determined cracks or discontinuities. The method is based on explicit numerical integration of equations of motion of rigid blocks over time. Blocks can have an arbitrary displacement, and the method also includes mutual interaction of blocks. In addition to dynamic calculations, the method offers the possibility of obtaining static solutions using viscous dumping, just as in methods of dynamic relaxation. An increasing number of numerical models presenting features of the discrete element method has been developed over time. They have found their use in the analysis of masonry structures [88][89][90]. The main feature of the discrete element method, which enables its use in the analysis of masonry structures, is presentation of a structure as a group of individual blocks mutually connected with contact elements. This approach has enabled simulation of collapse of a structure as a result of rotational sliding among joints and impact load. A wide palette of numerical models based on the discrete-element method has been developed so far. All these models differ from one another by the discrete element shape, by the way of calculation of contact forces among discrete elements, by the method used to identify a contact, by the way of calculation of equations of motion over time, etc.
In base of a discrete element shape, it is possible to distinguish block models, in which the blocks are presented with polygonal elements [87,91], and grain models of discrete elements, in which the blocks are presented like a group of circular discs in 2D, or spheres in 3D [92]. The latter ones are appropriate for micro-modelling of soil and other granular materials. Grain models are very effective for calculation due to the fact that it suffices to calculate the distance of the centre of two discs or spheres to enable identification and interaction of contacts, while with block models that part of the calculation is much more complex. Lemos [93] was among the first ones to use grain models when analysing irregular stone structures. Here stone blocks are modelled with bigger and mortar with smaller particles (Figure 14). Various hardnesses were attributed to connections between these two types of particles depending on different hardnesses of materials. Petrinic [94] has developed a model that permits interaction of grain particles and polygonal blocks. He used this model to analyse a stone bridge and model its stone blocks with four node block elements. He modelled the filling between stone blocks with grain elements. As to calculation of contact forces between discrete elements, Cundall and Hart [86] classify contacts as rigid and soft. Soft contacts [86,92] that are mostly part of discrete element methods permit analysis between two discrete elements in contact. The value of contact force is calculated based on the size of overlapping that is regulated with penalty coefficient.
In literature, this type of contact formulation is also called smooth contact or the force-displacement formulation. Conversely, the possibility of overlapping of discrete elements is excluded with rigid contacts [91,95]. In literature, this GRAĐEVINAR 65 (2013) 7, 603-618 Overview of the methods for the modelling of historical masonry structures contact formulation is also known as the non-smooth contact formulation, and is most frequently effectuated with numerical iteration in every period [91].
Due to the presentation of a contact among discrete elements, with soft contacts it is possible to talk about concentrated and divided contact forces ( Figure 15). Concentrated contact forces that are present in many models of the discrete element method are simulated with a series of contact points. At every contact point, the force can be obtained based on the constitutive law of behaviour in contact that is most frequently written in relation to the strain-relative displacement. A sufficient number of contact points is needed to properly describe the strain along contact. In linear or surface contacts [94,96,97], the strain along the contact is described with continuous function so that numerical problems that can lead to strain concentration are avoided, which is very important in situations where the occurrence and development of cracks is simulated.

Figure 15. Distributed and concentrated contact forces
The very nature of the discrete-element method serves to describe behaviour of blocks that can have arbitrary shifts. Possible interaction among blocks, and explicit nonlinearity in contact elements, makes the solving techniques using matrix presentation less attractive and inappropriate. For that reason, most numerical models that are based on the discrete-element method use explicit numerical integration of equations of motion in time, as taken over from the molecular dynamics model. In an absolutely rigid presentation of blocks, the motion of every block is described in 2D with two translations and one rotation, while it is described in 3D with three translations and three rotations. With numerical models that appropriated soft presentation of contacts that is based on penalty method, the explicit approach leads to the need for very small time period so as to provide for numerical stability. Besides explicit approach, some numerical models use implicit approach and matrix techniques to solve equation systems [91,95]. This approach enables selection of a larger time period, but the calculation time within one time period is longer, and problems connected to convergence of a solution often appear along the way. The discrete-element method is appropriate for modelling masonry structures on the simplified micro level where the blocks are presented like discrete elements that are mutually connected, with contact elements that simulate the presence of mortar, or on the real micro level where mortar and blocks are discretized with a series of smaller elements in such a way that contact elements in block have one, contact elements in mortar other, and contact elements between mortar and block third features. The rocking motion of block on rigid base [99][100][101], the static and dynamic analysis of masonry cap wall ( Figure 16) [89,98,102], the analysis of stone bridges [103,104], stability of pillars with architrave [105][106][107], the analysis of stone arches [88,89,108], the dynamic analysis of stone bell-towers and basilicas ( Figure 17) [109,110], are only some of the examples of use of the discrete-element method in the analysis of masonry structures. . Seismic behaviour and collapse mode for two different masonry structures exposed to seismic action [109] In most of the previously mentioned numerical models that are based on the discrete-element method, blocks are mostly treated as being rigid, which makes them inappropriate for the analysis of the type of structures in which the state of strain and deformations inside a discrete element can not be ignored. The assumption of absolutely rigid blocks is appropriate for modelling those types of masonry structures which break down mostly because of the loss of stability, due to creation of a mechanism inside the structure, which often happens with old stone structures that do not have a big pre-compression stress. For this type of issues, elastic features of blocks can be concentrated in contact elements if they are soft, or they can be ignored. The state of strain and deformations inside discrete element can be considered in such a way that every discrete element gets GRAĐEVINAR 65 (2013) 7, 603-618 Hrvoje Smoljanović, Nikolina Živaljić, Željana Nikolić discretized with its own finite element mesh. In this approach, the finite-element method is used to calculate the field of strain and deformations within a discrete element, while the discrete element method is used to calculate the contact forces. Many recently presented numerical models have tried to make best use of advantages of the finite and discrete element method [91,94,96,[111][112][113]. Thus Cundall [111] and Hart [112] use deformable blocks with their own finite element mesh: triangles in 2D and pyramids in 3D. Both numerical codes have algorithms for automatic identification and contact interaction. Barbosa [96] proposes a discrete-finite element model in which deformable blocks are presented with quadrilateral isoparametric finite elements. Petrinic [94] has developed a 2D model using polygonal blocks discretized with the triangular finite element mesh and rigid discs. The method developed by Mamaghani [113], called Discrete Finite Elements, is also based on the presentation of blocks with an internal finite element mesh. Shi and Goodman [91] have developed a method called Discontinuous Deformation Analysis (DDA) in which it is assumed that the condition for strain and deformation in deformable blocks is homogeneous. An improved model of deformability has been made within this method by using basic functions of higher order with which is possible to take into consideration non homogeneous condition of strain and deformation inside a block, or using a concept of under blocks in which every single block is divided into under blocks between which the splitting is enabled [114].
One of the approaches that make use of advantages of the finite and discrete element method is the Combined Finite-Discrete Elements Method (FEM/DEM), developed by Munjiza [115,116]. The FEM/DEM method is in the first place destined to simulate a fragmentation process [97,117,118] taking into consideration deformable blocks that can crack so that, as a consequence, several blocks can originate from one block during the analysis. Blocks are discretized with their own triangular finite element mesh between which contact elements can be inserted. The material nonlinearity is modelled and, on this basis, the initiation and propagation of cracks is described. Contact forces are calculated using potential contact forces, taking into consideration the Coulomb model of dry friction. The method uses an explicit numerical integration of equations of motion in time. The FEM/DEM method has proven to be quite efficient during seismic analysis of old dry built stone structures [119].

Conclusion
The description of mechanical behaviour of masonry structures represents a real challenge because of the very nature of masonry structures which are characterized by complex and extremely nonlinear behaviour, due to presence of joints between blocks that can but do not need to be filled with mortar.
One of important factors that affect numerical modelling of masonry structures is estimation of mechanical characteristics Overview of the methods for the modelling of historical masonry structures of material. In older masonry, the mortar between the joints can deteriorate under the influence of weather conditions and due to various chemical reactions, with the resulting loss of its bonding properties. That means that the behaviour of these structures is very similar to the behaviour of old dry stone masonry structures. Changes of mechanical characteristics can happen even in stone blocks due to change in temperature, freezing and similar occurrences, and they are manifested in creation of initial cracks that are very hard to completely embrace with numerical models. The second important problem in the analysis of masonry structures is the description of their geometry. Numerous Middle Age cathedrals and churches have very complex geometry that consists of the combination of curved 1D elements (arches), 2D elements (vaults), and 3D elements (domes). Describing such very complex geometries requires 3D numerical models with a great quantity of finite elements which, combined with complex constitutive behaviour of materials, can significantly extend the calculation time.
The problem of describing geometry is further aggravated by the fact that many historical masonry structures experience deformations over time, as a consequence of seismic activity or due to unequal deformation of the foundation soil, which can only be detected by means of very accurate space surveying measurements. Due to high complexity in the modelling of masonry structures, a considerable number of numerical models has been developed in recent times. They greatly differ in the level of accuracy, amount of input data, results that need to be obtained, effects that appear in a structure during a certain action, and time of calculation. The use of a particular method is also greatly dependent on the engineer's experience. In some cases very simple methods largely dependent on the engineer's knowledge can produce required results without much effort while, on the other hand, the use of very sophisticated methods with wrongly presumed input parameters can lead us in the wrong direction. Limit analysis methods are appropriate as a means to determine failure mechanisms for structures subjected to an assumed load. Their imperfection is the impossibility to analyze dynamic response of a structure over time. To analyze dynamic response of simple structures such as standing pillars, dynamic analytic methods are appropriate but, with an increase in the number of structural elements, these methods become inappropriate because of a large number of analytic equations. In such cases, it is best to use discrete element methods that have been developed to analyse dynamic response of structures where the complete collapse is due to the loss of stability. Unlike finite element methods, these methods are not capable of describing the deformability and state of strain within the structure that can cause full structural collapse due to force that exceeds the hardness of material. Advantages of discrete and finite element methods are combined in finite-discrete element methods.
In recent times, a lot of attention is given to the development of numerical models based on micro approach that requires detailed knowledge of mechanical characteristics of materials. This type of modelling involves a highly expensive experimental research. For the time being, their use is limited only to parts of a structure, due to big calculation demands. It may reasonably be expected that in the near future, with the development of computers, these models will give the most accurate results and will be used for calculation of entire structures. Until that time, a significant place in the analysis of masonry structures will be taken by homogenization techniques that represent the connection between modelling on the micro and macro levels.