From dd964959a3d70b4b4b56209c3565dc0caa782e5c Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Thu, 27 Jun 2013 13:05:07 +0000
Subject: [PATCH] Enable prisms and pyramids in homology solver
---
Geo/Cell.cpp | 157 +++++++++++++++++++++++++++++++++++++++----
Geo/CellComplex.cpp | 5 +-
Geo/Chain.cpp | 16 ++---
Geo/ChainComplex.cpp | 2 +-
Geo/MPyramid.cpp | 6 ++
Geo/MPyramid.h | 1 +
6 files changed, 162 insertions(+), 25 deletions(-)
diff --git a/Geo/Cell.cpp b/Geo/Cell.cpp
index d8c773c4f8..85a0163ea8 100644
--- a/Geo/Cell.cpp
+++ b/Geo/Cell.cpp
@@ -144,16 +144,12 @@ void Cell::findBdElement(int i, std::vector<MVertex*>& vertices) const
vertices.push_back(_v[MTetrahedron::faces_tetra(i, j)]);
return;
case 5:
- if(i < 3) {
- for(int j = 0; j < 3; j++)
- vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
- }
- else {
- vertices.push_back(_v[0]);
- vertices.push_back(_v[3]);
- vertices.push_back(_v[2]);
- vertices.push_back(_v[1]);
- }
+ if(i < 4)
+ for(int j = 0; j < 3; j++)
+ vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
+ else
+ for(int j = 0; j < 4; j++)
+ vertices.push_back(_v[MPyramid::faces_pyramid(i, j)]);
return;
case 6:
if(i < 2)
@@ -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, 2)]->getNum() == v[0]->getNum())
return -1;
- case 5: return 0;
- case 6: return 0;
+ case 5:
+ 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:
if(_v[MHexahedron::faces_hexa(i, 0)]->getNum() == v[0]->getNum() &&
_v[MHexahedron::faces_hexa(i, 1)]->getNum() == v[1]->getNum() &&
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 6638f0fd74..85990e92fc 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -82,11 +82,12 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
MElement* element = elements.at(i);
int dim = element->getDim();
int type = element->getType();
- if(type == TYPE_PYR || type == TYPE_PRI ||
+ if(//type == TYPE_PYR || type == TYPE_PRI ||
type == TYPE_POLYG || type == TYPE_POLYH) {
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;
std::pair<Cell*, bool> maybeCell = Cell::createCell(element, domain);
if(!maybeCell.second) {
diff --git a/Geo/Chain.cpp b/Geo/Chain.cpp
index 3c7b04b3f6..13ba5623c3 100644
--- a/Geo/Chain.cpp
+++ b/Geo/Chain.cpp
@@ -222,16 +222,12 @@ void ElemChain::getBoundaryVertices(int i, int dim, int numVertices,
vertices.push_back(v[MTetrahedron::faces_tetra(i, j)]);
return;
case 5:
- if(i < 3) {
- for(int j = 0; j < 3; j++)
- vertices.push_back(v[MPyramid::faces_pyramid(i, j)]);
- }
- else {
- vertices.push_back(v[0]);
- vertices.push_back(v[3]);
- vertices.push_back(v[2]);
- vertices.push_back(v[1]);
- }
+ if(i < 3)
+ for(int j = 0; j < 3; j++)
+ vertices.push_back(v[MPyramid::faces_pyramid(i, j)]);
+ else
+ for(int j = 0; j < 4; j++)
+ vertices.push_back(v[MPyramid::faces_pyramid(i, j)]);
return;
case 6:
if(i < 2)
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 114f90766e..d0da0319e7 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -337,7 +337,7 @@ void ChainComplex::computeHomology(bool dual)
}
// 2) this dimension is empty
- else if(getHMatrix(lowDim) == NULL){
+ else if(getHMatrix(setDim) == NULL){
setHbasis(setDim, NULL);
}
// 3) No higher dimension cells -> none of the cycles are boundaries
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index 650cb5115b..bd553a6e60 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -27,6 +27,12 @@ int MPyramid::getVolumeSign()
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
{
int order = (o == -1) ? getPolynomialOrder() : o;
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index bad11f2bea..2403d674ea 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -149,6 +149,7 @@ class MPyramid : public MElement {
return false;
return true;
}
+ void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
static int edges_pyramid(const int edge, const int vert)
{
static const int e[8][2] = {
--
GitLab