MueLu: Adding preliminary el2node dump capability

Build/Test Cases Summary Enabled Packages: MueLu 0) MPI_RELEASE_DEBUG_SHARED_PT => passed: passed=51,notpassed=0 (13.88 min)
ibaned · Mar 7, 2017 · c0d0d38 · c0d0d38
1 parent fd6c8c1
commit c0d0d38
Show file tree

Hide file tree

Showing 5 changed files with 81 additions and 8 deletions.
diff --git a/packages/muelu/src/Interface/MueLu_HierarchyManager.hpp b/packages/muelu/src/Interface/MueLu_HierarchyManager.hpp
@@ -63,6 +63,14 @@
 #include "MueLu_MasterList.hpp"
 #include "MueLu_PerfUtils.hpp"
 
+#ifdef HAVE_MUELU_INTREPID2
+#ifdef HAVE_MUELU_INTREPID2_REFACTOR
+#include "Kokkos_DynRankView.hpp"
+#else
+#include "Intrepid2_FieldContainer.hpp"
+#endif
+#endif
+
 namespace MueLu {
 
   // This class stores the configuration of a Hierarchy.
@@ -235,6 +243,15 @@ namespace MueLu {
       WriteData<Matrix>(H, restrictorsToPrint_,  "R");
       WriteData<MultiVector>(H, nullspaceToPrint_,  "Nullspace");
       WriteData<MultiVector>(H, coordinatesToPrint_,  "Coordinates");
+#ifdef HAVE_MUELU_INTREPID2
+#ifdef HAVE_MUELU_INTREPID2_REFACTOR
+      typedef Kokkos::DynRankView<LocalOrdinal,typename Node::device_type> FCi;
+#else
+      typedef Intrepid2::FieldContainer<LocalOrdinal> FCi;
+#endif
+      WriteDataFC<FCi>(H,elementToNodeMapsToPrint_, "pcoarsen: element to node map","el2node");
+#endif    
+
 
     } //SetupHierarchy
 
@@ -283,13 +300,14 @@ namespace MueLu {
     Teuchos::Array<int>   restrictorsToPrint_;
     Teuchos::Array<int>   nullspaceToPrint_;
     Teuchos::Array<int>   coordinatesToPrint_;
+    Teuchos::Array<int>   elementToNodeMapsToPrint_;
 
     std::map<int, std::vector<keep_pair> > keep_;
 
   private:
 
     template<class T>
-    void WriteData(Hierarchy& H, const Teuchos::Array<int>& data, const std::string& name) const {
+    void WriteData(Hierarchy& H, const Teuchos::Array<int>& data, const std::string& name, bool isFieldContainer=false) const {
       for (int i = 0; i < data.size(); ++i) {
         std::string fileName = name + "_" + Teuchos::toString(data[i]) + ".m";
 
@@ -298,13 +316,63 @@ namespace MueLu {
 
           if (L->IsAvailable(name)) {
             RCP<T> M = L->template Get< RCP<T> >(name);
-            if (!M.is_null())
-              Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(fileName,* M);
-          }
-        }
+            if (!M.is_null()) {
+	      Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(fileName,* M);              		
+	    }
+	  }
+	}
       }
     }
 
+
+    template<class T>
+    void WriteDataFC(Hierarchy& H, const Teuchos::Array<int>& data, const std::string& name, const std::string & ofname, bool isFieldContainer=false) const {
+      for (int i = 0; i < data.size(); ++i) {
+        const std::string fileName = ofname + "_" + Teuchos::toString(data[i]) + ".m";
+
+        if (data[i] < H.GetNumLevels()) {
+          RCP<Level> L = H.GetLevel(data[i]);
+
+          if (L->IsAvailable(name)) {
+            RCP<T> M = L->template Get< RCP<T> >(name);
+            if (!M.is_null()) {
+	      RCP<Matrix> A = L->template Get<RCP<Matrix> >("A");
+	      RCP<const CrsGraph> AG = A->getCrsGraph();
+	      WriteFieldContainer<T>(fileName,*M,*AG->getColMap());
+	    }
+	  }
+	}
+      }
+    }
+
+    // For dumping an IntrepidPCoarsening element-to-node map to disk     
+    template<class T>
+    void WriteFieldContainer(const std::string& fileName, T & fcont,const Map &colMap) const {
+      typedef LocalOrdinal LO;
+      typedef GlobalOrdinal GO;
+      typedef Node NO;
+      typedef Xpetra::MultiVector<GO,LO,GO,NO> GOMultiVector;  
+
+      size_t num_els = (size_t) fcont.dimension(0);
+      size_t num_vecs =(size_t) fcont.dimension(1);
+
+      // Generate rowMap
+      Teuchos::RCP<const Map> rowMap = Xpetra::MapFactory<LO,GO,NO>::Build(colMap.lib(),Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),fcont.dimension(0),colMap.getIndexBase(),colMap.getComm());
+
+      // Fill multivector to use *petra dump routines
+      RCP<GOMultiVector> vec = Xpetra::MultiVectorFactory<GO, LO, GO, NO>::Build(rowMap,num_vecs);
+
+      for(size_t j=0; j<num_vecs; j++)  {
+	Teuchos::ArrayRCP<GO> v = vec->getDataNonConst(j);
+	for(size_t i=0; i<num_els; i++) 
+	  v[i] = colMap.getGlobalElement(fcont(i,j));
+      }
+
+      Xpetra::IO<GO,LO,GO,NO>::Write(fileName,*vec);
+    }
+
+
+
     // Levels
     Array<RCP<FactoryManagerBase> > levelManagers_;        // one FactoryManager per level (the last levelManager is used for all the remaining levels)
 

diff --git a/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp b/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp
@@ -276,6 +276,8 @@ namespace MueLu {
         this->nullspaceToPrint_  = Teuchos::getArrayFromStringParameter<int>(printList, "Nullspace");
       if (printList.isParameter("Coordinates"))
         this->coordinatesToPrint_  = Teuchos::getArrayFromStringParameter<int>(printList, "Coordinates");
+      if (printList.isParameter("pcoarsen: element to node map"))
+        this->elementToNodeMapsToPrint_  = Teuchos::getArrayFromStringParameter<int>(printList, "pcoarsen: element to node map");
     }
 
     // Set verbosity parameter

diff --git a/packages/muelu/src/Transfers/PCoarsen/MueLu_IntrepidPCoarsenFactory_decl.hpp b/packages/muelu/src/Transfers/PCoarsen/MueLu_IntrepidPCoarsenFactory_decl.hpp
@@ -258,11 +258,10 @@ namespace MueLu {
     // ! will be stored in seeds.
     template<class Basis, class LOFieldContainer, class LocalOrdinal, class GlobalOrdinal, class Node>
     void FindGeometricSeedOrdinals(Teuchos::RCP<Basis> basis, const LOFieldContainer &elementToNodeMap,
-                                   std::vector<std::vector<LocalOrdinal>> &seeds,
+                                   std::vector<std::vector<LocalOrdinal> > &seeds,
                                    const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &rowMap,
                                    const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &columnMap);
 
-
   }//namespace MueLuIntrepid
 } //namespace MueLu
 

diff --git a/packages/muelu/src/Transfers/PCoarsen/MueLu_IntrepidPCoarsenFactory_def.hpp b/packages/muelu/src/Transfers/PCoarsen/MueLu_IntrepidPCoarsenFactory_def.hpp
@@ -47,6 +47,7 @@
 #define MUELU_IPCFACTORY_DEF_HPP
 
 #include <Xpetra_Matrix.hpp>
+#include <Xpetra_IO.hpp>
 #include <sstream>
 #include <algorithm>
 
@@ -107,9 +108,11 @@ inline std::string tolower(const std::string & str) {
   return data;
 }
 
+
+/*********************************************************************************************************/
   template<class Basis, class LOFieldContainer, class LocalOrdinal, class GlobalOrdinal, class Node>
   void FindGeometricSeedOrdinals(Teuchos::RCP<Basis> basis, const LOFieldContainer &elementToNodeMap,
-                                 std::vector<std::vector<LocalOrdinal>> &seeds,
+                                 std::vector<std::vector<LocalOrdinal> > &seeds,
                                  const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &rowMap,
                                  const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &columnMap)
   {

diff --git a/packages/muelu/test/unit_tests/IntrepidPCoarsenFactory.cpp b/packages/muelu/test/unit_tests/IntrepidPCoarsenFactory.cpp
@@ -3139,6 +3139,7 @@ TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(IntrepidPCoarsenFactory,BuildP_PseudoPoisson_L
     dump.set("A","{0,1,2}");
     dump.set("P","{0,1}");
     dump.set("R","{0,1}");
+    dump.set("pcoarsen: element to node map","{0,1,2}");
     Params.set("export data",dump);
 #endif