#  Seismic Wave Scattering Through a Compressed Hybrid BEM/FEM Method.

Nicolás Guarín-Zapata, Juan Gomez and Juan Jaramillo
Departamento de Ingeniería Civil
Universidad EAFIT
Medellín,
Colombia
<EMAIL_ADDRESS>
###### Abstract

Approximate numerical techniques, for the solution of the elastic wave
scattering problem over semi-infinite domains are reviewed. The approximations
involve the representation of the half-space by a boundary condition described
in terms of 2D boundary element discretizations. The classical BEM matrices
are initially re-written into the form of a dense dynamic stiffness matrix and
later approximated to a banded matrix. The resulting final banded matrix is
then used like a standard finite element to solve the wave scattering problem
at lower memory requirements. The accuracy of the reviewed methods is
benchmarked against the classical problems of a semi-circular and a
rectangular canyon. Results are presented in the time and frequency domain, as
well as in terms of relative errors in the considered approximations. The main
goal of the paper is to give the analyst a method that can be used at the
practising level where an approximate solution is enough in order to support
engineering decisions.

Keywords: Wave scattering, boundary element method, finite element method,
hybrid BEM/FEM.

##  1 Introduction

One of the main challenges in the numerical simulation of wave scattering
problems, over infinite or semi-infinite domains, is the proper imposition of
radiation boundary conditions: When the analysis is performed using a full
domain method, e.g., the finite differences method (FDM) or the finite element
method (FEM), one has to render finite the computational domain, and use
artificial boundary elements to approximate radiation damping. A wide variety
of absorbing boundary elements are available in the literature, (Lysmer &
Kuhlemeyer, 1969; Engquist & Majda, 1977; Clayton & Engquist, 1977). Most of
these elements suffer deficiencies, like dependence of the performance on the
angle of the incident wave or a limited frequency band of efficient operation.
On the other hand, among the many techniques available to approximate
Sommerfeld’s radiation boundary conditions, those based on integral
formulations are the most accurate (Pao & Varatharajulu, 1976; Banerjee,
1994a; Sánchez-Sesma et al., 1993); although their coupling to existing finite
element codes introduces inconvenient features, since their resulting
coefficient matrices are asymmetric and fully populated and destroy many of
the appealing features of finite element methods. In this work we explore
approximations to the radiation boundary condition, that combine the accuracy
of integral-based formulations with the advantages of the finite element
algorithm. Our approximate representation is intended to be useful in the
solution of scattering problems of moderate size using existing finite element
codes and the resources available in current personal computers.

In general, the complexity of a formulation to approximate the radiation
boundary condition is a function of the degree of space-time non-locality
introduced in the boundary model. In the case of classical FDs and FEs, non-
locality is an unnatural condition: Consequently, in order to obtain effective
solutions, the analyst must sacrifice accuracy at the expense of non-locality;
with accurate results obtained after a trade-off between non-locality and the
use of Saint-Venant-end-effects leading to large computational models. One of
the first contributions in absorbing boundaries, but still one of the most
popular ones, is the fully-local, viscous boundary element by Lysmer &
Kuhlemeyer (1969). This boundary provides the conceptual basis of most of the
currently used elements, which exhibit various levels of efficiency with
increasing levels of non-locality, (Berenger, 1994; Soudkhah & Pak, 2012; Lee
et al., 2012; Zheng et al., 2013; Duru & Kreiss, 2012). Moreover, many of
these absorbing conditions have been used in recent FDs and FEs simulations of
seismic regions, as large as the state of California and for frequencies up to
$4.0$ Hz, (Frankel, 1993; Komatitsch & Tromp, 1999; Min et al., 2003; Bielak
et al., 2003; Komatitsch et al., 2004; Givoli, 2004; Frehner et al., 2008;
Ichimura et al., 2009; Lee et al., 2009a, b, 2010; Bielak et al., 2010;
Chaljub et al., 2010; Lan & Zhang, 2011). In a second class of methods, which
are inherently non-local in space and time, (Pao & Varatharajulu, 1976;
Manolis & Beskos, 1988; Sánchez-Sesma & Campillo, 1991; Banerjee, 1994a;
Sánchez-Sesma & Luzón, 1995; Janod & Coutant, 2000; Iturrarán-Viveros et al.,
2005), the radiation condition is imposed through fully coupled boundary
integral formulations: In those cases the half-space condition is analytically
satisfied by the Green’s function, but the resulting algorithm, although
highly accurate is also computationally demanding with coefficient matrices
which are non-symmetric and fully populated.

A third solution strategy that has been popular during many years, is a
combination that approximates radiation damping via the BEM technique, while
the scatterer, (including heterogeneities and/or nonlinearities), is treated
with an FEM algorithm (Brebbia & Georgiou, 1979; Beer & Meek, 1981; Banerjee,
1994b; Zienkiewicz et al., 1977; Zienkiewicz & Taylor, 2005; Boutchicha et
al., 2007; Helldorfer et al., 2008; Seghir et al., 2009; Bielak et al., 2009).
That approach becomes especially attractive when the BEM discretization is
written in the form of a displacements-based finite element, representing the
semi-infinite domain as a half-space-super-element (HSSE) and favouring its
straightforward coupling to a finite element model for the scatterer. Although
the coupling of BEM and FEM algorithms can be accomplished in a wide majority
of commercially available finite element codes that assimilate user element
subroutines, e.g., (Abaqus, 1989; Taylor, 2011), the asymmetric and dense
character of the HSSE-coefficient-matrix still restrains the advantages of the
finite element method. Therefore, additional modifications must be imposed to
the resulting radiation boundary approximation, in order to obtain an
efficient and accurate formulation that preserves the nice properties of the
finite element method. In this work we enforce the HSSE-coefficient-matrix to
remain banded and symmetric, and asses the effect upon the solution accuracy
of such compression scheme.

Compression schemes for the boundary element matrices have been used in the
past, but all of them performed the matrix compression at the boundary element
level,i.e. in the discrete versions of the displacements and traction kernels.
For instance, Bouchon et al. (1995) used the criteria of making zero each term
of the traction and displacement matrices, where an original absolute value
below a preselected threshold was identified; Ostrowski et al. (2006) used
hierarchical matrices and adaptive cross approximation, which is a purely
algebraic method, to compress the discrete BEM integral operators; Abe et al.
(2001) and Eppler & Harbrecht (2005) used a wavelet BEM technique to achieve
the packing of the resulting matrices. The fast multi-poles boundary element
methods (FM-BEM) introduced by Liu (2009), have gained popularity due to its
low computational cost in both, storage requirements and CPU time; these
methods are based on a far field expansion of the kernels of the integral
equations. In this work, we follow a slightly different approach since we
directly operate in the final resulting BEM dynamic stiffness matrix; this
operation is a mandatory requirement if the goal is to preserve the advantages
of the classical displacement-based finite element method. Such requirement is
clarified if one considers the fact that the inverse of a banded matrix is not
generally banded, and although the product of banded matrices is still banded,
the resulting half-band-width is larger than in the original single matrices.

In the current work we first review the well-known direct and indirect BEM
algorithms, and establish the connection between the final HSSE matrices
obtained with both methods. The resulting matrices are then compressed using
(i) a threshold criteria and (ii) a method based on the distance to the
diagonal of the different terms; this method is termed herein a half-band-
width-method. The resulting HSSEs are coupled into an existing finite element
code and tested in the solution of two well documented scattering problems,
namely a semicircular and a rectangular canyon. Both problems are rich in
scattered motions and contain different sources of diffraction: As shown by
different authors (Bielak & Christiano, 1984; Jaramillo et al., 2013) the
scattered and diffracted parts of the total field contain the motions that
must be absorbed by the HSSE. The article is organized as follows: The first
part summarizes the scattering problem in terms of integral equations, leading
to boundary element algorithms and the formulation of the HSSE coefficient
matrices, together with its coupling to existing finite element codes. We then
discuss and apply the compression criteria to the solution of both canyons.
The results, which are shown in the frequency and time domain are then
compared in terms of transfer functions and synthetic seismograms along the
canyon surfaces.

##  2 Solution of the Scattering Problem

###  2.1 Direct and Indirect BEM Formulation and coupling to the finite
element method

The general scattering problem is schematized in fig. 1. Since the scatterer
is treated with a classical FEM algorithm, we only discuss the formulation of
the boundary value problem for the half-space part of the domain. In this work
we use a method due to Bielak & Christiano (1984) where the only unknown in
the half-space corresponds to the scattered field. In this way, the external
plane wave excitation is converted into an equivalent system of internal
sources, located along the contact interface between the two media. The
problem is then formulated in terms of the Green’s function using either an
integral representation theorem or a plain superposition of sources.

Figure 1: Schematic representation of the scattering problem.

Regardless of whether a direct or indirect BEM approach is used, both methods
find their basic foundation in the concept of the fundamental solution or the
specific problem Green’s functions $G_{ij}(\vec{x},\vec{\xi})$ and
$H_{ij}(\vec{x},\vec{\xi},\hat{n}_{x})$. These functions allow us to express
respectively the $i-th$ displacement $u_{i}$ and traction components $t_{i}$
at a field point $\vec{x}$ due to a load $P_{j}$ applied in the $j-th$
direction at a point $\vec{\xi}$ according to eq. 1;

$\displaystyle u_{i}(\vec{x})$
$\displaystyle=G_{ij}\left(\vec{x},\vec{\xi}\right)P_{j}\left(\vec{\xi}\right)$
(1) $\displaystyle t_{i}(\vec{x},\hat{n}_{x})$
$\displaystyle=H_{ij}\left(\vec{x},\vec{\xi},\hat{n}_{x}\right)P_{j}\left(\vec{\xi}\right).$

and where it is clear that $H_{ij}(\vec{x},\vec{\xi},\hat{n}_{x})$ is the
tractions counterpart of the displacement Green’s function
$G_{ij}(\vec{x},\vec{\xi})$. In the above $\hat{n}_{x}$ denotes the normal
vector outward to the surface at a field point $\vec{x}$. Direct use of eq. 1
into Betti’s reciprocity theorem for the scattered and the fundamental states,
yields Somigliana identity which is the basis of the integral equation to be
solved in a direct boundary element method as given in eq. 2

$C_{ij}(\vec{\xi})u_{i}(\vec{\xi})=\int_{S}G_{ij}(\vec{x},\vec{\xi})t_{i}(\vec{x},\hat{n}_{x})\,\mathrm{d}S(\vec{x})-\int_{S}H_{ij}(\vec{x},\hat{n}_{x};\vec{\xi},\hat{n}_{\xi})u_{i}(\vec{x})\,\mathrm{d}S(\vec{x}).$
(2)

where again $\hat{n}_{x}$ is the surface outward normal unit vector; $C_{ij}$
is a coefficient that depends on the smoothness of the boundary; $u_{i}$ and
$t_{i}$ are the boundary displacements and tractions vectors respectively.
Similarly, direct application of eq. 1 after assuming that the solution field
is produced by an arbitrary distribution of source densities
$\phi_{j}(\vec{\xi})$, along the boundary of the considered domain, yields the
fundamental equations to be solved in an indirect boundary element method, eq.
3;

$\displaystyle u_{i}(\vec{x})$
$\displaystyle=\int_{S}G_{ij}(\vec{x},\vec{\xi})\phi_{j}(\vec{\xi})\,\mathrm{d}S(\vec{\xi})$
(3) $\displaystyle t_{i}(\vec{x},\hat{n}_{x})$
$\displaystyle=\dfrac{1}{2}\phi_{j}(\vec{x})+\int_{S}H_{ij}(\vec{x},\vec{\xi},\hat{n}_{\xi})\phi_{j}(\vec{\xi})\,\mathrm{d}S(\vec{\xi})$

and where the problem is now solved in a two step algorithm, where one has to
find the distribution of source densities $\phi_{j}(\vec{\xi})$ that solves
the specified boundary conditions and then finds the solution field anywhere
inside the domain.

The corresponding system of algebraic equations in either case can be obtained
after discretization of the field and the boundary into $N$ (constant)
elements. In the case of the direct BEM approach an algebraic equation is
generated after selecting an observation point $\vec{\xi}^{l}$ and performing
collocation along the $N$ nodal points $\vec{x}^{k}$ leading to;

$C_{ij}U^{l}_{i}=G_{ij}^{kl}t_{i}^{k}-H_{ij}^{kl}U_{i}^{k}.$ (4)

In eq. 4 the terms $G_{ij}^{kl}$ and $H_{ij}^{kl}$ correspond to the integrals
of the Green’s functions over the $\vec{x}^{k}$-element and evaluated at the
$\vec{\xi}^{l}$ collocation point. Similarly, in the case of the indirect
algorithm we generate an algebraic equation after selecting an observation
point $\vec{x}^{k}$ and performing collocation along the $\vec{\xi}^{l}$ nodal
points leading to;

$\displaystyle U_{i}^{k}$ $\displaystyle=\tilde{G}_{ij}^{kl}\phi_{j}^{l}$ (5)
$\displaystyle t_{i}^{k}$ $\displaystyle=\tilde{H}_{ij}^{kl}\phi_{j}^{l}$

In eq. 5 the terms $\tilde{G}_{ij}^{kl}$ and $\tilde{H}_{ij}^{kl}$ have the
same meaning as in the direct algorithm, however, it should be noticed that
the observation and field points are now reversed. Observe also that the
displacement matrices $G_{ij}^{kl}$ and $\tilde{G}_{ij}^{kl}$ are the same in
both methods because of the symmetries in the Green’s function. On the
contrary, the traction matrices have both directional and evaluation arguments
which are transposed. The corresponding collocation scheme in each case is
described in fig. 2.

Figure 2: Definition of the observation and collocation points in the direct
(left) and indirect (right) boundary element algorithms.

The problem of coupling the BEM and FEM discretizations, reduces now to the
simple exercise of writing the BEM equations–relating unknown nodal tractions
and displacements to prescribed nodal tractions and displacements–into the
equivalent form of a nodal forces-displacements relationship like in a
standard displacements based finite element method formulation. Moreover, in
the case of the discussed wave propagation problem, formulated in terms of the
scattered field over the half-space and when the used Green’s functions are
those of a semi-infinite medium, the problem simplifies even more, since there
are no remote forces originating at the far surface of the BEM domain, labeled
$S_{M}$ in fig. 3, to be carried over to the FEM domain.

Figure 3: Schematic description of the coupling between a BEM and a FEM
discretization.

In the case of plane wave incidence, the stiffness relation to be coupled and
directly implemented into an existing FEM code, reduces to the Half-Space-
Super-Element (HSSE) stiffness matrix and its associated loads vector
containing the incoming field, but expressed in terms of consistent nodal
forces , see Bielak and Christiano (1984). Referring again to fig. 3, we let
the finite element part of the domain be $\Omega_{F}$. This domain is assumed
to have been discretized into internal degrees of freedom $U_{I}$ and the
contact surface into degrees of freedom $U_{S}$. In general, we will assume
that the discretizations of the BEM and the FEM meshes have different number
of nodes and therefore compatibility must be enforced through displacement and
coupling matrices $R_{u}$ and $R_{t}$ (Tassoulas, 1988), as specified in eq.
6;

$\displaystyle U_{N}$ $\displaystyle=R_{u}U_{S}$ (6) $\displaystyle
F_{N}(\hat{n}^{*})$ $\displaystyle=R_{t}t_{S}(\hat{n})$

and where the normal vectors satisfy $\hat{n}=-\hat{n}^{*}$. Now, we can write
the discrete finite element equations corresponding to the scatterer like;

$\begin{bmatrix}K_{II}&K_{IS}\\\
K_{SI}&K_{SS}\end{bmatrix}\begin{bmatrix}U_{I}\\\
U_{S}\end{bmatrix}=\begin{bmatrix}F_{I}\\\ F_{S}(\hat{n})\end{bmatrix}$ (7)

where $F_{S}(\hat{n})$ are the consistent nodal forces along the contact
surface. Similarly the equations for the half-space, after having converted
the BEM system into a generalized force displacement relationship are written
like;

$F_{N}^{S}(\hat{n}^{*})=K_{HS}U_{S}^{S}$ (8)

and where the forces along the contact surface have been denoted like
$F_{N}^{S}(\hat{n}^{*})$. Coupling of eq. 7 and eq. 8 using the conditions
along the contact surface between displacements and forces involving the
incoming field represented by the discrete terms $U^{0}_{S}$ and
$F^{0}_{S}(\hat{n}^{*})$ and given by;

$\displaystyle U_{S}$ $\displaystyle={U^{S}_{S}}+{U^{0}_{S}}$ (9)
$\displaystyle F_{S}(\hat{n})$
$\displaystyle+{F^{S}_{S}(\hat{n}^{*})}+{F^{0}_{S}(\hat{n}^{*})}=0$

yields;

$\begin{bmatrix}K_{II}&K_{IS}\\\
K_{SI}&K_{SS}+K_{HS}\end{bmatrix}\begin{bmatrix}U_{I}\\\
U_{S}\end{bmatrix}=\begin{bmatrix}F_{I}\\\
-F_{S}^{0}(\hat{n}^{*})+K_{HS}U_{S}^{0}\end{bmatrix}$ (10)

In eq. 10 the terms $U_{S}^{0}$ and $F_{S}^{0}(\hat{n}^{*})$ represent the
incoming field in terms of displacements and forces, evaluated at the contact
surface and $K_{HS}$ describes the half-space contribution to the global
equilibrium equations. Denoting the Green tractions matrix for the direct
boundary element and indirect boundary element method like $H^{BEM}$ and
$H^{IBEM}$ respectively, yields after direct comparison of the resulting
stiffness matrices the following connection between both algorithms

$\displaystyle K_{HS}^{BEM}$ $\displaystyle=R_{t}H^{BEM}G^{-1}R_{u}$ (11)
$\displaystyle K_{HS}^{IBEM}$ $\displaystyle=R_{t}G^{-1}H^{IBEM}R_{u}$

###  2.2 Compression of the stiffness matrix

The non-local nature of the BEM matrices $G$ and $H$ in eq. 11, results in a
finite element system of equations with a dense coefficient matrix. This
result, associated to the non-local nature of the half-space superelement, is
disadvantageous for its coupling in a finite element algorithm since the HSSE
now controls the memory requirements and avoids the use of available and
robust iterative solvers based on banded stiffness matrices. In this work we
used two compression methods to impose a banded condition in the $HSSE$
stiffness matrix, so we can increase computer resources and take advantage of
iterative solvers available in most commercial codes. Although several works
have been performed related to compression of BEM equations, we directly
operate here on the resulting stiffness matrix. To the best of our knowledge a
study of how this manipulation affects accuracy has not been developed. Two
compression methods were used and defined as follows:

  * •

A threshold method, where all the terms smaller than a preselected percentage
of the maximum absolute value present in the original matrix are made zero

  * •

A half-band-width method, where all the values located outside of a
preselected distance from the diagonal are eliminated.

Both methods are based on the tendency of the BEM matrices to be diagonally
dominant due to the singularities in the Green’s functions. The threshold
criteria however, is convenient as it can be applied with a sound physical
basis, since the terms in the BEM matrices can be understood like influence
coefficients indicating how strong is the effect of a given source in the
response at a given nodal point; therefore the preselected threshold value is
in fact a source intensity coefficient. The trend of the matrices to be
diagonally dominant can be observed in fig. 4, where we plotted the appearance
of the reduced BEM stiffness matrices for the model of a semi-circular canyon
after applying the threshold criteria. The BEM stiffness matrix can be
computed and compressed inside a user element subroutine, coupled to an
existing finite element code. It suffices to specify the nodal points
conforming the contact surface $S_{M}$ and the subroutine must return the
elemental contribution (see eq. 10 and eq. 11). If the problem is directly
solved in the time domain, it suffices to implement a user element subroutine
UEL where the element contribution to the global system is either one of the
$K$ ’s specified in eq. 11 and the contribution to the loads vector is the one
appearing in eq. 10. On the other hand, if the analysis is performed in the
frequency domain, where the $K$’s are complex valued, one needs to assume that
the model has real and imaginary degrees of freedom at each node, thus
doubling the number of degrees of freedom with respect to the ones used in a
direct time domain analysis. With this assumption at hand, the implementation
of the finite element proceeds exactly like in the time domain. For instance,
such a generalized framework has been proposed by Mosquera (2013).

Figure 4: Appearance of the stiffness matrices for a semicircular canyon after
applying the threshold criteria with values (a) 0.01 (less accurate), (b)
0.001 and (c) 0.0001 (more accurate).

##  3 Results

###  3.1 Compressed Stiffness Matrix Formulation

The well known benchmark problems of a semi-circular and a rectangular canyon
were selected in order to test and evaluate the loss of accuracy of the
stiffness matrix approximation. These problems constitute a worst case
scenario for the approximation if the resulting matrices are interpreted like
absorbing boundaries. The response of both canyons to $P$ and $SV$ waves
incident at $0^{\circ}$ and $30^{\circ}$ with respect to the vertical was
calculated. For the excitation we applied an incoming field corresponding to a
Ricker pulse. As a threshold value we applied a variation from 0.0 to 0.01,
with 0.0 corresponding to the original complete matrix and with a larger value
representing a less accurate approximation. In the case of the band-width
criteria the size of the relative half-band width was varied from 0.03 to 1.0
with the complete matrix corresponding to a value of 1.0 and a null value to a
matrix with only the diagonal terms. The values of the relative half-band
width for matrices approximated with the threshold criteria are shown in table
1.

Threshold | Relative half-bandwidth | Relative storage
---|---|---
Semicircular | Rectangular | Semicircular | Rectangular
$0$ | $1$ | $1$ | $100\%$ | $100\%$
$10^{-5}$ | $0.26$ | $0.25$ | $45.2\%$ | $43.8\%$
$10^{-4}$ | $0.13$ | $0.09$ | $24.3\%$ | $17.2\%$
$10^{-3}$ | $0.06$ | $0.05$ | $11.6\%$ | $9.8\%$
$10^{-2}$ | $0.03$ | $0.02$ | $5.9\%$ | $4.0\%$
Table 1: Relative half-band-width for the dynamic matrices compressed by the
threshold criteria. Relative storage values are also presented –storing just
the terms in the band.

After imposing the banded condition we are interested in determining the
storage requirements for the new resultant matrix. The original needs
corresponding to the full matrix requires $n^{2}$ memory words. For the new
approximated matrix and storing just the terms in the band, the relative
needed storage requirement, $Rst$, can be approximated after neglecting the
words required to store the diagonal itself like;

$Rst(Rhbw)=2Rhbw-Rhbw^{2}$ (12)

and where $Rhbw$ is the relative half-band width. A storage requirement
corresponding to $50\%$ of the full matrix is equivalent to a relative half-
band width of 0.3, while 0.13 would be equivalent to $25\%$ of the original
storage requirement.

Figures 5 to 12 show the transfer functions and synthetic displacement time-
histories (synthetic seismograms) for a dimensionless frequency $\eta=\omega
L/\pi\beta=1.0$ where $\omega$ is the circular frequency, $\beta$ the shear
wave velocity and $L$ is the characteristic canyon dimension for different
cases. On the other hand we show in fig. 13 the relative errors with respect
to the square-norm for different levels of approximation and for all
considered incident waves. In general, the best results were obtained for the
vertical displacements under $P-$wave incidence and for the horizontal
displacements in the case of $SV-$wave incidence. At the same time, as was
originally expected better results are obtained for the semi-circular canyon
which contains less diffraction sources as compared to the rectangular canyon.
This behaviour is observed for both the threshold and half-band-width
criteria.

Figure 5: Transfer function for a semicircular canyon for a vertically
incident P wave at a dimensionless frequency of 1.

(a) Complete matrix (b) Threshold value=0.00001

(c) Threshold value=0.0001 (d) Threshold value=0.01

Figure 6: Synthetic seismograms for horizontal displacement for a semicircular
canyon for a vertically incident P wave for different threshold values. Figure
7: Transfer function for a rectangular canyon for a $30^{\circ}$ incident SV
wave at a dimensionless frequency of 1.

(a) Complete matrix (b) Threshold value=0.00001

(c) Threshold value=0.0001 (d) Threshold value=0.01

Figure 8: Synthetic seismograms for horizontal displacement for a rectangular
canyon for a vertically incident SV wave for different threshold values.
Figure 9: Transfer function for a semicircular canyon for a vertically
incident SV wave at a dimensionless frequency of 1.

(a) Complete matrix (b) Relative half-bandwidth=0.26

(c) Relative half-bandwidth=0.13 (d) Relative half-bandwidth=0.03

Figure 10: Synthetic seismograms for horizontal displacement for a
semicircular canyon for vertically incident SV wave different half-bandwidth
values. Figure 11: Transfer function for a rectangular canyon for a
$30^{\circ}$ incident P wave at a dimensionless frequency of 1.

(a) Complete matrix (b) Relative half-bandwidth=0.25

(c) Relative half-bandwidth=0.09 (d) Relative half-bandwidth=0.02

Figure 12: Synthetic seismograms for horizontal displacement for a
semicircular canyon for a $30^{\circ}$ incident P wave for different half-
bandwidth values.

(a) Semicircular canyon (b) Rectangular canyon

(c) Semicircular canyon (d) Rectangular canyon

Figure 13: Comparison for different matrix compression criteria: a), b) for
Threshold criterion and c), d) for Half-bandwidth criterion. In all the cases
the geometry and the incidence type of wave SV or P and angle is varied for
well known cases. The error is computed for the transfer function at a
dimensionless frequency of 1.

##  4 Conclusions and Furtherwork

Two approximate numerical solution methods for the elastic wave scattering
problem, were discussed and evaluated in order to facilitate the study of
dispersion problems at the practising engineering level in available personal
computers and within existing finite element architectures. First, the loss of
accuracy associated with two different ways of converting fully dense
stiffness matrices into banded stiffness matrix was evaluated. The original
stiffness matrix appeared after the discretization of a half-space sub-domain
with two different boundary element algorithms, e.g., a direct boundary
element method and an indirect boundary element method. The proposed
techniques are oriented to the solution of the problem at the engineering
practising level since one arrives to a stiffness matrix that can be coupled
into existing commercial finite element codes, e.g., ABAQUS and FEAP. The
compression procedure results in savings in memory requirements. Although
similar procedures have been previously used in order to compress full BEM
matrices, to the best of our knowledge there are no published works for the
compression of coupled BEM/FEM algorithms. A good approximation is obtained
for vertical displacements under P-wave incidence and horizontal displacements
under SV-wave incidence. A larger error is obtained for the rectangular canyon
because of the presence of singularity points in the tractions field along the
corners. Errors as large as $50\%$ for the considered range of approximations
were found. For an error of the order of $10\%$ a relative half-band-width of
0.13 or equivalently, a relative storage requirement of 0.25 can be used. This
is a memory saving of $75\%$.

### Aknowledgements

This project was conducted with financial support from “Departamento
Administrativo de Ciencia, Tecnología e Innovación, COLCIENCIAS’ and from
Universidad EAFIT.

## References

  * Abaqus (1989) Abaqus, H., 1989.  Karlsson and sorensen, Inc.,” User’s Manual,” Version, 4.
  * Abe et al. (2001) Abe, K., Koro, K., & Itami, K., 2001.  An h-hierarchical Galerkin BEM using Haar wavelets, Engineering Analysis with Boundary Elements, 25, 581–591.
  * Banerjee (1994a) Banerjee, P., 1994a.  The boundary element methods in engineering, McGraw-Hill, London.
  * Banerjee (1994b) Banerjee, P. K., 1994b.  The Boundary Element Methods in Engineering, Mcgraw-Hill College, 2nd edn.
  * Beer & Meek (1981) Beer, G. & Meek, J., 1981.  The coupling of boundary and finite element methods for infinite domain problems in elasticity, Boundary element method, 1, 575–591.
  * Berenger (1994) Berenger, J.-P., 1994.  A perfectly matched layer for the absorption of electromagnetic waves, Journal of computational physics, 114(2), 185–200.
  * Bielak & Christiano (1984) Bielak, J. & Christiano, P., 1984.  On the effective seismic input for non-linear soil-structure interaction systems, Earthquake engineering & structural dynamics, 12(1), 107–119.
  * Bielak et al. (2003) Bielak, J., Loukakis, K., Hisada, Y., & Yoshimura, C., 2003.  Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions, Part I: Theory, Bulletin of the Seismological Society of America, 93, 817–824.
  * Bielak et al. (2009) Bielak, J., MacCamy, R., McGhee, D., & Barry, A., 2009.  Unified Symmetric BEM-FEM for Site Effects on Ground Motion – SH Waves, Journal of Engineering Mechanics, 117, 2265–2285.
  * Bielak et al. (2010) Bielak, J., Graves, R., Olsen, K., Taborda, R., Ramírez-Guzmán, L., Day, S., Ely, G., Roten, D., Jordan, T., Maechling, P., Urbanic, J., Cui, Y., & Juve, G., 2010.  The ShakeOut earthquake scenario: Verification of three simulation sets, Geophysical Journal International, 180, 375–404.
  * Bouchon et al. (1995) Bouchon, M., Schultz, C. A., & Tokstz, M. N., 1995.  A Fast Implementation of Boundary Integral Equation Methods to Calculate the Propagation of Seismic Waves in Laterally Varying Layered Media, Bulletin of the Seismological Society of America, 85(6), 1679–1687.
  * Boutchicha et al. (2007) Boutchicha, D., Rahmani, O., Abdelkader, M., Sahli, A., & Belarbi, A., 2007.  Merging of the indirect discrete boundary elements to the finite element analysis and its application to two dimensional elastostatics problems, International Journal of Applied Engineering Research, 2, 441–452.
  * Brebbia & Georgiou (1979) Brebbia, C. & Georgiou, P., 1979.  Combination of boundary and finite elements for elastostatics, Applied Mathematical Modelling, 3, 212–220.
  * Chaljub et al. (2010) Chaljub, E., Moczo, P., Tsuno, S., Bard, P.-Y., Kristek, J., Käser, M., Stupazzini, M., & Kristekova, M., 2010.  Quantitative Comparison of Four Numerical Predictions of 3D Ground Motion in the Grenoble Valley, France, Bulletin of the Seismological Society of America, 100(4), 1427–1455.
  * Clayton & Engquist (1977) Clayton, R. & Engquist, B., 1977.  Absorbing boundary conditions for acoustic and elastic wave equations, Bulletin of the Seismological Society of America, 67(6), 1529–1540.
  * Duru & Kreiss (2012) Duru, K. & Kreiss, G., 2012.  A well-posed and discretely stable perfectly matched layer for elastic wave equations in second order formulation, Communications in Computational Physics, 11(5), 1643.
  * Engquist & Majda (1977) Engquist, B. & Majda, A., 1977.  Absorbing boundary conditions for numerical simulation of waves, Proceedings of the National Academy of Sciences, 74(5), 1765–1766.
  * Eppler & Harbrecht (2005) Eppler, K. & Harbrecht, H., 2005.  Fast wavelet BEM for 3D electromagnetic shaping, Applied Numerical Mathematics, 54, 537–554.
  * Frankel (1993) Frankel, A., 1993.  Three-dimensional simulations of ground motions in the san bernardino valley, california, for hypothetical earthquakes on the san andreas fault, Bulletin of the Seismological Society of America, 83(4), 1020–1041.
  * Frehner et al. (2008) Frehner, M., Schmalholz, S., Saenger, E., & Steeb, H., 2008.  Comparison of finite difference and finite element methods for simulating two-dimensional scattering of elastic waves, Physics of the Earth and Planetary Interiors, 171, 112–121.
  * Givoli (2004) Givoli, D., 2004.  High-order local non-reflecting boundary conditions: a review, Wave Motion, 39(4), 319–326.
  * Guarín-Zapata (2012) Guarín-Zapata, N., 2012.  Simulación Numérica de Problemas de Propagación de Ondas: Dominios Infinitos y Semi-infinitos, Master’s thesis, Universidad EAFIT.
  * Helldorfer et al. (2008) Helldorfer, B., Haas, M., & Kuhn, G., 2008.  Automatic coupling of a boundary element code with a commercial finite element system, Advances in Engineering Software, 39, 699–709.
  * Ichimura et al. (2009) Ichimura, T., Hori, M., & Bielak, J., 2009.  A hybrid multiresolution meshing technique for finite element three-dimensional earthquake ground motion modelling in basins including topography, Geophysical Journal International, 177, 1221–1232.
  * Iturrarán-Viveros et al. (2005) Iturrarán-Viveros, U., Vai, R., & Sánchez-Sesma, F., 2005.  Scattering of elastic waves by a 2-D crack using the Indirect Boundary Element Method (IBEM)., Geophysical Journal International,, 162, 927–934.
  * Janod & Coutant (2000) Janod, F. & Coutant, O., 2000.  Seismic response of three-dimensional topographies using a time-domain boundary element method, Geophysical Journal International, 142(2), 603–614.
  * Jaramillo et al. (2013) Jaramillo, J., Gomez, J., Saenz, M., & Vergara, J., 2013.  Analytic approximation to the scattering of antiplane shear waves by free surfaces of arbitrary shape via superposition of incident, reflected and diffracted rays, Geophysical Journal International, 192(3), 1132–1143.
  * Komatitsch & Tromp (1999) Komatitsch, D. & Tromp, J., 1999.  Introduction to the spectral element method for three-dimensional seismic wave propagation, Geophysical Journal International, 139(3), 806–822.
  * Komatitsch et al. (2004) Komatitsch, D., Liu, Q., Tromp, J., Süss, P., Stidham, C., & Shaw, J. H., 2004\.  Simulations of Ground Motion in the Los Angeles Basin Based upon the Spectral-Element Method, Bulletin of the Seismological Society of America,, 94(1), 187–206.
  * Lan & Zhang (2011) Lan, H. & Zhang, Z., 2011.  Three-Dimensional Wave-Field Simulation in Heterogeneous Transversely Isotropic Medium with Irregular Free Surface, Bulletin of the Seismological Society of America, 101(3), 1354–1370.
  * Lee et al. (2010) Lee, C.-H., Kihou, K., Horigane, K., Tsutsui, S., Fukuda, T., Eisaki, H., Iyo, A., Yamaguchi, H., Baron, A. Q. R., Braden, M., & Yamada, K., 2010.  Effect of K Doping on Phonons in Ba1-x Kx Fe2 As2, Journal of the Physical Society of Japan, 79, 014714.
  * Lee et al. (2012) Lee, J. H., Kim, J. K., & Tassoulas, J. L., 2012.  Dynamic analysis of foundations in a layered half-space using a consistent transmitting boundary, Earthquake and Structures, 3(3-4), 203–230.
  * Lee et al. (2009a) Lee, S.-J., Chan, Y.-C., Komatitsch, D., & Tromp, B.-S. H. J., 2009a.  Effects of Realistic Surface Topography on Seismic Ground Motion in the Yangminshan Region of Taiwan Based Upon the Spectral-Element Method and LiDAR DTM, Bulletin of the Seismological Society of America, 99(2), 681–693.
  * Lee et al. (2009b) Lee, S.-J., Komatitsch, D., Huang, B.-S., & Tromp, J., 2009b.  Effects of Topography on Seismic-Wave Propagation: An Example from Northern Taiwan, Bulletin of the Seismological Society of America, 99(1), 314–325.
  * Liu (2009) Liu, Y., 2009.  Fast Multipole Boundary Element Method: Theory and Applications in Engineering, Cambridge University Press, 1st edn.
  * Lysmer & Kuhlemeyer (1969) Lysmer, J. & Kuhlemeyer, R., 1969.  Finite dynamic model for infinite media, JOURNAL OF ENGINEERING MECHANICS-ASCE (JOURNAL OF THE ENGINEERING MECHANICS DIVISION), 859, 877\.
  * Manolis & Beskos (1988) Manolis, G. & Beskos, D., 1988.  Boundary Element Methods in Elastodynamics, Spon Press, 1st edn.
  * Min et al. (2003) Min, D.-J., Shin, C., Pratt, R. G., & Yoo, H. S., 2003.  Weighted-Averaging Finite-Element Method for 2D Elastic Wave Equations in the Frequency Domain, Bulletin of the Seismological Society of America, 93(2), 904–921.
  * Mosquera (2013) Mosquera, F., 2013.  Implementation of User Element Subroutines for Frequency Domain Analysis of Wave Scattering Problems with Commercial Finite Element Codes, Master’s thesis, Universidad EAFIT.
  * Ostrowski et al. (2006) Ostrowski, J., Andjelic, Z., Bebendorf, M., Cranganu-Cretu, B., & Smajic, J., 2006\.  Fast BEM-Solution of Laplace Problems With $\cal{H}$-Matrices and ACA, IEEE TRANSACTIONS ON MAGNETICS, 42(4), 627–630.
  * Pao & Varatharajulu (1976) Pao, Y. H. & Varatharajulu, V., 1976.  Huygens principle, radiation conditions, and integral formulas for the scattering of elastic waves, The Journal of the Acoustical Society of America, 59(6), 1361–1371.
  * Sánchez-Sesma & Campillo (1991) Sánchez-Sesma, F. & Campillo, M., 1991.  Diffraction of P, SV, and Rayleigh waves by topographic features: A boundary integral formulation, Bulletin of the Seismological Society of America, 81(6), 2234–2253.
  * Sánchez-Sesma & Luzón (1995) Sánchez-Sesma, F. & Luzón, F., 1995.  Seismic response of three-dimensional alluvial valleys for incident P, S and Rayleigh waves, Bulletin of the Seismological Society of America, 85, 269–284.
  * Sánchez-Sesma et al. (1993) Sánchez-Sesma, F., Ramos-Martínez, J., & Campillo, M., 1993.  An indirect boundary element method applied to simulate the seismic response of alluvial valleys for incident p, s and rayleigh waves, Earthquake engineering & structural dynamics, 22(4), 279–295.
  * Seghir et al. (2009) Seghir, A., Tahakourt, A., & Bonnet, G., 2009.  Coupling FEM and symmetric BEM for dynamic interaction of dam-reservoir systems, Engineering Analysis with Boundary Elements, 33, 1201–1210.
  * Soudkhah & Pak (2012) Soudkhah, M. & Pak, R., 2012.  Wave absorbing-boundary method in seismic centrifuge simulation of vertical free-field ground motion, Computers and Geotechnics, 43, 155–164.
  * Tassoulas (1988) Tassoulas, J. L., 1988.  Dynamic soil-structure interaction, in Boundary element methods in structural analysis, pp. 273–308, ASCE.
  * Taylor (2011) Taylor, R., 2011.  FEAP/Programmer Manual version 8.3, University of Berkeley.
  * Zheng et al. (2013) Zheng, Y., Huang, X., et al., 2013.  Anisotropic perfectly matched layers for elastic waves in cartesian and curvilinear coordinates, Tech. rep., Massachusetts Institute of Technology. Earth Resources Laboratory.
  * Zienkiewicz & Taylor (2005) Zienkiewicz, O. & Taylor, R., 2005.  The Finite Element Method for Solid and Structural Mechanics, vol. 2, Butterworth-Heinemann, 6th edn.
  * Zienkiewicz et al. (1977) Zienkiewicz, O. C., Kelly, D. W., & Bettess, P., 1977.  The Coupling of the Finite Element Method and Boundary Solution Procedures, International Journal of Numerical Methods in Engineering, 11(2), 355–375.