This notebook solves a problem with a buried fault and a traction free surface. There's an explanation of the continuous problem and the discretization process included.

Setup

"%pylab inline" sets matplotlib to show plots in the browser rather than opening them on the server's desktop.

Import the standard numpy, matplotlib libraries,

Import the codim1 functions.

Define the necessary parameters -- elastic parameters, number of elements, element degree, and a few parameters controlling the order of the quadrature methods.

Create a bunch of kernel objects. These will be passed to the matrix and right hand side assembly routines in order to specify what sort of terms to produce. The DisplacementKernel and TractionKernel are the classical kernels appearing in Somigliana's formula. The AdjointTractionKernel results from taking the source surface-normal derivative of the DisplacementKernel and converting to tractions using the elastic moduli. The HypersingularKernel results in the same way from the TractionKernel. The two regularized hypersingular kernels are only equivalent to the HypersingularKernel when used underneath a double integral with the appropriate basis function derivatives.

The continuous problem

Boundary element methods begin from Somigliana's formula: $$\vec{u}(x) + \int_{S} T^*(x, y)\vec{u}(y) dS = \int_{S} U^*(x,y)\vec{t}(y) dS$$ where $u$ is the displacement, $t$ is the traction, $T^*$ is the TractionKernel, $U^*$ is the DisplacementKernel. When the point $x$ is placed on the boundary, the kernels become singular. The integral of $T^*$ is no longer integrable. But, it can be interpreted as a Cauchy Principal Value (CPV) integral, where the integral is written: $$\int_{S} T^*(x, y)\vec{u}(y) dS = \lim_{\epsilon \to 0}\int_{S - B_{\epsilon}} T^*(x, y)\vec{u}(y)dS + \lim_{\epsilon \to 0} \left(\int_{B_{\epsilon}} T^*(x, y)\vec{u}(y)dS \right)$$ $$\int_{S} T^*(x, y)\vec{u}(y) dS = PV\int_{S} T^*(x, y)\vec{u}(y)dS + \lim_{\epsilon \to 0} \left(\int_{B_{\epsilon}} T^*(x, y)\vec{u}(y)dS \right)$$ where $B_{\epsilon}$ is the the ball (circle, sphere) of radius $\epsilon$ around the singular point.

It turns out that the second integral has a non-zero value in the limit: $$\lim_{\epsilon \to 0} \int_{B_{\epsilon}} T^*(x, y)\vec{u}(y)dS = -\frac{\vec{u}(x)}{2}$$

So, for a boundary point Somigliana's formula becomes: $$\frac{\vec{u}(x)}{2} + PV\int_{S} T^*(x, y)\vec{u}(y) dS = \int_{S} U^*(x,y)\vec{t}(y) dS$$ I like to call this the displacement boundary integral equation, to emphasize that the only term outside of the integrals is a displacement. Take the surface normal derivatives and multiply by the elastic moduli to get the traction boundary integral equation. This is: $$\frac{\vec{t}(x)}{2} + PV\int_{S} W^*(x, y)\vec{u}(y) dS = \int_{S} \hat{T}^*(x,y)\vec{t}(y) dS$$ where $\hat{T}^*$ is the AdjointTractionKernel, $W^*$ is the HypersingularKernel. The AdjointTractionKernel is given that name because it represents a source point traction whereas TractionKernel represents an observation point traction. In terms of their respective formulae, this just means that the AdjointTractionKernel contains terms multiplied by the source normal vector whereas TractionKernel has terms multiplied by the observation normal vector.

In this problem, I want to apply zero traction boundary conditions at the surface and displacement discontinuity boundary conditions on the fault. I use the traction boundary integral equation and separate out each surface in the problem. I also use a simplified notation... $$\frac{t}{2} + \int_{fault+} W^*u + \int_{fault-} W^*u + \int_{surface} W^*u = \int_{fault+} \hat{T}^*t + \int_{fault-} \hat{T}^*t + \int_{surface} \hat{T}^*t$$ Most of the terms end up dropping out of this equation. The $\hat{T}^*$ terms are equal and opposite for the fault and because tractions are continuous: $$\int_{fault+} \hat{T}^*t + \int_{fault-} \hat{T}^*t = \int_{fault+}\hat{T}^*(t^+ - t^-) = 0$$ The $W^*$ are also equal and opposite for the fault, but $u^+ = -u^-$, so: $$\int_{fault+} W^*u + \int_{fault-} W^*u = \int_{fault+} W^*(u^+ - u^-) = \int_{fault+} W^*\Delta u$$ The terms including surface tractions are zero by virtue of the boundary condition $$\int_{surface} \hat{T}^*t = \frac{t}{2} = 0$$ So, I am left with a simple equation to solve: $$\int_{fault+} W^*\Delta u = -\int_{surface} W^*u$$

However, there are still two problems left. First, the kernels $W^{*}$ are highly singular ($1/(r^2)$). Second, this equation is still continuous... I'm going to run a few more lines of code before solving these issues.

Meshing!

Define the fault geometry, including the angle with the surface (delta), the fault normal direction and fault tangential direction.

Create the free surface mesh. I just use a flat surface with n_elements_surface separate elements. The surface extends for about 10 fault lengths in both directions.

Define the standard objects required by the matrix and right hand side assembly routines. The basis functions are a set of Lagrange interpolating polynomials with equally spaced nodes. The quadrature strategy provides singular quadrature formulae for the singular integrations. It also provides nonsingular Gauss quadrature methods that decrease in number of points with greater distance between the source and observation points. Finally, the DOFHandler organizes the degrees of freedom on the mesh. In particular, it ensures continuity of the solution on adjacent elements.

The discretized problem

The basis functions just described above are used to discretize the continuous problem. Instead of considering arbitrary functions from some space (specifically, a certain type of Sobolev space), I take piece-wise polynomial approximations of those functions. The "pieces" in "piece-wise" are the elements in the boundary mesh. I can consider the coefficients of these polynomials in order to get a discrete representation of the solution function. But, there is still a continuous part of the problem. From the previous section, the equation I am solving: $$-\int_{surface} W^*(x,y)u(y)dy = -\sum_{j}U_j\int_{surface}W^*(x,y)\phi_j(y) dy = \int_{fault} W^*(x,y)\Delta u(y)dy$$ where $U_j$ is the coefficient of the basis function $\phi_j(y)$.

By discretizing the solution function, I have a discrete representation of $u(y)$ on the free surface. Also, $\Delta u(y)$ on the fault is assumed to be known exactly because it is a boundary condition. But, the equation still has to be satisfied for all values of $x$ on the free surface. A computer can only solve a finite set of equations... There are two approaches to solving this problem. The traditional "collocation" boundary element method chooses a finite set of values of x -- normally, these are the interpolation nodes for the solution function. However, there is no theoretical basis for this decision and it can lead to some convergence problems. Also, exactly solving the equation at a subset of points results in higher error away from those points.

A more accurate method is to find the minimal weighted residual of the equation. But, weighted by what? First, I multiply by the weighting function $v$ and integrate over $x$. $$\int_{surface}v(x)\left(-\sum_{j}U_j\int_{surface}W^*(x,y)\phi_j(y) dy\right) dx = -\int_{surface}v(x)\int_{fault} W^*(x,y)\Delta u(y)dy$$ In order to satisfy the original integral equation, this new integral equation must be true for any function $v(x)$. But, this is still continuous. So, what did I gain by doing this? The trick is to choose $v(x)$ equal to each of the basis functions of the polynomial approximation. This results in a so-called "Galerkin" method. Then, I have $N$ different solution coefficients to solve (one coefficient for each basis function) for and $N$ different equations (one equation for each choice of a basis function as the weighting function). So, the problem to solve becomes: $$ \int_{surface}\phi_i(x)\left(-\sum_{j}U_j\int_{surface}W^*(x,y)\phi_j(y) dy\right) dx = -\int_{surface}\phi_i(x) \int_{fault} W^*(x,y)\Delta u(y)dy~~~~ \forall i \in 1 ... N$$ Or: $$ \sum_{j}\left(\int_{surface}\int_{surface}\phi_i(x)W^*(x,y)\phi_j(y)dydx\right) U_j = -\int_{surface}\int_{fault}\phi_i(x) W^*(x,y)\Delta u(y)dy~~~~ \forall i \in 1 ... N$$ There are (degree + 1) polynomial basis functions per element and n_elements elements. So, there are a total of $N$ = (nelements * (degree + 1)) basis functions. The problem can be rewritten as a linear system: $$Ax = b$$ $$A{ij} = \int{surface}\int{surface}\phi_i(x)W^(x,y)\phi_j(y)dydx$$ $$x_j = U_j$$ $$bi = -\int{surface}\int_{fault} \phi_i(x) W^(x,y)\Delta u(y)dydx$$ Computers are good at solving linear systems! This weighted residual approach results in lower error than the collocation method. Most people get between 1/2th and 1/100th the error for the same number of degrees of freedom. The gain in error versus a collocation method is higher for higher order polynomials.

Regularization

The last remaining problem is that the $W^*$ kernels are very very singular. This results in integrals that are extremely difficult to compute numerically. But, the double integrals provide a simple way around this problem. A simple, conceptual way to explain the solution: The kernel represents the derivative of a traction. A traction is proportional to the derivative of a displacement. Hence, the hypersingular kernel can be written as the second derivative of another kernel. $$W^*(x,y) = \frac{d^2}{dsdt}B^*(x,y)$$ Then, the terms in the problem can be integrated by parts. It turns out that the boundary terms in the integration are zero for a continuous solution. So, $$A_{ij} = \int_{surface}\int_{surface}\frac{d \phi_i(x)}{dx}B^*(x,y)\frac{d \phi_j(y)}{dy}dydx$$ $$b_i = -\int_{surface}\int_{fault} \frac{\phi_i(x)}{dx} B^*(x,y)\frac{\Delta u(y)}{dy}dydx$$ The $B^*$ kernel is $\log(r)$ singular whereas $W^*$ is $1/(r^2)$ singular. Thus, $B^*$ is much easier to compute. This regularization process is another advantage of the Galerkin Boundary Element Method. It would not be possible without the second integral of the Galerkin method.

Actually assembling and solving the problem!

These lines assemble the vector corresponding to the term: $$b_i = -\int_{surface}\int_{fault} \frac{\phi_i(x)}{dx} B^*(x,y)\frac{\Delta u(y)}{dy}dydx$$ Because $\Delta u(y)$ is constant on the fault, the only contribution comes from jumps at the fault tips, so the integrals reduce to $$b_i = -s\int_{surface}\frac{\phi_i(x)}{dx} B^*(x,y)dx$$ where $s$ is the slip on the fault.

The integrals are computed using Gaussian quadrature.

Assemble the matrix corresponding to the term: $$A_{ij} = \int_{surface}\int_{surface}\frac{d \phi_i(x)}{dx}B^*(x,y)\frac{d \phi_j(y)}{dy}dydx$$ There are (degree + 1) basis functions per element. Thus, there are $N$ = n_elements * (degree + 1) basis functions in total. So, this matrix will be a dense matrix with $N$ rows and columns.

The current problem is actually ill-posed so far. Adding an arbitrary rigid body displacement would still satisfy the equations. This is because I only applied traction boundary conditions. The tractions would be the same whether the fault was in Kazakhstan or on Mars. This translates into a singular matrix in the discretized problem -- the null space consists of rigid body motions. To enforce a reasonable solution, I just specify that the average x and y displacements are 0. Alternatively, I could specify the exact displacement for one point. However, the mathematics suggest that a point-wise constraint is a bad idea. The problem is posed in a Sobolev space, which means that the point-wise value of the solution is irrelevant -- only integrals of the solution (and its derivatives) are important. This means that as I refine the mesh, a point-wise constraint will result in an increasing amount of error -- convergence will be destroyed.

Solve the linear system and grab the value of the solution along the boundary.

This is just the analytical solution to a half space edge dislocation. I took it from Chapter 3 of the Segall book. The boundary element solution will not converge to this solution, but it should be close because the free surface is taken to be very long.

Plots and Verification

YAY! Now, I get to actually plot the results. The exact result is plotted with asterisks. The estimated solution is plotted with solid lines. The solution gets worse further from the fault. I expected this would be the case, because the free surface is not infinite.

At some point, I could add some interior displacement and stress plots...