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

Enable prisms and pyramids in homology solver

parent c4c6d08c
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/Cell.cpp
+ 145
− 12
View file @ dd964959
...@@ -144,16 +144,12 @@ void Cell::findBdElement(int i, std::vector<MVertex*>& vertices) const ...@@ -144,16 +144,12 @@ void Cell::findBdElement(int i, std::vector<MVertex*>& vertices) const
vertices.push_back(_v[MTetrahedron::faces_tetra(i, j)]); vertices.push_back(_v[MTetrahedron::faces_tetra(i, j)]);
return; return;
case 5: case 5:
if(i < 3) { if(i < 4)
for(int j = 0; j < 3; j++) for(int j = 0; j < 3; j++)
vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]); vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
} else
else { for(int j = 0; j < 4; j++)
vertices.push_back(_v[0]); vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
vertices.push_back(_v[3]);
vertices.push_back(_v[2]);
vertices.push_back(_v[1]);
}
return; return;
case 6: case 6:
if(i < 2) if(i < 2)
...@@ -270,8 +266,145 @@ int Cell::findBdCellOrientation(Cell* cell, int i) const ...@@ -270,8 +266,145 @@ int Cell::findBdCellOrientation(Cell* cell, int i) const
_v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[1]->getNum() && _v[MTetrahedron::faces_tetra(i, 1)]->getNum() == v[1]->getNum() &&
_v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[0]->getNum()) _v[MTetrahedron::faces_tetra(i, 2)]->getNum() == v[0]->getNum())
return -1; return -1;
case 5: return 0; case 5:
case 6: return 0; if(i < 4) {
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum())
return 1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum())
return 1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum())
return 1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum())
return -1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum())
return -1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum())
return -1;
}
else {
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum() &&
_v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[3]->getNum())
return 1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[3]->getNum() &&
_v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[0]->getNum())
return 1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[3]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum() &&
_v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[1]->getNum())
return 1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[3]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum() &&
_v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[2]->getNum())
return 1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[0]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[3]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[2]->getNum() &&
_v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[1]->getNum())
return -1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[3]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[2]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[1]->getNum() &&
_v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[0]->getNum())
return -1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[2]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[1]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[0]->getNum() &&
_v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[3]->getNum())
return -1;
if(_v[MPyramid::faces_pyramid(i, 0)]->getNum() == v[1]->getNum() &&
_v[MPyramid::faces_pyramid(i, 1)]->getNum() == v[0]->getNum() &&
_v[MPyramid::faces_pyramid(i, 2)]->getNum() == v[3]->getNum() &&
_v[MPyramid::faces_pyramid(i, 3)]->getNum() == v[2]->getNum())
return -1;
}
return 0;
case 6:
if(i < 2) {
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum())
return 1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum())
return 1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum())
return 1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum())
return -1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum())
return -1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum())
return -1;
}
else {
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum() &&
_v[MPrism::faces_prism(i, 3)]->getNum() == v[3]->getNum())
return 1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[3]->getNum() &&
_v[MPrism::faces_prism(i, 3)]->getNum() == v[0]->getNum())
return 1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[3]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum() &&
_v[MPrism::faces_prism(i, 3)]->getNum() == v[1]->getNum())
return 1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[3]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum() &&
_v[MPrism::faces_prism(i, 3)]->getNum() == v[2]->getNum())
return 1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[0]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[3]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[2]->getNum() &&
_v[MPrism::faces_prism(i, 3)]->getNum() == v[1]->getNum())
return -1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[3]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[2]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[1]->getNum() &&
_v[MPrism::faces_prism(i, 3)]->getNum() == v[0]->getNum())
return -1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[2]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[1]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[0]->getNum() &&
_v[MPrism::faces_prism(i, 3)]->getNum() == v[3]->getNum())
return -1;
if(_v[MPrism::faces_prism(i, 0)]->getNum() == v[1]->getNum() &&
_v[MPrism::faces_prism(i, 1)]->getNum() == v[0]->getNum() &&
_v[MPrism::faces_prism(i, 2)]->getNum() == v[3]->getNum() &&
_v[MPrism::faces_prism(i, 3)]->getNum() == v[2]->getNum())
return -1;
}
case 8: case 8:
if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[0]->getNum() && if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[0]->getNum() &&
_v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[1]->getNum() && _v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[1]->getNum() &&
......
Geo/CellComplex.cpp
+ 3
− 2
View file @ dd964959
...@@ -82,11 +82,12 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements, ...@@ -82,11 +82,12 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
MElement* element = elements.at(i); MElement* element = elements.at(i);
int dim = element->getDim(); int dim = element->getDim();
int type = element->getType(); int type = element->getType();
if(type == TYPE_PYR || type == TYPE_PRI || if(//type == TYPE_PYR || type == TYPE_PRI ||
type == TYPE_POLYG || type == TYPE_POLYH) { type == TYPE_POLYG || type == TYPE_POLYH) {
Msg::Error("Mesh element type %d not implemented in homology solver", type); Msg::Error("Mesh element type %d not implemented in homology solver", type);
} }
if(type == TYPE_QUA || type == TYPE_HEX) if(type == TYPE_QUA || type == TYPE_HEX ||
type == TYPE_PYR || type == TYPE_PRI)
_simplicial = false; _simplicial = false;
std::pair<Cell*, bool> maybeCell = Cell::createCell(element, domain); std::pair<Cell*, bool> maybeCell = Cell::createCell(element, domain);
if(!maybeCell.second) { if(!maybeCell.second) {
......
Geo/Chain.cpp
+ 6
− 10
View file @ dd964959
...@@ -222,16 +222,12 @@ void ElemChain::getBoundaryVertices(int i, int dim, int numVertices, ...@@ -222,16 +222,12 @@ void ElemChain::getBoundaryVertices(int i, int dim, int numVertices,
vertices.push_back(v[MTetrahedron::faces_tetra(i, j)]); vertices.push_back(v[MTetrahedron::faces_tetra(i, j)]);
return; return;
case 5: case 5:
if(i < 3) { if(i < 3)
for(int j = 0; j < 3; j++) for(int j = 0; j < 3; j++)
vertices.push_back(v[MPyramid::faces_pyramid(i, j)]); vertices.push_back(v[MPyramid::faces_pyramid(i, j)]);
} else
else { for(int j = 0; j < 4; j++)
vertices.push_back(v[0]); vertices.push_back(v[MPyramid::faces_pyramid(i, j)]);
vertices.push_back(v[3]);
vertices.push_back(v[2]);
vertices.push_back(v[1]);
}
return; return;
case 6: case 6:
if(i < 2) if(i < 2)
......
Geo/ChainComplex.cpp
+ 1
− 1
View file @ dd964959
...@@ -337,7 +337,7 @@ void ChainComplex::computeHomology(bool dual) ...@@ -337,7 +337,7 @@ void ChainComplex::computeHomology(bool dual)
} }
// 2) this dimension is empty // 2) this dimension is empty
else if(getHMatrix(lowDim) == NULL){ else if(getHMatrix(setDim) == NULL){
setHbasis(setDim, NULL); setHbasis(setDim, NULL);
} }
// 3) No higher dimension cells -> none of the cycles are boundaries // 3) No higher dimension cells -> none of the cycles are boundaries
......
Geo/MPyramid.cpp
+ 6
− 0
View file @ dd964959
...@@ -27,6 +27,12 @@ int MPyramid::getVolumeSign() ...@@ -27,6 +27,12 @@ int MPyramid::getVolumeSign()
else return 0; else return 0;
} }
void MPyramid::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = getNGQPyrPts(pOrder);
*pts = getGQPyrPts(pOrder);
}
const nodalBasis* MPyramid::getFunctionSpace(int o) const const nodalBasis* MPyramid::getFunctionSpace(int o) const
{ {
int order = (o == -1) ? getPolynomialOrder() : o; int order = (o == -1) ? getPolynomialOrder() : o;
......
Geo/MPyramid.h
+ 1
− 0
View file @ dd964959
...@@ -149,6 +149,7 @@ class MPyramid : public MElement { ...@@ -149,6 +149,7 @@ class MPyramid : public MElement {
return false; return false;
return true; return true;
} }
void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
static int edges_pyramid(const int edge, const int vert) static int edges_pyramid(const int edge, const int vert)
{ {
static const int e[8][2] = { static const int e[8][2] = {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment