From 9c71d28a3ea21938bcbdd939011eb0732db2f132 Mon Sep 17 00:00:00 2001 From: Amaury Johnan <amjohnen@gmail.com> Date: Wed, 2 Oct 2013 14:26:42 +0000 Subject: [PATCH] cleanup --- Mesh/HighOrder.cpp | 137 +++++---------------------------------------- 1 file changed, 13 insertions(+), 124 deletions(-) diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp index 20f96999cf..03ee51ccff 100644 --- a/Mesh/HighOrder.cpp +++ b/Mesh/HighOrder.cpp @@ -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), -- GitLab