Closure...

FunctionSpace/FunctionSpace.cpp

@@ -8,6 +8,15 @@ FunctionSpace::FunctionSpace(void){
 }
 
 FunctionSpace::~FunctionSpace(void){
+  map<const MElement*, vector<bool>*>::iterator it 
+    = edgeClosure->begin();
+  map<const MElement*, vector<bool>*>::iterator stop 
+    = edgeClosure->end();
+  
+  for(; it != stop; it++)
+    delete it->second;
+  delete edgeClosure;
+
   delete basis;
 }
 
@@ -48,6 +57,66 @@ void FunctionSpace::build(const GroupOfElement& goe,
     fPerFace = 0;  
   
   fPerCell = basis->getNCellBased(); // We always got 1 cell 
+
+  // Build Closure
+  edgeClosure = new map<const MElement*, vector<bool>*>;
+  closure();
+}
+
+void FunctionSpace::closure(void){
+  // Get Elements //
+  const vector<const MElement*>& element = goe->getAll();
+  const unsigned int             size    = element.size();    
+
+  // Iterate on elements //
+  for(unsigned int i = 0; i < size; i++){
+    // Get Element data
+    MElement& myElement =
+      const_cast<MElement&>(*element[i]);
+
+    const unsigned int nVertex = myElement.getNumVertices();
+    const unsigned int nEdge   = myElement.getNumEdges();
+    const unsigned int nFace   = myElement.getNumFaces();
+
+    const unsigned int nTotVertex = nVertex * fPerVertex;
+    const unsigned int nTotEdge   = nEdge   * fPerEdge;
+    const unsigned int nTotFace   = nFace   * fPerFace;
+    
+    const unsigned int nTot    = nTotVertex + nTotEdge + nTotFace + fPerCell;
+
+    // Closure
+    vector<bool>* closure = new vector<bool>(nTot);
+    unsigned int it = 0;
+
+    // Closure for vertices
+    for(unsigned int j = 0; j < nTotVertex; j++, it++)
+      (*closure)[it] = true;
+
+    // Closure for edges 
+    for(unsigned int j = 0; j < fPerEdge; j++){
+      for(unsigned int k = 0; k < nEdge; k++, it++){
+	// Orientation 
+	int orientation = mesh->getOrientation(myElement.getEdge(k));
+	
+	if(orientation == 1)
+	  (*closure)[it] = true;
+	
+	else
+	  (*closure)[it] = false;
+      }
+    }
+
+    // Closure for faces
+    // TODO
+
+    // Closure for cells
+    for(unsigned int j = 0; j < fPerCell; j++, it++)
+      (*closure)[it] = true;    
+
+    // Add Closure
+    edgeClosure->insert
+      (pair<const MElement*, vector<bool>*>(element[i], closure));
+  }
 }
 
 vector<Dof> FunctionSpace::getKeys(const MElement& elem) const{ 

FunctionSpace/FunctionSpace.h

 #ifndef _FUNCTIONSPACE_H_
 #define _FUNCTIONSPACE_H_
 
+#include <map>
 #include <vector>
 
 #include "Basis.h"
@@ -40,6 +41,8 @@ class FunctionSpace{
   const GroupOfElement* goe;
   const Basis*          basis;
 
+  std::map<const MElement*, std::vector<bool>*>* edgeClosure;
+
   unsigned int fPerVertex;
   unsigned int fPerEdge;
   unsigned int fPerFace;
@@ -71,6 +74,8 @@ class FunctionSpace{
 
   void build(const GroupOfElement& goe,
 	     int basisType, int order);
+
+  void closure(void);
 };
 
 

FunctionSpace/FunctionSpaceNode.cpp

@@ -25,8 +25,8 @@ interpolate(const MElement& element,
     const_cast<MElement&>(element);
   
   // Get Basis Functions //
-  const vector<const Polynomial*>& fun = getBasis(element).getFunctions();
-  const unsigned int nFun              = fun.size();
+  const vector<const Polynomial*> fun = getLocalFunctions(element);
+  const unsigned int nFun             = fun.size();
   
   // Get Reference coordinate //
   double phys[3] = {xyz(0), xyz(1), xyz(2)};

FunctionSpace/FunctionSpaceScalar.cpp

 #include "FunctionSpaceScalar.h"
+#include "Exception.h"
+
+using namespace std;
 
 FunctionSpaceScalar::~FunctionSpaceScalar(void){
 }
+
+const vector<const Polynomial*> FunctionSpaceScalar::
+getLocalFunctions(const MElement& element) const{
+
+  // Get Closure //
+  map<const MElement*, vector<bool>*>::iterator it = 
+    edgeClosure->find(&element);
+
+  if(it == edgeClosure->end())
+    throw Exception("Element not found for closure");
+
+  vector<bool>* closure   = it->second;
+  const unsigned int size = closure->size();
+
+  // Get Basis //
+  const vector<const Polynomial*>&    basis = getBasis(element).getFunctions(0);
+  const vector<const Polynomial*>& revBasis = getBasis(element).getFunctions(1);
+
+  // Get Functions //
+  vector<const Polynomial*> fun(size);
+  
+  for(unsigned int i = 0; i < size; i++)
+    if((*closure)[i])
+      fun[i] = basis[i];
+
+    else
+      fun[i] = revBasis[i];
+
+  // Return 
+  return fun;
+}

FunctionSpace/FunctionSpaceScalar.h

@@ -28,7 +28,10 @@ class FunctionSpaceScalar : public FunctionSpace{
     interpolate(const MElement& element, 
 		const std::vector<double>& coef,
 		const fullVector<double>& xyz) const = 0;
-
+  
+  const std::vector<const Polynomial*> 
+    getLocalFunctions(const MElement& element) const;
+  
   const BasisScalar& getBasis(const MElement& element) const;
 };