From 5dfd1e2b1ea7d49a74cc0a0261686c115a1bd065 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Mon, 28 Jul 2014 13:50:36 +0000
Subject: [PATCH] Precomputing of FunctionSpace Dof + Unique ID for
 GroupOfElement + Correct Condition in Haroche (Be compy-past cursedsvn st
 --quiet)

---
 FunctionSpace/FunctionSpace.cpp       | 130 ++++++++++++++------------
 FunctionSpace/FunctionSpace.h         |  34 ++-----
 FunctionSpace/FunctionSpaceScalar.cpp |  61 ++++++++----
 FunctionSpace/FunctionSpaceScalar.h   |  31 ++++++
 FunctionSpace/FunctionSpaceVector.cpp |  57 ++++++++---
 FunctionSpace/FunctionSpaceVector.h   |  31 ++++++
 6 files changed, 232 insertions(+), 112 deletions(-)

diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index 6ffdb4289b..780b654338 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -9,10 +9,18 @@
 
 using namespace std;
 
-const size_t FunctionSpace::nGeoType = 9;
-size_t FunctionSpace::nxtOffset = 0;
+const size_t FunctionSpace::nGeoType  = 9;
+      size_t FunctionSpace::nxtOffset = 0;
 
 FunctionSpace::FunctionSpace(void){
+  // Alloc Basis Vector for all possible geomtrical types //
+  basis.resize(nGeoType, NULL);
+
+  // Alloc Function per Entity //
+  fPerVertex.resize(nGeoType, 0);
+  fPerEdge.resize(nGeoType  , 0);
+  fPerFace.resize(nGeoType  , 0);
+  fPerCell.resize(nGeoType  , 0);
 }
 
 FunctionSpace::~FunctionSpace(void){
@@ -25,32 +33,20 @@ void FunctionSpace::build(const GroupOfElement& goe, string family){
   // Save Dof type offset //
   offset = nxtOffset;
 
-  // Save GroupOfElement & Mesh //
-  this->goe  = &goe;
+  // Save Mesh //
   this->mesh = &(goe.getMesh());
 
-  // Alloc Basis Vector for all possible geomtrical types //
-  basis.resize(nGeoType, NULL);
-
   // Generate Bases //
-
-  // Get geomtrical type statistics
   const vector<size_t>& geoTypeStat = goe.getTypeStats();
-  const size_t nGeoType = geoTypeStat.size();
+  const size_t             nGeoType = geoTypeStat.size();
 
-  // Buils basis for existing geomtrical type
   for(size_t i = 0; i < nGeoType; i++)
-    if(geoTypeStat[i])
+    if(geoTypeStat[i] != 0 && basis[i] == NULL)
       basis[i] = BasisGenerator::generate(i, form, order, family);
 
   // Get Number of Function per Entity //
-  fPerVertex.resize(nGeoType);
-  fPerEdge.resize(nGeoType);
-  fPerFace.resize(nGeoType);
-  fPerCell.resize(nGeoType);
-
   for(size_t i = 0; i < nGeoType; i++){
-    if(geoTypeStat[i]){
+    if(geoTypeStat[i] != 0 && fPerVertex[i] == 0){
       int nVertex = ReferenceSpaceManager::getNVertex(i);
       int nEdge   = ReferenceSpaceManager::getNEdge(i);
       int nFace   = ReferenceSpaceManager::getNFace(i);
@@ -69,55 +65,70 @@ void FunctionSpace::build(const GroupOfElement& goe, string family){
 
       fPerCell[i] = this->basis[i]->getNCellBased();
     }
-
-    else{
-      fPerVertex[i] = 0;
-      fPerEdge[i]   = 0;
-      fPerFace[i]   = 0;
-      fPerCell[i]   = 0;
-    }
   }
 
   // Build Dof //
-   buildDof();
-
-  // Find next offset //
-  nxtOffset = findMaxType() + 1;
+  buildDof(goe);
 }
 
-void FunctionSpace::buildDof(void){
+void FunctionSpace::buildDof(const GroupOfElement& goe){
   // Get Elements //
-  const size_t nElement = goe->getNumber();
-  const vector<const MElement*>& element = goe->getAll();
+  const size_t                   nElement = goe.getNumber();
+  const vector<const MElement*>&  element = goe.getAll();
 
-  vector<Dof> myDof;
-  size_t nDof;
+  // Push GroupOfElement into map //
+  pair<size_t, vector<vector<Dof> > >                      toInsert;
+  pair<map<size_t, vector<vector<Dof> > >::iterator, bool> isInserted;
+
+  toInsert.first  = goe.getId();
+  toInsert.second = vector<vector<Dof> >(0);
+  isInserted      = dof.insert(toInsert);
+
+  if(!isInserted.second)
+    throw Exception("FunctionSpace: cannot computed Dofs for GroupOfElement %d",
+                    goe.getId());
+
+  // Reference & Allocate //
+  vector<vector<Dof> >& myDof = isInserted.first->second;
+  myDof.resize(nElement);
 
   // Create Dofs //
-  for(size_t i = 0; i < nElement; i++){
-    // Get Dof for this Element
-    getKeys(*(element[i]), myDof);
-    nDof = myDof.size();
-
-    // Add Dofs
-    for(size_t j = 0; j < nDof; j++)
-      dof.insert(myDof[j]);
-  }
+  for(size_t i = 0; i < nElement; i++)
+    getKeys(*(element[i]), myDof[i]);
 }
 
 size_t FunctionSpace::findMaxType(void){
   // Maximum type //
   size_t maxType = 0;
 
-  // Iterate for dof //
-  const set<Dof>::iterator end = dof.end();
-        set<Dof>::iterator it  = dof.begin();
+  // Iterate on GroupOfElement Id //
+  map<size_t, vector<vector<Dof> > >::iterator it  = dof.begin();
+  map<size_t, vector<vector<Dof> > >::iterator end = dof.end();
+
+  size_t nElement;
+  size_t nDof;
+  size_t type;
+
+  for(; it != end; it++){
+    // Iterate on Elements of this GroupOfElement //
+    nElement = it->second.size();
 
-  for(; it != end; it++)
-    // If this type is bigger, it becomes the new 'maxType'
-    if(it->getType() > maxType)
-      maxType = it->getType();
+    for(size_t e = 0; e < nElement; e++){
+      // Iterate on Dofs of this Element //
+      nDof = it->second[e].size();
 
+      for(size_t d = 0; d < nDof; d++){
+        // This Dof Type
+        type = it->second[e][d].getType();
+
+        // If this Dof type is bigger, it becomes the new 'maxType'
+        if(type > maxType)
+          maxType = type;
+      }
+    }
+  }
+
+  // Return maxType //
   return maxType;
 }
 
@@ -243,16 +254,15 @@ void FunctionSpace::getKeys(const GroupOfElement& goe,
   }
 }
 
-void FunctionSpace::getKeys(const GroupOfElement& goe,
-                            std::vector<std::vector<Dof> >& dof) const{
-  // Get Elements //
-  const size_t nElement = goe.getNumber();
-  const vector<const MElement*>& element = goe.getAll();
+const std::vector<std::vector<Dof> >&
+FunctionSpace::getKeys(const GroupOfElement& goe) const{
+  // Find vector of Dof from map //
+  map<size_t, vector<vector<Dof> > >::const_iterator it = dof.find(goe.getId());
 
-  // Init Struct //
-  dof.resize(nElement);
+  if(it == dof.end())
+    throw Exception("FunctionSpace: cannot find Dofs of GroupOfElement %d",
+                    goe.getId());
 
-  // Create Dofs //
-  for(size_t i = 0; i < nElement; i++)
-    getKeys(*(element[i]), dof[i]);
+  // Return vector //
+  return it->second;
 }
diff --git a/FunctionSpace/FunctionSpace.h b/FunctionSpace/FunctionSpace.h
index 0bb99c8f4d..49248be93d 100644
--- a/FunctionSpace/FunctionSpace.h
+++ b/FunctionSpace/FunctionSpace.h
@@ -32,15 +32,14 @@ class FunctionSpace{
  protected:
   // Number of possible geomtrical topologies & Dof Type offset //
   static const size_t nGeoType;
-  static size_t nxtOffset;
+  static       size_t nxtOffset;
 
  protected:
   // Offset //
   size_t offset;
 
-  // Geometry //
-  const Mesh*     mesh;
-  const GroupOfElement* goe;
+  // Mesh //
+  const Mesh* mesh;
 
   // Basis //
   std::vector<const Basis*> basis;
@@ -56,7 +55,7 @@ class FunctionSpace{
   size_t order;
 
   // Dofs //
-  std::set<Dof> dof;
+  std::map<size_t, std::vector<std::vector<Dof> > > dof;
 
  public:
   virtual ~FunctionSpace(void);
@@ -65,20 +64,17 @@ class FunctionSpace{
   size_t getForm(void)  const;
   size_t getOrder(void) const;
 
-  const Basis&          getBasis(const MElement& element) const;
-  const Basis&          getBasis(size_t eType)            const;
-  const GroupOfElement& getSupport(void)                  const;
+  const Basis& getBasis(const MElement& element) const;
+  const Basis& getBasis(size_t eType)            const;
 
   void getKeys(const MElement& element, std::vector<Dof>& dof) const;
   void getKeys(const GroupOfElement& goe, std::set<Dof>& dof)  const;
-  void getKeys(const GroupOfElement& goe,
-               std::vector<std::vector<Dof> >& dof)            const;
+  const std::vector<std::vector<Dof> >& getKeys(const GroupOfElement& goe)const;
 
  protected:
   FunctionSpace(void);
-
   void   build(const GroupOfElement& goe, std::string family);
-  void   buildDof(void);
+  void   buildDof(const GroupOfElement& goe);
   size_t findMaxType(void);
 
   void getUnorderedKeys(const MElement& element, std::vector<Dof>& dof) const;
@@ -120,10 +116,6 @@ class FunctionSpace{
    @return Returns the Basis associated to the given geomtrical element type tag
    **
 
-   @fn FunctionSpace::getSupport
-   @return Returns the support of this FunctionSpace
-   **
-
    @fn void FunctionSpace::getKeys(const MElement&,std::vector<Dof>&) const
    @param element A MElement
    @param dof A vector of Dof%s
@@ -139,11 +131,9 @@ class FunctionSpace{
    of the given GroupOfElement
    **
 
-   @fn void FunctionSpace::getKeys(const GroupOfElement&, std::vector<std::vector<Dof> >&) const
+   @fn const std::vector<std::vector<Dof> >& FunctionSpace::getKeys(const GroupOfElement&) const
    @param goe A GroupOfElement
-   @param dof A vector of vector of Dof%s
-
-   Populates the given vector such that:
+   @return Returns a vector of vector of Dof such that:
    dof[i][j] is the jth Dof of the ith element of the given GroupOfElement
 */
 
@@ -172,8 +162,4 @@ inline const Basis& FunctionSpace::getBasis(size_t eType) const{
   return *basis[eType];
 }
 
-inline const GroupOfElement& FunctionSpace::getSupport(void) const{
-  return *goe;
-}
-
 #endif
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index 4472dc188b..635af0a8e8 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -2,35 +2,64 @@
 #include "Exception.h"
 #include "FunctionSpaceScalar.h"
 
-FunctionSpaceScalar::FunctionSpaceScalar(const GroupOfElement& goe,
-                                         size_t order){
-  if(order == 0)
-    throw Exception("%s: %s",
-                    "FunctionSpaceScalar",
-                    "Cannot have a order 0 scalar function space");
+FunctionSpaceScalar::
+FunctionSpaceScalar(const GroupOfElement& goe, size_t order){
+  // Temp vector
+  std::vector<const GroupOfElement*> tmp(1);
+  tmp[0] = &goe;
+
+  // Init
+  init(tmp, order, "hierarchical");
+}
 
-  this->scalar = true;
-  this->form   = 0;
-  this->order  = order;
+FunctionSpaceScalar::
+FunctionSpaceScalar(const std::vector<const GroupOfElement*>& goe,
+                    size_t order){
+  // Init
+  init(goe, order, "hierarchical");
+}
+
+FunctionSpaceScalar::
+FunctionSpaceScalar(const GroupOfElement& goe,
+                    size_t order, std::string family){
+  // Temp vector
+  std::vector<const GroupOfElement*> tmp(1);
+  tmp[0] = &goe;
 
-  build(goe, "hierarchical");
+  // Init
+  init(tmp, order, family);
 }
 
-FunctionSpaceScalar::FunctionSpaceScalar(const GroupOfElement& goe,
-                                         size_t order, std::string family){
+FunctionSpaceScalar::
+FunctionSpaceScalar(const std::vector<const GroupOfElement*>& goe,
+                    size_t order, std::string family){
+  // Init
+  init(goe, order, family);
+}
+
+FunctionSpaceScalar::~FunctionSpaceScalar(void){
+  // Done by FunctionSpace
+}
+
+void FunctionSpaceScalar::init(const std::vector<const GroupOfElement*>& goe,
+                               size_t order, std::string family){
+  // Check
   if(order == 0)
     throw Exception("%s: %s",
                     "FunctionSpaceScalar",
                     "Cannot have a order 0 scalar function space");
+  // Init
   this->scalar = true;
   this->form   = 0;
   this->order  = order;
 
-  build(goe, family);
-}
+  // Build FunctionSpace
+  const size_t nGoe = goe.size();
+  for(size_t i = 0; i < nGoe; i++)
+    build(*goe[i], family);
 
-FunctionSpaceScalar::~FunctionSpaceScalar(void){
-  // Done by FunctionSpace
+  // Next Offset for next FunctionSpace
+  nxtOffset = findMaxType() + 1;
 }
 
 double FunctionSpaceScalar::interpolateInABC(const MElement& element,
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 0cae40d3ba..d0a9380241 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -17,8 +17,13 @@
 class FunctionSpaceScalar : public FunctionSpace{
  public:
   FunctionSpaceScalar(const GroupOfElement& goe, size_t order);
+  FunctionSpaceScalar(const std::vector<const GroupOfElement*>& goe,
+                      size_t order);
+
   FunctionSpaceScalar(const GroupOfElement& goe, size_t order,
                       std::string family);
+  FunctionSpaceScalar(const std::vector<const GroupOfElement*>& goe,
+                      size_t order, std::string family);
 
   virtual ~FunctionSpaceScalar(void);
 
@@ -38,6 +43,9 @@ class FunctionSpaceScalar : public FunctionSpace{
                           const fullVector<double>& xyz) const;
 
  private:
+  void init(const std::vector<const GroupOfElement*>& goe,
+            size_t order, std::string family);
+
   double interpolateInABC(const MElement& element,
                           const std::vector<double>& coef,
                           double abc[3]) const;
@@ -58,6 +66,15 @@ class FunctionSpaceScalar : public FunctionSpace{
    The instanciated FunctionSpace will use a hierarchical Basis
    **
 
+   @fn FunctionSpaceScalar::FunctionSpaceScalar(const std::vector<const GroupOfElement*>&,size_t)
+   @param goe A vector of GroupOfElement
+   @param order A natural number
+   Instanciates a new FunctionSpaceScalar
+   on the given GroupOfElement%s and with the given order
+
+   The instanciated FunctionSpace will use a hierarchical Basis
+   **
+
    @fn FunctionSpaceScalar::FunctionSpaceScalar(const GroupOfElement&,size_t,std::string)
    @param goe A GroupOfElement
    @param order A natural number
@@ -72,6 +89,20 @@ class FunctionSpaceScalar : public FunctionSpace{
    @see See BasisGenerator::generate()
    **
 
+   @fn FunctionSpaceScalar::FunctionSpaceScalar(const std::vector<const GroupOfElement*>&,size_t,std::string)
+   @param goe A vector of GroupOfElement
+   @param order A natural number
+   @param family A stringr
+   Instanciates a new FunctionSpaceScalar
+   on the given GroupOfElement%s and with the given order
+
+   The instanciated FunctionSpace will use the requested Basis family:
+   @li If family is equal to 'lagrange' a Lagrange Basis will be used
+   @li If family is equal to 'hierarchical' a hierarchical Basis will be used
+
+   @see See BasisGenerator::generate()
+   **
+
    @fn FunctionSpaceScalar::~FunctionSpaceScalar
    Deletes this FunctionSpaceScalar
    **
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index 3303886db1..5a9385ed12 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -1,28 +1,61 @@
 #include "Mapper.h"
 #include "FunctionSpaceVector.h"
 
-FunctionSpaceVector::FunctionSpaceVector(const GroupOfElement& goe,
-                                         size_t order){
-  this->scalar = false;
-  this->form   = 1;
-  this->order  = order;
+FunctionSpaceVector::
+FunctionSpaceVector(const GroupOfElement& goe, size_t order){
+  // Temp vector
+  std::vector<const GroupOfElement*> tmp(1);
+  tmp[0] = &goe;
+
+  // Init
+  init(tmp, order, "hierarchical");
+}
 
-  build(goe, "hierarchical");
+FunctionSpaceVector::
+FunctionSpaceVector(const std::vector<const GroupOfElement*>& goe,
+                    size_t order){
+  // Init
+  init(goe, order, "hierarchical");
 }
 
-FunctionSpaceVector::FunctionSpaceVector(const GroupOfElement& goe,
-                                         size_t order, std::string family){
-  this->scalar = false;
-  this->form   = 1;
-  this->order  = order;
+FunctionSpaceVector::
+FunctionSpaceVector(const GroupOfElement& goe,
+                    size_t order, std::string family){
+  // Temp vector
+  std::vector<const GroupOfElement*> tmp(1);
+  tmp[0] = &goe;
+
+  // Init
+  init(tmp, order, family);
+}
 
-  build(goe, family);
+FunctionSpaceVector::
+FunctionSpaceVector(const std::vector<const GroupOfElement*>& goe,
+                    size_t order, std::string family){
+  // Init
+  init(goe, order, family);
 }
 
 FunctionSpaceVector::~FunctionSpaceVector(void){
   // Done by FunctionSpace
 }
 
+void FunctionSpaceVector::init(const std::vector<const GroupOfElement*>& goe,
+                               size_t order, std::string family){
+  // Init
+  this->scalar = false;
+  this->form   = 1;
+  this->order  = order;
+
+  // Build FunctionSpace
+  const size_t nGoe = goe.size();
+  for(size_t i = 0; i < nGoe; i++)
+    build(*goe[i], family);
+
+  // Next Offset for next FunctionSpace
+  nxtOffset = findMaxType() + 1;
+}
+
 fullVector<double> FunctionSpaceVector::
 interpolateInABC(const MElement& element,
                  const std::vector<double>& coef,
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index 5f2a027760..33f07bb34e 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -18,8 +18,13 @@
 class FunctionSpaceVector : public FunctionSpace{
  public:
   FunctionSpaceVector(const GroupOfElement& goe, size_t order);
+  FunctionSpaceVector(const std::vector<const GroupOfElement*>& goe,
+                      size_t order);
+
   FunctionSpaceVector(const GroupOfElement& goe, size_t order,
                       std::string family);
+  FunctionSpaceVector(const std::vector<const GroupOfElement*>& goe,
+                      size_t order, std::string family);
 
   virtual ~FunctionSpaceVector(void);
 
@@ -43,6 +48,9 @@ class FunctionSpaceVector : public FunctionSpace{
                                     const std::vector<double>& coef,
                                     const fullVector<double>& uvw) const;
  private:
+  void init(const std::vector<const GroupOfElement*>& goe,
+            size_t order, std::string family);
+
   fullVector<double>
     interpolateInABC(const MElement& element,
                      const std::vector<double>& coef,
@@ -65,6 +73,15 @@ class FunctionSpaceVector : public FunctionSpace{
    The instanciated FunctionSpace will use a hierarchical Basis
    **
 
+   @fn FunctionSpaceVector::FunctionSpaceVector(const std::vector<const GroupOfElement*>&,size_t)
+   @param goe A vector of GroupOfElement
+   @param order A natural number
+   Instanciates a new FunctionSpaceVector
+   on the given GroupOfElement%s and with the given order
+
+   The instanciated FunctionSpace will use a hierarchical Basis
+   **
+
    @fn FunctionSpaceVector::FunctionSpaceVector(const GroupOfElement&,size_t,std::string)
    @param goe A GroupOfElement
    @param order A natural number
@@ -79,6 +96,20 @@ class FunctionSpaceVector : public FunctionSpace{
    @see See BasisGenerator::generate()
    **
 
+   @fn FunctionSpaceVector::FunctionSpaceVector(const std::vector<const GroupOfElement*>&,size_t,std::string)
+   @param goe A vector of GroupOfElement
+   @param order A natural number
+   @param family A stringr
+   Instanciates a new FunctionSpaceVector
+   on the given GroupOfElement%s and with the given order
+
+   The instanciated FunctionSpace will use the requested Basis family:
+   @li If family is equal to 'lagrange' a Lagrange Basis will be used
+   @li If family is equal to 'hierarchical' a hierarchical Basis will be used
+
+   @see See BasisGenerator::generate()
+   **
+
    @fn FunctionSpaceVector::~FunctionSpaceVector
    Deletes this FunctionSpaceVector
    **
-- 
GitLab