Change the way top and base surfaces are made from the way layers are…

… extracted. When surfaces are made we look for active pillars. When layers are made we look for active columns.

LICENCE.md

@@ -0,0 +1,2 @@
+## License
+[GNU LGPL-3.0](https://www.gnu.org/licenses/lgpl-3.0.txt)

README

@@ -22,9 +22,6 @@ To have Git remember the credentials for you, do this
 
   Intel MKL -> Must be installed
 
-               (PAAL: Compilers and libraries for version 2016.2.181 installed under /opt/software/intel)
-
-
 2. To make an executable using CLion or another IDE
 ---------------------------------------------------
 

README.md

@@ -0,0 +1,72 @@
+## GitHub
+
+To be able to check out the repository to your file system, you must have
+generated a personal access token in GitHub. This token is used in place
+of your personal password when accessing the repo
+
+  https://github.com/settings/tokens
+
+To have Git remember the credentials for you, do this
+
+  $ git config --global credential.helper cache                  # (activate memory cache)
+  $ git config --global credential.helper 'cache --timeout=3600' # (timeout in seconds)
+
+
+## Prerequisites
+
+  Boost     -> version 1.65.1 is included under 3rd-party
+               NB! If another version of Boost is installed in the operating
+               system cmake may become confused and stop working.
+
+  Intel MKL -> Must be installed
+
+## To make an executable using CLion or another IDE
+
+Go to the directory containing the SeismicForward git repository
+
+```
+clion CMakeLists.txt
+```
+
+## To make an executable using cmake
+
+Go to the directory above the directory containing the SeismicForward git repository
+
+```
+cd dir-above-SeismicForward
+```
+Make directory where you want the project and executable to be build
+
+```
+mkdir my-proj-dir
+```
+
+Go to this directory and run run_cmake.sh to set up compiler and library dependencies
+
+```
+cd my-proj-dir
+../SeismicForward/run_cmake.sh
+```
+
+Generate the executable
+
+```
+make
+```
+
+## Run the tests
+
+To run all tests do
+
+```
+../SeismicForward/TestScript.pl
+```
+
+To run a single test do
+
+```
+../SeismicForward/TestScript.pl case=1
+```
+
+NBNB! If more than one angle is involved, I think that results for the first
+angle only is checked. This is the way NRLib currently works.

geo2seis/seismic_parameters.cpp

@@ -294,13 +294,14 @@ void SeismicParameters::FindTopAndBaseSurfaces(NRLib::RegularSurface<double> & t
   topeclipse = NRLib::RegularSurface<double>(x0, y0, lx, ly, nx, ny, missing);
   boteclipse = NRLib::RegularSurface<double>(x0, y0, lx, ly, nx, ny, missing);
 
-  double etdx     = topeclipse.GetDX();
-  double etdy     = topeclipse.GetDY();
-  double ebdx     = boteclipse.GetDX(); // equals etdx
-  double ebdy     = boteclipse.GetDY(); // equals etdy
-  double angle    = 0.0;
-  bool   cornerpt = model_settings->GetUseCornerpointInterpol();
-  bool   bilinear = false;
+  double etdx       = topeclipse.GetDX();
+  double etdy       = topeclipse.GetDY();
+  double ebdx       = boteclipse.GetDX(); // equals etdx
+  double ebdy       = boteclipse.GetDY(); // equals etdy
+  double angle      = 0.0;
+  bool   cornerpt   = model_settings->GetUseCornerpointInterpol();
+  bool   bilinear   = false;
+  bool   is_surface = true;
 
   NRLib::Grid2D<bool> mask(nx, ny, true);
   FindExtrapolationRegion(mask, *seismic_geometry, x0, y0, etdx, etdy);  // Setup extrapolation mask grid
@@ -313,18 +314,20 @@ void SeismicParameters::FindTopAndBaseSurfaces(NRLib::RegularSurface<double> & t
   NRLib::Grid2D<double> tvalues(nx, ny, missing);
   NRLib::Grid2D<double> bvalues(nx, ny, missing);
 
-  eclipse_geometry.FindLayer(tvalues, mask, top_k, 0, etdx, etdy, x0, y0, 0.0, cornerpt, bilinear, missing);
-  eclipse_geometry.FindLayer(bvalues, mask, bot_k, 1, ebdx, ebdy, x0, y0, 0.0, cornerpt, bilinear, missing);
+  eclipse_geometry.FindLayer(tvalues, mask, top_k, 0, etdx, etdy, x0, y0, 0.0, cornerpt, bilinear, is_surface, missing);
+  eclipse_geometry.FindLayer(bvalues, mask, bot_k, 1, ebdx, ebdy, x0, y0, 0.0, cornerpt, bilinear, is_surface, missing);
 
   ExtrapolateLayer(tvalues, mask, etdx, etdy, missing);
   ExtrapolateLayer(bvalues, mask, ebdx, ebdy, missing);
 
   double avg_top = tvalues.FindAvg(missing);
   double avg_bot = bvalues.FindAvg(missing);
 
+  assert(avg_top != missing);
+  assert(avg_bot != missing);
+
   for (size_t i = 0; i < topeclipse.GetNI(); i++) {
     for (size_t j = 0; j < topeclipse.GetNJ(); j++) {
-
       if (tvalues(i, j) == missing && bvalues(i, j) == missing) {
         topeclipse(i, j) = avg_top;
         boteclipse(i, j) = avg_bot;

geo2seis/seismic_regridding.cpp

@@ -181,19 +181,9 @@ void SeismicRegridding::FindZValues(SeismicParameters & seismic_parameters,
   const bool                            bilinear               = false;
   const double                          missing                = -999.0;
 
-
   NRLib::RegularSurface<double>         orig_top(topeclipse);
   NRLib::RegularSurface<double>         orig_bot(boteclipse);
 
-  //const double                          z_top_shift      = model_settings->GetZWaveletTop();
-  //const double                          z_bot_shift      = model_settings->GetZWaveletBot();
-
-  //orig_top.Add(     z_top_shift);
-  //orig_bot.Add(-1 * z_bot_shift);
-
-  //
-  // XXXxxxx: NBNB-PAL Dette her er kanskje litt juks????
-  //
   for (size_t i = 0 ; i < orig_top.GetNI() ; ++i)
     for (size_t j = 0 ; j < orig_top.GetNJ() ; ++j)
       if (orig_bot(i,j) < orig_top(i,j))

nr/nrlib/eclipsegrid/eclipsegeometry.cpp

@@ -1400,6 +1400,8 @@ void EclipseGeometry::FindRegularGridOfZValues(NRLib::StormContGrid
 
   NRLib::Grid2D<bool> dummy_mask(ni, nj, false);
 
+  bool         is_surface = false; // We generate layers and not top and bot Eclipse grid surfaces in this case
+
   float        monitor_size;
   float        next_monitor;
   float        count;
@@ -1433,6 +1435,7 @@ void EclipseGeometry::FindRegularGridOfZValues(NRLib::StormContGrid
             angle,
             cornerpoint_interpolation,
             bilinear_else_triangles,
+            is_surface,
             missingValue);
 
   //
@@ -1456,6 +1459,7 @@ void EclipseGeometry::FindRegularGridOfZValues(NRLib::StormContGrid
               angle,
               cornerpoint_interpolation,
               bilinear_else_triangles,
+              is_surface,
               missingValue);
     Monitor(nk - 1, nk - 2 - k,  monitor_size, next_monitor, count, carets);
   }
@@ -1659,6 +1663,7 @@ void EclipseGeometry::FindLayer(NRLib::Grid2D<double>     & z_grid,
                                 const double                angle,
                                 const bool                  cornerpoint_interpolation,
                                 const bool                  bilinear_else_triangles,
+                                const bool                  is_surface,
                                 const double                missingValue) const
 {
   if (cornerpoint_interpolation)
@@ -1682,6 +1687,7 @@ void EclipseGeometry::FindLayer(NRLib::Grid2D<double>     & z_grid,
                                       y0,
                                       angle,
                                       bilinear_else_triangles,
+                                      is_surface,
                                       missingValue);
 }
 
@@ -1843,6 +1849,7 @@ void EclipseGeometry::FindLayerCenterPointInterpolation(NRLib::Grid2D<double>
                                                         const double                y0,
                                                         const double                angle,
                                                         const bool                  bilinear_else_triangles,
+                                                        const bool                  is_surface,
                                                         const double                missingValue) const
 {
   double                    cosA = cos(angle);
@@ -1855,16 +1862,31 @@ void EclipseGeometry::FindLayerCenterPointInterpolation(NRLib::Grid2D<double>
 
   for (size_t j = 0 ; j < nj_ - 2 ; j++) { // Loops over each i,j trace in Eclipse grid
     for (size_t i = 0 ; i < ni_ - 2 ; i++) {
-      if (IsColumnActive(i  , j  ) &&   //  o-----o-----o  (i+2, j+2)
+      bool active;
+      if (is_surface)
+        active =
+          IsPillarActive(i  , j  ) &&
+          IsPillarActive(i+1, j  ) &&
+          IsPillarActive(i  , j+1) &&
+          IsPillarActive(i+1, j+1) &&
+          IsPillarActive(i+2,   j) &&
+          IsPillarActive(i+2, j+1) &&
+          IsPillarActive(i  , j+2) &&
+          IsPillarActive(i+1, j+2) &&
+          IsPillarActive(i+2, j+2);
+      else
+        active =
+          IsColumnActive(i  , j  ) &&   //  o-----o-----o  (i+2, j+2)
           IsColumnActive(i+1, j  ) &&   //  | C2  | C3  |
           IsColumnActive(i  , j+1) &&   //  o-----o-----o
           IsColumnActive(i+1, j+1) &&   //  | C0  | C1  |
           IsColumnActive(i+2,   j) &&   //  o-----o-----o  (i+2, j)
           IsColumnActive(i+2, j+1) &&
           IsColumnActive(i  , j+2) &&
           IsColumnActive(i+1, j+2) &&
-          IsColumnActive(i+2, j+2)) {
+          IsColumnActive(i+2, j+2);
 
+      if (active) {
         //NRLib::LogKit::LogFormatted(NRLib::LogKit::Low, "i,j =  %d, %d\n",i,j);
 
         C0 = FindPointCellSurface(i  , j  , k, lower_or_upper, 0.5, 0.5); // Find centre of eclipse grid cell (i  , j  )

nr/nrlib/eclipsegrid/eclipsegeometry.hpp

@@ -150,6 +150,7 @@ class EclipseGeometry {
                  const double                angle,
                  const bool                  cornerpoint_interpolation,
                  const bool                  bilinear_else_triangles,
+                 const bool                  is_surface,
                  const double                missingValue) const;
 
   void TranslateAndRotate(NRLib::Point       & corners,
@@ -284,6 +285,7 @@ class EclipseGeometry {
                                           const double            y0,
                                           const double            angle,
                                           const bool              bilinear_else_triangles,
+                                          const bool              is_surface,
                                           const double            missingValue) const;
 
   ///Function used by FindLayerSurface to fill in values to z_grid in the area inside the (NB) four corners (listed clockwise)