diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp index 7b46f6ebdce65d3e989eff98956844cc57a02635..a87c799e6f8ecd59113f2ba8795f7d68c4f88d36 100644 --- a/Geo/MElement.cpp +++ b/Geo/MElement.cpp @@ -107,6 +107,19 @@ double MElement::rhoShapeMeasure() return 0.; } +void MElement::getNode(int num, double &u, double &v, double &w) +{ + // only for MElements that don't have a lookup table for this + // (currently only 1st order elements have) + double uvw[3]; + MVertex* ver = getVertex(num); + double xyz[3] = {ver->x(), ver->y(), ver->z()}; + xyz2uvw(xyz, uvw); + u = uvw[0]; + v = uvw[1]; + w = uvw[2]; +} + void MElement::getShapeFunctions(double u, double v, double w, double s[], int o) { const polynomialBasis* fs = getFunctionSpace(o); @@ -565,12 +578,21 @@ double MElement::integrateFlux(double val[], int face, int pOrder, int order) getFaceVertices(face, v); MElementFactory f; int type = 0; - if(getType() == TYPE_TRI || getType() == TYPE_TET) type = getTriangleType(getPolynomialOrder()); - else if(getType() == TYPE_QUA || getType() == TYPE_HEX) type = getQuadType(getPolynomialOrder()); - else if(getType() == TYPE_PYR && face < 4) type = getTriangleType(getPolynomialOrder()); - else if(getType() == TYPE_PYR && face >= 4) type = getQuadType(getPolynomialOrder()); - else if(getType() == TYPE_PRI && face < 2) type = getTriangleType(getPolynomialOrder()); - else if(getType() == TYPE_PRI && face >= 2) type = getQuadType(getPolynomialOrder()); + switch(getType()) { + case TYPE_TRI : type = getTriangleType(getPolynomialOrder()); break; + case TYPE_TET : type = getTriangleType(getPolynomialOrder()); break; + case TYPE_QUA : type = getQuadType(getPolynomialOrder()); break; + case TYPE_HEX : type = getQuadType(getPolynomialOrder()); break; + case TYPE_PYR : + if(face < 4) type = getTriangleType(getPolynomialOrder()); + else type = getQuadType(getPolynomialOrder()); + break; + case TYPE_PRI : + if(face < 2) type = getTriangleType(getPolynomialOrder()); + else type = getQuadType(getPolynomialOrder()); + break; + default: type = 0; break; + } MElement* fe = f.create(type, v); double intv[3]; diff --git a/Geo/MElement.h b/Geo/MElement.h index aa649f97d17747bd6b9d849febf4ddcd644621e4..b4836a4dc12105a22b58aa16e6ebe2d1bca95da7 100644 --- a/Geo/MElement.h +++ b/Geo/MElement.h @@ -220,12 +220,8 @@ class MElement // get the function space for the jacobian of the element virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const { return 0; } - // return parametric coordinates (u, v, w) of a vertex - virtual void getNode(int num, double &u, double &v, double &w) - { - u = v = w = 0.; - Msg::Error("Node get not available for this element"); - } + // return parametric coordinates (u,v,w) of a vertex + virtual void getNode(int num, double &u, double &v, double &w); // return the interpolating nodal shape functions evaluated at point // (u,v,w) in parametric coordinates (if order == -1, use the diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h index f015c66e3f0cb3b1e0d5d5436a145d2651426bd2..4336efe7cc25329525d5b80cc1be48111b472c7a 100644 --- a/Geo/MHexahedron.h +++ b/Geo/MHexahedron.h @@ -161,7 +161,7 @@ class MHexahedron : public MElement { s[7][1] = 0.125 * (1. - u) * (1. + w); s[7][2] = 0.125 * (1. - u) * (1. + v); } - void getNode(int num, double &u, double &v, double &w) + virtual void getNode(int num, double &u, double &v, double &w) { switch(num) { case 0 : u = -1.; v = -1.; w = -1.; break; @@ -352,6 +352,10 @@ class MHexahedron20 : public MHexahedron { _vs[8] = old[10]; _vs[10] = old[8]; _vs[9] = old[11]; _vs[11] = old[9]; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; /* @@ -493,6 +497,10 @@ class MHexahedron27 : public MHexahedron { _vs[13] = old[15]; _vs[15] = old[13]; _vs[14] = old[16]; _vs[16] = old[14]; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; /* @@ -632,6 +640,10 @@ class MHexahedronN : public MHexahedron { { Msg::Error("not done yet reverse hexN"); } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; #endif diff --git a/Geo/MLine.h b/Geo/MLine.h index 207aaceae6beee3b08ef9c82eff16a523ff61e45..bfbcde1eacadecf7afcdf6a42ccbd5654f333b8e 100644 --- a/Geo/MLine.h +++ b/Geo/MLine.h @@ -78,7 +78,7 @@ class MLine : public MElement { return false; return true; } - void getNode(int num, double &u, double &v, double &w) + virtual void getNode(int num, double &u, double &v, double &w) { v = w = 0.; switch(num) { @@ -142,6 +142,10 @@ class MLine3 : public MLine { //virtual int getTypeForVTK() const { return 21; } virtual const char *getStringForPOS() const { return "SL2"; } virtual const char *getStringForINP() const { return "C1D3"; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 2 ? MLine::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; /* @@ -198,6 +202,10 @@ class MLineN : public MLine { if(_vs.size() == 9) return MSH_LIN_11; return 0; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 2 ? MLine::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; struct compareMLinePtr { diff --git a/Geo/MPoint.h b/Geo/MPoint.h index 1d0e201f82513066bc6668d8ad0e17b7bd7ecb1c..01402a55536b8a5066f9e18a6fa75c8de6bf5a6b 100644 --- a/Geo/MPoint.h +++ b/Geo/MPoint.h @@ -42,7 +42,7 @@ class MPoint : public MElement { virtual int getTypeForMSH() const { return MSH_PNT; } virtual int getTypeForVTK() const { return 1; } virtual const char *getStringForPOS() const { return "SP"; } - void getNode(int num, double &u, double &v, double &w) + virtual void getNode(int num, double &u, double &v, double &w) { u = v = w = 0.; } diff --git a/Geo/MPrism.h b/Geo/MPrism.h index cd3b22c2f21d73965842d33507389c0fb65459c6..8702faa1c16a0d770f6b16f2ffe18ced4de107aa 100644 --- a/Geo/MPrism.h +++ b/Geo/MPrism.h @@ -159,7 +159,7 @@ class MPrism : public MElement { s[5][1] = 0.5 * (1. + w) ; s[5][2] = 0.5 * v ; }*/ - void getNode(int num, double &u, double &v, double &w) + virtual void getNode(int num, double &u, double &v, double &w) { switch(num) { case 0 : u = 0.; v = 0.; w = -1.; break; @@ -329,6 +329,10 @@ class MPrism15 : public MPrism { tmp = _vs[2]; _vs[2] = _vs[4]; _vs[4] = tmp; tmp = _vs[7]; _vs[7] = _vs[8]; _vs[8] = tmp; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 6 ? MPrism::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; /* @@ -451,6 +455,10 @@ class MPrism18 : public MPrism { // quad face vertices tmp = _vs[10]; _vs[10] = _vs[11]; _vs[11] = tmp; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 6 ? MPrism::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; #endif diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h index 423352f07eb7178d041465c1b1deadb9a668db1a..6ee6e36ea525c001dce033621a2152befcc5ddde 100644 --- a/Geo/MPyramid.h +++ b/Geo/MPyramid.h @@ -165,7 +165,7 @@ class MPyramid : public MElement { s[4][1] = 0.; s[4][2] = 1.; } - void getNode(int num, double &u, double &v, double &w) + virtual void getNode(int num, double &u, double &v, double &w) { switch(num) { case 0 : u = -1.; v = -1.; w = 0.; break; @@ -321,6 +321,10 @@ class MPyramid13 : public MPyramid { tmp = _vs[1]; _vs[1] = _vs[5]; _vs[5] = tmp; tmp = _vs[2]; _vs[2] = _vs[6]; _vs[6] = tmp; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 5 ? MPyramid::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; /* @@ -437,6 +441,10 @@ class MPyramid14 : public MPyramid { tmp = _vs[1]; _vs[1] = _vs[5]; _vs[5] = tmp; tmp = _vs[2]; _vs[2] = _vs[6]; _vs[6] = tmp; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 5 ? MPyramid::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; #endif diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h index 28c2fa5cb1fc3fcaf4e99dc6cbe3d5849fa628d7..11e1bc18dcf29d6a8ab66aaf4af72e6d3be3f695 100644 --- a/Geo/MQuadrangle.h +++ b/Geo/MQuadrangle.h @@ -117,7 +117,7 @@ void projectInMeanPlane(double *xn, double *yn); virtual const char *getStringForINP() const { return "C2D4"; } virtual const polynomialBasis* getFunctionSpace(int o=-1) const; virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const; - void getNode(int num, double &u, double &v, double &w) + virtual void getNode(int num, double &u, double &v, double &w) { w = 0.; switch(num) { @@ -237,6 +237,10 @@ class MQuadrangle8 : public MQuadrangle { tmp = _vs[0]; _vs[0] = _vs[3]; _vs[3] = tmp; tmp = _vs[1]; _vs[1] = _vs[2]; _vs[2] = tmp; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; /* @@ -309,6 +313,10 @@ class MQuadrangle9 : public MQuadrangle { tmp = _vs[0]; _vs[0] = _vs[3]; _vs[3] = tmp; tmp = _vs[1]; _vs[1] = _vs[2]; _vs[2] = tmp; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; /* @@ -399,6 +407,10 @@ class MQuadrangleN : public MQuadrangle { inv.insert(inv.begin(), _vs.rbegin(), _vs.rend()); _vs = inv; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; template <class T> diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h index 6d4b761bedde063c6d9eab92d45d183c5f1b8af3..f666db4a18903c99ed0caf04791614c964573044 100644 --- a/Geo/MTetrahedron.h +++ b/Geo/MTetrahedron.h @@ -130,7 +130,7 @@ class MTetrahedron : public MElement { void xyz2uvw(double xyz[3], double uvw[3]); virtual const polynomialBasis* getFunctionSpace(int o=-1) const; virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const; - void getNode(int num, double &u, double &v, double &w) + virtual void getNode(int num, double &u, double &v, double &w) { switch(num) { case 0 : u = 0.; v = 0.; w = 0.; break; @@ -269,6 +269,10 @@ class MTetrahedron10 : public MTetrahedron { tmp = _vs[1]; _vs[1] = _vs[2]; _vs[2] = tmp; tmp = _vs[5]; _vs[5] = _vs[3]; _vs[3] = tmp; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 4 ? MTetrahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; /* @@ -394,6 +398,10 @@ class MTetrahedronN : public MTetrahedron { virtual int getNumEdgesRep(); virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n); virtual int getNumFacesRep(); + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 4 ? MTetrahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; #endif diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h index 0e64bb3e276fbf2421c403a8e6dc5c6891f2eb6b..6b675fe4ac9a717a4c4b18dc784de39953162ed5 100644 --- a/Geo/MTriangle.h +++ b/Geo/MTriangle.h @@ -125,7 +125,7 @@ class MTriangle : public MElement { } virtual const polynomialBasis* getFunctionSpace(int o=-1) const; virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const; - void getNode(int num, double &u, double &v, double &w) + virtual void getNode(int num, double &u, double &v, double &w) { w = 0.; switch(num) { @@ -226,6 +226,10 @@ class MTriangle6 : public MTriangle { tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp; tmp = _vs[0]; _vs[0] = _vs[2]; _vs[2] = tmp; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 3 ? MTriangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; /* @@ -328,6 +332,10 @@ class MTriangleN : public MTriangle { inv.insert(inv.begin(), _vs.rbegin(), _vs.rend()); _vs = inv; } + virtual void getNode(int num, double &u, double &v, double &w) + { + num < 3 ? MTriangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w); + } }; template <class T>