Feature/exp force fun

[krober10nd]

@nschloe Here's a throwaway PR just to show what I observed: a test comparing the two forcing functions for the rectangle with refinement case with everything else the same. Kind of jerry-rigged this so it's not very clean..

For the Persson-Strang (PS) forcing function I normally would use a psuedo-timestep of 0.10 but this goes unstable with the Bossen-Heckbert one so I used a psuedo-timestep of 0.10 here.

 import matplotlib.pyplot as plt
 import matplotlib.tri as tri
 import numpy

 import SeismicMesh

 def rect_with_refinement(h, force_function="PS"):
     domain = SeismicMesh.geometry.Rectangle((-1.0, +1.0, -1.0, +1.0))

     def f(x):
         return h + 0.1 * numpy.sqrt(x[:, 0] ** 2 + x[:, 1] ** 2)

     points, cells = SeismicMesh.generate_mesh(
         domain=domain,
         h0=h,
         edge_length=f,
         verbose=0,
         force_function=force_function,
         mesh_improvement=False,
     )
     return points, cells

 p_ps, c_ps = rect_with_refinement(0.001, force_function="PS")
 p_bh, c_bh = rect_with_refinement(0.001, force_function="BH")

 ax = plt.gca()
 plt.xlabel("Iteration count")
 plt.ylabel("mean mesh quality")
 ax.legend()
 plt.show()

 triang = tri.Triangulation(p_bh[:, 0], p_bh[:, 1], c_bh)
 plt.triplot(triang)
 plt.gca().set_aspect("equal", adjustable="box")
 plt.show()

 triang = tri.Triangulation(p_ps[:, 0], p_ps[:, 1], c_ps)
 plt.triplot(triang)
 plt.gca().set_aspect("equal", adjustable="box")
 plt.show()

The mean quality per meshing iteration.

[krober10nd]

So in summary: PS can run with a higher pseudo dt than compared to BH at least in 2D. For BH, higher than 0.15, goes nuclear.

PostScript: how do you copy and paste code from a terminal into GitHub? I never get this right.

import matplotlib.pyplot as plt
import numpy

import SeismicMesh

def rect_with_refinement(h, force_function="PS", delta_t=0.10):
      domain = SeismicMesh.geometry.Rectangle((-1.0, +1.0, -1.0, +1.0))

     def f(x):
           return h + 0.1 * numpy.sqrt(x[:, 0] ** 2 + x[:, 1] ** 2)

       points, cells = SeismicMesh.generate_mesh(
           domain=domain,
           h0=h,
           edge_length=f,
           verbose=0,
           force_function=force_function,
           mesh_improvement=False,
           delta_t=delta_t,
       )
       return points, cells

   delta_t = [0.05, 0.10, 0.15, 0.20, 0.25, 0.30]
   for dt in delta_t:
       p_ps, c_ps = rect_with_refinement(0.001, force_function="PS", delta_t=dt)
       try:
           p_bh, c_bh = rect_with_refinement(0.001, force_function="BH", delta_t=dt)
      except:
           pass

   ax = plt.gca()
   plt.xlabel("Iteration count")
   plt.ylabel("mean mesh quality")
   ax.legend()
   plt.show()