From 5acb53136fc026ec505701be17463f404c078657 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Mon, 15 Mar 2010 15:21:35 +0000
Subject: [PATCH] Interface shift.

---
 Geo/Cell.cpp        | 51 ++++++++++++++++++++++------
 Geo/Cell.h          | 60 ++++++++++++++-------------------
 Geo/CellComplex.cpp | 82 ++++++++++++++++-----------------------------
 Geo/CellComplex.h   |  2 +-
 4 files changed, 96 insertions(+), 99 deletions(-)

diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index 0b30b209c2..379b5bdff5 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -14,14 +14,14 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const
 {  
   //cells with fever vertices first
   
-  if(c1->getNumVertices() != c2->getNumVertices()){
-    return (c1->getNumVertices() < c2->getNumVertices());
+  if(c1->getNumSortedVertices() != c2->getNumSortedVertices()){
+    return (c1->getNumSortedVertices() < c2->getNumSortedVertices());
   }
 
   // "natural number" -like ordering 
   // (the number of a vertice corresponds a digit)
   
-  for(int i=0; i < c1->getNumVertices();i++){
+  for(int i=0; i < c1->getNumSortedVertices();i++){
     if(c1->getSortedVertex(i) < c2->getSortedVertex(i)) return true;
     else if (c1->getSortedVertex(i) > c2->getSortedVertex(i)) return false;
   }
@@ -46,6 +46,37 @@ Cell::~Cell()
   if(_delimage) delete _image; 
 }
 
+bool Cell::findBoundaryCells(std::vector<Cell*>& bdCells)
+{
+  bdCells.clear();
+  MElementFactory factory;
+  for(int i = 0; i < getNumFacets(); i++){
+    std::vector<MVertex*> vertices;
+    getFacetVertices(i, vertices);
+    int type = _image->getType();
+    int newtype = 0;
+    if(_dim == 3){
+      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) newtype = MSH_LIN_2;
+    else if(_dim == 1) newtype = MSH_PNT;
+    if(newtype == 0){
+      printf("Error: mesh element %d not implemented yet! \n", type);
+      return false;
+    }
+    MElement* element = factory.create(newtype, vertices, 0, 
+				       _image->getPartition());
+    Cell* cell = new Cell(element);
+    bdCells.push_back(cell);
+  }
+  return true;
+}
+
 int Cell::getNumFacets() const 
 {
   if(_image == NULL){ 
@@ -121,11 +152,11 @@ void Cell::printCell()
 {
   printf("%d-cell: \n" , getDim());
   printf("Vertices: ");
-  for(int i = 0; i < this->getNumVertices(); i++){
+  for(int i = 0; i < this->getNumSortedVertices(); i++){
     printf("%d ", this->getSortedVertex(i));
   }
-  printf(", in subdomain: %d\n", inSubdomain());
-  printf("Combined: %d \n" , isCombined() );
+  printf(", in subdomain: %d, ", inSubdomain());
+  printf("combined: %d. \n" , isCombined() );
 };
 
 void Cell::restoreCell(){
@@ -252,7 +283,7 @@ void Cell::printCoboundary(bool orig)
 CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() 
 {  
   // use "smaller" cell as c2
-  if(c1->getNumVertices() < c2->getNumVertices()){
+  if(c1->getNumSortedVertices() < c2->getNumSortedVertices()){
     Cell* temp = c1;
     c1 = c2;
     c2 = temp;
@@ -265,12 +296,12 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
   _image = NULL;
 
   // vertices
-  _vs.reserve(c1->getNumVertices() + c2->getNumVertices());
-  for(int i = 0; i < c1->getNumVertices(); i++){
+  _vs.reserve(c1->getNumSortedVertices() + c2->getNumSortedVertices());
+  for(int i = 0; i < c1->getNumSortedVertices(); i++){
     _vs.push_back(c1->getSortedVertex(i));
   }
   std::vector<int> v;
-  for(int i = 0; i < c2->getNumVertices(); i++){
+  for(int i = 0; i < c2->getNumSortedVertices(); i++){
     if(!this->hasVertex(c2->getSortedVertex(i))){
       v.push_back(c2->getSortedVertex(i)); 
     }
diff --git a/Geo/Cell.h b/Geo/Cell.h
index 101de646bc..d79a7e10a5 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -18,10 +18,9 @@
 #include <set>
 #include <list>
 #include <map>
-#include <queue>
 #include "CellComplex.h"
 #include "MElement.h"
-//#include "GModel.h"
+#include "MVertex.h"
 
 class Cell;
 
@@ -69,38 +68,41 @@ class Cell
   // sorted vertices of this cell (used for ordering of the cells)
   std::vector<int> _vs;
 
-  virtual MVertex* getVertex(int vertex) const { 
+  // get vertex 
+  MVertex* getVertex(int vertex) const { 
     if(_image == NULL) printf("ERROR: No image mesh element for cell. \n");
     return _image->getVertex(vertex); }  
- 
+   // get the number of vertices this cell has
+  int getNumVertices() const { 
+    if(_image == NULL) printf("ERROR: No image mesh element for cell. \n");
+    return _image->getNumVertices(); }
+  // get the number of facets of this cell
+  int getNumFacets() const;
+  // get the vertices on a facet of this cell
+  void getFacetVertices(const int num, std::vector<MVertex*> &v) const;
+  
  public:
 
  Cell() : _combined(false), _index(0), _immune(false), _image(NULL), 
     _delimage(false), _subdomain(false) {}
   Cell(MElement* image);
-  virtual ~Cell();
-  
+  ~Cell();
+
   // the mesh element this cell is represented by
-  virtual MElement* getImageMElement() const { 
+  MElement* getImageMElement() const { 
     if(_image == NULL) printf("ERROR: No image mesh element for cell. \n");
     return _image; }
-  // get the number of vertices this cell has
-  virtual int getNumVertices() const { 
-    if(_image == NULL) return _vs.size();
-    else return _image->getNumVertices(); }
-  // get the number of facets of this cell
-  virtual int getNumFacets() const;
-  // 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(Cell* cell);
   // get/set whether the image mesh element should be deleted when
   // this cell gets deleted 
   // (was the mesh element in the original mesh,
   // is it needed after homology computation for visualization?) 
-  virtual void setDeleteImage(bool delimage) { _delimage = delimage; }
-  virtual bool getDeleteImage() const { return _delimage; }
-  
+  void setDeleteImage(bool delimage) { _delimage = delimage; }
+  bool getDeleteImage() const { return _delimage; }
+
+  // find the cells on the boundary of this cell 
+  bool findBoundaryCells(std::vector<Cell*>& bdCells);
+  // get boundary cell orientation
+  int getFacetOri(Cell* cell);
 
   virtual int getDim() const { return _dim; };
   virtual int getIndex() const { return _index; };
@@ -109,6 +111,7 @@ class Cell
   virtual bool getImmune() const { return _immune; };
   virtual bool inSubdomain() const { return _subdomain; }
   virtual void setInSubdomain(bool subdomain)  { _subdomain = subdomain; }
+  virtual int getNumSortedVertices() const { return _vs.size(); }
   virtual int getSortedVertex(int vertex) const { return _vs.at(vertex); }
   
   // restores the cell information to its original state before reduction
@@ -172,32 +175,22 @@ class Cell
   
   // equivalence
   bool operator==(const Cell& c2) const {  
-    if(this->getNumVertices() != c2.getNumVertices()){
+    if(this->getNumSortedVertices() != c2.getNumSortedVertices()){
       return false;
     }
-    for(int i=0; i < this->getNumVertices();i++){
+    for(int i=0; i < this->getNumSortedVertices();i++){
       if(this->getSortedVertex(i) != c2.getSortedVertex(i)){
 	return false;
       }
     }
     return true;
   }
-  /*
-  Cell operator=(const Cell& c2) {
-    Cell cell;
-    cell._ocbd = c2._ocbd;
-
-    return cell;
-  }
-  Cell(const Cell& c2){ *this = c2; }*/
 };
 
 // A cell that is a combination of cells of same dimension
 class CombinedCell : public Cell{
   
  private:
-  // sorted list of vertices
-  std::vector<int> _vs;
   // list of cells this cell is a combination of
   std::list< std::pair<int, Cell*> > _cells;
   
@@ -206,9 +199,6 @@ class CombinedCell : public Cell{
   CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false);
   ~CombinedCell() {}
   
-  int getNumVertices() const { return _vs.size(); } 
-  int getSortedVertex(int vertex) const { return _vs.at(vertex); }
-
   void getCells(std::list< std::pair<int, Cell*> >& cells) { cells = _cells; }
   int getNumCells() {return _cells.size();}
   
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index a938889832..2698f839f1 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -37,6 +37,7 @@ void CellComplex::panic_exit(){
   }
   for(unsigned int i = 0; i < _newcells.size(); i++) delete _newcells.at(i);
   _newcells.clear();
+  printf("ERROR: No proper cell complex could be constructed!\n");
 }
 
 bool CellComplex::insert_cells(std::vector<MElement*>& elements,
@@ -45,65 +46,42 @@ bool CellComplex::insert_cells(std::vector<MElement*>& elements,
   // add highest dimensional cells
   for(unsigned int i=0; i < elements.size(); i++){
     MElement* element = elements.at(i);    
-    int dim = element->getDim();
-    int type = element->getType();
     
-    // simplex types
-    if( !(type == TYPE_PNT || type == TYPE_LIN || type == TYPE_TRI
-	  || type == TYPE_TET )){
+    int type = element->getType();   
+    if(_simplicial && !(type == TYPE_PNT || type == TYPE_LIN 
+		     || type == TYPE_TRI || type == TYPE_TET) ){
       _simplicial = false;
     }
     //FIXME: no getFaceInfo methods for these MElements
     if(type == TYPE_PYR){
-      printf("Error: mesh element %d not implemented yet! \n", type);
+      printf("ERROR: mesh element %d not implemented! \n", type);
       return false;
     }
     Cell* cell = new Cell(element);
     cell->setInSubdomain(subdomain);
-    std::pair<citer, bool> insertInfo = _cells[dim].insert(cell);
+    std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell);
     if(!insertInfo.second) delete cell;  
   }
   
   // add lower dimensional cells recursively
-  MElementFactory factory;
   for (int dim = 3; dim > 0; dim--){
     for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
       Cell* cell = *cit;
-      std::vector<MVertex*> vertices;
-      for(int i = 0; i < cell->getNumFacets(); i++){ 
-        cell->getFacetVertices(i, vertices);
-        int type = cell->getImageMElement()->getType();
-        int newtype = 0;
-        //FIXME: high order meshes don't create high order cells
-        if(dim == 3){
-          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) newtype = MSH_LIN_2;
-	else if(dim == 1)  newtype = MSH_PNT;
-	if(newtype == 0){
-	  printf("Error: mesh element %d not implemented yet! \n", type);
-	  return false;
-	}
-	
-	MElement* element = factory.create(newtype, vertices, 0, 
-					   cell->getImageMElement()->
-					   getPartition());
-        Cell* newCell = new Cell(element);
+      std::vector<Cell*> bdCells;
+      if(!cell->findBoundaryCells(bdCells)) return false;
+      for(unsigned int i = 0; i < bdCells.size(); i++){
+	Cell* newCell = bdCells.at(i);
 	newCell->setInSubdomain(subdomain);
-        newCell->setDeleteImage(true);
-        std::pair<citer, bool> insertInfo = _cells[dim-1].insert(newCell);
-        if(!insertInfo.second){
-          delete newCell; 
-          newCell = *(insertInfo.first);
-        }
-        if(!subdomain) {
-          int ori = cell->getFacetOri(newCell);
-          cell->addBoundaryCell( ori, newCell, true);
+	newCell->setDeleteImage(true);
+	std::pair<citer, bool> insertInfo = 
+	  _cells[newCell->getDim()].insert(newCell);
+	if(!insertInfo.second){ // the cell was already in the cell complex
+	  delete newCell; 
+	  newCell = *(insertInfo.first); 
+	}
+	if(!subdomain) {
+	  int ori = cell->getFacetOri(newCell);
+	  cell->addBoundaryCell( ori, newCell, true);
 	  newCell->addCoboundaryCell( ori, cell, true);
 	}
       }
@@ -132,9 +110,7 @@ void CellComplex::insertCell(Cell* cell)
 {
   _newcells.push_back(cell);      
   std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell);
-  if(!insertInfo.second){
-    //printf("Warning: Cell not inserted! \n");
-  }
+  if(!insertInfo.second) printf("Warning: Cell not inserted! \n");
 }
 
 void CellComplex::removeCell(Cell* cell, bool other)
@@ -178,15 +154,15 @@ void CellComplex::enqueueCells(std::map<Cell*, int, Less_Cell>& cells,
   }
 }
 
-int CellComplex::coreduction(Cell* generator, int omitted)
+int CellComplex::coreduction(Cell* startCell, int omitted)
 {  
   int coreductions = 0;
   
   std::queue<Cell*> Q;
   std::set<Cell*, Less_Cell> Qset;
   
-  Q.push(generator);
-  Qset.insert(generator);
+  Q.push(startCell);
+  Qset.insert(startCell);
   
   std::map<Cell*, int, Less_Cell > bd_s;
   std::map<Cell*, int, Less_Cell > cbd_c;
@@ -304,7 +280,7 @@ int CellComplex::reduceComplex(bool omit)
       
       citer cit = firstCell(getDim());
       Cell* cell = *cit;
-      removeCell(cell);
+      removeCell(cell, false);
       
       omitted++;
       std::set< Cell*, Less_Cell > omittedCells;
@@ -420,9 +396,9 @@ int CellComplex::cocombine(int dim)
 
         if(!(*c1 == *c2) && abs(or1) == abs(or2)
            && inSameDomain(s, c1) && inSameDomain(s, c2)
-           && c1->getNumVertices() < getSize(dim) 
+           && c1->getNumSortedVertices() < getSize(dim) 
 	   // heuristics for mammoth cell birth control
-           && c2->getNumVertices() < getSize(dim)){
+           && c2->getNumSortedVertices() < getSize(dim)){
 
 	  removeCell(s);
           
@@ -482,9 +458,9 @@ int CellComplex::combine(int dim)
 
         if(!(*c1 == *c2) && abs(or1) == abs(or2)
            && inSameDomain(s, c1) && inSameDomain(s, c2)
-           && c1->getNumVertices() < getSize(dim) 
+           && c1->getNumSortedVertices() < getSize(dim) 
 	   // heuristics for mammoth cell birth control
-           && c2->getNumVertices() < getSize(dim)){
+           && c2->getNumSortedVertices() < getSize(dim)){
 
           removeCell(s);
           
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index 246eff7bb4..a801d1f928 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -61,7 +61,7 @@ class CellComplex
   void insertCell(Cell* cell);
   
   // queued coreduction presented in Mrozek's paper
-  int coreduction(Cell* generator, int omitted=0);
+  int coreduction(Cell* startCell, int omitted=0);
   
  public: 
   
-- 
GitLab