Skip to content
Snippets Groups Projects
Commit 159cf407 authored by Matti Pellika's avatar Matti Pellika
Browse files

Enabled dummy mpz interface when GMP library is not present. Cleaning.

parent ecc5aa6f
No related branches found
No related tags found
No related merge requests found
...@@ -425,12 +425,12 @@ if(ENABLE_GMM) ...@@ -425,12 +425,12 @@ if(ENABLE_GMM)
endif(ENABLE_GMM) endif(ENABLE_GMM)
if(ENABLE_KBIPACK) if(ENABLE_KBIPACK)
find_library(GMP_LIB NAMES libgmp.a libgmp.dll.a)
if(GMP_LIB)
message("-- Found GMP")
set(HAVE_KBIPACK TRUE) set(HAVE_KBIPACK TRUE)
list(APPEND CONFIG_OPTIONS "Kbipack") list(APPEND CONFIG_OPTIONS "Kbipack")
add_subdirectory(contrib/kbipack) add_subdirectory(contrib/kbipack)
find_library(GMP_LIB NAMES libgmp.a libgmp.dll.a)
if(GMP_LIB)
message("-- Found GMP")
set(HAVE_GMP TRUE) set(HAVE_GMP TRUE)
list(APPEND EXTERNAL_LIBRARIES ${GMP_LIB}) list(APPEND EXTERNAL_LIBRARIES ${GMP_LIB})
endif(GMP_LIB) endif(GMP_LIB)
......
...@@ -10,8 +10,8 @@ ...@@ -10,8 +10,8 @@
#if defined(HAVE_KBIPACK) #if defined(HAVE_KBIPACK)
bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const { bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const
{
//cells with fever vertices first //cells with fever vertices first
if(c1->getNumVertices() != c2->getNumVertices()){ if(c1->getNumVertices() != c2->getNumVertices()){
...@@ -32,7 +32,8 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const { ...@@ -32,7 +32,8 @@ bool Less_Cell::operator()(const Cell* c1, const Cell* c2) const {
Cell::Cell(MElement* image, bool subdomain, bool boundary) : Cell::Cell(MElement* image, bool subdomain, bool boundary) :
_combined(false), _index(0), _immune(false), _image(NULL), _combined(false), _index(0), _immune(false), _image(NULL),
_deleteImage(false) { _deleteImage(false)
{
_onDomainBoundary = boundary; _onDomainBoundary = boundary;
_inSubdomain = subdomain; _inSubdomain = subdomain;
_dim = image->getDim(); _dim = image->getDim();
...@@ -42,25 +43,30 @@ Cell::Cell(MElement* image, bool subdomain, bool boundary) : ...@@ -42,25 +43,30 @@ Cell::Cell(MElement* image, bool subdomain, bool boundary) :
std::sort(_vs.begin(), _vs.end()); std::sort(_vs.begin(), _vs.end());
} }
Cell::~Cell() { Cell::~Cell()
{
if(_deleteImage) delete _image; if(_deleteImage) delete _image;
} }
bool Cell::hasVertex(int vertex) const { bool Cell::hasVertex(int vertex) const
{
std::vector<int>::const_iterator it = std::find(_vs.begin(), _vs.end(), std::vector<int>::const_iterator it = std::find(_vs.begin(), _vs.end(),
vertex); vertex);
if (it != _vs.end()) return true; if (it != _vs.end()) return true;
else return false; else return false;
} }
int Cell::getNumFacets() const { int Cell::getNumFacets() const
{
if(getDim() == 0) return 0; if(getDim() == 0) return 0;
else if(getDim() == 1) return 2; else if(getDim() == 1) return 2;
else if(getDim() == 2) return _image->getNumEdges(); else if(getDim() == 2) return _image->getNumEdges();
else if(getDim() == 3) return _image->getNumFaces(); else if(getDim() == 3) return _image->getNumFaces();
else return 0; else return 0;
} }
void Cell::getFacetVertices(const int num, std::vector<MVertex*> &v) const {
void Cell::getFacetVertices(const int num, std::vector<MVertex*> &v) const
{
if(getDim() == 0) return; if(getDim() == 0) return;
else if(getDim() == 1) { v.resize(1); v[0] = getVertex(num); } else if(getDim() == 1) { v.resize(1); v[0] = getVertex(num); }
else if(getDim() == 2) _image->getEdgeVertices(num, v); else if(getDim() == 2) _image->getEdgeVertices(num, v);
...@@ -68,8 +74,8 @@ void Cell::getFacetVertices(const int num, std::vector<MVertex*> &v) const { ...@@ -68,8 +74,8 @@ void Cell::getFacetVertices(const int num, std::vector<MVertex*> &v) const {
return; return;
} }
int Cell::getFacetOri(std::vector<MVertex*> &v)
int Cell::getFacetOri(std::vector<MVertex*> &v) { {
if(getDim() == 0) return 0; if(getDim() == 0) return 0;
else if(getDim() == 1){ else if(getDim() == 1){
if(v.size() != 1) return 0; if(v.size() != 1) return 0;
...@@ -97,7 +103,8 @@ int Cell::getFacetOri(std::vector<MVertex*> &v) { ...@@ -97,7 +103,8 @@ int Cell::getFacetOri(std::vector<MVertex*> &v) {
else return 0; else return 0;
} }
void Cell::printCell() { void Cell::printCell()
{
printf("%d-cell %d: \n" , getDim(), getNum()); printf("%d-cell %d: \n" , getDim(), getNum());
printf("Vertices: "); printf("Vertices: ");
for(int i = 0; i < this->getNumVertices(); i++){ for(int i = 0; i < this->getNumVertices(); i++){
...@@ -115,7 +122,8 @@ void Cell::restoreCell(){ ...@@ -115,7 +122,8 @@ void Cell::restoreCell(){
_immune = false; _immune = false;
} }
bool Cell::addBoundaryCell(int orientation, Cell* cell, bool org) { bool Cell::addBoundaryCell(int orientation, Cell* cell, bool org)
{
if(org) _obd.insert( std::make_pair(cell, orientation ) ); if(org) _obd.insert( std::make_pair(cell, orientation ) );
biter it = _boundary.find(cell); biter it = _boundary.find(cell);
if(it != _boundary.end()){ if(it != _boundary.end()){
...@@ -131,7 +139,8 @@ bool Cell::addBoundaryCell(int orientation, Cell* cell, bool org) { ...@@ -131,7 +139,8 @@ bool Cell::addBoundaryCell(int orientation, Cell* cell, bool org) {
return true; return true;
} }
bool Cell::addCoboundaryCell(int orientation, Cell* cell, bool org) { bool Cell::addCoboundaryCell(int orientation, Cell* cell, bool org)
{
if(org) _ocbd.insert( std::make_pair(cell, orientation ) ); if(org) _ocbd.insert( std::make_pair(cell, orientation ) );
biter it = _coboundary.find(cell); biter it = _coboundary.find(cell);
if(it != _coboundary.end()){ if(it != _coboundary.end()){
...@@ -147,7 +156,8 @@ bool Cell::addCoboundaryCell(int orientation, Cell* cell, bool org) { ...@@ -147,7 +156,8 @@ bool Cell::addCoboundaryCell(int orientation, Cell* cell, bool org) {
return true; return true;
} }
int Cell::removeBoundaryCell(Cell* cell, bool other) { int Cell::removeBoundaryCell(Cell* cell, bool other)
{
biter it = _boundary.find(cell); biter it = _boundary.find(cell);
if(it != _boundary.end()){ if(it != _boundary.end()){
_boundary.erase(it); _boundary.erase(it);
...@@ -157,7 +167,9 @@ int Cell::removeBoundaryCell(Cell* cell, bool other) { ...@@ -157,7 +167,9 @@ int Cell::removeBoundaryCell(Cell* cell, bool other) {
return 0; return 0;
} }
int Cell::removeCoboundaryCell(Cell* cell, bool other) {
int Cell::removeCoboundaryCell(Cell* cell, bool other)
{
biter it = _coboundary.find(cell); biter it = _coboundary.find(cell);
if(it != _coboundary.end()){ if(it != _coboundary.end()){
_coboundary.erase(it); _coboundary.erase(it);
...@@ -167,7 +179,8 @@ int Cell::removeCoboundaryCell(Cell* cell, bool other) { ...@@ -167,7 +179,8 @@ int Cell::removeCoboundaryCell(Cell* cell, bool other) {
return 0; return 0;
} }
bool Cell::hasBoundary(Cell* cell, bool org){ bool Cell::hasBoundary(Cell* cell, bool org)
{
if(!org){ if(!org){
biter it = _boundary.find(cell); biter it = _boundary.find(cell);
if(it != _boundary.end()) return true; if(it != _boundary.end()) return true;
...@@ -180,7 +193,8 @@ bool Cell::hasBoundary(Cell* cell, bool org){ ...@@ -180,7 +193,8 @@ bool Cell::hasBoundary(Cell* cell, bool org){
} }
} }
bool Cell::hasCoboundary(Cell* cell, bool org){ bool Cell::hasCoboundary(Cell* cell, bool org)
{
if(!org){ if(!org){
biter it = _coboundary.find(cell); biter it = _coboundary.find(cell);
if(it != _coboundary.end()) return true; if(it != _coboundary.end()) return true;
...@@ -193,14 +207,16 @@ bool Cell::hasCoboundary(Cell* cell, bool org){ ...@@ -193,14 +207,16 @@ bool Cell::hasCoboundary(Cell* cell, bool org){
} }
} }
void Cell::makeDualCell(){ void Cell::makeDualCell()
{
std::map<Cell*, int, Less_Cell > temp = _boundary; std::map<Cell*, int, Less_Cell > temp = _boundary;
_boundary = _coboundary; _boundary = _coboundary;
_coboundary = temp; _coboundary = temp;
_dim = 3-_dim; _dim = 3-_dim;
} }
void Cell::printBoundary(bool org) { void Cell::printBoundary(bool org)
{
for(biter it = firstBoundary(org); it != lastBoundary(org); it++){ for(biter it = firstBoundary(org); it != lastBoundary(org); it++){
printf("Boundary cell orientation: %d ", (*it).second); printf("Boundary cell orientation: %d ", (*it).second);
Cell* cell2 = (*it).first; Cell* cell2 = (*it).first;
...@@ -210,7 +226,9 @@ void Cell::printBoundary(bool org) { ...@@ -210,7 +226,9 @@ void Cell::printBoundary(bool org) {
printf("Cell boundary is empty. \n"); printf("Cell boundary is empty. \n");
} }
} }
void Cell::printCoboundary(bool org) {
void Cell::printCoboundary(bool org)
{
for(biter it = firstCoboundary(org); it != lastCoboundary(org); it++){ for(biter it = firstCoboundary(org); it != lastCoboundary(org); it++){
printf("Coboundary cell orientation: %d, ", (*it).second); printf("Coboundary cell orientation: %d, ", (*it).second);
Cell* cell2 = (*it).first; Cell* cell2 = (*it).first;
...@@ -221,7 +239,8 @@ void Cell::printCoboundary(bool org) { ...@@ -221,7 +239,8 @@ void Cell::printCoboundary(bool org) {
} }
} }
CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell() { CombinedCell::CombinedCell(Cell* c1, Cell* c2, bool orMatch, bool co) : Cell()
{
// use "smaller" cell as c2 // use "smaller" cell as c2
if(c1->getNumVertices() < c2->getNumVertices()){ if(c1->getNumVertices() < c2->getNumVertices()){
Cell* temp = c1; Cell* temp = c1;
......
...@@ -10,7 +10,9 @@ ...@@ -10,7 +10,9 @@
#if defined(HAVE_KBIPACK) #if defined(HAVE_KBIPACK)
CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain ){ CellComplex::CellComplex( std::vector<GEntity*> domain,
std::vector<GEntity*> subdomain )
{
_domain = domain; _domain = domain;
_subdomain = subdomain; _subdomain = subdomain;
...@@ -48,8 +50,8 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su ...@@ -48,8 +50,8 @@ CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> su
} }
void CellComplex::find_boundary(std::vector<GEntity*>& domain, void CellComplex::find_boundary(std::vector<GEntity*>& domain,
std::vector<GEntity*>& subdomain){ std::vector<GEntity*>& subdomain)
{
// determine mesh entities on boundary of the domain // determine mesh entities on boundary of the domain
bool duplicate = false; bool duplicate = false;
for(unsigned int j=0; j < _domain.size(); j++){ for(unsigned int j=0; j < _domain.size(); j++){
...@@ -113,8 +115,8 @@ void CellComplex::find_boundary(std::vector<GEntity*>& domain, ...@@ -113,8 +115,8 @@ void CellComplex::find_boundary(std::vector<GEntity*>& domain,
} }
} }
void CellComplex::insert_cells(bool subdomain, bool boundary){ void CellComplex::insert_cells(bool subdomain, bool boundary)
{
std::vector<GEntity*> domain; std::vector<GEntity*> domain;
if(subdomain) domain = _subdomain; if(subdomain) domain = _subdomain;
...@@ -231,7 +233,8 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){ ...@@ -231,7 +233,8 @@ void CellComplex::insert_cells(bool subdomain, bool boundary){
} }
} }
CellComplex::~CellComplex(){ CellComplex::~CellComplex()
{
for(int i = 0; i < 4; i++){ for(int i = 0; i < 4; i++){
_cells[i].clear(); _cells[i].clear();
...@@ -247,7 +250,8 @@ CellComplex::~CellComplex(){ ...@@ -247,7 +250,8 @@ CellComplex::~CellComplex(){
} }
/* /*
CellComplex::CellComplex(CellComplex* cellComplex){ CellComplex::CellComplex(CellComplex* cellComplex)
{
_domain = cellComplex->getDomain(); _domain = cellComplex->getDomain();
_subdomain = cellComplex->getSubdomain(); _subdomain = cellComplex->getSubdomain();
_boundary = cellComplex->getBoundary(); _boundary = cellComplex->getBoundary();
...@@ -258,16 +262,17 @@ CellComplex::CellComplex(CellComplex* cellComplex){ ...@@ -258,16 +262,17 @@ CellComplex::CellComplex(CellComplex* cellComplex){
} }
_dim = cellComplex->getDim(); _dim = cellComplex->getDim();
restoreComplex(); restoreComplex();
} }*/
*/
void CellComplex::insertCell(Cell* cell){ void CellComplex::insertCell(Cell* cell)
{
_newcells.push_back(cell); _newcells.push_back(cell);
std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell); std::pair<citer, bool> insertInfo = _cells[cell->getDim()].insert(cell);
if(!insertInfo.second) Msg::Debug("Warning: Cell not inserted! \n"); if(!insertInfo.second) Msg::Debug("Warning: Cell not inserted! \n");
} }
void CellComplex::removeCell(Cell* cell, bool other){ void CellComplex::removeCell(Cell* cell, bool other)
{
std::map<Cell*, int, Less_Cell > coboundary; std::map<Cell*, int, Less_Cell > coboundary;
cell->getCoboundary(coboundary); cell->getCoboundary(coboundary);
std::map<Cell*, int, Less_Cell > boundary; std::map<Cell*, int, Less_Cell > boundary;
...@@ -284,31 +289,31 @@ void CellComplex::removeCell(Cell* cell, bool other){ ...@@ -284,31 +289,31 @@ void CellComplex::removeCell(Cell* cell, bool other){
} }
_cells[cell->getDim()].erase(cell); _cells[cell->getDim()].erase(cell);
} }
void CellComplex::removeCellQset(Cell* cell, void CellComplex::removeCellQset(Cell* cell,
std::set<Cell*, Less_Cell>& Qset){ std::set<Cell*, Less_Cell>& Qset)
{
Qset.erase(cell); Qset.erase(cell);
} }
void CellComplex::enqueueCells(std::map<Cell*, int, Less_Cell>& cells, void CellComplex::enqueueCells(std::map<Cell*, int, Less_Cell>& cells,
std::queue<Cell*>& Q, std::queue<Cell*>& Q,
std::set<Cell*, Less_Cell>& Qset){ std::set<Cell*, Less_Cell>& Qset)
{
for(std::map<Cell*, int, Less_Cell>::iterator cit = cells.begin(); for(std::map<Cell*, int, Less_Cell>::iterator cit = cells.begin();
cit != cells.end(); cit++){ cit != cells.end(); cit++){
Cell* cell = (*cit).first; Cell* cell = (*cit).first;
citer it = Qset.find(cell); citer it = Qset.find(cell);
//citer it2 = _cells[cell->getDim()].find(cell); if(it == Qset.end()){
if(it == Qset.end()){// && it2 != _cells[cell->getDim()].end()){
Qset.insert(cell); Qset.insert(cell);
Q.push(cell); Q.push(cell);
} }
} }
} }
int CellComplex::coreduction(Cell* generator, int omitted){ int CellComplex::coreduction(Cell* generator, int omitted)
{
int coreductions = 0; int coreductions = 0;
std::queue<Cell*> Q; std::queue<Cell*> Q;
...@@ -346,7 +351,8 @@ int CellComplex::coreduction(Cell* generator, int omitted){ ...@@ -346,7 +351,8 @@ int CellComplex::coreduction(Cell* generator, int omitted){
return coreductions; return coreductions;
} }
int CellComplex::reduction(int dim, int omitted){ int CellComplex::reduction(int dim, int omitted)
{
if(dim < 1 || dim > 3) return 0; if(dim < 1 || dim > 3) return 0;
int count = 0; int count = 0;
...@@ -378,7 +384,8 @@ int CellComplex::reduction(int dim, int omitted){ ...@@ -378,7 +384,8 @@ int CellComplex::reduction(int dim, int omitted){
return count; return count;
} }
int CellComplex::coreduction(int dim, int omitted){ int CellComplex::coreduction(int dim, int omitted)
{
if(dim < 1 || dim > 3) return 0; if(dim < 1 || dim > 3) return 0;
int count = 0; int count = 0;
...@@ -410,8 +417,8 @@ int CellComplex::coreduction(int dim, int omitted){ ...@@ -410,8 +417,8 @@ int CellComplex::coreduction(int dim, int omitted){
return count; return count;
} }
int CellComplex::reduceComplex(){ int CellComplex::reduceComplex()
{
double t1 = Cpu(); double t1 = Cpu();
Msg::Debug("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n", Msg::Debug("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
...@@ -498,8 +505,8 @@ int CellComplex::coreduceComplex() ...@@ -498,8 +505,8 @@ int CellComplex::coreduceComplex()
return omitted; return omitted;
} }
void CellComplex::computeBettiNumbers(){ void CellComplex::computeBettiNumbers()
{
removeSubdomain(); removeSubdomain();
for(int i = 0; i < 4; i++){ for(int i = 0; i < 4; i++){
...@@ -526,8 +533,8 @@ void CellComplex::computeBettiNumbers(){ ...@@ -526,8 +533,8 @@ void CellComplex::computeBettiNumbers(){
} }
int CellComplex::cocombine(int dim){ int CellComplex::cocombine(int dim)
{
double t1 = Cpu(); double t1 = Cpu();
Msg::Debug("Cell complex before cocombining: %d volumes, %d faces, %d edges and %d vertices.\n", Msg::Debug("Cell complex before cocombining: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0)); getSize(3), getSize(2), getSize(1), getSize(0));
...@@ -591,7 +598,8 @@ int CellComplex::cocombine(int dim){ ...@@ -591,7 +598,8 @@ int CellComplex::cocombine(int dim){
return count; return count;
} }
int CellComplex::combine(int dim){ int CellComplex::combine(int dim)
{
double t1 = Cpu(); double t1 = Cpu();
Msg::Debug("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n", Msg::Debug("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n",
getSize(3), getSize(2), getSize(1), getSize(0)); getSize(3), getSize(2), getSize(1), getSize(0));
...@@ -656,7 +664,8 @@ int CellComplex::combine(int dim){ ...@@ -656,7 +664,8 @@ int CellComplex::combine(int dim){
return count; return count;
} }
void CellComplex::printComplex(int dim){ void CellComplex::printComplex(int dim)
{
for (citer cit = firstCell(dim); cit != lastCell(dim); cit++){ for (citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit; Cell* cell = *cit;
cell->printCell(); cell->printCell();
...@@ -666,7 +675,8 @@ void CellComplex::printComplex(int dim){ ...@@ -666,7 +675,8 @@ void CellComplex::printComplex(int dim){
} }
} }
bool CellComplex::checkCoherence(){ bool CellComplex::checkCoherence()
{
bool coherent = true; bool coherent = true;
for(int i = 0; i < 4; i++){ for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
...@@ -727,7 +737,8 @@ bool CellComplex::checkCoherence(){ ...@@ -727,7 +737,8 @@ bool CellComplex::checkCoherence(){
return coherent; return coherent;
} }
void CellComplex::restoreComplex(){ void CellComplex::restoreComplex()
{
for(int i = 0; i < 4; i++){ for(int i = 0; i < 4; i++){
_betti[i] = 0; _betti[i] = 0;
_cells[i] = _ocells[i]; _cells[i] = _ocells[i];
...@@ -739,7 +750,8 @@ void CellComplex::restoreComplex(){ ...@@ -739,7 +750,8 @@ void CellComplex::restoreComplex(){
_store.clear(); _store.clear();
} }
bool CellComplex::hasCell(Cell* cell, bool org){ bool CellComplex::hasCell(Cell* cell, bool org)
{
if(!org){ if(!org){
citer cit = _cells[cell->getDim()].find(cell); citer cit = _cells[cell->getDim()].find(cell);
if( cit == lastCell(cell->getDim()) ) return false; if( cit == lastCell(cell->getDim()) ) return false;
...@@ -752,7 +764,8 @@ bool CellComplex::hasCell(Cell* cell, bool org){ ...@@ -752,7 +764,8 @@ bool CellComplex::hasCell(Cell* cell, bool org){
} }
} }
bool CellComplex::swapSubdomain(){ bool CellComplex::swapSubdomain()
{
if(_multidim) return false; if(_multidim) return false;
for(int i = 0; i < 4; i++){ for(int i = 0; i < 4; i++){
for(citer cit = firstCell(i); cit != lastCell(i); cit++){ for(citer cit = firstCell(i); cit != lastCell(i); cit++){
...@@ -766,7 +779,8 @@ bool CellComplex::swapSubdomain(){ ...@@ -766,7 +779,8 @@ bool CellComplex::swapSubdomain(){
return true; return true;
} }
void CellComplex::makeDualComplex(){ void CellComplex::makeDualComplex()
{
std::set<Cell*, Less_Cell> temp = _cells[0]; std::set<Cell*, Less_Cell> temp = _cells[0];
_cells[0] = _cells[3]; _cells[0] = _cells[3];
_cells[3] = temp; _cells[3] = temp;
......
...@@ -13,8 +13,8 @@ ...@@ -13,8 +13,8 @@
#if defined(HAVE_KBIPACK) #if defined(HAVE_KBIPACK)
ChainComplex::ChainComplex(CellComplex* cellComplex){ ChainComplex::ChainComplex(CellComplex* cellComplex)
{
_dim = cellComplex->getDim(); _dim = cellComplex->getDim();
_cellComplex = cellComplex; _cellComplex = cellComplex;
...@@ -25,13 +25,11 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){ ...@@ -25,13 +25,11 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
_QMatrix[4] = NULL; _QMatrix[4] = NULL;
_Hbasis[4] = NULL; _Hbasis[4] = NULL;
for(int dim = 0; dim < 4; dim++){ for(int dim = 0; dim < 4; dim++){
unsigned int cols = cellComplex->getSize(dim); unsigned int cols = cellComplex->getSize(dim);
unsigned int rows = 0; unsigned int rows = 0;
if(dim > 0) rows = cellComplex->getSize(dim-1); if(dim > 0) rows = cellComplex->getSize(dim-1);
int index = 1; int index = 1;
// ignore subdomain cells // ignore subdomain cells
for(CellComplex::citer cit = cellComplex->firstCell(dim); for(CellComplex::citer cit = cellComplex->firstCell(dim);
...@@ -118,8 +116,8 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){ ...@@ -118,8 +116,8 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
} }
void ChainComplex::KerCod(int dim){ void ChainComplex::KerCod(int dim)
{
if(dim < 0 || dim > 3 || _HMatrix[dim] == NULL) return; if(dim < 0 || dim > 3 || _HMatrix[dim] == NULL) return;
gmp_matrix* HMatrix gmp_matrix* HMatrix
...@@ -167,8 +165,8 @@ void ChainComplex::KerCod(int dim){ ...@@ -167,8 +165,8 @@ void ChainComplex::KerCod(int dim){
} }
//j:B_(k+1)->Z_k //j:B_(k+1)->Z_k
void ChainComplex::Inclusion(int lowDim, int highDim){ void ChainComplex::Inclusion(int lowDim, int highDim)
{
if(getKerHMatrix(lowDim) == NULL if(getKerHMatrix(lowDim) == NULL
|| getCodHMatrix(highDim) == NULL || getCodHMatrix(highDim) == NULL
|| abs(lowDim-highDim) != 1) return; || abs(lowDim-highDim) != 1) return;
...@@ -191,7 +189,6 @@ void ChainComplex::Inclusion(int lowDim, int highDim){ ...@@ -191,7 +189,6 @@ void ChainComplex::Inclusion(int lowDim, int highDim){
cols = gmp_matrix_cols(Zbasis); cols = gmp_matrix_cols(Zbasis);
if(rows < cols) return; if(rows < cols) return;
// A*inv(V) = U*S // A*inv(V) = U*S
gmp_normal_form* normalForm gmp_normal_form* normalForm
= create_gmp_Smith_normal_form(Zbasis, INVERTED, INVERTED); = create_gmp_Smith_normal_form(Zbasis, INVERTED, INVERTED);
...@@ -199,7 +196,6 @@ void ChainComplex::Inclusion(int lowDim, int highDim){ ...@@ -199,7 +196,6 @@ void ChainComplex::Inclusion(int lowDim, int highDim){
mpz_t elem; mpz_t elem;
mpz_init(elem); mpz_init(elem);
for(int i = 1; i <= cols; i++){ for(int i = 1; i <= cols; i++){
gmp_matrix_get_elem(elem, i, i, normalForm->canonical); gmp_matrix_get_elem(elem, i, i, normalForm->canonical);
...@@ -211,8 +207,6 @@ void ChainComplex::Inclusion(int lowDim, int highDim){ ...@@ -211,8 +207,6 @@ void ChainComplex::Inclusion(int lowDim, int highDim){
gmp_matrix_left_mult(normalForm->left, Bbasis); gmp_matrix_left_mult(normalForm->left, Bbasis);
gmp_matrix* LB = copy_gmp_matrix(Bbasis, 1, 1, gmp_matrix* LB = copy_gmp_matrix(Bbasis, 1, 1,
gmp_matrix_cols(Zbasis), gmp_matrix_cols(Zbasis),
gmp_matrix_cols(Bbasis)); gmp_matrix_cols(Bbasis));
...@@ -253,8 +247,8 @@ void ChainComplex::Inclusion(int lowDim, int highDim){ ...@@ -253,8 +247,8 @@ void ChainComplex::Inclusion(int lowDim, int highDim){
} }
void ChainComplex::Quotient(int dim){ void ChainComplex::Quotient(int dim)
{
if(dim < 0 || dim > 4 || _JMatrix[dim] == NULL) return; if(dim < 0 || dim > 4 || _JMatrix[dim] == NULL) return;
gmp_matrix* JMatrix = gmp_matrix* JMatrix =
...@@ -296,13 +290,12 @@ void ChainComplex::Quotient(int dim){ ...@@ -296,13 +290,12 @@ void ChainComplex::Quotient(int dim){
return; return;
} }
void ChainComplex::computeHomology(bool dual){ void ChainComplex::computeHomology(bool dual)
{
int lowDim = 0; int lowDim = 0;
int highDim = 0; int highDim = 0;
int setDim = 0; int setDim = 0;
for(int i=-1; i < 4; i++){ for(int i=-1; i < 4; i++){
if(dual){ if(dual){
......
...@@ -25,10 +25,18 @@ ...@@ -25,10 +25,18 @@
#include "GVertex.h" #include "GVertex.h"
#include "CellComplex.h" #include "CellComplex.h"
#include "mpz.h"
#if defined(HAVE_GMP)
#include "gmp.h"
extern "C" { extern "C" {
#include "gmp_normal_form.h" #include "gmp_normal_form.h"
} }
#else
extern "C" {
#include "mpz.h"
#include "gmp_normal_form.h"
}
#endif
// A class representing a chain complex of a cell complex. // A class representing a chain complex of a cell complex.
// This should only be constructed for a reduced cell complex because of // This should only be constructed for a reduced cell complex because of
......
...@@ -11,11 +11,11 @@ ...@@ -11,11 +11,11 @@
#include "OS.h" #include "OS.h"
#if defined(HAVE_KBIPACK) #if defined(HAVE_KBIPACK)
Homology::Homology(GModel* model, std::vector<int> physicalDomain,
std::vector<int> physicalSubdomain){
Homology::Homology(GModel* model, std::vector<int> physicalDomain,
std::vector<int> physicalSubdomain)
{
_model = model; _model = model;
_domain = physicalDomain; _domain = physicalDomain;
_subdomain = physicalSubdomain; _subdomain = physicalSubdomain;
...@@ -27,7 +27,6 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, ...@@ -27,7 +27,6 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain,
std::map<int, std::vector<GEntity*> > groups[4]; std::map<int, std::vector<GEntity*> > groups[4];
model->getPhysicalGroups(groups); model->getPhysicalGroups(groups);
std::map<int, std::vector<GEntity*> >::iterator it; std::map<int, std::vector<GEntity*> >::iterator it;
std::vector<GEntity*> domainEntities; std::vector<GEntity*> domainEntities;
std::vector<GEntity*> subdomainEntities; std::vector<GEntity*> subdomainEntities;
...@@ -77,6 +76,7 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain, ...@@ -77,6 +76,7 @@ Homology::Homology(GModel* model, std::vector<int> physicalDomain,
_cellComplex->getSize(1), _cellComplex->getSize(0)); _cellComplex->getSize(1), _cellComplex->getSize(0));
} }
Homology::~Homology(){ Homology::~Homology(){
delete _cellComplex; delete _cellComplex;
for(int i = 0; i < 4; i++) { for(int i = 0; i < 4; i++) {
...@@ -88,8 +88,8 @@ Homology::~Homology(){ ...@@ -88,8 +88,8 @@ Homology::~Homology(){
} }
} }
void Homology::findGenerators(std::string fileName){ void Homology::findGenerators(std::string fileName)
{
Msg::Info("Reducing the Cell Complex..."); Msg::Info("Reducing the Cell Complex...");
Msg::StatusBar(1, false, "Reducing..."); Msg::StatusBar(1, false, "Reducing...");
double t1 = Cpu(); double t1 = Cpu();
...@@ -163,8 +163,6 @@ void Homology::findGenerators(std::string fileName){ ...@@ -163,8 +163,6 @@ void Homology::findGenerators(std::string fileName){
_generators[j].push_back(chain); _generators[j].push_back(chain);
} }
} }
} }
createPViews(); createPViews();
...@@ -191,8 +189,8 @@ void Homology::findGenerators(std::string fileName){ ...@@ -191,8 +189,8 @@ void Homology::findGenerators(std::string fileName){
return; return;
} }
void Homology::findDualGenerators(std::string fileName){ void Homology::findDualGenerators(std::string fileName)
{
Msg::Info("Reducing Cell Complex..."); Msg::Info("Reducing Cell Complex...");
Msg::StatusBar(1, false, "Reducing..."); Msg::StatusBar(1, false, "Reducing...");
double t1 = Cpu(); double t1 = Cpu();
...@@ -215,7 +213,6 @@ void Homology::findDualGenerators(std::string fileName){ ...@@ -215,7 +213,6 @@ void Homology::findDualGenerators(std::string fileName){
_cellComplex->getSize(3), _cellComplex->getSize(2), _cellComplex->getSize(3), _cellComplex->getSize(2),
_cellComplex->getSize(1), _cellComplex->getSize(0)); _cellComplex->getSize(1), _cellComplex->getSize(0));
Msg::Info("Computing homology spaces..."); Msg::Info("Computing homology spaces...");
Msg::StatusBar(1, false, "Computing..."); Msg::StatusBar(1, false, "Computing...");
t1 = Cpu(); t1 = Cpu();
...@@ -227,8 +224,6 @@ void Homology::findDualGenerators(std::string fileName){ ...@@ -227,8 +224,6 @@ void Homology::findDualGenerators(std::string fileName){
int dim = _cellComplex->getDim(); int dim = _cellComplex->getDim();
int HRank[4]; int HRank[4];
for(int i = 0; i < 4; i++) HRank[i] = 0; for(int i = 0; i < 4; i++) HRank[i] = 0;
for(int j = 3; j > -1; j--){ for(int j = 3; j > -1; j--){
...@@ -275,7 +270,7 @@ void Homology::findDualGenerators(std::string fileName){ ...@@ -275,7 +270,7 @@ void Homology::findDualGenerators(std::string fileName){
createPViews(); createPViews();
if(fileName != "") writeGeneratorsMSH(fileName); if(fileName != "") writeGeneratorsMSH(fileName);
Msg::Info("Ranks of homology spaces for dual cell complex:"); Msg::Info("Ranks of homology spaces for the dual cell complex:");
Msg::Info("H0* = %d", HRank[0]); Msg::Info("H0* = %d", HRank[0]);
Msg::Info("H1* = %d", HRank[1]); Msg::Info("H1* = %d", HRank[1]);
Msg::Info("H2* = %d", HRank[2]); Msg::Info("H2* = %d", HRank[2]);
...@@ -296,8 +291,8 @@ void Homology::findDualGenerators(std::string fileName){ ...@@ -296,8 +291,8 @@ void Homology::findDualGenerators(std::string fileName){
return; return;
} }
void Homology::computeBettiNumbers(){ void Homology::computeBettiNumbers()
{
Msg::Info("Running coreduction..."); Msg::Info("Running coreduction...");
Msg::StatusBar(1, false, "Computing..."); Msg::StatusBar(1, false, "Computing...");
double t1 = Cpu(); double t1 = Cpu();
...@@ -319,12 +314,14 @@ void Homology::computeBettiNumbers(){ ...@@ -319,12 +314,14 @@ void Homology::computeBettiNumbers(){
return; return;
} }
void Homology::restoreHomology() { void Homology::restoreHomology()
{
_cellComplex->restoreComplex(); _cellComplex->restoreComplex();
for(int i = 0; i < 4; i++) _generators[i].clear(); for(int i = 0; i < 4; i++) _generators[i].clear();
} }
std::string Homology::getDomainString() { std::string Homology::getDomainString()
{
std::string domainString = "({"; std::string domainString = "({";
for(unsigned int i = 0; i < _domain.size(); i++){ for(unsigned int i = 0; i < _domain.size(); i++){
std::string temp = ""; std::string temp = "";
...@@ -354,7 +351,8 @@ std::string Homology::getDomainString() { ...@@ -354,7 +351,8 @@ std::string Homology::getDomainString() {
return domainString; return domainString;
} }
void Homology::createPViews(){ void Homology::createPViews()
{
for(int i = 0; i < 4; i++){ for(int i = 0; i < 4; i++){
for(int j = 0; j < _generators[i].size(); j++){ for(int j = 0; j < _generators[i].size(); j++){
Chain* chain = _generators[i].at(j); Chain* chain = _generators[i].at(j);
...@@ -363,7 +361,8 @@ void Homology::createPViews(){ ...@@ -363,7 +361,8 @@ void Homology::createPViews(){
} }
} }
bool Homology::writeGeneratorsMSH(std::string fileName, bool binary){ bool Homology::writeGeneratorsMSH(std::string fileName, bool binary)
{
if(!_model->writeMSH(fileName, 2.0, binary)) return false; if(!_model->writeMSH(fileName, 2.0, binary)) return false;
Msg::Info("Wrote homology computation results to %s.", fileName.c_str()); Msg::Info("Wrote homology computation results to %s.", fileName.c_str());
Msg::Debug("Wrote homology computation results to %s. \n", Msg::Debug("Wrote homology computation results to %s. \n",
......
...@@ -15,9 +15,9 @@ ...@@ -15,9 +15,9 @@
*/ */
#include "mpz.h" #include "mpz.h"
#include "GmshConfig.h"
// disabled for now #if ! defined(HAVE_GMP)
#if 0
#include "limits.h" #include "limits.h"
...@@ -27,7 +27,6 @@ void overflow() ...@@ -27,7 +27,6 @@ void overflow()
} }
long int addcheck(long int a, long int b){ long int addcheck(long int a, long int b){
long int c = a + b; long int c = a + b;
if (b >= 0){ if (b >= 0){
if (c < a) overflow(); if (c < a) overflow();
......
...@@ -17,10 +17,12 @@ ...@@ -17,10 +17,12 @@
#ifndef __MPZ_H__ #ifndef __MPZ_H__
#define __MPZ_H__ #define __MPZ_H__
#include "GmshConfig.h"
#if defined(HAVE_GMP)
#include "gmp.h" #include "gmp.h"
#else
// disabled for now
#if 0
#include <stdlib.h> #include <stdlib.h>
#include <stdio.h> #include <stdio.h>
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment