<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi
Step1:  Triangular mesh generation
In the last lesson we learned how to create a quad mesh by Transfinite Interpolation to accurately approximate the strong topography of a sea dike. We can use this mesh for spectral element modelling. But what should we do, if we need a triangular mesh for example for Finite Element modelling?
Yigma Tepe - another problem with strong topography
You might think that the sea dike topography was already quite complex. Well, here is another problem we are currently working on at the "Applied Geophysics" group at Christian-Albrechts University Kiel
Step2:  This quad mesh is already able to accurately describe the free-surface topography. 
Triangular mesh generation
If we need a triangular mesh, for example for finite element or finite volume modelling,  we could apply Delaunay triangulation to the node point distribution of the Yigma Tepe TFI mesh. For further details related to Delaunay triangulation, I refer to Computational Geometry in Python by Francisco Blanco-Silva. 
In a first step, we assemble the x- and z-vectors in a (NX*NZ x 2) matrix.
Step3:  Next, we compute the Voronoi diagram for the mesh points. This describes the partitioning of a plane with n points into convex polygons such that each polygon contains exactly one generating point and every point in a given polygon is closer to its generating point than to any other.
Step4:  The Delaunay triangulation creates triangles by connecting the points in neighbouring Voronoi cells.
Step5:  Let's take a look at the final mesh for the Yigma Tepe model
Step6:  The regular triangulation within the tumulus looks reasonable. However, the Delaunay triangulation also added unwanted triangles above the topography. To solve this problem we have to use constrained Delaunay triangulation in order to restrict the triangulation to the model below the free-surface topography. Unfortunately, constrained Delaunay triangulation is not available within SciPy. 
However, there exists a python wrapper to the mesh generator Triangle by Richard Shewchuck. This wrapper can be installed by 
git clone --depth=1 https
Step7:  In order to use the constrained Delaunay triangulation, we obviously have to define the constraining vertex points lying on the boundaries of our model. In this case it is quite easy, because the TFI mesh is regular. 
OK, perhaps not so easy, because we have to be sure that no redundant points are in the final list and the points are sequentially defined in clockwise direction.
Step8:  The above code looks a little bit chaotic, but you can check that the points in the resulting array model_bound are correctly sorted and contains no redundant points.
Step9:  Good, now we have defined the model boundary points. Time for some constrained Delaunay triangulation ...
Step10:  Very good, compared to the SciPy Delaunay triangulation, no triangles are added above the topography. However, most triangles have very small minimum angles, which would lead to serious numerical issues in later finite element modelling runs. So in the next step we restrict the minimum angle to 20° using the option q20.
Step11:  Finally, we want a more evenly distribution of the triangle sizes. This can be achieved by imposing a maximum area to the triangles with the option a20.
<ASSISTANT_TASK:>
Python Code:
# Execute this cell to load the notebook's style sheet, then ignore it
from IPython.core.display import HTML
css_file = '../style/custom.css'
HTML(open(css_file, "r").read())
Explanation: Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi
End of explanation
# Import Libraries
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# Here, I introduce a new library, which is useful 
# to define the fonts and size of a figure in a notebook
from pylab import rcParams
# Get rid of a Matplotlib deprecation warning
import warnings
warnings.filterwarnings("ignore")
# Load Yigma Tepe quad mesh created by TFI
X = np.loadtxt('data/yigma_tepe_TFI_mesh_X.dat', delimiter=' ', skiprows=0, unpack='True')
Z = np.loadtxt('data/yigma_tepe_TFI_mesh_Z.dat', delimiter=' ', skiprows=0, unpack='True')
# number of grid points in each spatial direction
NZ, NX = np.shape(X)
print("NX = ", NX)
print("NZ = ", NZ)
# Define figure size
rcParams['figure.figsize'] = 10, 7
# Plot Yigma Tepe TFI mesh
plt.plot(X, Z, 'k')
plt.plot(X.T, Z.T, 'k')
plt.plot(X, Z, 'bo', markersize=4)
plt.title("Yigma Tepe TFI mesh" )
plt.xlabel("x [m]")
plt.ylabel("z [m]")
plt.axes().set_aspect('equal')
#plt.savefig('yigma_tepe_TFI.pdf', bbox_inches='tight', format='pdf')
plt.show()
Explanation: Triangular mesh generation
In the last lesson we learned how to create a quad mesh by Transfinite Interpolation to accurately approximate the strong topography of a sea dike. We can use this mesh for spectral element modelling. But what should we do, if we need a triangular mesh for example for Finite Element modelling?
Yigma Tepe - another problem with strong topography
You might think that the sea dike topography was already quite complex. Well, here is another problem we are currently working on at the "Applied Geophysics" group at Christian-Albrechts University Kiel:
<img src="images/yigmatepe_1.jpg" width="100%">
This is Yigma Tepe, a tumulus which might be a burial of an Attilad King located in Pergamon (Turkey). Extensive geophysical prospecting and archaeological excavations revealed small scale near-surface structures, which makes Yigma Tepe an ideal target to further develop 2D SH and 3D seismic full waveform inversion in the field of archaeological prospection. A critical part is the correct discretization of the complex free-surface topography, crucial for accurate modelling of surface waves. 
Let's take a look at a 2D quad mesh created by Transfinite-Interpolation for a SH-profile along the tumulus.
End of explanation
# Reshape X and Z vector
x = X.flatten()
z = Z.flatten()
# Assemble x and z vector into NX*NZ x 2 matrix
points = np.vstack([x,z]).T
Explanation: This quad mesh is already able to accurately describe the free-surface topography. 
Triangular mesh generation
If we need a triangular mesh, for example for finite element or finite volume modelling,  we could apply Delaunay triangulation to the node point distribution of the Yigma Tepe TFI mesh. For further details related to Delaunay triangulation, I refer to Computational Geometry in Python by Francisco Blanco-Silva. 
In a first step, we assemble the x- and z-vectors in a (NX*NZ x 2) matrix.
End of explanation
# calculate and plot Voronoi diagram for mesh points
from scipy.spatial import Voronoi, voronoi_plot_2d
vor = Voronoi(points)
plt.figure(figsize=(12,6))
ax = plt.subplot(111, aspect='equal')
voronoi_plot_2d(vor, ax=ax)
plt.title("Part of Yigma Tepe (Voronoi diagram)" )
plt.xlabel("x [m]")
plt.ylabel("z [m]")
plt.xlim( 25,  75)
plt.ylim(10, 35)
plt.show()
Explanation: Next, we compute the Voronoi diagram for the mesh points. This describes the partitioning of a plane with n points into convex polygons such that each polygon contains exactly one generating point and every point in a given polygon is closer to its generating point than to any other.
End of explanation
# Apply Delaunay triangulation to the quad mesh node points
from scipy.spatial import Delaunay
tri = Delaunay(points)
plt.figure(figsize=(12,6))
ax = plt.subplot(111, aspect='equal')
voronoi_plot_2d(vor, ax=ax)
plt.triplot(points[:,0], points[:,1], tri.simplices.copy(), linewidth=3, color='b')
plt.title("Part of Yigma Tepe (Voronoi diagram & Delaunay triangulation)" )
plt.xlabel("x [m]")
plt.ylabel("z [m]")
plt.xlim( 25,  75)
plt.ylim(10, 35)
plt.show()
Explanation: The Delaunay triangulation creates triangles by connecting the points in neighbouring Voronoi cells.
End of explanation
# Plot triangular mesh
plt.triplot(points[:,0], points[:,1], tri.simplices.copy())
plt.title("Yigma Tepe Delaunay mesh" )
plt.xlabel("x [m]")
plt.ylabel("z [m]")
plt.axes().set_aspect('equal')
plt.show()
Explanation: Let's take a look at the final mesh for the Yigma Tepe model
End of explanation
# import triangulate library
from triangle import triangulate, show_data, plot as tplot
import triangle
Explanation: The regular triangulation within the tumulus looks reasonable. However, the Delaunay triangulation also added unwanted triangles above the topography. To solve this problem we have to use constrained Delaunay triangulation in order to restrict the triangulation to the model below the free-surface topography. Unfortunately, constrained Delaunay triangulation is not available within SciPy. 
However, there exists a python wrapper to the mesh generator Triangle by Richard Shewchuck. This wrapper can be installed by 
git clone --depth=1 https://github.com/drufat/triangle.git
cd triangle
python setup.py install
A detailed documentation of the Triangle Python wrapper by Dzhelil Rufat can be found here.
End of explanation
# Estimate boundary points 
# surface topography
surf = np.vstack([X[9,:-2],Z[9,:-2]]).T
# right model boundary
right = np.vstack([X[1:,69],Z[1:,69]]).T
# bottom model boundary
bottom = np.vstack([X[0,1:],Z[0,1:]]).T
# left model boundary
left = np.vstack([X[:-2,0],Z[:-2,0]]).T
# assemble model boundary
model_stack = np.vstack([surf,np.flipud(right)])
model_stack1 = np.vstack([model_stack,np.flipud(bottom)])
model_bound = np.vstack([model_stack1,left])
Explanation: In order to use the constrained Delaunay triangulation, we obviously have to define the constraining vertex points lying on the boundaries of our model. In this case it is quite easy, because the TFI mesh is regular. 
OK, perhaps not so easy, because we have to be sure that no redundant points are in the final list and the points are sequentially defined in clockwise direction.
End of explanation
plt.plot(model_bound[:,0],model_bound[:,1],'bo')
plt.title("Yigma Tepe model boundary" )
plt.xlabel("x [m]")
plt.ylabel("z [m]")
plt.axes().set_aspect('equal')
plt.show()
Explanation: The above code looks a little bit chaotic, but you can check that the points in the resulting array model_bound are correctly sorted and contains no redundant points.
End of explanation
# define vertices (no redundant points)
vert = model_bound
# apply Delaunay triangulation to vertices
tri = triangle.delaunay(vert)
# define vertex markers
vertm = np.array(np.zeros((len(vert),1)),dtype='int32')
# define how the vertices are connected, e.g. point 0 is connected to point 1, 
# point 1 to point 2 and so on ... 
points1 = np.arange(len(vert))
points2 = np.arange(len(vert))+1
# last point is connected to the first point
points2[-1] = 0
# define connectivity of boundary polygon
seg = np.array(np.vstack([points1,points2]).T,dtype='int32')
# define marker for boundary polygon
segm = np.array(np.ones((len(seg),1)),dtype='int32')
# assemble dictionary for triangle optimisation
A = dict(vertices=vert, vertex_markers=vertm, segments=seg, segment_markers=segm,triangles=tri)
# Optimise initial triangulation
cndt = triangle.triangulate(A,'pD')
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax,**cndt)
Explanation: Good, now we have defined the model boundary points. Time for some constrained Delaunay triangulation ...
End of explanation
cncfq20dt = triangulate(A,'pq20D')
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax,**cncfq20dt)
Explanation: Very good, compared to the SciPy Delaunay triangulation, no triangles are added above the topography. However, most triangles have very small minimum angles, which would lead to serious numerical issues in later finite element modelling runs. So in the next step we restrict the minimum angle to 20° using the option q20.
End of explanation
cncfq20adt = triangulate(A,'pq20a20D')
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax,**cncfq20adt)
Explanation: Finally, we want a more evenly distribution of the triangle sizes. This can be achieved by imposing a maximum area to the triangles with the option a20.
End of explanation
<END_TASK>