diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index c19622eb2463c7641e2ba1f109e612eaac27e6df..cc54acc90f203393b192a35ee7a3a5e84c191efe 100755
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -32,7 +32,6 @@ Cell::Cell(MElement* image) :
_combined(false), _index(0), _immune(false), _image(NULL),
_delimage(false), _subdomain(false)
{
- _dim = image->getDim();
_image = image;
for(int i = 0; i < getNumVertices(); i++)
_vs.push_back(getVertex(i)->getNum());
@@ -46,6 +45,7 @@ Cell::~Cell()
bool Cell::findBoundaryCells(std::vector<Cell*>& bdCells)
{
+ if(_combined) return false;
bdCells.clear();
MElementFactory factory;
for(int i = 0; i < getNumFacets(); i++){
@@ -53,7 +53,7 @@ bool Cell::findBoundaryCells(std::vector<Cell*>& bdCells)
getFacetVertices(i, vertices);
int type = _image->getType();
int newtype = 0;
- if(_dim == 3){
+ if(getDim() == 3){
if(type == TYPE_TET) newtype = MSH_TRI_3;
else if(type == TYPE_HEX) newtype = MSH_QUA_4;
else if(type == TYPE_PRI) {
@@ -61,8 +61,8 @@ bool Cell::findBoundaryCells(std::vector<Cell*>& bdCells)
else if(vertices.size() == 4) newtype = MSH_QUA_4;
}
}
- else if(_dim == 2) newtype = MSH_LIN_2;
- else if(_dim == 1) newtype = MSH_PNT;
+ else if(getDim() == 2) newtype = MSH_LIN_2;
+ else if(getDim() == 1) newtype = MSH_PNT;
if(newtype == 0){
printf("Error: mesh element %d not implemented yet! \n", type);
return false;
@@ -77,7 +77,7 @@ bool Cell::findBoundaryCells(std::vector<Cell*>& bdCells)
int Cell::getNumFacets() const
{
- if(_image == NULL){
+ if(_image == NULL || _combined){
printf("ERROR: No image mesh element for cell. \n");
return 0;
}
@@ -90,7 +90,7 @@ int Cell::getNumFacets() const
void Cell::getFacetVertices(const int num, std::vector<MVertex*> &v) const
{
- if(_image == NULL){
+ if(_image == NULL || _combined){
printf("ERROR: No image mesh element for cell. \n");
return;
}
@@ -103,7 +103,7 @@ void Cell::getFacetVertices(const int num, std::vector<MVertex*> &v) const
int Cell::findBoundaryCellOrientation(Cell* cell)
{
- if(_image == NULL){
+ if(_image == NULL || _combined){
printf("ERROR: No image mesh element for cell. \n");
return 0;
}
@@ -146,6 +146,15 @@ bool Cell::hasVertex(int vertex) const
else return false;
}
+bool CombinedCell::hasVertex(int vertex) const
+{
+ for(std::map<Cell*, int, Less_Cell>::const_iterator cit = _cells.begin();
+ cit != _cells.end(); cit++){
+ if(cit->first->hasVertex(vertex)) return true;
+ }
+ return false;
+}
+
void Cell::printCell()
{
printf("%d-cell: \n" , getDim());
@@ -282,12 +291,9 @@ CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
}
_num = ++_globalNum;
-
_index = c1->getIndex();
- _dim = c1->getDim();
_subdomain = c1->inSubdomain();
_combined = true;
- _image = NULL;
// cells
c1->getCells(_cells);
diff --git a/Geo/Cell.h b/Geo/Cell.h
index 0d0dd3e6caf3312c295e3852176ea06a334d8a9f..ae0ed442ca44a425c140fd7c39e47ae6feced0c7 100644
--- a/Geo/Cell.h
+++ b/Geo/Cell.h
@@ -35,9 +35,6 @@ class Cell
{
protected:
- // cell dimension
- int _dim;
-
// whether this cell belongs to a subdomain
// used in relative homology computation
bool _subdomain;
@@ -58,11 +55,11 @@ class Cell
// cell complex
std::map<Cell*, int, Less_Cell > _obd;
std::map<Cell*, int, Less_Cell > _ocbd;
+
+ private:
// sorted vertices of this cell (used for ordering of the cells)
std::vector<int> _vs;
-
-
// the mesh element that is the image of this cell
MElement* _image;
@@ -71,13 +68,9 @@ class Cell
bool _delimage;
// get vertex
- MVertex* getVertex(int vertex) const {
- if(_image == NULL) printf("ERROR: No image mesh element for cell. \n");
- return _image->getVertex(vertex); }
+ MVertex* getVertex(int vertex) const { 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->getNumPrimaryVertices(); }
+ int getNumVertices() const { return _image->getNumPrimaryVertices(); }
// get the number of facets of this cell
int getNumFacets() const;
// get the vertices on a facet of this cell
@@ -92,22 +85,25 @@ class Cell
// the mesh element this cell is represented by
MElement* getImageMElement() const {
- if(_image == NULL) printf("ERROR: No image mesh element for cell. \n");
+ if(_image == NULL || _combined){
+ printf("ERROR: No image mesh element for cell. \n"); }
return _image; }
+ // get the cell dimension
+ virtual int getDim() const { return _image->getDim(); };
// 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?)
- void setDeleteImage(bool delimage) { _delimage = delimage; }
- bool getDeleteImage() const { return _delimage; }
+ void setDeleteImage(bool delimage) { if(!_combined) _delimage = delimage; }
+ bool getDeleteImage() const {
+ if(!_combined) return _delimage;
+ else return false; }
// find the cells on the boundary of this cell
bool findBoundaryCells(std::vector<Cell*>& bdCells);
// get boundary cell orientation
int findBoundaryCellOrientation(Cell* cell);
-
- int getDim() const { return _dim; };
int getIndex() const { return _index; };
void setIndex(int index) { _index = index; };
void setImmune(bool immune) { _immune = immune; };
@@ -213,10 +209,14 @@ class CombinedCell : public Cell{
CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co=false);
~CombinedCell() {}
-
+
+ int getDim() const { return _cells.begin()->first->getDim(); }
+ int getNumSortedVertices() const { return 0; }
+ int getSortedVertex(int vertex) const { return 0; }
void getCells(std::map< Cell*, int, Less_Cell >& cells) { cells = _cells; }
int getNumCells() const {return _cells.size();}
int getNum() const { return _num; }
+ bool hasVertex(int vertex) const;
bool operator==(const Cell& c2) const {
if(this->isCombined() != c2.isCombined()) return false;
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index e15c8a0af6d057d52b09109d6db0657155e975a1..3dcb47e8201831e03bc85d05907457595cba90d0 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -691,13 +691,10 @@ HomologySequence::HomologySequence(ChainComplex* subcomplex,
mpz_t elem;
mpz_init_set_si(elem, -1);
-
for(int i = 0; i < 4; i++){
-
- //printf("Dimension %d. \n", i);
_Ic_sub[i] = NULL;
_Ic_rel[i] = NULL;
-
+
_Dh[i] = NULL;
_invDh[i] = NULL;
@@ -705,7 +702,38 @@ HomologySequence::HomologySequence(ChainComplex* subcomplex,
_Ih[i] = NULL;
_invJh[i] = NULL;
_invIh[i] = NULL;
+
+ }
+
+ findIcMaps();
+
+ /*findDhMap(1);
+ findInvIhMap(0);
+ blockHBasis(_Dh[1], _invIh[0], _subcomplex, 0);*/
+
+ /*findJhMap(1);
+ findInvDhMap(1);
+ blockHBasis(_Jh[1], _invDh[1], _relcomplex, 1);*/
+
+}
+
+HomologySequence::~HomologySequence()
+{
+ for(int i = 0; i < 4; i++){
+ destroy_gmp_matrix(_Ic_sub[i]);
+ destroy_gmp_matrix(_Ic_rel[i]);
+ destroy_gmp_matrix(_Ih[i]);
+ destroy_gmp_matrix(_Jh[i]);
+ destroy_gmp_matrix(_invIh[i]);
+ destroy_gmp_matrix(_invJh[i]);
+ destroy_gmp_matrix(_Dh[i]);
+ destroy_gmp_matrix(_invDh[i]);
+ }
+}
+void HomologySequence::findIcMaps()
+{
+ for(int i = 0; i < 4; i++){
mpz_t one;
mpz_init_set_si(one, 1);
if(_complex->getBasisSize(i, 0) > 0
@@ -739,104 +767,95 @@ HomologySequence::HomologySequence(ChainComplex* subcomplex,
}
}
mpz_clear(one);
+ }
+}
- if(_Ic_sub[i] != NULL
- && _complex->getBasisSize(i, 3) > 0
- && _subcomplex->getBasisSize(i, 3) > 0){
- gmp_matrix* IH = copy_gmp_matrix(_Ic_sub[i], 1, 1,
- gmp_matrix_rows(_Ic_sub[i]),
- gmp_matrix_cols(_Ic_sub[i]));
- gmp_matrix_right_mult(IH, _subcomplex->getBasis(i, 3));
- _Ih[i] = createIncMap(IH, _complex->getBasis(i, 3));
- }
- if(_Ic_sub[i] != NULL
- && _complex->getBasisSize(i, 3) > 0
- && _subcomplex->getBasisSize(i, 3) > 0){
- gmp_matrix* IH = copy_gmp_matrix(_Ic_sub[i], 1, 1,
- gmp_matrix_rows(_Ic_sub[i]),
- gmp_matrix_cols(_Ic_sub[i]));
- gmp_matrix_transp(IH);
- gmp_matrix_right_mult(IH, _complex->getBasis(i, 3));
- _invIh[i] = createIncMap(IH, _subcomplex->getBasis(i, 3));
- }
-
- if(_Ic_rel[i] != NULL
- && _complex->getBasisSize(i, 3) > 0
- && _relcomplex->getBasisSize(i, 3) > 0){
- gmp_matrix* JH = copy_gmp_matrix(_Ic_rel[i], 1, 1,
- gmp_matrix_rows(_Ic_rel[i]),
- gmp_matrix_cols(_Ic_rel[i]));
- gmp_matrix_transp(JH);
- gmp_matrix_right_mult(JH, _complex->getBasis(i, 3));
- _Jh[i] = createIncMap(JH, _relcomplex->getBasis(i, 3));
- }
- if(_Ic_rel[i] != NULL
- && _complex->getBasisSize(i, 3) > 0
- && _relcomplex->getBasisSize(i, 3) > 0){
- gmp_matrix* JH = copy_gmp_matrix(_Ic_rel[i], 1, 1,
- gmp_matrix_rows(_Ic_rel[i]),
- gmp_matrix_cols(_Ic_rel[i]));
- gmp_matrix_right_mult(JH, _relcomplex->getBasis(i, 3));
- _invJh[i] = createIncMap(JH, _complex->getBasis(i, 3));
- }
+void HomologySequence::findIhMap(int i)
+{
+ if(_Ic_sub[i] != NULL
+ && _complex->getBasisSize(i, 3) > 0
+ && _subcomplex->getBasisSize(i, 3) > 0){
+ gmp_matrix* IH = copy_gmp_matrix(_Ic_sub[i], 1, 1,
+ gmp_matrix_rows(_Ic_sub[i]),
+ gmp_matrix_cols(_Ic_sub[i]));
+ gmp_matrix_right_mult(IH, _subcomplex->getBasis(i, 3));
+ _Ih[i] = createIncMap(IH, _complex->getBasis(i, 3));
+ }
+}
- //printMatrix(_Ih[i]);
- //printMatrix(_invIh[i]);
- //printMatrix(_Jh[i]);
- //printMatrix(_invJh[i]);
-
- if(i > 0 && _relcomplex->getBasisSize(i, 3) > 0
- && _subcomplex->getBasisSize(i-1, 3) > 0
- && _complex->getBoundaryOp(i) != NULL){
- gmp_matrix* JDIH = copy_gmp_matrix(_Ic_sub[i-1], 1, 1,
- gmp_matrix_rows(_Ic_sub[i-1]),
- gmp_matrix_cols(_Ic_sub[i-1]));
- gmp_matrix_transp(JDIH);
- gmp_matrix_right_mult(JDIH, _complex->getBoundaryOp(i));
- gmp_matrix_right_mult(JDIH, _Ic_rel[i]);
- gmp_matrix_right_mult(JDIH, _relcomplex->getBasis(i, 3));
- _Dh[i] = createIncMap(JDIH, _subcomplex->getBasis(i-1, 3));
- }
+void HomologySequence::findInvIhMap(int i)
+{
+ if(_Ic_sub[i] != NULL
+ && _complex->getBasisSize(i, 3) > 0
+ && _subcomplex->getBasisSize(i, 3) > 0){
+ gmp_matrix* IH = copy_gmp_matrix(_Ic_sub[i], 1, 1,
+ gmp_matrix_rows(_Ic_sub[i]),
+ gmp_matrix_cols(_Ic_sub[i]));
+ gmp_matrix_transp(IH);
+ gmp_matrix_right_mult(IH, _complex->getBasis(i, 3));
+ _invIh[i] = createIncMap(IH, _subcomplex->getBasis(i, 3));
+ }
+}
- if(i > 0 && _relcomplex->getBasisSize(i, 3) > 0
- && _subcomplex->getBasisSize(i-1, 3) > 0
- && _complex->getBoundaryOp(i) != NULL){
- gmp_matrix* JDIH = copy_gmp_matrix(_Ic_rel[i], 1, 1,
- gmp_matrix_rows(_Ic_rel[i]),
- gmp_matrix_cols(_Ic_rel[i]));
- gmp_matrix_transp(JDIH);
- gmp_matrix* bd = _complex->getBoundaryOp(i);
- gmp_matrix_transp(bd);
- gmp_matrix_right_mult(JDIH, bd);
- gmp_matrix_transp(bd);
- gmp_matrix_right_mult(JDIH, _Ic_sub[i-1]);
- gmp_matrix_right_mult(JDIH, _subcomplex->getBasis(i-1, 3));
- _invDh[i] = createIncMap(JDIH, _relcomplex->getBasis(i, 3));
- }
-
- //printMatrix(_Dh[i]);
- //printMatrix(_invDh[i]);
+void HomologySequence::findJhMap(int i)
+{
+ if(_Ic_rel[i] != NULL
+ && _complex->getBasisSize(i, 3) > 0
+ && _relcomplex->getBasisSize(i, 3) > 0){
+ gmp_matrix* JH = copy_gmp_matrix(_Ic_rel[i], 1, 1,
+ gmp_matrix_rows(_Ic_rel[i]),
+ gmp_matrix_cols(_Ic_rel[i]));
+ gmp_matrix_transp(JH);
+ gmp_matrix_right_mult(JH, _complex->getBasis(i, 3));
+ _Jh[i] = createIncMap(JH, _relcomplex->getBasis(i, 3));
+ }
+}
+void HomologySequence::findInvJhMap(int i)
+{
+ if(_Ic_rel[i] != NULL
+ && _complex->getBasisSize(i, 3) > 0
+ && _relcomplex->getBasisSize(i, 3) > 0){
+ gmp_matrix* JH = copy_gmp_matrix(_Ic_rel[i], 1, 1,
+ gmp_matrix_rows(_Ic_rel[i]),
+ gmp_matrix_cols(_Ic_rel[i]));
+ gmp_matrix_right_mult(JH, _relcomplex->getBasis(i, 3));
+ _invJh[i] = createIncMap(JH, _complex->getBasis(i, 3));
}
- for(int i = 0; i < 3; i++){
- blockHBasis(_Dh[i+1], _invIh[i], _subcomplex, i);
- blockHBasis(_Ih[i], _invJh[i], _complex, i);
- blockHBasis(_Jh[i], _invDh[i], _relcomplex, i);
+}
+
+void HomologySequence::findDhMap(int i)
+{
+ if(i > 0 && _relcomplex->getBasisSize(i, 3) > 0
+ && _subcomplex->getBasisSize(i-1, 3) > 0
+ && _complex->getBoundaryOp(i) != NULL){
+ gmp_matrix* JDIH = copy_gmp_matrix(_Ic_sub[i-1], 1, 1,
+ gmp_matrix_rows(_Ic_sub[i-1]),
+ gmp_matrix_cols(_Ic_sub[i-1]));
+ gmp_matrix_transp(JDIH);
+ gmp_matrix_right_mult(JDIH, _complex->getBoundaryOp(i));
+ gmp_matrix_right_mult(JDIH, _Ic_rel[i]);
+ gmp_matrix_right_mult(JDIH, _relcomplex->getBasis(i, 3));
+ _Dh[i] = createIncMap(JDIH, _subcomplex->getBasis(i-1, 3));
}
-
}
-HomologySequence::~HomologySequence()
+void HomologySequence::findInvDhMap(int i)
{
- for(int i = 0; i < 4; i++){
- destroy_gmp_matrix(_Ic_sub[i]);
- destroy_gmp_matrix(_Ic_rel[i]);
- destroy_gmp_matrix(_Ih[i]);
- destroy_gmp_matrix(_Jh[i]);
- destroy_gmp_matrix(_invIh[i]);
- destroy_gmp_matrix(_invJh[i]);
- destroy_gmp_matrix(_Dh[i]);
- destroy_gmp_matrix(_invDh[i]);
+ if(i > 0 && _relcomplex->getBasisSize(i, 3) > 0
+ && _subcomplex->getBasisSize(i-1, 3) > 0
+ && _complex->getBoundaryOp(i) != NULL){
+ gmp_matrix* JDIH = copy_gmp_matrix(_Ic_rel[i], 1, 1,
+ gmp_matrix_rows(_Ic_rel[i]),
+ gmp_matrix_cols(_Ic_rel[i]));
+ gmp_matrix_transp(JDIH);
+ gmp_matrix* bd = _complex->getBoundaryOp(i);
+ gmp_matrix_transp(bd);
+ gmp_matrix_right_mult(JDIH, bd);
+ gmp_matrix_transp(bd);
+ gmp_matrix_right_mult(JDIH, _Ic_sub[i-1]);
+ gmp_matrix_right_mult(JDIH, _subcomplex->getBasis(i-1, 3));
+ _invDh[i] = createIncMap(JDIH, _relcomplex->getBasis(i, 3));
}
}
@@ -913,15 +932,14 @@ gmp_matrix* HomologySequence::createIncMap(gmp_matrix* domBasis,
return LB;
}
-
-gmp_matrix* HomologySequence::removeZeroCols(gmp_matrix* matrix) // FIXME
+gmp_matrix* HomologySequence::removeZeroCols(gmp_matrix* matrix)
{
mpz_t elem;
mpz_init(elem);
int rows = gmp_matrix_rows(matrix);
int cols = gmp_matrix_cols(matrix);
-
+ //printMatrix(matrix);
std::vector<int> zcols;
for(int j = 1; j <= cols; j++){
@@ -936,7 +954,7 @@ gmp_matrix* HomologySequence::removeZeroCols(gmp_matrix* matrix) // FIXME
if(zcol) zcols.push_back(j);
}
if(zcols.empty()) return matrix;
-
+
gmp_matrix* newMatrix = create_gmp_matrix_zero(rows, cols-zcols.size());
if(cols-zcols.size() == 0) return newMatrix;
@@ -949,26 +967,32 @@ gmp_matrix* HomologySequence::removeZeroCols(gmp_matrix* matrix) // FIXME
gmp_matrix_set_elem(elem, i, j-k, newMatrix);
}
}
-
+ //printMatrix(newMatrix);
destroy_gmp_matrix(matrix);
mpz_clear(elem);
return newMatrix;
}
-
+
void HomologySequence::blockHBasis(gmp_matrix* block1T,
gmp_matrix* block2T,
ChainComplex* complex, int dim)
{
+
+ printMatrix(block1T);
+ printMatrix(block2T);
+
if(block1T == NULL && block2T == NULL) return;
gmp_matrix* Hbasis = complex->getBasis(dim, 3);
if(block1T == NULL && block2T != NULL){
gmp_matrix_right_mult(Hbasis, block2T);
+ printMatrix(Hbasis);
return;
}
if(block1T != NULL && block2T == NULL){
gmp_matrix_right_mult(Hbasis, block1T);
+ printMatrix(Hbasis);
return;
}
@@ -977,14 +1001,18 @@ void HomologySequence::blockHBasis(gmp_matrix* block1T,
gmp_matrix* temp1 = copy_gmp_matrix(Hbasis, 1, 1, rows, cols);
gmp_matrix* temp2 = copy_gmp_matrix(Hbasis, 1, 1, rows, cols);
+ printMatrix(temp1);
+ printMatrix(temp2);
+
gmp_matrix_right_mult(temp1, block1T);
gmp_matrix_right_mult(temp2, block2T);
-
+ printMatrix(temp1);
+ printMatrix(temp2);
temp1 = removeZeroCols(temp1);
temp2 = removeZeroCols(temp2);
- //printMatrix(temp1);
- //printMatrix(temp2);
+ printMatrix(temp1);
+ printMatrix(temp2);
int bcol = gmp_matrix_cols(temp1);
@@ -998,7 +1026,7 @@ void HomologySequence::blockHBasis(gmp_matrix* block1T,
}
}
- //printMatrix(Hbasis);
+ printMatrix(Hbasis);
mpz_clear(elem);
destroy_gmp_matrix(temp1);
destroy_gmp_matrix(temp2);
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index ea1c7fe210fd4fac49cdee7382f0a4f4b4e0b314..7092e3b11ff8ecf12ef88c6ecc381efc22924176 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -204,6 +204,14 @@ class HomologySequence
gmp_matrix* _Dh[4];
gmp_matrix* _invDh[4];
+ void findIcMaps();
+ void findIhMap(int i);
+ void findInvIhMap(int i);
+ void findJhMap(int i);
+ void findInvJhMap(int i);
+ void findDhMap(int i);
+ void findInvDhMap(int i);
+
public:
HomologySequence(ChainComplex* subcomplex, ChainComplex* complex,
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 3603054fc02ca7ba611f4128d79d79c76d56db4c..9bec179e16ec03c7e07c4f0e0bf69183f09b9561 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -122,7 +122,8 @@ void Homology::findGenerators(CellComplex* cellComplex)
Msg::StatusBar(1, false, "Reducing...");
double t1 = Cpu();
- int omitted = cellComplex->reduceComplex();
+ int omitted = cellComplex->reduceComplex();
+
double t2 = Cpu();
Msg::Info("Cell Complex reduction complete (%g s).", t2 - t1);
Msg::Info("%d volumes, %d faces, %d edges and %d vertices.",
@@ -374,7 +375,7 @@ void Homology::findHomSequence(){
std::string name = "H" + dimension + domainString + generator;
std::map<Cell*, int, Less_Cell> protoChain;
- chains->getBasisChain(protoChain, i, j, 3);
+ chains->getBasisChain(protoChain, i, j, 3, true);
Chain* chain = new Chain(protoChain, cellComplex, _model, name,
chains->getTorsion(j,i));
if(chain->getSize() == 0) {
@@ -612,19 +613,19 @@ int Chain::createPGroup()
physicalMap[entityNum] = physicalInfo;
// hide mesh
- opt_mesh_points(0, GMSH_SET, 0);
+ /*opt_mesh_points(0, GMSH_SET, 0);
opt_mesh_lines(0, GMSH_SET, 0);
opt_mesh_triangles(0, GMSH_SET, 0);
opt_mesh_quadrangles(0, GMSH_SET, 0);
opt_mesh_tetrahedra(0, GMSH_SET, 0);
opt_mesh_hexahedra(0, GMSH_SET, 0);
opt_mesh_prisms(0, GMSH_SET, 0);
- opt_mesh_pyramids(0, GMSH_SET, 0);
+ opt_mesh_pyramids(0, GMSH_SET, 0);*/
// show post-processing normals, tangents and element boundaries
//opt_view_normals(0, GMSH_SET, 20);
//opt_view_tangents(0, GMSH_SET, 20);
- opt_view_show_element(0, GMSH_SET, 1);
+ //opt_view_show_element(0, GMSH_SET, 1);
if(!data.empty()){
_model->storeChain(getDim(), entityMap, physicalMap);
diff --git a/benchmarks/homology/trefoil.geo b/benchmarks/homology/trefoil.geo
index 52aaf789b479a03014dc9ebaa840dee4e06737da..d5a75d9b4ffe982f05f299162008592526b3f0f3 100644
--- a/benchmarks/homology/trefoil.geo
+++ b/benchmarks/homology/trefoil.geo
@@ -282,4 +282,4 @@ Mesh 3;
HomGen("trefoil.msh") = {{179},{}};
HomGen("trefoil.msh") = {{179},{180}};
-//HomCut("trefoil.msh") = {{179},{}};
+HomCut("trefoil.msh") = {{179},{}};
diff --git a/contrib/kbipack/gmp_matrix.h b/contrib/kbipack/gmp_matrix.h
index 318dd85657c06d91b0bea17a99be6563d220af2d..931f3828036a15ed73221d9d54a43cdddabeebe5 100644
--- a/contrib/kbipack/gmp_matrix.h
+++ b/contrib/kbipack/gmp_matrix.h
@@ -127,10 +127,6 @@ gmp_matrix_row_inz (size_t r, size_t c1, size_t c2,
int
gmp_matrix_transp(gmp_matrix * M);
-/* A <- A+B */
-int
-gmp_matrix_sum(gmp_matrix * A, const gmp_matrix * B);
-
/* A <- A*B */
int
gmp_matrix_right_mult(gmp_matrix * A, const gmp_matrix * B);