Bugfix & cleaning.

Geo/Cell.cpp

@@ -46,14 +46,6 @@ Cell::~Cell()
   if(_deleteImage) delete _image; 
 }
 
-bool Cell::hasVertex(int vertex) const 
-{
-  std::vector<int>::const_iterator it = std::find(_vs.begin(), _vs.end(), 
-						  vertex);
-  if (it != _vs.end()) return true;
-  else return false;
-}
-
 int Cell::getNumFacets() const 
 { 
   if(getDim() == 0) return 0;
@@ -72,8 +64,12 @@ void Cell::getFacetVertices(const int num, std::vector<MVertex*> &v) const
   return;
 }
 
-int Cell::getFacetOri(std::vector<MVertex*> &v) 
+int Cell::getFacetOri(Cell* cell) 
 {
+  std::vector<MVertex*> v; 
+  for(int i = 0; i < cell->getNumVertices(); i++) {
+    v.push_back(cell->getVertex(i));
+  }
   if(getDim() == 0) return 0;
   else if(getDim() == 1){
     if(v.size() != 1) return 0;
@@ -101,6 +97,14 @@ int Cell::getFacetOri(std::vector<MVertex*> &v)
   else return 0;
 }  
 
+bool Cell::hasVertex(int vertex) const 
+{
+  std::vector<int>::const_iterator it = std::find(_vs.begin(), _vs.end(), 
+						  vertex);
+  if (it != _vs.end()) return true;
+  else return false;
+}
+
 void Cell::printCell() 
 {
   printf("%d-cell: \n" , getDim());

Geo/Cell.h

@@ -92,14 +92,7 @@ class Cell
   // get the vertices on a facet of this cell
   virtual void getFacetVertices(const int num, std::vector<MVertex*> &v) const;
   // get boundary cell orientation
-  virtual int getFacetOri(std::vector<MVertex*> &v); 
-  virtual int getFacetOri(Cell* cell) { 
-    std::vector<MVertex*> v; 
-    for(int i = 0; i < cell->getNumVertices(); i++) {
-      v.push_back(cell->getVertex(i));
-    }
-    return getFacetOri(v);
-  }
+  virtual int getFacetOri(Cell* cell);
 
   virtual int getDim() const { return _dim; };
   virtual int getIndex() const { return _index; };

Geo/CellComplex.cpp

@@ -42,42 +42,25 @@ void CellComplex::panic_exit(){
 bool CellComplex::insert_cells(std::vector<MElement*>& elements,
 			       bool subdomain)
 {
-  std::vector<MVertex*> vertices;
-  std::pair<citer, bool> insertInfo;
-
   // add highest dimensional cells
   for(unsigned int i=0; i < elements.size(); i++){
-    vertices.clear();
-    MElement* element = elements.at(i);
-    
-    for(int k=0; k < element->getNumVertices(); k++){
-      MVertex* vertex = element->getVertex(k);
-      vertices.push_back(vertex);
-    }
-      
+    MElement* element = elements.at(i);    
     int dim = element->getDim();
-    int type = element->getTypeForMSH();
+    int type = element->getType();
     
-    Cell* cell = new Cell(element);
     // simplex types
-    if( !(type == MSH_PNT
-	  || type == MSH_LIN_2 || type == MSH_TRI_3 || type == MSH_TET_4
-	  || type == MSH_LIN_3 || type == MSH_TRI_6 || type == MSH_TET_10 
-	  || type == MSH_PNT || type == MSH_TRI_9 || type == MSH_TRI_10 
-	  || type == MSH_TRI_12 || type == MSH_TRI_15 || type == MSH_TRI_15I 
-	  || type == MSH_TRI_21 || type == MSH_LIN_4 || type == MSH_LIN_5 
-	  || type == MSH_LIN_6 || type == MSH_TET_20 || type == MSH_TET_35 
-	  || type == MSH_TET_56 || type == MSH_TET_34 || type == MSH_TET_52) ){
+    if( !(type == TYPE_PNT || type == TYPE_LIN || type == TYPE_TRI
+	  || type == TYPE_TET )){
       _simplicial = false;
     }
-    /* FIXME: no getFaceInfo methods for these MElements */
-    if(type == MSH_PYR_5 || type == MSH_PYR_14 || type == MSH_PYR_13){
-      //printf("Error: mesh element %d not implemented yet! \n", type);
+    //FIXME: no getFaceInfo methods for these MElements
+    if(type == TYPE_PYR){
+      printf("Error: mesh element %d not implemented yet! \n", type);
       return false;
     }
-    cell->setImmune(false);
+    Cell* cell = new Cell(element);
     cell->setInSubdomain(subdomain);
-    insertInfo = _cells[dim].insert(cell);
+    std::pair<citer, bool> insertInfo = _cells[dim].insert(cell);
     if(!insertInfo.second) delete cell;  
   }
   
@@ -89,36 +72,21 @@ bool CellComplex::insert_cells(std::vector<MElement*>& elements,
       std::vector<MVertex*> vertices;
       for(int i = 0; i < cell->getNumFacets(); i++){ 
         cell->getFacetVertices(i, vertices);
-        int type = cell->getImageMElement()->getTypeForMSH();
+        int type = cell->getImageMElement()->getType();
         int newtype = 0;
-        //FIXME: add missing boundary cell type relations
-        //FIXME: high order meshes don't work
+        //FIXME: high order meshes don't create high order cells
         if(dim == 3){
-          if(type == MSH_TET_4) newtype = MSH_TRI_3;
-          /*else if(type == MSH_TET_10) newtype = MSH_TRI_6;
-          else if(type == MSH_TET_20) newtype = MSH_TRI_9;*/
-	  /*else if(type == MSH_HEX_8) newtype = MSH_QUA_4;*/
-          /*else if(type == MSH_HEX_20) newtype = MSH_QUA_8;
-          else if(type == MSH_HEX_27) newtype = MSH_QUA_9;*/
-	  else if(type == MSH_PRI_6 
-		  && vertices.size() == 3) newtype = MSH_TRI_3;
-	  else if(type == MSH_PRI_6 
-		  && vertices.size() == 4) newtype = MSH_QUA_4;
+          if(type == TYPE_TET) newtype = MSH_TRI_3;
+	  else if(type == TYPE_HEX) newtype = MSH_QUA_4;
+	  else if(type == TYPE_PRI) {
+	    if(vertices.size() == 3) newtype = MSH_TRI_3;
+	    else if(vertices.size() == 4) newtype = MSH_QUA_4;
+	  }
         }
-        else if(dim == 2){
-	  if(type == MSH_TRI_3 || type == MSH_QUA_4) newtype = MSH_LIN_2;
-	  /*else if(type == MSH_TRI_6 
-	    || type == MSH_QUA_8) newtype = MSH_LIN_3;
-	    else if(type == MSH_TRI_9) newtype = MSH_LIN_4;
-	    else if(type == MSH_QUA_9) newtype = MSH_LIN_3;*/
-	}
-	else if(dim == 1){
-	  if(type == MSH_LIN_2) newtype = MSH_PNT;
-	  /*else if(type == MSH_LIN_3 || type == MSH_LIN_4 ||
-	    type == MSH_LIN_5 || type == MSH_LIN_6) newtype = MSH_PNT;*/
-	}  
+        else if(dim == 2) newtype = MSH_LIN_2;
+	else if(dim == 1)  newtype = MSH_PNT;
 	if(newtype == 0){
-	  //printf("Error: mesh element %d not implemented yet! \n", type);
+	  printf("Error: mesh element %d not implemented yet! \n", type);
 	  return false;
 	}
 	
@@ -126,21 +94,15 @@ bool CellComplex::insert_cells(std::vector<MElement*>& elements,
 					   cell->getImageMElement()->
 					   getPartition());
         Cell* newCell = new Cell(element);
-        newCell->setImmune(cell->getImmune());
 	newCell->setInSubdomain(subdomain);
         newCell->setDeleteImage(true);
-        insertInfo = _cells[dim-1].insert(newCell);
+        std::pair<citer, bool> insertInfo = _cells[dim-1].insert(newCell);
         if(!insertInfo.second){
           delete newCell; 
-          Cell* oldCell = *(insertInfo.first);
-          if(!subdomain){
-            int ori = cell->getFacetOri(oldCell);
-            oldCell->addCoboundaryCell( ori, cell, true);
-            cell->addBoundaryCell( ori, oldCell, true);
-          }
+          newCell = *(insertInfo.first);
         }
-        else if(!subdomain) {
-          int ori = cell->getFacetOri(vertices);
+        if(!subdomain) {
+          int ori = cell->getFacetOri(newCell);
           cell->addBoundaryCell( ori, newCell, true);
 	  newCell->addCoboundaryCell( ori, cell, true);
 	}

Geo/Homology.cpp

@@ -112,8 +112,8 @@ void Homology::findGenerators(CellComplex* cellComplex)
 {
   bool ownComplex = false;
   if(cellComplex==NULL){
-    CellComplex* cellComplex = createCellComplex(_domainEntities, 
-						 _subdomainEntities);
+    cellComplex = createCellComplex(_domainEntities, 
+						                        _subdomainEntities);
     ownComplex = true;
   }
   std::string domainString = getDomainString(_domain, _subdomain);
@@ -229,8 +229,8 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
 { 
   bool ownComplex = false;
   if(cellComplex==NULL){
-    CellComplex* cellComplex = createCellComplex(_domainEntities, 
-						 _subdomainEntities);
+    cellComplex = createCellComplex(_domainEntities, 
+						                        _subdomainEntities);
     ownComplex = true;
   }