diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 2a0c624583cf60e26c53da26973fbd1ab6a4fab9..f9b7572dc208a769608b38765d1d071c18a1c380 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -440,7 +440,7 @@ gmp_matrix* ChainComplex::getBasis(int dim, int basis)
}
void ChainComplex::getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
- int num, int dim, int basis)
+ int num, int dim, int basis, bool deform)
{
gmp_matrix* basisMatrix;
if(basis == 0) basisMatrix = getBasis(dim, 0);
@@ -461,6 +461,9 @@ void ChainComplex::getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
mpz_t elem;
mpz_init(elem);
+ int torsion = 1;
+ if(basis == 3) torsion = getTorsion(dim, num);
+
for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = cit->first;
int index = cit->second;
@@ -473,11 +476,13 @@ void ChainComplex::getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
for(Cell::citer it = subCells.begin(); it != subCells.end(); it++){
Cell* subCell = (*it).first;
int coeff = (*it).second;
- chain[subCell] = coeff*elemi;
+ chain[subCell] = coeff*elemi*torsion;
}
}
}
mpz_clear(elem);
+
+ if(deform) smoothenChain(chain);
}
int ChainComplex::getBasisSize(int dim, int basis)
@@ -506,6 +511,175 @@ int ChainComplex::getTorsion(int dim, int num)
else return _torsion[dim].at(num-1);
}
+bool ChainComplex::deform(std::map<Cell*, int, Less_Cell>& cells,
+ std::map<Cell*, int, Less_Cell>& cellsInChain,
+ std::map<Cell*, int, Less_Cell>& cellsNotInChain)
+{
+ std::vector<int> cc;
+ std::vector<int> bc;
+
+ for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++){
+ Cell* c = (*cit).first;
+ c->setImmune(false);
+ if(!c->inSubdomain()) {
+ int coeff = 0;
+ citer it = cells.find(c);
+ if(it != cells.end()) coeff = it->second;
+ cc.push_back(coeff);
+ bc.push_back((*cit).second);
+ }
+ }
+
+ if(cc.empty() || (getDim() == 2 && cc.size() < 2) ) return false;
+ int inout = cc[0]*bc[0];
+ for(unsigned int i = 0; i < cc.size(); i++){
+ if(cc[i]*bc[i] != inout) return false;
+ }
+
+ for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++){
+ Cell* cell = cit->first;
+ citer it = cells.find(cell);
+ if(it != cells.end()) cells[cell] = 0;
+ }
+
+ int n = 1;
+ for(citer cit = cellsNotInChain.begin(); cit != cellsNotInChain.end();
+ cit++){
+ Cell* cell = (*cit).first;
+ if(n == 2) cell->setImmune(true);
+ else cell->setImmune(false);
+ int coeff = -1*inout*(*cit).second;
+
+ std::pair<citer,bool> insert = cells.insert( std::make_pair( cell, coeff));
+ if(!insert.second && (*insert.first).second == 0){
+ (*insert.first).second = coeff;
+ }
+ else if (!insert.second && (*insert.first).second != 0){
+ printf("Error: invalid chain smoothening add! \n");
+ }
+ n++;
+ }
+ return true;
+}
+
+bool ChainComplex::deformChain(std::map<Cell*, int, Less_Cell>& cells,
+ std::pair<Cell*, int> cell, bool bend)
+{
+ Cell* c1 = cell.first;
+ int dim = c1->getDim();
+ for(citer cit = c1->firstCoboundary(true); cit != c1->lastCoboundary(true);
+ cit++){
+
+ std::map<Cell*, int, Less_Cell> cellsInChain;
+ std::map<Cell*, int, Less_Cell> cellsNotInChain;
+ Cell* c1CbdCell = (*cit).first;
+
+ for(citer cit2 = c1CbdCell->firstBoundary(true);
+ cit2 != c1CbdCell->lastBoundary(true); cit2++){
+ Cell* c1CbdBdCell = (*cit2).first;
+ int coeff = (*cit2).second;
+ if( (cells.find(c1CbdBdCell) != cells.end() && cells[c1CbdBdCell] != 0)
+ || c1CbdBdCell->inSubdomain()){
+ cellsInChain.insert(std::make_pair(c1CbdBdCell, coeff));
+ }
+ else cellsNotInChain.insert(std::make_pair(c1CbdBdCell, coeff));
+ }
+
+ bool next = false;
+
+ for(citer cit2 = cellsInChain.begin(); cit2 != cellsInChain.end(); cit2++){
+ Cell* c = (*cit2).first;
+ int coeff = 0;
+ citer it = cells.find(c);
+ if(it != cells.end()) coeff = it->second;
+ if(c->getImmune()) next = true;
+ if(c->inSubdomain()) bend = false;
+ if(coeff > 1 || coeff < -1) next = true;
+ }
+
+ for(citer cit2 = cellsNotInChain.begin(); cit2 != cellsNotInChain.end();
+ cit2++){
+ Cell* c = (*cit2).first;
+ if(c->inSubdomain()) next = true;
+ }
+ if(next) continue;
+
+ if( (dim == 1 && cellsInChain.size() == 2
+ && cellsNotInChain.size() == 1) ||
+ (dim == 2 && cellsInChain.size() == 3
+ && cellsNotInChain.size() == 1)){
+ //printf("straighten \n");
+ return deform(cells, cellsInChain, cellsNotInChain);
+ }
+ else if ( (dim == 1 && cellsInChain.size() == 1
+ && cellsNotInChain.size() == 2 && bend) ||
+ (dim == 2 && cellsInChain.size() == 2
+ && cellsNotInChain.size() == 2 && bend)){
+ //printf("bend \n");
+ return deform(cells, cellsInChain, cellsNotInChain);
+ }
+ else if ((dim == 1 && cellsInChain.size() == 3
+ && cellsNotInChain.size() == 0) ||
+ (dim == 2 && cellsInChain.size() == 4
+ && cellsNotInChain.size() == 0)){
+ //printf("remove boundary \n");
+ return deform(cells, cellsInChain, cellsNotInChain);
+ }
+ }
+
+ return false;
+}
+
+void ChainComplex::smoothenChain(std::map<Cell*, int, Less_Cell>& cells)
+{
+ if(!_cellComplex->simplicial() || cells.empty()) return;
+ int dim = cells.begin()->first->getDim();
+ int start = cells.size();
+
+ int useless = 0;
+ for(int i = 0; i < 20; i++){
+ int size = cells.size();
+ for(citer cit = cells.begin(); cit != cells.end(); cit++){
+ //if(!deformChain(*cit, false) && getDim() == 2) deformChain(*cit, true);
+ if(dim == 2) deformChain(cells, *cit, true);
+ deformChain(cells, *cit, false);
+
+ }
+
+ deImmuneCells(cells);
+ eraseNullCells(cells);
+
+ if (size >= cells.size()) useless++;
+ else useless = 0;
+ if (useless > 5) break;
+ }
+
+ deImmuneCells(cells);
+ for(citer cit = cells.begin(); cit != cells.end(); cit++){
+ deformChain(cells, *cit, false);
+ }
+ eraseNullCells(cells);
+ int size = cells.size();
+ printf("Smoothened a %d-chain from %d cells to %d cells.\n",
+ dim, start, size);
+}
+
+void ChainComplex::eraseNullCells(std::map<Cell*, int, Less_Cell>& cells)
+{
+ std::vector<Cell*> toRemove;
+ for(citer cit = cells.begin(); cit != cells.end(); cit++){
+ if(cit->second == 0) toRemove.push_back(cit->first);
+ }
+ for(unsigned int i = 0; i < toRemove.size(); i++) cells.erase(toRemove[i]);
+}
+
+void ChainComplex::deImmuneCells(std::map<Cell*, int, Less_Cell>& cells)
+{
+ for(citer cit = cells.begin(); cit != cells.end(); cit++){
+ Cell* cell = (*cit).first;
+ cell->setImmune(false);
+ }
+}
HomologySequence::HomologySequence(ChainComplex* subcomplex,
ChainComplex* complex,
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index dcf9b2fe3923202566cb2eb076864be3166aa186..ea1c7fe210fd4fac49cdee7382f0a4f4b4e0b314 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -92,7 +92,16 @@ class ChainComplex
gmp_matrix* getHbasis(int dim) const {
if(dim > -1 && dim < 5) return _Hbasis[dim]; else return NULL;}
-
+ // local deformation tools for chains
+ bool deformChain(std::map<Cell*, int, Less_Cell>& cells,
+ std::pair<Cell*, int> cell, bool bend);
+ bool deform(std::map<Cell*, int, Less_Cell>& cells,
+ std::map<Cell*, int, Less_Cell>& cellsInChain,
+ std::map<Cell*, int, Less_Cell>& cellsNotInChain);
+ void smoothenChain(std::map<Cell*, int, Less_Cell>& cells);
+ void eraseNullCells(std::map<Cell*, int, Less_Cell>& cells);
+ void deImmuneCells(std::map<Cell*, int, Less_Cell>& cells);
+
public:
// domain = 0 : relative chain space
// domain = 1 : absolute chain space of all cells in cellComplex
@@ -147,8 +156,11 @@ class ChainComplex
// get coefficient vector for dim-dimensional Hbasis chain chainNumber
std::vector<int> getCoeffVector(int dim, int chainNumber);
// get basis chain from a basis matrix
+ // (deform: with local deformations to make chain smoother and to have
+ // smaller support, deformed chain is homologous to the old one,
+ // only works for chains of the primary chain complex)
void getBasisChain(std::map<Cell*, int, Less_Cell>& chain,
- int num, int dim, int basis);
+ int num, int dim, int basis, bool deform=false);
// get rank of a basis
int getBasisSize(int dim, int basis);
// homology torsion coefficient for dim-dimensional chain num
@@ -160,7 +172,7 @@ class ChainComplex
gmp_matrix_right_mult(_Hbasis[dim], T);
}
}
-
+
// debugging aid
int printMatrix(gmp_matrix* matrix){
printf("%d rows and %d columns\n",
@@ -213,7 +225,7 @@ class HomologySequence
printf("%d rows and %d columns\n",
(int)gmp_matrix_rows(matrix), (int)gmp_matrix_cols(matrix));
return gmp_matrix_printf(matrix); }
-
+
};
#endif
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index eef7aca4cb1be41eb8711ba1416128cab67b88e4..3603054fc02ca7ba611f4128d79d79c76d56db4c 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -152,26 +152,17 @@ void Homology::findGenerators(CellComplex* cellComplex)
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));
- /* std::set<Cell*, Less_Cell> cells;
- cellComplex->getCells(cells, j);
- Chain* chain = new Chain(cells,
- chains->getCoeffVector(j,i), cellComplex,
- _model, name, chains->getTorsion(j,i));*/
- t1 = Cpu();
- int start = chain->getSize();
- chain->smoothenChain();
- t2 = Cpu();
- Msg::Info("Smoothened H%d %d from %d cells to %d cells (%g s).",
- j, i, start, chain->getSize(), t2 - t1);
- if(chain->getSize() != 0) {
- HRank[j] = HRank[j] + 1;
- if(chain->getTorsion() != 1){
- Msg::Warning("H%d %d has torsion coefficient %d!",
- j, i, chain->getTorsion());
- }
+ if(chain->getSize() == 0) {
+ delete chain;
+ continue;
+ }
+ HRank[j] = HRank[j] + 1;
+ if(chain->getTorsion() != 1){
+ Msg::Warning("H%d %d has torsion coefficient %d!",
+ j, i, chain->getTorsion());
}
_generators.push_back(chain->createPGroup());
delete chain;
@@ -184,7 +175,11 @@ void Homology::findGenerators(CellComplex* cellComplex)
std::vector<int> coeffs (cellComplex->getOmitted(i).size(),1);
Chain* chain = new Chain(cellComplex->getOmitted(i), coeffs,
cellComplex, _model, name, 1);
- if(chain->getSize() != 0) HRank[j] = HRank[j] + 1;
+ if(chain->getSize() == 0){
+ delete chain;
+ continue;
+ }
+ HRank[j] = HRank[j] + 1;
_generators.push_back(chain->createPGroup());
delete chain;
}
@@ -260,23 +255,18 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
chains->getBasisChain(protoChain, i, j, 3);
Chain* chain = new Chain(protoChain, cellComplex, _model, name,
chains->getTorsion(j,i));
- /*std::set<Cell*, Less_Cell> cells;
- cellComplex->getCells(cells, j);
- Chain* chain = new Chain(cells,
- chains->getCoeffVector(j,i), cellComplex,
- _model, name, chains->getTorsion(j,i));*/
+ if(chain->getSize() == 0) {
+ delete chain;
+ continue;
+ }
_generators.push_back(chain->createPGroup());
- delete chain;
- if(chain->getSize() != 0){
- HRank[dim-j] = HRank[dim-j] + 1;
- if(chain->getTorsion() != 1){
- Msg::Warning("H%d* %d has torsion coefficient %d!",
- dim-j, i, chain->getTorsion());
- }
+ HRank[dim-j] = HRank[dim-j] + 1;
+ if(chain->getTorsion() != 1){
+ Msg::Warning("H%d* %d has torsion coefficient %d!",
+ dim-j, i, chain->getTorsion());
}
}
-
if(j == 0 && cellComplex->getNumOmitted() > 0){
for(int i = 0; i < cellComplex->getNumOmitted(); i++){
std::string generator;
@@ -287,7 +277,11 @@ void Homology::findDualGenerators(CellComplex* cellComplex)
std::vector<int> coeffs (cellComplex->getOmitted(i).size(),1);
Chain* chain = new Chain(cellComplex->getOmitted(i), coeffs,
cellComplex, _model, name, 1);
- if(chain->getSize() != 0) HRank[dim-j] = HRank[dim-j] + 1;
+ if(chain->getSize() == 0) {
+ delete chain;
+ continue;
+ }
+ HRank[dim-j] = HRank[dim-j] + 1;
_generators.push_back(chain->createPGroup());
delete chain;
}
@@ -383,13 +377,10 @@ void Homology::findHomSequence(){
chains->getBasisChain(protoChain, i, j, 3);
Chain* chain = new Chain(protoChain, cellComplex, _model, name,
chains->getTorsion(j,i));
- if(chain->getSize() == 0) continue;
- t1 = Cpu();
- int start = chain->getSize();
- //chain->smoothenChain();
- t2 = Cpu();
- //Msg::Info("Smoothened H%d %d from %d cells to %d cells (%g s).",
- // j, i, start, chain->getSize(), t2 - t1);
+ if(chain->getSize() == 0) {
+ delete chain;
+ continue;
+ }
HRank[j] = HRank[j] + 1;
if(chain->getTorsion() != 1){
Msg::Warning("H%d %d has torsion coefficient %d!",
@@ -547,144 +538,6 @@ Chain::Chain(std::map<Cell*, int, Less_Cell>& chain,
eraseNullCells();
}
-bool Chain::deform(std::map<Cell*, int, Less_Cell>& cellsInChain,
- std::map<Cell*, int, Less_Cell>& cellsNotInChain)
-{
- std::vector<int> cc;
- std::vector<int> bc;
-
- for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++){
- Cell* c = (*cit).first;
- c->setImmune(false);
- if(!c->inSubdomain()) {
- cc.push_back(getCoeff(c));
- bc.push_back((*cit).second);
- }
- }
-
- if(cc.empty() || (getDim() == 2 && cc.size() < 2) ) return false;
- int inout = cc[0]*bc[0];
- for(unsigned int i = 0; i < cc.size(); i++){
- if(cc[i]*bc[i] != inout) return false;
- }
-
- for(citer cit = cellsInChain.begin(); cit != cellsInChain.end(); cit++){
- removeCell((*cit).first);
- }
-
- int n = 1;
- for(citer cit = cellsNotInChain.begin(); cit != cellsNotInChain.end();
- cit++){
- Cell* c = (*cit).first;
- if(n == 2) c->setImmune(true);
- else c->setImmune(false);
- int coeff = -1*inout*(*cit).second;
- addCell(c, coeff);
- n++;
- }
-
- return true;
-}
-
-bool Chain::deformChain(std::pair<Cell*, int> cell, bool bend)
-{
- Cell* c1 = cell.first;
- for(citer cit = c1->firstCoboundary(true); cit != c1->lastCoboundary(true);
- cit++){
-
- std::map<Cell*, int, Less_Cell> cellsInChain;
- std::map<Cell*, int, Less_Cell> cellsNotInChain;
- Cell* c1CbdCell = (*cit).first;
-
- for(citer cit2 = c1CbdCell->firstBoundary(true);
- cit2 != c1CbdCell->lastBoundary(true); cit2++){
- Cell* c1CbdBdCell = (*cit2).first;
- int coeff = (*cit2).second;
- if( (hasCell(c1CbdBdCell) && getCoeff(c1CbdBdCell) != 0)
- || c1CbdBdCell->inSubdomain()){
- cellsInChain.insert(std::make_pair(c1CbdBdCell, coeff));
- }
- else cellsNotInChain.insert(std::make_pair(c1CbdBdCell, coeff));
- }
-
- bool next = false;
-
- for(citer cit2 = cellsInChain.begin(); cit2 != cellsInChain.end(); cit2++){
- Cell* c = (*cit2).first;
- int coeff = getCoeff(c);
- if(c->getImmune()) next = true;
- if(c->inSubdomain()) bend = false;
- if(coeff > 1 || coeff < -1) next = true;
- }
-
- for(citer cit2 = cellsNotInChain.begin(); cit2 != cellsNotInChain.end();
- cit2++){
- Cell* c = (*cit2).first;
- if(c->inSubdomain()) next = true;
- }
- if(next) continue;
-
- if( (getDim() == 1 && cellsInChain.size() == 2
- && cellsNotInChain.size() == 1) ||
- (getDim() == 2 && cellsInChain.size() == 3
- && cellsNotInChain.size() == 1)){
- //printf("straighten \n");
- return deform(cellsInChain, cellsNotInChain);
- }
- else if ( (getDim() == 1 && cellsInChain.size() == 1
- && cellsNotInChain.size() == 2 && bend) ||
- (getDim() == 2 && cellsInChain.size() == 2
- && cellsNotInChain.size() == 2 && bend)){
- //printf("bend \n");
- return deform(cellsInChain, cellsNotInChain);
- }
- else if ((getDim() == 1 && cellsInChain.size() == 3
- && cellsNotInChain.size() == 0) ||
- (getDim() == 2 && cellsInChain.size() == 4
- && cellsNotInChain.size() == 0)){
- //printf("remove boundary \n");
- return deform(cellsInChain, cellsNotInChain);
- }
- }
-
- return false;
-}
-
-void Chain::smoothenChain()
-{
- if(!_cellComplex->simplicial()) return;
-
- int start = getSize();
- double t1 = Cpu();
-
- int useless = 0;
- for(int i = 0; i < 20; i++){
- int size = getSize();
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- //if(!deformChain(*cit, false) && getDim() == 2) deformChain(*cit, true);
- if(getDim() == 2) deformChain(*cit, true);
- deformChain(*cit, false);
- }
- deImmuneCells();
- eraseNullCells();
- if (size >= getSize()) useless++;
- else useless = 0;
- if (useless > 5) break;
- }
-
- deImmuneCells();
- for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
- deformChain(*cit, false);
- }
-
- eraseNullCells();
- double t2 = Cpu();
- Msg::Debug("Smoothened a %d-chain from %d cells to %d cells (%g s).\n",
- getDim(), start, getSize(), t2-t1);
- return;
-}
-
-
int Chain::writeChainMSH(const std::string &name)
{
if(getSize() == 0) return 1;
@@ -831,10 +684,8 @@ int Chain::getCoeff(Cell* c)
void Chain::eraseNullCells()
{
std::vector<Cell*> toRemove;
- for(int i = 0; i < 4; i++){
- for(citer cit = _cells.begin(); cit != _cells.end(); ++cit){
- if((*cit).second == 0) toRemove.push_back((*cit).first);
- }
+ for(citer cit = _cells.begin(); cit != _cells.end(); cit++){
+ if(cit->second == 0) toRemove.push_back(cit->first);
}
for(unsigned int i = 0; i < toRemove.size(); i++) _cells.erase(toRemove[i]);
return;
diff --git a/Geo/Homology.h b/Geo/Homology.h
index 2f7298ea0cef021627cee8d4f0dd3c89e1045288..24d8963f0d53add3954a8bfc63c4b1b730f79395 100644
--- a/Geo/Homology.h
+++ b/Geo/Homology.h
@@ -105,11 +105,6 @@ class Chain{
int _dim;
- bool deform(std::map<Cell*, int, Less_Cell> &cellsInChain,
- std::map<Cell*, int, Less_Cell> &cellsNotInChain);
- bool deformChain(std::pair<Cell*, int> cell, bool bend);
-
-
public:
Chain() {}
Chain(std::set<Cell*, Less_Cell> cells, std::vector<int> coeffs,
@@ -166,11 +161,6 @@ class Chain{
int getNum() const { return _num; }
void setNum(int num) { _num=num; }
- // make local deformations to the chain to make it smoother and smaller
- // (for primary complex chains only, not for dual chains represented
- // by primary cells (yet).)
- void smoothenChain();
-
// append this chain to a MSH ASCII file as $ElementData
// for debugging only
int writeChainMSH(const std::string &name);
diff --git a/benchmarks/homology/torsion.geo b/benchmarks/homology/torsion.geo
new file mode 100644
index 0000000000000000000000000000000000000000..3fe925778b6dadb47c15056175d16b25588e3d0e
--- /dev/null
+++ b/benchmarks/homology/torsion.geo
@@ -0,0 +1,136 @@
+// Gmsh
+
+m=0.3;
+
+Point(1) = {0, 0, 0, m};
+Point(2) = {10, 0, 0, m};
+Point(3) = {10, 10, 0, m};
+Point(4) = {0, 10, 0, m};
+Point(5) = {4, 4, 0, m};
+Point(6) = {6, 4, 0, m};
+Point(7) = {6, 6, 0, m};
+Point(8) = {4, 6, 0, m};
+Point(newp) = {0, 0, 1, m};
+Point(newp) = {10, 0, 1, m};
+Point(newp) = {10, 10, 1, m};
+Point(newp) = {0, 10, 1, m};
+Point(newp) = {4, 4, 1, m};
+Point(newp) = {6, 4, 1, m};
+Point(newp) = {6, 6, 1, m};
+Point(newp) = {4, 6, 1, m};
+
+Point(newp) = {0, 1, 0, m};
+Point(newp) = {7, 1, 0, m};
+Point(newp) = {7, 0, 0, m};
+
+Line(5) = {17, 4};
+Line(6) = {4, 12};
+Line(7) = {9, 12};
+
+Point(newp) = {9, 10, 1, m};
+Point(newp) = {10, 3, 1, m};
+Point(newp) = {9, 3, 1, m};
+Point(newp) = {9, 4, 1, m};
+
+
+Point(newp) = {5, 4, 1, m};
+Point(newp) = {5, 4, 1, m};
+Point(newp) = {4, 3, 1, m};
+
+Line(8) = {20, 11};
+Line(9) = {11, 3};
+Line(10) = {20, 12};
+Line(11) = {4, 3};
+Line(12) = {9, 1};
+Line(13) = {20, 23};
+Line(14) = {11, 21};
+Line(15) = {22, 21};
+Line(16) = {23, 22};
+Line(17) = {14, 23};
+Line(18) = {14, 24};
+
+Line(19) = {24, 13};
+Line(20) = {13, 26};
+Line(21) = {26, 22};
+
+Point(newp) = {6, 5, 1, m};
+Point(newp) = {6, 5, 0, m};
+Point(newp) = {5, 4, 0, m};
+
+
+Line(22) = {8, 7};
+Line(23) = {7, 15};
+Line(24) = {15, 16};
+Line(25) = {16, 8};
+Line(26) = {13, 16};
+Line(27) = {5, 8};
+Line(28) = {5, 13};
+Line(29) = {28, 27};
+Line(30) = {28, 7};
+Line(31) = {15, 27};
+Line(32) = {5, 29};
+Line(33) = {29, 24};
+
+Point(newp) = {6, 1, 0, m};
+Line(34) = {17, 1};
+Line(35) = {1, 19};
+Line(36) = {19, 18};
+Line(37) = {18, 30};
+Line(38) = {30, 17};
+Line(39) = {28, 6};
+
+Point(newp) = {7, 6, 0, m};
+Line(40) = {7, 31};
+Line(41) = {31, 18};
+Line(42) = {6, 30};
+
+
+
+Line Loop(45) = {35, 36, 37, 38, 34};
+Plane Surface(46) = {45};
+Line Loop(47) = {41, 37, -42, -39, 30, 40};
+Plane Surface(48) = {47};
+Line Loop(49) = {34, -12, 7, -6, -5};
+Plane Surface(50) = {49};
+Line Loop(51) = {11, -9, -8, 10, -6};
+Plane Surface(52) = {51};
+Line Loop(53) = {13, 16, 15, -14, -8};
+Plane Surface(54) = {53};
+Line Loop(55) = {21, -16, -17, 18, 19, 20};
+Plane Surface(56) = {55};
+Line Loop(57) = {33, 19, -28, 32};
+Plane Surface(58) = {57};
+Line Loop(59) = {26, 25, -27, 28};
+Plane Surface(60) = {59};
+Line Loop(61) = {24, 25, 22, 23};
+Plane Surface(62) = {61};
+Line Loop(63) = {30, 23, 31, -29};
+Plane Surface(64) = {63};
+Line(65) = {2, 10};
+Line(66) = {2, 3};
+Line(67) = {21, 10};
+Line(68) = {2, 19};
+Line(69) = {10, 9};
+Line(70) = {27, 14};
+Line(71) = {14, 6};
+Line(72) = {6, 29};
+Line Loop(73) = {66, -9, 14, 67, -65};
+Plane Surface(74) = {73};
+Line Loop(75) = {69, 12, 35, -68, 65};
+Plane Surface(76) = {75};
+Line Loop(77) = {67, 69, 7, -10, 13, -17, -70, -31, 24, -26, 20, 21, 15};
+Plane Surface(78) = {77};
+Line Loop(79) = {68, 36, -41, -40, -22, -27, 32, -72, 42, 38, 5, 11, -66};
+Plane Surface(80) = {79};
+Line Loop(81) = {72, 33, -18, 71};
+Plane Surface(82) = {81};
+Line Loop(83) = {29, 70, 71, -39};
+Plane Surface(84) = {83};
+Surface Loop(85) = {76, 78, 74, 80, 46, 48, 84, 64, 62, 60, 58, 82, 56, 54, 52, 50};
+Volume(86) = {85};
+Physical Volume(87) = {86};
+Physical Surface(88) = {46, 48, 64, 62, 60, 58, 56, 54, 52, 50};
+
+Mesh 3;
+
+HomGen("torsion.msh") = {{87},{88}};
diff --git a/benchmarks/homology/trefoil.geo b/benchmarks/homology/trefoil.geo
new file mode 100644
index 0000000000000000000000000000000000000000..52aaf789b479a03014dc9ebaa840dee4e06737da
--- /dev/null
+++ b/benchmarks/homology/trefoil.geo
@@ -0,0 +1,285 @@
+// Gmsh
+
+m=1;
+
+Point(newp) = {-2, 0, 0, m};
+Point(newp) = {10, 0, 0, m};
+Point(newp) = {10, 10, 0, m};
+Point(newp) = {-2, 10, 0, m};
+Point(newp) = {-2, 0, 2, m};
+Point(newp) = {10, 0, 2, m};
+Point(newp) = {10, 10, 2, m};
+Point(newp) = {-2, 10, 2, m};
+
+Point(newp) = {0, 2, 0, m};
+Point(newp) = {8, 2, 0, m};
+Point(newp) = {8, 8, 0, m};
+Point(newp) = {0, 8, 0, m};
+Point(newp) = {0, 2, 2, m};
+Point(newp) = {8, 2, 2, m};
+Point(newp) = {8, 8, 2, m};
+Point(newp) = {0, 8, 2, m};
+Line(1) = {4, 3};
+Line(2) = {7, 8};
+Line(3) = {8, 4};
+Line(4) = {3, 7};
+Line(5) = {11, 12};
+Line(6) = {12, 16};
+Line(7) = {16, 15};
+Line(8) = {15, 11};
+Line(9) = {12, 9};
+Line(10) = {13, 9};
+Line(11) = {13, 16};
+Line(12) = {4, 1};
+Line(13) = {1, 5};
+Line(14) = {5, 8};
+
+Point(newp) = {4, 2, 0, m};
+Point(newp) = {6, 2, 0, m};
+Point(newp) = {4, 0, 0, m};
+Point(newp) = {6, 0, 0, m};
+
+Point(newp) = {4, 2, 2, m};
+Point(newp) = {6, 2, 2, m};
+Point(newp) = {4, 0, 2, m};
+Point(newp) = {6, 0, 2, m};
+
+Point(newp) = {4, 2, -8, m};
+Point(newp) = {6, 2, -8, m};
+Point(newp) = {4, 0, -10, m};
+Point(newp) = {6, 0, -10, m};
+
+Line(15) = {9, 17};
+Line(16) = {13, 21};
+Line(17) = {17, 25};
+Line(18) = {21, 22};
+Line(19) = {22, 18};
+Line(20) = {18, 26};
+Line(21) = {5, 23};
+Line(22) = {23, 24};
+Line(23) = {24, 20};
+Line(24) = {20, 28};
+Line(25) = {1, 19};
+Line(26) = {19, 27};
+
+Point(newp) = {4, 22, -10, m};
+Point(newp) = {6, 22, -10, m};
+Point(newp) = {4, 20, -8, m};
+Point(newp) = {6, 20, -8, m};
+
+Point(newp) = {4, 20, 8, m};
+Point(newp) = {6, 20, 8, m};
+Point(newp) = {4, 22, 10, m};
+Point(newp) = {6, 22, 10, m};
+Line(27) = {25, 31};
+Line(28) = {32, 26};
+Line(29) = {29, 27};
+Line(30) = {30, 28};
+Line(31) = {31, 33};
+Line(32) = {34, 32};
+Line(33) = {29, 35};
+Line(34) = {36, 30};
+
+Point(newp) = {4, 6, 8, m};
+Point(newp) = {6, 4, 8, m};
+Point(newp) = {4, 6, 10, m};
+Point(newp) = {6, 4, 10, m};
+Line(35) = {34, 38};
+Line(36) = {40, 36};
+Line(37) = {33, 37};
+Line(38) = {35, 39};
+
+Point(newp) = {1, 6, 8, m};
+Point(newp) = {3, 4, 8, m};
+Point(newp) = {1, 6, 10, m};
+Point(newp) = {3, 4, 10, m};
+Point(newp) = {1, 4, 8, m};
+Point(newp) = {1, 4, 10, m};
+
+Point(newp) = {3, 6, 8, m};
+Point(newp) = {3, 6, 10, m};
+
+Line(39) = {37, 47};
+Line(40) = {47, 41};
+Line(41) = {39, 48};
+Line(42) = {48, 43};
+Line(43) = {38, 42};
+Line(44) = {40, 44};
+Line(45) = {44, 46};
+Line(46) = {42, 45};
+Line(47) = {25, 26};
+Line(48) = {28, 27};
+Line(49) = {29, 30};
+Line(50) = {31, 32};
+Line(51) = {33, 34};
+Line(52) = {35, 36};
+
+Point(newp) = {3, 6, -18, m};
+Point(newp) = {1, 6, -18, m};
+Point(newp) = {3, 4, -18, m};
+Point(newp) = {1, 4, -18, m};
+Line(53) = {47, 49};
+Line(54) = {50, 41};
+Line(55) = {45, 41};
+Line(56) = {45, 46};
+Line(57) = {41, 43};
+Line(58) = {45, 52};
+Line(59) = {42, 51};
+Line(60) = {43, 46};
+Line(61) = {38, 40};
+Line(62) = {39, 37};
+
+Point(newp) = {3, 6, -20, m};
+Point(newp) = {1, 6, -20, m};
+Point(newp) = {3, 4, -20, m};
+Point(newp) = {1, 4, -20, m};
+Line(63) = {50, 54};
+Line(64) = {52, 56};
+Line(65) = {56, 55};
+Line(66) = {54, 56};
+Line(67) = {49, 53};
+Line(68) = {54, 53};
+
+Point(newp) = {8, 6, -18, m};
+Point(newp) = {10, 6, -20, m};
+Point(newp) = {8, 4, -18, m};
+Point(newp) = {10, 4, -20, m};
+Line(69) = {51, 59};
+Line(70) = {57, 49};
+Line(71) = {55, 60};
+Line(72) = {58, 53};
+
+Point(newp) = {8, 6, 0, m};
+Point(newp) = {10, 6, 2, m};
+Point(newp) = {8, 4, 0, m};
+Point(newp) = {10, 4, 2, m};
+Point(newp) = {8, 6, 2, m};
+Point(newp) = {10, 6, 0, m};
+Point(newp) = {8, 4, 2, m};
+Point(newp) = {10, 4, 0, m};
+Line(73) = {57, 61};
+Line(74) = {58, 66};
+Line(75) = {59, 63};
+Line(76) = {60, 68};
+Line(77) = {68, 64};
+Line(78) = {64, 62};
+Line(79) = {62, 7};
+Line(80) = {3, 66};
+Line(81) = {64, 67};
+Line(82) = {67, 65};
+Line(83) = {65, 15};
+Line(84) = {11, 61};
+Line(85) = {63, 67};
+Line(86) = {66, 61};
+Line(87) = {58, 60};
+Line(88) = {59, 57};
+Line(89) = {49, 51};
+Delete {
+ Point{2, 10, 14, 6};
+}
+Line(90) = {22, 24};
+Line(91) = {17, 19};
+Line Loop(92) = {1, 80, 86, -84, 5, 9, 15, 91, -25, -12};
+Plane Surface(93) = {92};
+Line Loop(94) = {79, -4, 80, -74, 87, 76, 77, 78};
+Plane Surface(95) = {94};
+Line Loop(96) = {73, -86, -74, 72, -67, -70};
+Plane Surface(97) = {96};
+Line Loop(98) = {71, -87, 72, -68, 66, 65};
+Plane Surface(99) = {98};
+Line Loop(100) = {70, 89, 69, 88};
+Plane Surface(101) = {100};
+Line Loop(102) = {73, -84, -8, -83, -82, -85, -75, 88};
+Plane Surface(103) = {102};
+Line Loop(104) = {76, 77, 81, -85, -75, -69, -59, 46, 58, 64, 65, 71};
+Plane Surface(105) = {104};
+Line Loop(106) = {32, 28, -20, -19, 90, 23, 24, -30, -34, -36, -61, -35};
+Plane Surface(107) = {106};
+Line Loop(108) = {32, -50, 31, 51};
+Plane Surface(109) = {108};
+Line(110) = {42, 47};
+Line Loop(111) = {35, 43, 110, -39, -37, 51};
+Plane Surface(112) = {111};
+Line Loop(113) = {39, 40, 57, -42, -41, 62};
+Plane Surface(114) = {113};
+Line Loop(115) = {36, -52, 38, 41, 42, 60, -45, -44};
+Plane Surface(116) = {115};
+Line Loop(117) = {55, 57, 60, -56};
+Plane Surface(118) = {117};
+Line Loop(119) = {44, 45, -56, -46, -43, 61};
+Plane Surface(120) = {119};
+Line Loop(121) = {2, -14, 21, 22, -90, -18, -16, 11, 7, -83, -82, -81, 78, 79};
+Plane Surface(122) = {121};
+Line Loop(123) = {25, 26, -48, -24, -23, -22, -21, -13};
+Plane Surface(124) = {123};
+Line Loop(125) = {30, 48, -29, 49};
+Plane Surface(126) = {125};
+Line Loop(127) = {53, 67, -68, -63, 54, -40};
+Plane Surface(128) = {127};
+Line Loop(129) = {59, -89, -53, -110};
+Plane Surface(130) = {129};
+Line Loop(131) = {34, -49, 33, 52};
+Plane Surface(132) = {131};
+Line Loop(133) = {27, 31, 37, -62, -38, -33, 29, -26, -91, 17};
+Plane Surface(134) = {133};
+Line Loop(135) = {1, 4, 2, 3};
+Plane Surface(136) = {135};
+Line Loop(137) = {12, 13, 14, 3};
+Plane Surface(138) = {137};
+Line Loop(139) = {27, 50, 28, -47};
+Plane Surface(140) = {139};
+Line Loop(141) = {16, 18, 19, 20, -47, -17, -15, -10};
+Plane Surface(142) = {141};
+Line Loop(143) = {9, -10, 11, -6};
+Plane Surface(144) = {143};
+Line Loop(145) = {58, 64, -66, -63, 54, -55};
+Plane Surface(146) = {145};
+Line Loop(147) = {7, 8, 5, 6};
+Plane Surface(148) = {147};
+Surface Loop(149) = {132, 107, 109, 140, 134, 112, 120, 116, 114, 128, 130, 105, 95, 122, 136, 93, 97, 103, 148, 144, 142, 101, 99, 146, 118, 124, 126, 138};
+Volume(150) = {149};
+Physical Volume(151) = {150};
+Physical Surface(152) = {132, 134, 109, 107, 126, 140, 95, 97, 103, 136, 148, 122, 146, 128, 99, 101, 105, 120, 116, 112, 114, 138, 93, 144, 124, 118, 130, 142};
+
+Point(newp) = {30, -30, 30, 10*m};
+Point(newp) = {-30, 30, 30, 10*m};
+Point(newp) = {30, 30, 30, 10*m};
+Point(newp) = {-30, -30, 30, 10*m};
+Point(newp) = {30, -30, -30, 10*m};
+Point(newp) = {-30, 30, -30, 10*m};
+Point(newp) = {30, 30, -30, 10*m};
+Point(newp) = {-30, -30, -30, 10*m};
+Line(153) = {70, 71};
+Line(154) = {71, 75};
+Line(155) = {75, 74};
+Line(156) = {74, 70};
+Line(157) = {69, 72};
+Line(158) = {72, 76};
+Line(159) = {76, 73};
+Line(160) = {73, 69};
+Line(161) = {69, 71};
+Line(162) = {70, 72};
+Line(163) = {74, 76};
+Line(164) = {75, 73};
+Line Loop(165) = {156, 162, 158, -163};
+Plane Surface(166) = {165};
+Line Loop(167) = {159, 160, 157, 158};
+Plane Surface(168) = {167};
+Line Loop(169) = {155, 163, 159, -164};
+Plane Surface(170) = {169};
+Line Loop(171) = {154, 164, 160, 161};
+Plane Surface(172) = {171};
+Line Loop(173) = {157, -162, 153, -161};
+Plane Surface(174) = {173};
+Line Loop(175) = {156, 153, 154, 155};
+Plane Surface(176) = {175};
+Surface Loop(177) = {166, 176, 174, 168, 170, 172};
+Volume(178) = {177, 149};
+Physical Volume(179) = {178};
+Physical Surface(180) = {166, 176, 174, 172, 170, 168, 93, 122, 128, 124, 140, 126, 109, 132, 134, 107, 112, 116, 136, 148, 130, 146, 142, 144, 118, 114, 120, 138, 99, 105, 95, 97, 103, 101};
+
+Mesh 3;
+
+HomGen("trefoil.msh") = {{179},{}};
+HomGen("trefoil.msh") = {{179},{180}};
+//HomCut("trefoil.msh") = {{179},{}};
diff --git a/benchmarks/homology/wirewound.geo b/benchmarks/homology/wirewound.geo
new file mode 100644
index 0000000000000000000000000000000000000000..6450e707a3932b4b8bdb831b00014123a33cac40
--- /dev/null
+++ b/benchmarks/homology/wirewound.geo
@@ -0,0 +1,246 @@
+// Gmsh
+m = 2;
+k = 1;
+Point(1) = {0, 0, 0, m};
+Point(2) = {1, 0, 0, m};
+Point(3) = {1, 1, 0, m};
+Point(4) = {0, 1, 0, m};
+Line(1) = {2, 3};
+Line(2) = {3, 4};
+Line(3) = {4, 1};
+Line(4) = {1, 2};
+Coherence;
+Line Loop(5) = {4, 1, 2, 3};
+Plane Surface(6) = {5};
+Extrude {0, 0, 5} {
+ Surface{6};
+}
+Point(15) = {-0.1, 1.1, 1};
+Point(16) = {-0.1, 1.2, 1};
+Point(17) = {-0.1, 1.2, 1.1};
+Point(18) = {-0.1, 1.1, 1.1};
+Line(29) = {16, 17};
+Line(30) = {18, 17};
+Line(31) = {18, 15};
+Line(32) = {16, 15};
+Line Loop(33) = {29, -30, 31, -32};
+Plane Surface(34) = {33};
+Extrude {1.2, 0, 0} {
+ Surface{34};
+}
+Translate {0, -1.3, 0.5} {
+ Duplicata { Volume{2}; }
+}
+Line(84) = {24, 46};
+Line(85) = {45, 28};
+Line Loop(86) = {84, -64, 85, -38};
+Plane Surface(87) = {86};
+Extrude {0.1, 0, 0} {
+ Surface{87};
+}
+Circle(110) = {50, 46, 56};
+Circle(111) = {54, 45, 60};
+Circle(112) = {55, 24, 20};
+Circle(113) = {64, 28, 19};
+Line Loop(114) = {37, 94, 112};
+Plane Surface(115) = {114};
+Line Loop(116) = {103, 113, -39};
+Plane Surface(117) = {116};
+Line Loop(118) = {95, -110, -65};
+Plane Surface(119) = {118};
+Line Loop(120) = {99, -111, 67};
+Plane Surface(121) = {120};
+Line Loop(122) = {110, 90, -111, -66};
+Ruled Surface(123) = {122};
+Line Loop(124) = {112, -36, -113, 92};
+Ruled Surface(125) = {124};
+Surface Loop(126) = {123, 119, 121, 63, 100};
+Volume(127) = {126};
+Surface Loop(128) = {125, 115, 117, 56, 108};
+Volume(129) = {128};
+
+
+Translate {0, 0, 0.5} {
+ Duplicata { Volume{2, 129, 58, 127, 57}; }
+}
+Line(267) = {29, 74};
+Line(268) = {70, 30};
+Line Loop(269) = {134, -267, 59, -268};
+Plane Surface(270) = {269};
+Extrude {-0.1, 0, 0} {
+ Surface{270};
+}
+Circle(293) = {420, 74, 65};
+Circle(294) = {419, 70, 66};
+Circle(295) = {38, 29, 424};
+Circle(296) = {34, 30, 428};
+Line Loop(297) = {282, -295, 62};
+Plane Surface(298) = {297};
+Line Loop(299) = {296, -286, 60};
+Plane Surface(300) = {299};
+Line Loop(301) = {135, -293, -278};
+Plane Surface(302) = {301};
+Line Loop(303) = {294, 133, 277};
+Plane Surface(304) = {303};
+Line Loop(305) = {293, 132, -294, 272};
+Ruled Surface(306) = {305};
+Line Loop(307) = {295, 274, -296, 61};
+Ruled Surface(308) = {307};
+Surface Loop(309) = {58, 287, 308, 298, 300};
+Volume(310) = {309};
+Surface Loop(311) = {131, 306, 302, 304, 279};
+Volume(312) = {311};
+Translate {5, 0, 0} {
+ Duplicata { Volume{1}; }
+}
+
+Line(340) = {445, 6};
+Line(341) = {10, 454};
+Line(342) = {438, 3};
+Line(343) = {2, 429};
+Line Loop(344) = {323, 340, 9, 341};
+Plane Surface(345) = {344};
+Line Loop(346) = {342, -1, 343, -318};
+Plane Surface(347) = {346};
+Extrude {0, 0, 1} {
+ Surface{345};
+}
+Extrude {0, 0, -1} {
+ Surface{347};
+}
+Circle(392) = {456, 445, 446};
+Circle(393) = {455, 454, 450};
+Circle(394) = {14, 10, 464};
+Circle(395) = {5, 6, 460};
+Circle(396) = {1, 2, 470};
+Circle(397) = {4, 3, 466};
+Circle(398) = {465, 438, 434};
+Circle(399) = {474, 429, 430};
+Line Loop(400) = {393, 322, 354};
+Plane Surface(401) = {400};
+Line Loop(402) = {363, -394, -10};
+Plane Surface(403) = {402};
+Line Loop(404) = {2, 397, -377};
+Plane Surface(405) = {404};
+Line Loop(406) = {4, 381, -396};
+Plane Surface(407) = {406};
+Line Loop(408) = {385, 399, -315};
+Plane Surface(409) = {408};
+Line Loop(410) = {392, -320, 355};
+Plane Surface(411) = {410};
+Line Loop(412) = {395, -359, -8};
+Plane Surface(413) = {412};
+Line Loop(414) = {394, -351, -395, -11};
+Ruled Surface(415) = {414};
+Line Loop(416) = {397, 372, -396, -3};
+Ruled Surface(417) = {416};
+Line Loop(418) = {399, 316, -398, -374};
+Ruled Surface(419) = {418};
+Line Loop(420) = {392, 321, -393, 349};
+Ruled Surface(421) = {420};
+Surface Loop(422) = {421, 411, 401, 356, 319};
+Volume(423) = {422};
+Line Loop(424) = {317, 376, 398};
+Plane Surface(425) = {424};
+Surface Loop(426) = {419, 409, 425, 314, 390};
+Volume(427) = {426};
+Surface Loop(428) = {382, 417, 405, 407, 6};
+Volume(429) = {428};
+Surface Loop(430) = {403, 415, 413, 364, 28};
+Volume(431) = {430};
+Translate {0, 0, 0.5} {
+ Duplicata { Volume{312, 241, 130, 161, 185, 216, 240}; }
+}
+
+Circle(624) = {414, 410, 576};
+Circle(625) = {418, 409, 572};
+Line Loop(626) = {243, 624, -479};
+Plane Surface(627) = {626};
+Line Loop(628) = {474, -625, 245};
+Plane Surface(629) = {628};
+Line Loop(630) = {624, -465, -625, -244};
+Ruled Surface(631) = {630};
+Surface Loop(632) = {631, 627, 629, 477, 241};
+Volume(633) = {632};
+Extrude {-5, 0, 0} {
+ Surface{598};
+}
+Extrude {-5, 0, 0} {
+ Surface{34};
+}
+
+Point(newp) = {-5.1, -5, 10, k*m};
+Point(newp) = {-5.1, 5, 10, k*m};
+Point(newp) = {-5.1, -5, -5, k*m};
+Point(newp) = {-5.1, 5, -5, k*m};
+
+Line(681) = {1022, 1021};
+Line(682) = {1023, 1021};
+Line(683) = {1022, 1024};
+Line(684) = {1024, 1023};
+
+Line Loop(685) = {683, 684, 682, -681};
+Line Loop(686) = {657, 658, 659, 660};
+Line Loop(687) = {635, 636, 637, 638};
+Plane Surface(688) = {685, 686, 687};
+
+Point(newp) = {10, -5, 10, k*m};
+Point(newp) = {10, 5, 10, k*m};
+Point(newp) = {10, -5, -5, k*m};
+Point(newp) = {10, 5, -5, k*m};
+Line(689) = {1028, 1027};
+Line(690) = {1028, 1026};
+Line(691) = {1026, 1025};
+Line(692) = {1025, 1027};
+Line Loop(693) = {692, -689, 690, 691};
+Plane Surface(694) = {693};
+Line(695) = {1027, 1023};
+Line(696) = {1024, 1028};
+Line(697) = {1025, 1021};
+Line(698) = {1026, 1022};
+Line Loop(699) = {695, -684, 696, 689};
+Plane Surface(700) = {699};
+Line Loop(701) = {697, -681, -698, 691};
+Plane Surface(702) = {701};
+Line Loop(703) = {690, 698, 683, 696};
+Plane Surface(704) = {703};
+Line Loop(705) = {697, -682, -695, -692};
+Plane Surface(706) = {705};
+
+Surface Loop(709) = {706, 702, 688, 704, 694, 700, 664, 668, 672, 676, 642, 646, 650, 654, 618, 623, 608, 613, 553, 543, 563, 548, 498, 503, 508, 513, 583, 574, 579, 528, 519, 524, 226, 217, 222, 302, 306, 304, 156, 141, 146, 151, 256, 261, 266, 251, 283, 270, 291, 292, 206, 186, 196, 191, 171, 162, 167, 443, 438, 447, 472, 457, 482, 462, 629, 631, 627, 47, 51, 55, 43, 115, 125, 117, 104, 87, 96, 109, 121, 123, 119, 83, 68, 73, 78, 298, 308, 300};
+Surface Loop(710) = {329, 334, 339, 324, 425, 419, 409, 391, 378, 347, 386, 417, 405, 407, 19, 23, 27, 15, 403, 415, 413, 369, 360, 345, 368, 401, 421, 411};
+Volume(711) = {709, 710};
+
+// coil
+Physical Volume(1) = {634, 635, 597, 573, 542, 518, 487, 432, 456, 633, 240, 216, 185, 161, 130, 312, 241, 310, 57, 127, 58, 129, 2};
+
+// coil terminals
+Physical Surface(2) = {655, 677};
+
+// air
+Physical Volume(3) = {711};
+
+// core
+Physical Volume(4) = {315, 429, 427, 313, 423, 314, 431, 1};
+
+// air surface
+Physical Surface(5) = {704, 700, 688, 694, 706, 702, 664, 668, 676, 672, 43, 55, 51, 47, 117, 125, 115, 96, 104, 87, 109, 121, 123, 119, 68, 83, 78, 73, 300, 308, 292, 298, 283, 270, 291, 302, 306, 304, 146, 151, 141, 156, 171, 162, 167, 196, 206, 186, 191, 226, 217, 222, 256, 261, 266, 251, 627, 631, 462, 457, 482, 472, 629, 443, 447, 438, 498, 513, 503, 508, 528, 519, 553, 563, 524, 548, 543, 579, 574, 583, 618, 613, 608, 623, 642, 654, 646, 650, 23, 19, 27, 15, 405, 417, 407, 347, 386, 378, 391, 425, 419, 409, 324, 339, 329, 334, 421, 411, 401, 368, 345, 369, 360, 413, 415, 403};
+
+// air surface & coil surface
+Physical Surface(6) = {704, 700, 688, 694, 706, 702, 664, 668, 676, 672, 43, 55, 51, 47, 117, 125, 115, 96, 104, 87, 109, 121, 123, 119, 68, 83, 78, 73, 300, 308, 292, 298, 283, 270, 291, 302, 306, 304, 146, 151, 141, 156, 171, 162, 167, 196, 206, 186, 191, 226, 217, 222, 256, 261, 266, 251, 627, 631, 462, 457, 482, 472, 629, 443, 447, 438, 498, 513, 503, 508, 528, 519, 553, 563, 524, 548, 543, 579, 574, 583, 618, 613, 608, 623, 642, 654, 646, 650};
+
+// air & core
+Physical Volume(7) = {711, 429, 315, 427, 313, 423, 314, 431, 1};
+
+// domain surface
+Physical Surface(8) = {688, 704, 700, 706, 694, 702, 655, 677};
+
+
+Mesh 3;
+
+HomGen("wirewound.msh") = {{3},{}};
+
+HomGen("wirewound.msh") = {{3},{5}};
+//HomCut("wirewound.msh") = {{3},{}};
+
+