Merge branch 'feature' into 'master'

[363-add-mixedlagrangetransfer] removed class Transfers, added new getTransfer...

See merge request mpp/mpp!372

.gitlab-ci.yml

@@ -40,6 +40,8 @@ build-mpp:
       --no-cache -t ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_MPP}
       -f docker/${IMAGE_NAME_MPP}.baseimage .
     - docker push ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_MPP}
+  after_script:
+    - docker rmi ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_MPP}
   dependencies: [  ]
   tags: [ docker ]
 
@@ -55,6 +57,8 @@ ctest-mpp:
   script:
     - cd /mpp/build
     - ctest
+  after_script:
+    - docker rmi ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_MPP}
   dependencies: [ "build-mpp" ]
   tags: [ docker ]
 
@@ -70,6 +74,8 @@ mpitest-mpp:
   script:
     - cd /mpp/build
     - python3 mppyrun.py --mpi_tests=1 --mute=1
+  after_script:
+    - docker rmi ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_MPP}
   dependencies: [ "build-mpp" ]
   tags: [ docker ]
 
@@ -84,6 +90,8 @@ mppytest-mpp:
   script:
     - cd /mpp/python/mppy/tests/
     - python3 test_log_utilities.py
+  after_script:
+    - docker rmi ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_MPP}
   dependencies: [ "build-mpp" ]
   tags: [ docker ]
 
@@ -177,6 +185,9 @@ deploy-mpp:
     - docker tag ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_MPP}
       ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_RELEASE}/mpp-${OS}${OS_VERSION}:${CI_COMMIT_TAG}
     - docker push ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_RELEASE}/mpp-${OS}${OS_VERSION}:${CI_COMMIT_TAG}
+  after_script:
+    - docker rmi ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_MPP}
+    - docker rmi ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_RELEASE}/mpp-${OS}${OS_VERSION}:${CI_COMMIT_TAG}
   only: [ tags ]
   dependencies: [ "ctest-mpp", "mpitest-mpp", "trigger-tutorial", "trigger-mluq",
    "trigger-cardmech", "trigger-navierstokes", "trigger-spacetime", "trigger-dgwave", trigger-fwi]

python/mppy/mppy.py

@@ -81,10 +81,12 @@ class Mpp:
                 for cid in reversed(sorted(kernels.keys())):
                     f.write("ma-pde{} slots={}\n".format(cid, kernels[cid]))
             kernel_count = str(sum(kernels.values()))
-            run_parameters = ['mpirun', '-hostfile', 'hosts', '-np', kernel_count,
+            run_parameters = ['mpirun', '--use-hwthread-cpus',
+                              '-hostfile', 'hosts', '-np', kernel_count,
                               executable]
         else:
-            run_parameters = ['mpirun', '-np', str(kernels), executable]
+            run_parameters = ['mpirun', '--use-hwthread-cpus',
+                              '-np', str(kernels), executable]
         if config:
             run_parameters.append(config)
         for key, arg in kwargs.items():

python/mppy/utilities/log_utilities.py

@@ -94,14 +94,18 @@ class LogParser:
             return None
         key, value = [k.strip() for k in splitted]
 
-        if len(vec := self.vec_find_pattern.findall(value)) > 0:
+        vec = self.vec_find_pattern.findall(value)
+        if len(vec) > 0:
             value = [float(v.strip()) for v in vec if v.strip()]
             self.add_to_dic(key, value, data)
-
-        elif len(floats := self.float_pattern.findall(value)) > 0:
-            self.add_to_dic(key=key, value=float(floats[0]), data_dict=data)
-        elif len(string := self.string_pattern.findall(value)) > 0:
-            self.add_to_dic(key=key, value=string[0], data_dict=data)
+        else:
+            floats = self.float_pattern.findall(value)
+            if len(floats) > 0:
+                self.add_to_dic(key=key, value=float(floats[0]), data_dict=data)
+            else:
+                string = self.string_pattern.findall(value)
+                if len(string) > 0:
+                    self.add_to_dic(key=key, value=string[0], data_dict=data)
 
     @staticmethod
     def add_to_dic(key, value, data_dict):

src/lib0_basic/Point.cpp

@@ -11,8 +11,8 @@ namespace mpp_geometry {
     return tolerance;
   }
 
-  bool isNear(double a, double b) {
-    return abs(a - b) < tolerance;
+  bool isNear(double a, double b, double tol) {
+    return abs(a - b) < tol;
   }
 
 #ifdef BUILD_IA

src/lib0_basic/Point.hpp

@@ -22,7 +22,7 @@ namespace mpp_geometry {
 
   double GetTolerance();
 
-  bool isNear(double a, double b);
+  bool isNear(double a, double b, double tol=GeometricTolerance);
 
 #ifdef BUILD_IA
 
@@ -183,6 +183,15 @@ public:
     return !(*this == y);
   }
 
+  constexpr bool isNear(const PointT &other, double tol) const {
+    for (int n = 0; n < sDim + tDim; ++n) {
+      if (!mpp_geometry::isNear(z[n], other.z[n], tol))
+        return false;
+    }
+
+    return true;
+  }
+
   constexpr bool operator<(const PointT &other) const {
     if constexpr(tDim > 0) {
       if (z[sDim] < other.z[sDim] - GeometricTolerance) return true;

src/lib0_basic/utility/Config.hpp

@@ -98,6 +98,10 @@ protected: //to use constructors in tests
   void print(const std::string &);
 
 public:
+  bool exists(const std::string &key) {
+    return configMap.find(key) != configMap.end();
+  }
+
   template<typename T>
   bool get(const std::string &key, T &out, bool mute = false) {
     if (configMap.find(key) == configMap.end()) {

src/lib0_basic/utility/Logging.hpp

@@ -8,6 +8,7 @@
 #include <unordered_map>
 #include <ostream>
 #include <complex>
+#include <iomanip>
 
 
 class Time;
@@ -29,7 +30,7 @@ struct _Setw;
 }
 
 
-//Use beginl to write in same line. Note that in output file there will be set an new line
+//Use beginl to write in same line. Note that in output file there will be set a new line
 namespace std {
 
 template<class CharT, class Traits>
@@ -84,6 +85,25 @@ public:
   int Lenght() const { return length; }
 };
 
+template<typename T>
+class PrintIterEntry {
+  const std::string msg;
+  const T value;
+  int verbose;
+  int setw;
+public:
+  PrintIterEntry(const std::string &msg, const T &value, int setw, int verbose = 1)
+    : msg(msg), value(value), verbose(verbose), setw(setw) {}
+
+  const std::string &Message() const { return msg; }
+
+  const T &Value() const { return value; }
+
+  int Verbose() const { return verbose; }
+
+  int SetW() const { return setw; }
+};
+
 typedef PrintInfoEntry<std::string> InfoEntry;
 
 typedef std::vector<InfoEntry> InfoEntries;
@@ -107,6 +127,7 @@ protected:
   bool doubleBlock;
 private:
   template<typename T>
+
   void printEntry(int classVerbose, const PrintInfoEntry<T> &entry) {
     if (classVerbose < entry.Verbose()) return;
     printEntryRaw(entry.Message(), entry.Value());
@@ -119,19 +140,36 @@ private:
     this->endl();
   }
 
+  template<typename T>
+  void printIterEntry(int classVerbose, const PrintIterEntry<T> &entry) {
+    if (classVerbose < entry.Verbose()) return;
+    printIterEntryRaw(entry.Message(), entry.Value(), entry.SetW());
+  }
+
   template<typename T>
   void printEntries(int classVerbose, const PrintInfoEntry<T> &entry) {
     printEntry(classVerbose, entry);
   }
 
   template<typename T, typename... Types>
-  void printEntries(int classVerbose,
-                    const PrintInfoEntry<T> &entry,
+  void printEntries(int classVerbose, const PrintInfoEntry<T> &entry,
                     const Types &... types) {
     printEntry(classVerbose, entry);
     printEntries(classVerbose, types...);
   }
 
+  template<typename T>
+  void printIterEntries(int classVerbose, const PrintIterEntry<T> &entry) {
+    printIterEntry(classVerbose, entry);
+  }
+
+  template<typename T, typename... Types>
+  void printIterEntries(int classVerbose, const PrintIterEntry<T> &entry,
+                        const Types &... types) {
+    printIterEntry(classVerbose, entry);
+    printIterEntries(classVerbose, types...);
+  }
+
 public:
   //File logging cannot be enabled again since file might not be opened
   static void DisableFileLogging();
@@ -204,37 +242,59 @@ public:
   template<class T>
   void printEntryRaw(const std::string &msg0, const T &val) {
     std::string msg(msg0 + ": ");
-    int l = 25 - int(msg0.size());
+    int l = 40 - int(msg0.size());
     if (l > 0)
       msg.append(std::string(l, '.')).append(" ");
     (*this) << "  " << msg << val;
     this->endl();
   }
 
+  template<class T>
+  void printIterEntryRaw(std::string msg, const T &val, int setw) {
+    (*this) << msg << "=" << std::setw(setw) << std::left << val << " ";
+  }
+
+  template<typename T>
+  bool hasPrintableEntries(int classVerbose,
+                           const PrintInfoEntry<T> &entry) {
+    return classVerbose >= entry.Verbose();
+  }
+
+  template<typename T, typename... Types>
+  bool hasPrintableEntries(int classVerbose,
+                           const PrintInfoEntry<T> &entry,
+                           const Types &... types) {
+    if (hasPrintableEntries(classVerbose, entry)) {
+      return true;
+    }
+    if (sizeof...(types) == 0) {
+      return false;
+    }
+    return hasPrintableEntries(classVerbose, types...);
+  }
+
   template<typename T>
   bool hasPrintableEntries(int classVerbose,
-                           const PrintInfoEntry<T> &entry){
+                           const PrintIterEntry<T> &entry) {
     return classVerbose >= entry.Verbose();
   }
 
   template<typename T, typename... Types>
   bool hasPrintableEntries(int classVerbose,
-                          const PrintInfoEntry<T> &entry,
-                          const Types &... types) {
-    if(hasPrintableEntries(classVerbose, entry)){
+                           const PrintIterEntry<T> &entry,
+                           const Types &... types) {
+    if (hasPrintableEntries(classVerbose, entry)) {
       return true;
     }
-    if(sizeof...(types) == 0){
+    if (sizeof...(types) == 0) {
       return false;
     }
     return hasPrintableEntries(classVerbose, types...);
   }
 
   template<typename... Types>
-  void PrintEntries(const std::string &title,
-                    int msgVerbose,
-                    const Types &... types) {
-    if(!hasPrintableEntries(msgVerbose, types...)){
+  void PrintEntries(const std::string &title, int msgVerbose, const Types &... types) {
+    if (!hasPrintableEntries(msgVerbose, types...)) {
       return;
     }
 
@@ -245,21 +305,27 @@ public:
   }
 
   template<typename... Types>
-  void PrintInfo(const std::string &className,
-                 int classVerbose,
-                 const Types &... types) {
+  void PrintIterEntries(int msgVerbose, const Types &... types) {
+    if (!hasPrintableEntries(msgVerbose, types...)) {
+      return;
+    }
+    printIterEntries(msgVerbose, types...);
+    this->endl();
+  }
+
+  template<typename... Types>
+  void PrintInfo(const std::string &className, int classVerbose, const Types &... types) {
     PrintEntries(className + " Info:", classVerbose, types...);
   }
 
   template<typename T>
-  void PrintInfo(const std::string &className,
-                 int classVerbose,
+  void PrintInfo(const std::string &className, int classVerbose,
                  const std::vector<PrintInfoEntry<T>> &entries) {
     bool hasPrintableEntry = false;
-    for (int i = 0; i < entries.size(); ++i){
-      hasPrintableEntry |= classVerbose >= entries[i].Verbose() ;
+    for (int i = 0; i < entries.size(); ++i) {
+      hasPrintableEntry |= classVerbose >= entries[i].Verbose();
     }
-    if(!hasPrintableEntry){
+    if (!hasPrintableEntry) {
       return;
     }
     (*this) << className << " Info:";
@@ -269,6 +335,11 @@ public:
     this->endl();
   }
 
+  template<typename... Types>
+  void PrintIteration(int classVerbose, const Types &... types) {
+    PrintIterEntries(classVerbose, types...);
+  }
+
   virtual void printFromPLogging(const std::string &line_p, int proc) {}
 
 private:
@@ -280,14 +351,14 @@ private:
 
 template<typename T>
 Logging &operator<<(Logging &s, const std::vector<T> &t) {
-  for (const auto &p : t)
+  for (const auto &p: t)
     s << p << std::endl;
   return s;
 }
 
 template<typename T>
 Logging &operator<<(Logging &s, const std::list<T> &t) {
-  for (const auto &p : t)
+  for (const auto &p: t)
     s << p << std::endl;
   return s;
 }

src/lib1_math/Tensor.hpp

@@ -507,12 +507,19 @@ constexpr inline TensorT<T, dim> operator-(const TensorT<T, dim> &S) {
   return R *= T(-1.0);
 }
 
+template<typename T, int dim>
+constexpr inline VectorFieldT<T, dim> operator*(const VectorFieldT<T,dim> &v, const TensorT<T, dim> &S) {
+  TensorT<T, dim> vS = transpose(S);
+  return vS*v;
+}
+
 template<typename T, typename TT, int dim>
 constexpr inline TensorT<T, dim> operator*(const TT &b, const TensorT<T, dim> &S) {
   TensorT<T, dim> bS(S);
   return bS *= b;
 }
 
+
 template<typename T, typename TT, int dim>
 constexpr inline TensorT<T, dim> operator*(const TensorT<T, dim> &S, const TT &b) {
   TensorT<T, dim> bS(S);

src/lib2_mesh/mesh/Mesh.cpp

@@ -8,6 +8,12 @@
 #include <iostream>
 #include <algorithm>
 
+std::string std::to_string(const LevelPair &level) {
+  if (level.time < 0)
+    return "[ " + std::to_string(level.space) + " ]";
+  return "[ " + std::to_string(level.space) + " | " + std::to_string(level.time) + " ]";
+}
+
 
 using namespace std;
 
@@ -107,6 +113,9 @@ bool Mesh::onBnd(const Point &x) const {
 
 bool Mesh::onBndDG(const cell &c, int f) const {
   face ff = find_face(c.Face(f));
+#if USE_SPACETIME
+  return ((ff.Right() == Infty) || (ff.Left() == Infty));
+#endif
   return (find_bnd_face(c.Face(f)) != bnd_faces_end())
     && ((ff.Right() == Infty) || (ff.Left() == Infty));
 }
@@ -287,10 +296,6 @@ void Mesh::Fill(const shared_ptr<CoarseGeometry> &cGeo) {
     }
   }
 
-  if (_cGeo->faceIds.empty())
-    for (face f = faces(); f != faces_end(); ++f)
-      if (f.Right() == Infty)
-        InsertBoundaryFace(f(), 1);
   for (const auto &i: _cGeo->faceIds) {
     int n = i.size();
     Assert(n == i.type());
@@ -300,6 +305,17 @@ void Mesh::Fill(const shared_ptr<CoarseGeometry> &cGeo) {
     if (find_face(x) != faces_end())
       InsertBoundaryFace(x, i.flag());
   }
+
+  int bndFaceCount = 0;
+  for (face f = faces(); f != faces_end(); ++f) {
+    if (f.Right() == Infty) {
+      bndFaceCount++;
+      InsertBoundaryFace(f(), 0);
+    }
+  }
+  /*if(bndFaceCount>_cGeo->faceIds.size()){
+    Warning("Not all boundary faces on geometry defined.\n Inserting missing faces with default value \"0\".")
+  }*/
 }
 
 int find_other_proc(const procset& ps){
@@ -802,6 +818,14 @@ cell Mesh::find_cell_or_overlap_cell(const Point &z) const {
   THROW("Error: no cell found")
 }
 
+cell Mesh::find_previous_cell(const cell &c) const {
+  face f = find_face(c.Face(c.Faces() - 2));
+  cell c_prev = find_cell_or_overlap_cell(f.Left());
+  if (c_prev == overlap_end())
+    c_prev = c;
+  return c_prev;
+}
+
 cell Mesh::find_neighbour_cell(const cell &c, int f) const {
   face ff = find_face(c.Face(f));
   if (c() != ff.Right())

src/lib2_mesh/mesh/Mesh.hpp

@@ -16,13 +16,21 @@ struct LevelPair {
 
   LevelPair NextInSpace() const { return LevelPair{space + 1, time}; }
 
+  LevelPair CoarserInSpace() const { return LevelPair{space - 1, time}; }
+
   LevelPair NextInTime() const { return LevelPair{space, time + 1}; }
+
+  LevelPair CoarserInTime() const { return LevelPair{space, time - 1}; }
 };
 
 inline bool operator==(const LevelPair &l1, const LevelPair &l2) {
   return (l1.space == l2.space) && (l1.time == l2.time);
 }
 
+inline bool operator!=(const LevelPair &l1, const LevelPair &l2) {
+  return !(l1 == l2);
+}
+
 namespace std {
   template<>
   struct hash<LevelPair> {
@@ -30,8 +38,17 @@ namespace std {
       return size_t(levels.space * 10000 + levels.time);
     }
   };
+
+  std::string to_string(const LevelPair &level);
 } // namespace std
 
+inline std::ostream &operator<<(std::ostream &os, LevelPair level) {
+  return os << std::to_string(level);
+}
+
+inline Logging &operator<<(Logging &os, LevelPair level) {
+  return os << std::to_string(level);
+}
 
 class Mesh {
 protected:
@@ -199,6 +216,8 @@ public:
 
   cell find_neighbour_cell(const cell &c, int f) const;
 
+  cell find_previous_cell(const cell &c) const;
+
   cell overlap() const;
 
   cell overlap_end() const;

src/lib2_mesh/mesh/Meshes.hpp

@@ -7,7 +7,7 @@
 
 class Meshes {
 protected:
-  typedef std::unordered_map<LevelPair , std::unique_ptr<Mesh>> MeshesContainer;
+  typedef std::unordered_map<LevelPair, std::unique_ptr<Mesh>> MeshesContainer;
 
   int verbose = 1;
 
@@ -89,15 +89,39 @@ public:
   // todo remove
   int dim() const { return fine().dim(); }
 
+  LevelPair FineLevel() const {
+    if (!IsStMeshes())
+      return {fineLevel, -1};
+    else
+      return {fineLevel, fineLevel + Scales};
+  }
+
+  LevelPair CoarseLevel() const {
+    if (!IsStMeshes())
+      return {coarseLevel, -1};
+    else
+      return {coarseLevel, fineLevel + Scales};
+  }
+
+  LevelPair PLevel() const {
+    if (!IsStMeshes())
+      return {plevel, -1};
+    else
+      return {plevel, fineLevel + Scales};
+  }
+
+  // Todo remove
   int Level() const { return fineLevel; }
 
+  // Todo remove
   int pLevel() const { return plevel; }
 
+  // Todo remove
   int cLevel() const { return coarseLevel; }
 
   int CommSplit() const { return commSplit; }
 
-  const Mesh &begin() const { return coarseMesh(); }
+  const Mesh &begin() const { return pMesh(); }
 
   const Mesh &end() const { return fineMesh(); }
 
@@ -105,46 +129,28 @@ public:
 
   const Mesh &fine() const { return fineMesh(); }
 
-  const Mesh &coarseMesh() const {
-    if (!IsStMeshes())
-      return *meshes.at({plevel, -1});
-    else
-      return *meshes.at({plevel, plevel + Scales});
-  }
+  const Mesh &coarseMesh() const { return (*this)[CoarseLevel()]; }
 
-  const Mesh &fineMesh() const {
-    if (!IsStMeshes())
-      return *meshes.at({fineLevel, -1});
-    else
-      return *meshes.at({fineLevel, fineLevel + Scales});
-  }
+  const Mesh &fineMesh() const { return (*this)[FineLevel()]; }
 
-  const std::unordered_map<LevelPair , std::unique_ptr<Mesh>> &GetMeshes() const {
+  const Mesh &pMesh() const { return (*this)[PLevel()]; }
+
+  const std::unordered_map<LevelPair, std::unique_ptr<Mesh>> &GetMeshes() const {
     return meshes;
   }
 
   const Mesh &operator[](int level) const {
     if (!IsStMeshes())
-      return *meshes.at({level, -1});
+      return (*this)[{level, -1}];
     else
-      return *meshes.at({level, level + Scales});
+      return (*this)[{level, level + Scales}];
   }
 
-//  const Mesh &operator[](int level) const {
-//    try {
-//      return *meshes.at({level, level});
-//    }
-//    catch (...) {
-//      Exit("No Mesh in Meshes with given level.")
-//    }
-//  };
-
-  const Mesh &operator[](LevelPair levels) const {
+  const Mesh &operator[](LevelPair level) const {
     try {
-      return *meshes.at(levels);
-    }
-    catch (...) {
-      Exit("No Mesh in Meshes with given level key.")
+      return *meshes.at(level);
+    } catch (...) {
+      THROW("No Mesh in Meshes with given level key " + std::to_string(level))
     }
   };
 
@@ -203,7 +209,7 @@ public:
   // Todo remove and replace
   std::shared_ptr<Mesh> getPlusPLevel(int i, int j) {
     std::shared_ptr<Mesh> ref =
-      std::move(meshes.at({i + plevel, j + plevel}));
+        std::move(meshes.at({i + plevel, j + plevel}));
     return ref;
   }