tests added for L matrix building

.github/workflows/wheels.yml

@@ -43,9 +43,12 @@ jobs:
           python -m pip install --upgrade pip
           pip install twine
 
-      - name: Build Python wheels
-        uses: pypa/[email protected]
-
+      - name: Install cibuildwheel
+        run: python -m pip install cibuildwheel==2.16.2
+
+      - name: Build wheels
+        run: python -m cibuildwheel --output-dir wheelhouse
+
       - uses: actions/upload-artifact@v3
         with:
           path: ./wheelhouse/*.whl

tests/test_grid3d.cpp

@@ -29,6 +29,9 @@
 #include <vtkXMLRectilinearGridReader.h>
 #pragma clang diagnostic pop
 
+#include <Eigen/Sparse>
+#include <unsupported/Eigen/SparseExtra>
+
 #include "Grid3D.h"
 #include "Rcv.h"
 #include "Src.h"
@@ -91,6 +94,10 @@ const char* models[] = {
     "./files/layers_medium.vtu",
     "./files/gradient_medium.vtu"
 };
+const char* models_L[] = {
+    "./files/layers_medium.vtr",
+    "./files/layers_medium.vtu",
+};
 const char* references[] = {
     "./files/sol_analytique_couches_tt.vtr",
     "./files/sol_analytique_gradient_tt.vtr",
@@ -220,6 +227,85 @@ BOOST_DATA_TEST_CASE(
     BOOST_TEST(error < 0.15);
 }
 
+BOOST_DATA_TEST_CASE(
+                     testGrid3D_L,
+                     bdata::make(models_L) * bdata::make(methods),
+                     model, method) {
+    Src<double> src("./files/src3d_in.dat");
+    src.init();
+    Rcv<double> rcv("./files/rcv3d_in.dat");
+    rcv.init(1);
+
+    input_parameters par;
+    par.method = method;
+    switch(method) {
+        case FAST_SWEEPING:
+            par.weno3 = 1;
+            break;
+        case SHORTEST_PATH:
+            par.nn[0] = 5;
+            par.nn[1] = 5;
+            par.nn[2] = 5;
+            break;
+        case DYNAMIC_SHORTEST_PATH:
+            par.radius_tertiary_nodes = 3.0;
+            par.nn[0] = 2;
+            par.nn[1] = 2;
+            par.nn[2] = 2;
+            break;
+        default:
+            // do nothing
+            break;
+    }
+    par.modelfile = model;
+    Grid3D<double,uint32_t> *g;
+    if (string(model).find("vtr") != string::npos) {
+        g = buildRectilinear3DfromVtr<double>(par, 1);
+    } else {
+        g = buildUnstructured3DfromVtu<double>(par, 1);
+    }
+    vector<vector<siv<double>>> l_data;
+    try {
+        g->raytrace(src.get_coord(), src.get_t0(), rcv.get_coord(), rcv.get_tt(0), l_data);
+    } catch (std::exception& e) {
+        std::cerr << e.what() << std::endl;
+        abort();
+    }
+
+    vector<double> slo;
+    g->getSlowness(slo);
+    Eigen::SparseMatrix<double> L(l_data.size(), slo.size());
+    vector<Eigen::Triplet<double>> coeff;
+    for ( auto n=0; n<l_data.size(); ++n ) {
+        for ( auto nn=0; nn<l_data[n].size(); ++nn ) {
+            coeff.push_back(Eigen::Triplet<double>(n, l_data[n][nn].i, l_data[n][nn].v));
+        }
+    }
+    L.setFromTriplets(coeff.begin(), coeff.end());
+
+    saveMarket(L, "./"+get_class_name(g)+"_L");
+
+    Eigen::VectorXd s(slo.size());
+    for ( auto n=0; n<slo.size(); ++n )
+        s[n] = slo[n];
+    Eigen::VectorXd tt1 = L * s;
+    auto tt2 = rcv.get_tt(0);
+
+    double error = 0.0;
+    size_t nn = 0;
+    for ( auto n=0; n<tt2.size(); ++n ) {
+        if ( tt2[n] == 0.0 )
+            continue;
+        error += std::abs(tt1[n] - tt2[n]) / tt2[n];
+        nn++;
+    }
+    error /= nn;
+
+    BOOST_TEST_MESSAGE( "\t\t" << get_class_name(g) << ", l_data - error = " << error );
+
+    BOOST_TEST(error < 0.05);
+}
+
 BOOST_DATA_TEST_CASE(
                      translate,
                      (bdata::make(models) ^ bdata::make(models_tr)) * bdata::make(methods_tr),

tests/test_tmesh3d.py

@@ -6,6 +6,7 @@
 import numpy as np
 import vtk
 from vtk.util.numpy_support import vtk_to_numpy
+from scipy.io import mmread
 
 import ttcrpy.tmesh as tm
 
@@ -50,7 +51,7 @@ def setUp(self):
         self.slowness = vtk_to_numpy(data.GetCellData().GetArray('Slowness'))
 
         self.src = np.loadtxt('./files/src.dat',skiprows=1)
-        self.src = self.src.reshape((1, 4))
+        self.src = np.roll(self.src, 1).reshape((1, 4))
         self.rcv = np.loadtxt('./files/rcv.dat',skiprows=1)
 
     def test_Mesh3Dfs(self):
@@ -83,6 +84,62 @@ def test_Mesh3Ddsp(self):
                         'DSPM accuracy failed (slowness in cells)')
 
 
+class TestMesh3Dc_L(unittest.TestCase):
+
+    def setUp(self):
+        reader = vtk.vtkXMLUnstructuredGridReader()
+        reader.SetFileName('./files/layers_medium.vtu')
+        reader.Update()
+
+        self.nodes = np.empty((reader.GetOutput().GetNumberOfPoints(), 3 ))
+        for n in range(reader.GetOutput().GetNumberOfPoints()):
+            x = reader.GetOutput().GetPoint(n)
+            self.nodes[n, 0] = x[0]
+            self.nodes[n, 1] = x[1]
+            self.nodes[n, 2] = x[2]
+
+        self.tet = np.empty((reader.GetOutput().GetNumberOfCells(), 4 ), dtype=int)
+        ind = vtk.vtkIdList()
+        for n in range(reader.GetOutput().GetNumberOfCells()):
+            reader.GetOutput().GetCellPoints(n, ind)
+            self.tet[n, 0] = ind.GetId(0)
+            self.tet[n, 1] = ind.GetId(1)
+            self.tet[n, 2] = ind.GetId(2)
+            self.tet[n, 3] = ind.GetId(3)
+
+        data = reader.GetOutput()
+        self.slowness = vtk_to_numpy(data.GetCellData().GetArray('Slowness'))
+
+        self.src = np.loadtxt('./files/src3d_in.dat', skiprows=1)
+        self.src = np.roll(self.src, 1).reshape((1, 4))
+        self.rcv = np.loadtxt('./files/rcv3d_in.dat', skiprows=1)
+
+    def test_Mesh3Dfs(self):
+        g = tm.Mesh3d(self.nodes, self.tet, method='FSM', tt_from_rp=0)
+        tt, L = g.raytrace(self.src, self.rcv, slowness=self.slowness, compute_L=True)
+        tt2 = L @ self.slowness
+        ind = tt2 != 0.0
+        err = np.sum(np.abs(tt[ind]-tt2[ind]) / tt2[ind]) / tt2[ind].size
+        self.assertLess(err, 0.01, 'FSM (L) accuracy failed (slowness in cells)')
+
+    def test_Mesh3Dsp(self):
+        g = tm.Mesh3d(self.nodes, self.tet, method='SPM', n_secondary=5, tt_from_rp=0)
+        tt, L = g.raytrace(self.src, self.rcv, slowness=self.slowness, compute_L=True)
+        tt2 = L @ self.slowness
+        ind = tt2 != 0.0
+        err = np.sum(np.abs(tt[ind]-tt2[ind]) / tt2[ind]) / tt2[ind].size
+        self.assertLess(err, 0.01, 'SPM (L) accuracy failed (slowness in cells)')
+
+    def test_Mesh3Ddsp(self):
+        g = tm.Mesh3d(self.nodes, self.tet, method='DSPM', n_secondary=2,
+                      n_tertiary=3, radius_factor_tertiary=3.0, tt_from_rp=0)
+        tt, L = g.raytrace(self.src, self.rcv, slowness=self.slowness, compute_L=True)
+        tt2 = L @ self.slowness
+        ind = tt2 != 0.0
+        err = np.sum(np.abs(tt[ind]-tt2[ind]) / tt2[ind]) / tt2[ind].size
+        self.assertLess(err, 0.01, 'DSPM (L) accuracy failed (slowness in cells)')
+
+
 class TestMesh3Dn(unittest.TestCase):
 
     def setUp(self):
@@ -110,7 +167,7 @@ def setUp(self):
         self.slowness = vtk_to_numpy(data.GetPointData().GetArray('Slowness'))
 
         self.src = np.loadtxt('./files/src.dat',skiprows=1)
-        self.src = self.src.reshape((1, 4))
+        self.src = np.roll(self.src, 1).reshape((1, 4))
         self.rcv = np.loadtxt('./files/rcv.dat',skiprows=1)
 
     def test_Mesh3Dfs(self):