Skip to content
Snippets Groups Projects
Commit 9c71d28a authored by Amaury Johnen's avatar Amaury Johnen
Browse files

cleanup

parent 7569fbee
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 13 additions and 124 deletions
Mesh/HighOrder.cpp
+ 13
− 124
View file @ 9c71d28a
......@@ -576,99 +576,6 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
}
}
static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
std::vector<MVertex*> &vr, bool linear, int nPts = 1)
{
fullMatrix<double> points;
int start = 0;
switch (incomplete->getType()){
case TYPE_TET :
switch (nPts){
case 0: return;
case 1: return;
case 2: points = BasisFactory::getNodalBasis(MSH_TET_20)->points; break;
case 3: points = BasisFactory::getNodalBasis(MSH_TET_35)->points; break;
case 4: points = BasisFactory::getNodalBasis(MSH_TET_56)->points; break;
case 5: points = BasisFactory::getNodalBasis(MSH_TET_84)->points; break;
case 6: points = BasisFactory::getNodalBasis(MSH_TET_120)->points; break;
case 7: points = BasisFactory::getNodalBasis(MSH_TET_165)->points; break;
case 8: points = BasisFactory::getNodalBasis(MSH_TET_220)->points; break;
case 9: points = BasisFactory::getNodalBasis(MSH_TET_286)->points; break;
default:
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
}
start = ((nPts+2)*(nPts+3)*(nPts+4)-(nPts-2)*(nPts-1)*(nPts))/6; // = 4+6*(p-1)+4*(p-2)*(p-1)/2 = 4*(p+1)*(p+2)/2-6*(p-1)-2*4 = 2*p*p+2
break;
case TYPE_HEX :
switch (nPts){
case 0: return;
case 1: points = BasisFactory::getNodalBasis(MSH_HEX_27)->points; break;
case 2: points = BasisFactory::getNodalBasis(MSH_HEX_64)->points; break;
case 3: points = BasisFactory::getNodalBasis(MSH_HEX_125)->points; break;
case 4: points = BasisFactory::getNodalBasis(MSH_HEX_216)->points; break;
case 5: points = BasisFactory::getNodalBasis(MSH_HEX_343)->points; break;
case 6: points = BasisFactory::getNodalBasis(MSH_HEX_512)->points; break;
case 7: points = BasisFactory::getNodalBasis(MSH_HEX_729)->points; break;
case 8: points = BasisFactory::getNodalBasis(MSH_HEX_1000)->points; break;
default :
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
}
start = (nPts+2)*(nPts+2)*(nPts+2) - (nPts)*(nPts)*(nPts) ; // = 6*(p-1)*(p-1)+12*(p-1)+8 = 6*(p+1)*(p+1)-12*(p-1)-2*8 = 6*p*p+2
break;
case TYPE_PRI :
switch (nPts){
case 0: return;
case 1: return;
case 2: points = BasisFactory::getNodalBasis(MSH_PRI_40)->points; break;
case 3: points = BasisFactory::getNodalBasis(MSH_PRI_75)->points; break;
case 4: points = BasisFactory::getNodalBasis(MSH_PRI_126)->points; break;
case 5: points = BasisFactory::getNodalBasis(MSH_PRI_196)->points; break;
case 6: points = BasisFactory::getNodalBasis(MSH_PRI_288)->points; break;
case 7: points = BasisFactory::getNodalBasis(MSH_PRI_405)->points; break;
case 8: points = BasisFactory::getNodalBasis(MSH_PRI_550)->points; break;
default:
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
}
start = 4*(nPts+1)*(nPts+1)+2; // = 4*p*p+2 = 6+9*(p-1)+2*(p-2)*(p-1)/2+3*(p-1)*(p-1) = 2*(p+1)*(p+2)/2+3*(p+1)*(p+1)-9*(p-1)-2*6
break;
case TYPE_PYR:
switch (nPts){
case 0:
case 1: return;
case 2: points = BasisFactory::getNodalBasis(MSH_PYR_30)->points; break;
case 3: points = BasisFactory::getNodalBasis(MSH_PYR_55)->points; break;
case 4: points = BasisFactory::getNodalBasis(MSH_PYR_91)->points; break;
case 5: points = BasisFactory::getNodalBasis(MSH_PYR_140)->points; break;
case 6: points = BasisFactory::getNodalBasis(MSH_PYR_204)->points; break;
case 7: points = BasisFactory::getNodalBasis(MSH_PYR_285)->points; break;
case 8: points = BasisFactory::getNodalBasis(MSH_PYR_385)->points; break;
default :
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
}
start = ( nPts+2) * ( (nPts+2) + 1) * (2*(nPts+2) + 1) / 6 -
(nPts-1) * ( (nPts-1) + 1) * (2*(nPts-1) + 1) / 6;
break;
}
for(int k = start; k < points.size1(); k++){
MVertex *v;
const double t1 = points(k, 0);
const double t2 = points(k, 1);
const double t3 = points(k, 2);
SPoint3 pos;
incomplete->pnt(t1, t2, t3, pos);
v = new MVertex(pos.x(), pos.y(), pos.z(), gr);
gr->mesh_vertices.push_back(v);
vr.push_back(v);
}
}
static void getFaceAndInteriorVertices(GRegion *gr, MElement *incomplete, MElement *ele,
std::vector<MVertex*> &vr, faceContainer &faceVertices,
bool linear, int nPts = 1)
......@@ -691,7 +598,7 @@ static void getFaceAndInteriorVertices(GRegion *gr, MElement *incomplete, MEleme
case 8: points = BasisFactory::getNodalBasis(MSH_TET_220)->points; break;
case 9: points = BasisFactory::getNodalBasis(MSH_TET_286)->points; break;
default:
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
Msg::Error("getFaceAndInteriorVertices not implemented for order %i", nPts+1);
break;
}
startFace = 4 + nPts*6;
......@@ -709,7 +616,7 @@ static void getFaceAndInteriorVertices(GRegion *gr, MElement *incomplete, MEleme
case 7: points = BasisFactory::getNodalBasis(MSH_HEX_729)->points; break;
case 8: points = BasisFactory::getNodalBasis(MSH_HEX_1000)->points; break;
default :
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
Msg::Error("getFaceAndInteriorVertices not implemented for order %i", nPts+1);
break;
}
startFace = 8 + nPts*12 ;
......@@ -727,7 +634,7 @@ static void getFaceAndInteriorVertices(GRegion *gr, MElement *incomplete, MEleme
case 7: points = BasisFactory::getNodalBasis(MSH_PRI_405)->points; break;
case 8: points = BasisFactory::getNodalBasis(MSH_PRI_550)->points; break;
default:
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
Msg::Error("getFaceAndInteriorVertices not implemented for order %i", nPts+1);
break;
}
startFace = 6 + nPts*9;
......@@ -745,7 +652,7 @@ static void getFaceAndInteriorVertices(GRegion *gr, MElement *incomplete, MEleme
case 7: points = BasisFactory::getNodalBasis(MSH_PYR_285)->points; break;
case 8: points = BasisFactory::getNodalBasis(MSH_PYR_385)->points; break;
default :
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
Msg::Error("getFaceAndInteriorVertices not implemented for order %i", nPts+1);
break;
}
startFace = 5 + nPts*8;
......@@ -933,10 +840,6 @@ static MTetrahedron *setHighOrder(MTetrahedron *t, GRegion *gr,
// it either way)
MTetrahedronN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2),
t->getVertex(3), ve, nPts + 1, 0, t->getPartition());
//getFaceVertices(gr, t, vf, faceVertices, edgeVertices, linear, nPts);
//ve.insert(ve.end(), vf.begin(), vf.end());
//getRegionVertices(gr, &incpl, t, vr, linear, nPts);
//ve.insert(ve.end(), vr.begin(), vr.end());
getFaceAndInteriorVertices(gr, &incpl, t, vf, faceVertices, linear, nPts);
ve.insert(ve.end(), vf.begin(), vf.end());
}
......@@ -979,8 +882,6 @@ static MHexahedron *setHighOrder(MHexahedron *h, GRegion *gr,
h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2],
ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10],
ve[11], 0, h->getPartition());
//getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts);
//getRegionVertices(gr, &incpl, h, vr, linear, nPts);
getFaceAndInteriorVertices(gr, &incpl, h, vf, faceVertices, linear, nPts);
return new MHexahedron27(h->getVertex(0), h->getVertex(1), h->getVertex(2),
h->getVertex(3), h->getVertex(4), h->getVertex(5),
......@@ -994,10 +895,6 @@ static MHexahedron *setHighOrder(MHexahedron *h, GRegion *gr,
h->getVertex(3), h->getVertex(4), h->getVertex(5),
h->getVertex(6), h->getVertex(7), ve, nPts + 1, 0,
h->getPartition());
//getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts);
//ve.insert(ve.end(), vf.begin(), vf.end());
//getRegionVertices(gr, &incpl, h, vr, linear, nPts);
//ve.insert(ve.end(), vr.begin(), vr.end());
getFaceAndInteriorVertices(gr, &incpl, h, vf, faceVertices, linear, nPts);
ve.insert(ve.end(), vf.begin(), vf.end());
return new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2),
......@@ -1029,12 +926,16 @@ static MPrism *setHighOrder(MPrism *p, GRegion *gr,
}
}
else{
// create serendipity element to place internal vertices (we used to
// saturate face vertices also, but the corresponding function spaces do
// not exist anymore, and there is no theoretical justification for doing
// it either way)
MPrismN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2),
p->getVertex(3), p->getVertex(4), p->getVertex(5),
ve, nPts + 1, 0, p->getPartition());
getFaceAndInteriorVertices(gr, &incpl, p, vf, faceVertices, linear, nPts);
if (nPts == 1) {
MPrismN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2),
p->getVertex(3), p->getVertex(4), p->getVertex(5),
ve, nPts + 1, 0, p->getPartition());
//getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
getFaceAndInteriorVertices(gr, &incpl, p, vf, faceVertices, linear, nPts);
return new MPrism18(p->getVertex(0), p->getVertex(1), p->getVertex(2),
p->getVertex(3), p->getVertex(4), p->getVertex(5),
ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8],
......@@ -1042,18 +943,6 @@ static MPrism *setHighOrder(MPrism *p, GRegion *gr,
0, p->getPartition());
}
else {
// create serendipity element to place internal vertices (we used to
// saturate face vertices also, but the corresponding function spaces do
// not exist anymore, and there is no theoretical justification for doing
// it either way)
MPrismN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2),
p->getVertex(3), p->getVertex(4), p->getVertex(5),
ve, nPts + 1, 0, p->getPartition());
//getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
//ve.insert(ve.end(), vf.begin(), vf.end());
//getRegionVertices(gr, &incpl, p, vr, linear, nPts);
//ve.insert(ve.end(), vr.begin(), vr.end());
getFaceAndInteriorVertices(gr, &incpl, p, vf, faceVertices, linear, nPts);
ve.insert(ve.end(), vf.begin(), vf.end());
return new MPrismN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
p->getVertex(3), p->getVertex(4), p->getVertex(5),
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment