diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp index f75c507b67b84b615ea0667ee1e4ed8a4473055e..4e76eb969336f83d975a84bc13467a807d12c5df 100644 --- a/Geo/MTetrahedron.cpp +++ b/Geo/MTetrahedron.cpp @@ -21,13 +21,13 @@ void MTetrahedron::getEdgeRep(bool curved, int num, double *x, double *y, double SVector3 *n) { // don't use MElement::_getEdgeRep: it's slow due to the creation of MFaces - MVertex *v0 = getVertex(edges_tetra(num, 0)); - MVertex *v1 = getVertex(edges_tetra(num, 1)); + MVertex *v0 = _v[edges_tetra(num, 0)]; + MVertex *v1 = _v[edges_tetra(num, 1)]; x[0] = v0->x(); y[0] = v0->y(); z[0] = v0->z(); x[1] = v1->x(); y[1] = v1->y(); z[1] = v1->z(); if(CTX::instance()->mesh.lightLines > 1){ static const int vv[6] = {2, 0, 1, 1, 0, 2}; - MVertex *v2 = getVertex(vv[num]); + MVertex *v2 = _v[vv[num]]; SVector3 t1(x[1] - x[0], y[1] - y[0], z[1] - z[0]); SVector3 t2(v2->x() - x[0], v2->y() - y[0], v2->z() - z[0]); SVector3 normal = crossprod(t1, t2); @@ -43,9 +43,9 @@ void MTetrahedron::getFaceRep(bool curved, int num, double *x, double *y, double SVector3 *n) { // don't use MElement::_getFaceRep: it's slow due to the creation of MFaces - MVertex *v0 = getVertex(faces_tetra(num, 0)); - MVertex *v1 = getVertex(faces_tetra(num, 1)); - MVertex *v2 = getVertex(faces_tetra(num, 2)); + MVertex *v0 = _v[faces_tetra(num, 0)]; + MVertex *v1 = _v[faces_tetra(num, 1)]; + MVertex *v2 = _v[faces_tetra(num, 2)]; x[0] = v0->x(); x[1] = v1->x(); x[2] = v2->x(); y[0] = v0->y(); y[1] = v1->y(); y[2] = v2->y(); z[0] = v0->z(); z[1] = v1->z(); z[2] = v2->z(); diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp index 48cd869369dc33f13ed7d7443d08d875cf2806cf..836aebd428790931e9d156b4a063c87e8ed85668 100644 --- a/Geo/MTriangle.cpp +++ b/Geo/MTriangle.cpp @@ -17,6 +17,42 @@ #define SQU(a) ((a)*(a)) +void MTriangle::getEdgeRep(bool curved, int num, double *x, double *y, double *z, + SVector3 *n) +{ + // don't use MElement::_getEdgeRep: it's slow due to the creation of MFace + MVertex *v0 = _v[edges_tri(num, 0)]; + MVertex *v1 = _v[edges_tri(num, 1)]; + x[0] = v0->x(); y[0] = v0->y(); z[0] = v0->z(); + x[1] = v1->x(); y[1] = v1->y(); z[1] = v1->z(); + if(CTX::instance()->mesh.lightLines > 1){ + static const int vv[3] = {2, 0, 1}; + MVertex *v2 = _v[vv[num]]; + SVector3 t1(x[1] - x[0], y[1] - y[0], z[1] - z[0]); + SVector3 t2(v2->x() - x[0], v2->y() - y[0], v2->z() - z[0]); + SVector3 normal = crossprod(t1, t2); + normal.normalize(); + n[0] = n[1] = normal; + } + else{ + n[0] = n[1] = SVector3(0., 0., 1.); + } +} + +void MTriangle::getFaceRep(bool curved, int num, double *x, double *y, double *z, + SVector3 *n) +{ + // don't use MElement::_getFaceRep: it's slow due to the creation of MFaces + x[0] = _v[0]->x(); x[1] = _v[1]->x(); x[2] = _v[2]->x(); + y[0] = _v[0]->y(); y[1] = _v[1]->y(); y[2] = _v[2]->y(); + z[0] = _v[0]->z(); z[1] = _v[1]->z(); z[2] = _v[2]->z(); + SVector3 t1(x[1] - x[0], y[1] - y[0], z[1] - z[0]); + SVector3 t2(x[2] - x[0], y[2] - y[0], z[2] - z[0]); + SVector3 normal = crossprod(t1, t2); + normal.normalize(); + for(int i = 0; i < 3; i++) n[i] = normal; +} + SPoint3 MTriangle::circumcenter() { double p1[3] = {_v[0]->x(), _v[0]->y(), _v[0]->z()}; @@ -250,6 +286,7 @@ void MTriangleN::getFaceRep(bool curved, int num, if (curved) _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges); else MTriangle::getFaceRep(false, num, x, y, z, n); } + void MTriangle6::getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n) { @@ -263,10 +300,10 @@ void MTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) *pts = getGQTPts(pOrder); } -void MTriangle::reorient(int rot,bool swap) { - +void MTriangle::reorient(int rot,bool swap) +{ if (rot == 0 && !swap) return; - + MVertex* tmp[3]; std::memcpy(tmp,_v,3*sizeof(MVertex*)); if (swap) for (int i=0;i<3;i++) _v[i] = tmp[(3-i+rot)%3]; @@ -275,8 +312,8 @@ void MTriangle::reorient(int rot,bool swap) { #include "HighOrder.h" -void MTriangle6::reorient(int rot, bool swap) { - +void MTriangle6::reorient(int rot, bool swap) +{ if (rot == 0 && !swap) return; MTriangle::reorient(rot,swap); @@ -286,18 +323,17 @@ void MTriangle6::reorient(int rot, bool swap) { else for (int i=0;i<3;i++) _vs[i] = tmp[(3-rot+i)%3]; } -void MTriangleN::reorient(int rot, bool swap) { - +void MTriangleN::reorient(int rot, bool swap) +{ if (rot == 0 && !swap) return; MTriangle::reorient(rot,swap); - + std::vector<MVertex*> tmp; int order = getPolynomialOrder(); int nbEdge = order - 1; - unsigned int idx = 0; - - + unsigned int idx = 0; + if (swap) { for (int iEdge=0;iEdge<3;iEdge++) { int edgeIdx = ((5-iEdge+rot)%3)*nbEdge; @@ -312,20 +348,15 @@ void MTriangleN::reorient(int rot, bool swap) { } idx += 3*nbEdge; - + if (_vs.size() > idx ) { if (order == 3) tmp.push_back(_vs[idx]); if (order == 4) { if (swap) for(int i=0;i<3;i++) tmp.push_back(_vs[idx+(3+rot-i)%3]); else for(int i=0;i<3;i++) tmp.push_back(_vs[idx+(3+i-rot)%3]); } - if (order >=5) + if (order >=5) Msg::Error("Reorientation of a triangle not supported above order 4"); } _vs = tmp; } - - - - - diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h index 31b856fb566b5d7eaa32aeda6171a646eabc15e5..a49ae4eb34e9a320df3dbc38e70a1a99bc2dcce2 100644 --- a/Geo/MTriangle.h +++ b/Geo/MTriangle.h @@ -87,11 +87,7 @@ class MTriangle : public MElement { Msg::Error("Could not get edge information for triangle %d", getNum()); } virtual int getNumEdgesRep(bool curved){ return 3; } - virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n) - { - MEdge e(getEdge(num)); - _getEdgeRep(e.getVertex(0), e.getVertex(1), x, y, z, n, 0); - } + virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n); virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const { v.resize(2); @@ -103,10 +99,7 @@ class MTriangle : public MElement { return MFace(_v[0], _v[1], _v[2]); } virtual int getNumFacesRep(bool curved){ return 1; } - virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n) - { - _getFaceRep(_v[0], _v[1], _v[2], x, y, z, n); - } + virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n); virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const { v.resize(3);