tests now include weakly anelliptical anisotorpy

.gitignore

@@ -52,3 +52,5 @@ ttcr_io.o
 tests/files/*_rp.vtp
 
 tests/files/*_tt_grid.vt*
+
+ambiguity-as_used_for_manuscript.py

tests/files/Grid2Drcsp_tt_grid_aniso.vtr → tests/files/Grid2Drcsp_tt_grid_elliptical.vtr

@@ -10,13 +10,13 @@
     <CellData>
     </CellData>
     <Coordinates>
-      <DataArray type="Float64" Name="Array 0x141e04080" format="binary" RangeMin="0" RangeMax="6000">
+      <DataArray type="Float64" Name="Array 0x10620e020" format="binary" RangeMin="0" RangeMax="6000">
         AQAAAACAAAAoAwAAFwEAAA==eJxNzi1MQmEUxvEbCQYiwUAkGAgG4g0GEiMYSI7AHIndOeacwd0x5phzehVUBJUroCL4cf2EQCAaDEQD4UaDgUgwGJ5/8C2/ved9zzmPYfw/CVOuSHuVe06O1qivS3OT9y3p2/wrSLfI/22ZLtG3I8O79O/Rv4+OHB8w71B6ZeZWpHPE/GNpnbCnKpOn7KvJaJ29ZzJ4zv4LOUW/IbMud0xdkgvjTfJhrEVOjLTJi6ErcmPgmvw4Q+tGfmO6I78weSs/0OzKAUZ7soPhO1nF4L0sofEgN3CK2Ufp44JHHdvo4/yTTGEFxzj3LONYxBH+YuxF5tHDH4y8ygy6OMHQm1xGBz8x8C6X0MYhznCxL62++QcopYc0
       </DataArray>
-      <DataArray type="Float64" Name="Array 0x141e16080" format="binary" RangeMin="0" RangeMax="0">
+      <DataArray type="Float64" Name="Array 0x106207550" format="binary" RangeMin="0" RangeMax="0">
         AQAAAACAAAAIAAAACwAAAA==eJxjYIAAAAAIAAE=
       </DataArray>
-      <DataArray type="Float64" Name="Array 0x141e04300" format="binary" RangeMin="0" RangeMax="3000">
+      <DataArray type="Float64" Name="Array 0x106208500" format="binary" RangeMin="0" RangeMax="3000">
         AQAAAACAAACYAQAAqQAAAA==eJxNzikOwkAYhuFKJLISiUQiR6JIJYqgCIo0CASCNIQQQggtW9lpS1ksR+AIHKFHQCIR3ysY82Rm/uWzrP9TNbIuvSb3tnx1eO9K0+O/LzOPuoGMh9SPZGNM30QWpvTP6PcxkO858xbyuWTuSgZr5ofS3bBnK50d+/aydGDvUeZP7D/LD2aRbMXcsZaQCysX8mE5JScWr+RF+0ZuzN3Jj190H+YHwvo/dQ==
       </DataArray>
     </Coordinates>

tests/files/prepare_tests_aniso.py

@@ -5,6 +5,7 @@
 import matplotlib.pyplot as plt
 import vtk
 
+import traveltime as wa_tt
 
 # %%
 
@@ -238,4 +239,136 @@
 # writer.SetInputData(ugrid)
 # writer.SetFileName('tests/files/sol_analytique_elliptical_2d_tt.vtu')
 # writer.SetDataModeToBinary()
-# writer.Write()
+# writer.Write()
+
+# %%
+
+VP0 = 1875.
+r2 = 0.1
+r4 = 0.15602005
+
+rrr = np.array([1., np.NAN, r2, np.NAN, r4])   # generic anisotropy parameters (here, for Dog Creek Shale)
+
+scs = wa_tt.precompute(v00=VP0, rrr=rrr)                                 # precompute energy velocity and construct an interpolation polynomial
+
+x = np.arange(0.0, 6000.01, 60.0)
+z = np.arange(0.0, 3000.01, 60.0)
+
+ncx = x.size - 1
+ncz = z.size - 1
+
+xx, zz = np.meshgrid(x, z)
+
+xx = xx.flatten()
+zz = zz.flatten()
+
+t = wa_tt.traveltime(xxx=xx, zzz=zz, scs=scs)
+
+t[0] = 0.0
+
+
+# %%
+xCoords = vtk.vtkDoubleArray()
+yCoords = vtk.vtkDoubleArray()
+zCoords = vtk.vtkDoubleArray()
+
+for i in np.arange(0.0, 6000.1, 60.0):
+    xCoords.InsertNextValue(i)
+for i in np.arange(0.0, 3000.1, 60.0):
+    zCoords.InsertNextValue(i)
+yCoords.InsertNextValue(0.0)
+
+rgrid = vtk.vtkRectilinearGrid()
+rgrid.SetDimensions(xCoords.GetNumberOfValues(),
+                    yCoords.GetNumberOfValues(),
+                    zCoords.GetNumberOfValues())
+rgrid.SetXCoordinates(xCoords)
+rgrid.SetYCoordinates(yCoords)
+rgrid.SetZCoordinates(zCoords)
+
+v = vtk.vtkDoubleArray()
+r2_v = vtk.vtkDoubleArray()
+r4_v = vtk.vtkDoubleArray()
+
+for i in range(ncx*ncz):
+    v.InsertNextValue(VP0)
+    r2_v.InsertNextValue(r2)
+    r4_v.InsertNextValue(r4)
+v.SetName("Velocity")
+r2_v.SetName("r2")
+r4_v.SetName("r4")
+rgrid.GetCellData().AddArray(v)
+rgrid.GetCellData().AddArray(r2_v)
+rgrid.GetCellData().AddArray(r4_v)
+
+writer = vtk.vtkXMLRectilinearGridWriter()
+writer.SetInputData(rgrid)
+writer.SetFileName('weakly_an_fine2d.vtr')
+writer.SetDataModeToBinary()
+writer.Write()
+
+
+rgrid = vtk.vtkRectilinearGrid()
+rgrid.SetDimensions(xCoords.GetNumberOfValues(),
+                    yCoords.GetNumberOfValues(),
+                    zCoords.GetNumberOfValues())
+rgrid.SetXCoordinates(xCoords)
+rgrid.SetYCoordinates(yCoords)
+rgrid.SetZCoordinates(zCoords)
+
+tt = vtk.vtkDoubleArray()
+for i in t.flatten():
+    tt.InsertNextValue(i)
+tt.SetName("Travel Time")
+rgrid.GetPointData().AddArray(tt)
+
+writer = vtk.vtkXMLRectilinearGridWriter()
+writer.SetInputData(rgrid)
+writer.SetFileName('sol_analytique_weakly_an_2d_tt.vtr')
+writer.SetDataModeToBinary()
+writer.Write()
+
+# %%
+
+npx = x.size
+
+tri = []
+for nz in range(ncz):
+    for nx in range(ncx):
+        tri.append([nz*npx + nx, nz*npx + nx+1, (nz+1)*npx + nx])
+        tri.append([nz*npx + nx+1, (nz+1)*npx + nx, (nz+1)*npx + nx + 1])
+
+ugrid = vtk.vtkUnstructuredGrid()
+pts = vtk.vtkPoints()
+
+for n in range(len(xx)):
+    pts.InsertNextPoint(xx[n], 0.0, zz[n])
+
+ugrid.SetPoints(pts)
+
+triangle = vtk.vtkTriangle()
+v = vtk.vtkDoubleArray()
+r2_v = vtk.vtkDoubleArray()
+r4_v = vtk.vtkDoubleArray()
+v.SetName("Velocity")
+r2_v.SetName("r2")
+r4_v.SetName("r4")
+
+for n in range(len(tri)):
+    triangle.GetPointIds().SetId(0, tri[n][0])
+    triangle.GetPointIds().SetId(1, tri[n][1])
+    triangle.GetPointIds().SetId(2, tri[n][2])
+    ugrid.InsertNextCell( triangle.GetCellType(), triangle.GetPointIds() )
+    v.InsertNextValue(VP0)
+    r2_v.InsertNextValue(r2)
+    r4_v.InsertNextValue(r4)
+
+ugrid.GetCellData().AddArray(v)
+ugrid.GetCellData().AddArray(r2_v)
+ugrid.GetCellData().AddArray(r4_v)
+
+writer = vtk.vtkXMLUnstructuredGridWriter()
+writer.SetInputData(ugrid)
+writer.SetFileName('weakly_an_fine2d.vtu')
+writer.SetDataModeToBinary()
+writer.Write()