diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 51969a9b653a65fc04f2ca0093d57ecc73be2a83..56a7378bac874abe5e2f9f786acbd9f23b2a61c9 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -446,8 +446,15 @@ double GFace::curvature(const SPoint2 ¶m) const
SVector3 dndu = 500 * (n2 - n1);
SVector3 dndv = 500 * (n4 - n3);
- double c = fabs(dot(dndu, du) + dot(dndv, dv)) / detJ;
+ // double c = fabs(dot(dndu, du) + dot(dndv, dv)) / detJ;
+
+ double ddu = dot(dndu,du);
+ double ddv = dot(dndv,dv);
+
+ double c = std::max(fabs(ddu),fabs(ddv))/detJ;
+
+
// Msg::Info("c = %g detJ %g", c, detJ);
return c;
@@ -729,6 +736,7 @@ SPoint2 GFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t)
}
}
+
// length of a curve drawn on a surface
// S = (X(u,v), Y(u,v), Z(u,v) );
// u = u(t) , v = v(t)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 62e51f709ebedebf5130a8bef3042b6b4420fc8f..d58c857965fd0a9986be6bad026e7143738b4b11 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -98,6 +98,23 @@ void MElement::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
Msg::Error("No integration points defined for this type of element");
}
+void MElement::getShapeFunction(int num,double u,double v,double w,double &s) {
+#if !defined(HAVE_GMSH_EMBEDDED)
+ double sf[256];
+ getFunctionSpace()->f(u,v,w,sf);
+ s = sf[num];
+#endif
+}
+
+void MElement::getGradShapeFunction(int num,double uu,double vv,double ww,double s[3]) {
+#if !defined(HAVE_GMSH_EMBEDDED)
+ double sf[256][3];
+ getFunctionSpace()->df(uu, vv, ww, sf);
+ for (int i=0;i<3;i++) s[i] = sf[num][i];
+#endif
+}
+
+
SPoint3 MTriangle::circumcenter()
{
#if defined(HAVE_GMSH_EMBEDDED)
@@ -152,6 +169,35 @@ double MTriangle::gammaShapeMeasure()
#endif
}
+const gmshFunctionSpace* MTriangle::getFunctionSpace(int o) const {
+
+ int order = (o == -1) ? getPolynomialOrder() : o;
+
+ int nf = getNumFaceVertices();
+
+ if ((nf ==0) && (o == -1)) {
+ switch (order) {
+ case 1: return &gmshFunctionSpaces::find(MSH_TRI_3); break;
+ case 2: return &gmshFunctionSpaces::find(MSH_TRI_6); break;
+ case 3: return &gmshFunctionSpaces::find(MSH_TRI_9); break;
+ case 4: return &gmshFunctionSpaces::find(MSH_TRI_12); break;
+ case 5: return &gmshFunctionSpaces::find(MSH_TRI_15I); break;
+ default: Msg::Error("Order %d triangle function space not implemented", order); break;
+ }
+ }
+ else {
+ switch (order) {
+ case 1: return &gmshFunctionSpaces::find(MSH_TRI_3); break;
+ case 2: return &gmshFunctionSpaces::find(MSH_TRI_6); break;
+ case 3: return &gmshFunctionSpaces::find(MSH_TRI_10); break;
+ case 4: return &gmshFunctionSpaces::find(MSH_TRI_15); break;
+ case 5: return &gmshFunctionSpaces::find(MSH_TRI_21); break;
+ default: Msg::Error("Order %d triangle function space implemented", order); break;
+ }
+ }
+ return NULL;
+}
+
double MTetrahedron::gammaShapeMeasure()
{
#if defined(HAVE_GMSH_EMBEDDED)
@@ -162,6 +208,7 @@ double MTetrahedron::gammaShapeMeasure()
#endif
}
+
double MTetrahedron::etaShapeMeasure()
{
#if defined(HAVE_GMSH_EMBEDDED)
@@ -189,6 +236,39 @@ void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3])
sys3x3(mat, b, uvw, &det);
}
+
+const gmshFunctionSpace* MTetrahedron::getFunctionSpace(int o) const {
+
+ int order = (o == -1) ? getPolynomialOrder() : o;
+
+ int nv = getNumVolumeVertices();
+
+
+ if ((nv == 0) && (o == -1)) {
+ switch (order) {
+ case 1: return &gmshFunctionSpaces::find(MSH_TET_4); break;
+ case 2: return &gmshFunctionSpaces::find(MSH_TET_10); break;
+ case 3: return &gmshFunctionSpaces::find(MSH_TET_20); break;
+ case 4: return &gmshFunctionSpaces::find(MSH_TET_34); break;
+ case 5: return &gmshFunctionSpaces::find(MSH_TET_52); break;
+ default: Msg::Error("Order %d tetrahedron function space not implemented", order); break;
+ }
+ }
+ else {
+ switch (order) {
+ case 1: return &gmshFunctionSpaces::find(MSH_TET_4); break;
+ case 2: return &gmshFunctionSpaces::find(MSH_TET_10); break;
+ case 3: return &gmshFunctionSpaces::find(MSH_TET_20); break;
+ case 4: return &gmshFunctionSpaces::find(MSH_TET_35); break;
+ case 5: return &gmshFunctionSpaces::find(MSH_TET_56); break;
+ default: Msg::Error("Order %d tetrahedron function space implemented", order); break;
+ }
+ }
+ return NULL;
+}
+
+
+
int MHexahedron::getVolumeSign()
{
double mat[3][3];
@@ -259,70 +339,160 @@ std::string MElement::getInfoString()
double MElement::getJacobian(double u, double v, double w, double jac[3][3])
{
+
+ const gmshFunctionSpace* fs = getFunctionSpace();
+
jac[0][0] = jac[0][1] = jac[0][2] = 0.;
jac[1][0] = jac[1][1] = jac[1][2] = 0.;
jac[2][0] = jac[2][1] = jac[2][2] = 0.;
- double s[3];
- switch(getDim()){
- case 3 :
- for(int i = 0; i < getNumVertices(); i++) {
- getGradShapeFunction(i, u, v, w, s);
- MVertex *p = getVertex(i);
- jac[0][0] += p->x() * s[0]; jac[0][1] += p->y() * s[0]; jac[0][2] += p->z() * s[0];
- jac[1][0] += p->x() * s[1]; jac[1][1] += p->y() * s[1]; jac[1][2] += p->z() * s[1];
- jac[2][0] += p->x() * s[2]; jac[2][1] += p->y() * s[2]; jac[2][2] += p->z() * s[2];
+
+ double gsf[256][3];
+ fs->df(u,v,w,gsf);
+ for (int i=0;i<getNumVertices();i++) {
+
+ const MVertex* v = getVertex(i);
+ double* gg = gsf[i];
+
+ for (int j=0;j<3;j++) {
+ jac[j][0] += v->x() * gg[j];
+ jac[j][1] += v->y() * gg[j];
+ jac[j][2] += v->z() * gg[j];
+ }
+ }
+
+ double dJ = 0;
+
+ switch (fs->monomials.size2()) {
+
+ case 1:
+ {
+ dJ = sqrt(jac[0][0] * jac[0][0] + jac[1][1] * jac[0][0] + jac[2][2] * jac[2][2]);
+ break;
}
- return fabs(jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
+ case 2:
+ {
+ dJ = sqrt(SQU(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) +
+ SQU(jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2]) +
+ SQU(jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1]));
+ break;
+ }
+ case 3:
+ {
+ dJ = fabs(jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2]);
- case 2 :
- for(int i = 0; i < getNumVertices(); i++) {
- getGradShapeFunction(i, u, v, w, s);
- MVertex *p = getVertex(i);
- jac[0][0] += p->x() * s[0]; jac[0][1] += p->y() * s[0]; jac[0][2] += p->z() * s[0];
- jac[1][0] += p->x() * s[1]; jac[1][1] += p->y() * s[1]; jac[1][2] += p->z() * s[1];
+ break;
+ }
+ }
+ return dJ;
+}
+
+
+
+double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+{
+
+ const gmshFunctionSpace* fs = getFunctionSpace(1);
+
+ jac[0][0] = jac[0][1] = jac[0][2] = 0.;
+ jac[1][0] = jac[1][1] = jac[1][2] = 0.;
+ jac[2][0] = jac[2][1] = jac[2][2] = 0.;
+
+ double gsf[256][3];
+ fs->df(u,v,w,gsf);
+
+ for (int i=0;i<getNumPrimaryVertices();i++) {
+
+ const MVertex* v = getVertex(i);
+ double* gg = gsf[i];
+
+ for (int j=0;j<3;j++) {
+ jac[j][0] += v->x() * gg[j];
+ jac[j][1] += v->y() * gg[j];
+ jac[j][2] += v->z() * gg[j];
}
+ }
+
+ double dJ = 0;
+
+ switch (fs->monomials.size2()) {
+
+ case 1:
{
- double a[3], b[3], c[3];
- a[0] = jac[0][0];
- a[1] = jac[0][1];
- a[2] = jac[0][2];
- b[0] = jac[1][0];
- b[1] = jac[1][1];
- b[2] = jac[1][2];
- prodve(a, b, c);
- norme(c);
- jac[2][0] = c[0]; jac[2][1] = c[1]; jac[2][2] = c[2];
+ dJ = sqrt(jac[0][0] * jac[0][0] + jac[0][0] * jac[0][0] + jac[0][0] * jac[0][0]);
+ break;
}
- return sqrt(SQU(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) +
+ case 2:
+ {
+ dJ = sqrt(SQU(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) +
SQU(jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2]) +
SQU(jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1]));
- case 1:
- for(int i = 0; i < getNumVertices(); i++) {
- getGradShapeFunction(i, u, v, w, s);
- MVertex *p = getVertex(i);
- jac[0][0] += p->x() * s[0]; jac[0][1] += p->y() * s[0]; jac[0][2] += p->z() * s[0];
+ break;
}
+ case 3:
{
- double a[3], b[3], c[3];
- a[0] = getVertex(1)->x() - getVertex(0)->x();
- a[1] = getVertex(1)->y() - getVertex(0)->y();
- a[2] = getVertex(1)->z() - getVertex(0)->z();
- if((fabs(a[0]) >= fabs(a[1]) && fabs(a[0]) >= fabs(a[2])) ||
- (fabs(a[1]) >= fabs(a[0]) && fabs(a[1]) >= fabs(a[2]))) {
- b[0] = a[1]; b[1] = -a[0]; b[2] = 0.;
- }
- else {
- b[0] = 0.; b[1] = a[2]; b[2] = -a[1];
- }
- prodve(a, b, c);
- jac[1][0] = b[0]; jac[1][1] = b[1]; jac[1][2] = b[2];
- jac[2][0] = c[0]; jac[2][1] = c[1]; jac[2][2] = c[2];
+ dJ = fabs(jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
+ jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
+ jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2]);
+ break;
}
- return sqrt(SQU(jac[0][0]) + SQU(jac[0][1]) + SQU(jac[0][2]));
- default:
- return 1.;
}
+ return dJ;
+}
+
+void MElement::pnt(double uu, double vv, double ww, SPoint3 &p) {
+
+
+ double x=0.;
+ double y=0.;
+ double z=0.;
+
+ const gmshFunctionSpace* fs = getFunctionSpace();
+
+ if (fs) {
+
+ double sf[256];
+ fs->f(uu,vv,ww,sf);
+
+ for (int j=0;j<getNumVertices();j++) {
+ const MVertex* v = getVertex(j);
+ x += sf[j] * v->x();
+ y += sf[j] * v->y();
+ z += sf[j] * v->z();
+ }
+ }
+ else Msg::Fatal("Could not find function space\n");
+
+ p = SPoint3(x,y,z);
+
+}
+
+void MElement::primaryPnt(double uu, double vv, double ww, SPoint3 &p) {
+
+ double x=0.;
+ double y=0.;
+ double z=0.;
+
+ const gmshFunctionSpace* fs = getFunctionSpace(1);
+
+ if (fs) {
+
+ double sf[256];
+ fs->f(uu,vv,ww,sf);
+ if (getNumPrimaryVertices() != 4) printf("Incorrect number of vertices %d\n",getNumPrimaryVertices()
+ );
+
+
+ for (int j=0;j<getNumPrimaryVertices();j++) {
+ const MVertex* v = getVertex(j);
+ x += sf[j] * v->x();
+ y += sf[j] * v->y();
+ z += sf[j] * v->z();
+ }
+ }
+
+ p = SPoint3(x,y,z);
+
}
void MElement::xyz2uvw(double xyz[3], double uvw[3])
@@ -424,21 +594,6 @@ double MElement::interpolateDiv(double val[], double u, double v, double w, int
return fx[0] + fy[1] + fz[2];
}
-// Expensive, generic version
-void MElement::pnt(double u, double v, double w, SPoint3 &p)
-{
- double x(0), y(0), z(0);
- for (int i = 0; i < getNumVertices(); i++) {
- double s;
- getShapeFunction(i, u, v, w, s);
- MVertex* v = getVertex(i);
- x += v->x() * s;
- y += v->y() * s;
- z += v->z() * s;
- }
- p = SPoint3(x, y, z);
-}
-
void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
int elementary, int physical)
{
@@ -737,339 +892,19 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
}
}
-void MLine::pnt(double u,double v,double w,SPoint3& p)
-{
- double f1, f2;
- getShapeFunction(0, u, v, w, f1);
- getShapeFunction(1, u, v, w, f2);
-
- double x = f1 * _v[0]->x() + f2 * _v[1]->x();
- double y = f1 * _v[0]->y() + f2 * _v[1]->y();
- double z = f1 * _v[0]->z() + f2 * _v[1]->z();
-
- p = SPoint3(x, y, z);
-}
-
-void MLine3::pnt(double u, double v, double w, SPoint3 &p)
-{
- double sf[3];
- gmshFunctionSpaces::find(MSH_LIN_3).f(u, v, w, sf);
- double x = sf[0] * _v[0]->x() + sf[1] * _v[1]->x() + sf[2] * _vs[0]->x();
- double y = sf[0] * _v[0]->y() + sf[1] * _v[1]->y() + sf[2] * _vs[0]->y();
- double z = sf[0] * _v[0]->z() + sf[1] * _v[1]->z() + sf[2] * _vs[0]->z();
-
- p = SPoint3(x,y,z);
-}
-
-void MLine3::getShapeFunction(int num,double uu, double vv,double ww,double& s)
-{
- if (num > 2) s = 0.;
- double sf[3];
- gmshFunctionSpaces::find(MSH_LIN_3).f(uu, vv, ww, sf);
- s = sf[num];
-}
-
-void MLine3::getGradShapeFunction(int num,double uu,double vv,double ww,double s[3])
-{
- if (num > 2) for (int i = 0; i < 3; i++) s[i] = 0.;
- double sf[3][3];
- gmshFunctionSpaces::find(MSH_LIN_3).df(uu, vv, ww, sf);
- for (int i = 0; i < 3; i++) s[i] = sf[num][i];
-}
-
-void MLineN::pnt(double uu, double vv, double ww, SPoint3 &p)
-{
- double sf[6];
-
- switch (getPolynomialOrder()) {
- case 1: gmshFunctionSpaces::find(MSH_LIN_2).f(uu, vv, ww, sf); break;
- case 2: gmshFunctionSpaces::find(MSH_LIN_3).f(uu, vv, ww, sf); break;
- case 3: gmshFunctionSpaces::find(MSH_LIN_4).f(uu, vv, ww, sf); break;
- case 4: gmshFunctionSpaces::find(MSH_LIN_5).f(uu, vv, ww, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_LIN_6).f(uu, vv, ww, sf); break;
- default:
- Msg::Error("Order %d line point interpolation not implemented",
- getPolynomialOrder());
- break;
- }
-
- double x = 0; for (int i = 0; i < 2; i++) x += sf[i] * _v[i]->x();
- double y = 0; for (int i = 0; i < 2; i++) y += sf[i] * _v[i]->y();
- double z = 0; for (int i = 0; i < 2; i++) z += sf[i] * _v[i]->z();
-
- for (int i = 0; i < _vs.size(); i++) x += sf[i+2] * _vs[i]->x();
- for (int i = 0; i < _vs.size(); i++) y += sf[i+2] * _vs[i]->y();
- for (int i = 0; i < _vs.size(); i++) z += sf[i+2] * _vs[i]->z();
-
- p = SPoint3(x, y, z);
-}
-
-void MLineN::getShapeFunction(int num, double uu, double vv, double ww, double &s)
-{
- double sf[6];
- switch (getPolynomialOrder()) {
- case 1: gmshFunctionSpaces::find(MSH_LIN_2).f(uu, vv, ww, sf); break;
- case 2: gmshFunctionSpaces::find(MSH_LIN_3).f(uu, vv, ww, sf); break;
- case 3: gmshFunctionSpaces::find(MSH_LIN_4).f(uu, vv, ww, sf); break;
- case 4: gmshFunctionSpaces::find(MSH_LIN_5).f(uu, vv, ww, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_LIN_6).f(uu, vv, ww, sf); break;
- default:
- Msg::Error("Order %d line shape function not provided",
- getPolynomialOrder());
- break;
- }
- s = sf[num];
-}
-
-void MLineN::getGradShapeFunction(int num, double uu, double vv, double ww, double s[3])
-{
- double sf[6][3];
- switch (getPolynomialOrder()) {
- case 1: gmshFunctionSpaces::find(MSH_LIN_2).df(uu, vv, ww, sf); break;
- case 2: gmshFunctionSpaces::find(MSH_LIN_3).df(uu, vv, ww, sf); break;
- case 3: gmshFunctionSpaces::find(MSH_LIN_4).df(uu, vv, ww, sf); break;
- case 4: gmshFunctionSpaces::find(MSH_LIN_5).df(uu, vv, ww, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_LIN_6).df(uu, vv, ww, sf); break;
- default:
- Msg::Error("Order %d line shape function gradient not implemented",
- getPolynomialOrder());
- break;
- }
- for (int i = 0; i < 3; i++) s[i] = sf[num][i];
-}
-
-void MTriangle::getShapeFunction(int num, double uu, double vv, double ww, double &s)
-{
- int nf = getNumFaceVertices();
- double fv[256];
-
- if (!nf){
- switch(getPolynomialOrder()){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, fv); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, fv); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu, vv, fv); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu, vv, fv); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu, vv, fv); break;
- default:
- Msg::Error("Order %d triangle shape function not implemented",
- getPolynomialOrder());
- break;
- }
- }
- else{
- switch(getPolynomialOrder()){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, fv); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, fv); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu, vv, fv); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu, vv, fv); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu, vv, fv); break;
- default:
- Msg::Error("Order %d triangle shape function not implemented",
- getPolynomialOrder());
- break;
- }
- }
- s = fv[num];
-}
-
-void MTriangle::getGradShapeFunction(int num, double uu, double vv, double ww, double s[3])
-{
- double grads[256][2];
- int nf = getNumFaceVertices();
-
- if (!nf){
- switch(getPolynomialOrder()){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, ww, grads); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, ww, grads); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_9).df(uu, vv, ww, grads); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_12).df(uu, vv, ww, grads); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_15I).df(uu, vv, ww, grads); break;
- default:
- Msg::Error("Order %d triangle shape function gradient not implemented",
- getPolynomialOrder());
- break;
- }
- }
- else{
- switch(getPolynomialOrder()){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, ww, grads); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, ww, grads); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_10).df(uu, vv, ww, grads); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_15).df(uu, vv, ww, grads); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_21).df(uu, vv, ww, grads); break;
- default:
- Msg::Error("Order %d triangle shape function gradient not implemented",
- getPolynomialOrder());
- break;
- }
- }
-
- for (int i = 0; i < 2; i++) s[i] = grads[num][i];
- s[2] = 0;
-}
-
-void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double ww,
- double j[2][3])
-{
- double grads[256][2];
- int nf = getNumFaceVertices();
-
- if (!nf){
- switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, ww, grads); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, ww, grads); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_9).df(uu, vv, ww, grads); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_12).df(uu, vv, ww, grads); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_15I).df(uu, vv, ww, grads); break;
- default: Msg::Error("Order %d triangle jac not implemented", ord); break;
- }
- }
- else{
- switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, ww, grads); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, ww, grads); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_10).df(uu, vv, ww, grads); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_15).df(uu, vv, ww, grads); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_21).df(uu, vv, ww, grads); break;
- default: Msg::Error("Order %d triangle jac not implemented", ord); break;
- }
- }
- j[0][0] = 0 ; for(int i = 0; i < 3; i++) j[0][0] += grads [i][0] * _v[i]->x();
- j[1][0] = 0 ; for(int i = 0; i < 3; i++) j[1][0] += grads [i][1] * _v[i]->x();
- j[0][1] = 0 ; for(int i = 0; i < 3; i++) j[0][1] += grads [i][0] * _v[i]->y();
- j[1][1] = 0 ; for(int i = 0; i < 3; i++) j[1][1] += grads [i][1] * _v[i]->y();
- j[0][2] = 0 ; for(int i = 0; i < 3; i++) j[0][2] += grads [i][0] * _v[i]->z();
- j[1][2] = 0 ; for(int i = 0; i < 3; i++) j[1][2] += grads [i][1] * _v[i]->z();
-
- if (ord == 1) return;
-
- for(int i = 3; i < 3 * ord + nf; i++) j[0][0] += grads[i][0] * vs[i - 3]->x();
- for(int i = 3; i < 3 * ord + nf; i++) j[1][0] += grads[i][1] * vs[i - 3]->x();
- for(int i = 3; i < 3 * ord + nf; i++) j[0][1] += grads[i][0] * vs[i - 3]->y();
- for(int i = 3; i < 3 * ord + nf; i++) j[1][1] += grads[i][1] * vs[i - 3]->y();
- for(int i = 3; i < 3 * ord + nf; i++) j[0][2] += grads[i][0] * vs[i - 3]->z();
- for(int i = 3; i < 3 * ord + nf; i++) j[1][2] += grads[i][1] * vs[i - 3]->z();
-}
-
-void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, double ww,
- SPoint3 &p)
-{
- double sf[256];
- int nf = getNumFaceVertices();
-
- if (!nf){
- switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu, vv, sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu, vv, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu, vv, sf); break;
- default: Msg::Error("Order %d triangle pnt not implemented", ord); break;
- }
- }
- else{
- switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu, vv, sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu, vv, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu, vv, sf); break;
- default: Msg::Error("Order %d triangle pnt not implemented", ord); break;
- }
- }
-
- double x = 0 ; for(int i = 0; i < 3; i++) x += sf[i] * _v[i]->x();
- double y = 0 ; for(int i = 0; i < 3; i++) y += sf[i] * _v[i]->y();
- double z = 0 ; for(int i = 0; i < 3; i++) z += sf[i] * _v[i]->z();
-
- for(int i = 3; i < 3 * ord + nf; i++) x += sf[i] * vs[i - 3]->x();
- for(int i = 3; i < 3 * ord + nf; i++) y += sf[i] * vs[i - 3]->y();
- for(int i = 3; i < 3 * ord + nf; i++) z += sf[i] * vs[i - 3]->z();
-
- p = SPoint3(x, y, z);
-}
+const gmshFunctionSpace* MLine::getFunctionSpace(int o) const {
-void MTetrahedron::jac(int ord, MVertex *vs[], double uu, double vv, double ww,
- double j[3][3])
-{
- double grads[256][3];
- switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TET_4) .df(uu, vv, ww, grads); break;
- case 2: gmshFunctionSpaces::find(MSH_TET_10).df(uu, vv, ww, grads); break;
- case 3: gmshFunctionSpaces::find(MSH_TET_20).df(uu, vv, ww, grads); break;
- case 4: gmshFunctionSpaces::find(MSH_TET_35).df(uu, vv, ww, grads); break;
- case 5: gmshFunctionSpaces::find(MSH_TET_56).df(uu, vv, ww, grads); break;
- default: Msg::Error("Order %d tetrahedron jac not implemented", ord); break;
- }
-
- j[0][0] = 0 ; for(int i = 0; i < 4; i++) j[0][0] += grads [i][0] * _v[i]->x();
- j[1][0] = 0 ; for(int i = 0; i < 4; i++) j[1][0] += grads [i][1] * _v[i]->x();
- j[2][0] = 0 ; for(int i = 0; i < 4; i++) j[2][0] += grads [i][2] * _v[i]->x();
-
- j[0][1] = 0 ; for(int i = 0; i < 4; i++) j[0][1] += grads [i][0] * _v[i]->y();
- j[1][1] = 0 ; for(int i = 0; i < 4; i++) j[1][1] += grads [i][1] * _v[i]->y();
- j[2][1] = 0 ; for(int i = 0; i < 4; i++) j[2][1] += grads [i][2] * _v[i]->y();
-
- j[0][2] = 0 ; for(int i = 0; i < 4; i++) j[0][2] += grads [i][0] * _v[i]->z();
- j[1][2] = 0 ; for(int i = 0; i < 4; i++) j[1][2] += grads [i][1] * _v[i]->z();
- j[2][2] = 0 ; for(int i = 0; i < 4; i++) j[2][2] += grads [i][2] * _v[i]->z();
+ int order = o == -1 ? getPolynomialOrder() : o;
- if (ord == 1) return;
-
- const int N = (ord+1)*(ord+2)*(ord+3)/6;
-
- for(int i = 4; i < N; i++) j[0][0] += grads[i][0] * vs[i - 4]->x();
- for(int i = 4; i < N; i++) j[1][0] += grads[i][1] * vs[i - 4]->x();
- for(int i = 4; i < N; i++) j[2][0] += grads[i][2] * vs[i - 4]->x();
-
- for(int i = 4; i < N; i++) j[0][1] += grads[i][0] * vs[i - 4]->y();
- for(int i = 4; i < N; i++) j[1][1] += grads[i][1] * vs[i - 4]->y();
- for(int i = 4; i < N; i++) j[2][1] += grads[i][2] * vs[i - 4]->y();
-
- for(int i = 4; i < N; i++) j[0][2] += grads[i][0] * vs[i - 4]->z();
- for(int i = 4; i < N; i++) j[1][2] += grads[i][1] * vs[i - 4]->z();
- for(int i = 4; i < N; i++) j[2][2] += grads[i][2] * vs[i - 4]->z();
-}
-
-void MTetrahedron::pnt(int ord, MVertex *vs[], double uu, double vv, double ww,
- SPoint3 &p)
-{
- double sf[256];
-
- int nv = getNumVolumeVertices();
-
- if (!nv){
- switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TET_4).f(uu, vv, ww, sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TET_10).f(uu, vv, ww, sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TET_20).f(uu, vv, ww, sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TET_34).f(uu, vv, ww, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TET_52).f(uu, vv, ww, sf); break;
- default: Msg::Error("Order %d tetrahedron pnt not implemented", ord); break;
- }
+ switch (order) {
+ case 1: return &gmshFunctionSpaces::find(MSH_LIN_2); break;
+ case 2: return &gmshFunctionSpaces::find(MSH_LIN_3); break;
+ case 3: return &gmshFunctionSpaces::find(MSH_LIN_4); break;
+ case 4: return &gmshFunctionSpaces::find(MSH_LIN_5); break;
+ case 5: return &gmshFunctionSpaces::find(MSH_LIN_6); break;
+ default: Msg::Error("Order %d line point interpolation not implemented", getPolynomialOrder()); break;
}
- else{
- switch(ord){
- case 4: gmshFunctionSpaces::find(MSH_TET_35).f(uu, vv, ww, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TET_56).f(uu, vv, ww, sf); break;
- default: Msg::Error("Order %d tetrahedron pnt not implemented", ord); break;
- }
- }
-
- double x = 0 ; for(int i = 0; i < 4; i++) x += sf[i] * _v[i]->x();
- double y = 0 ; for(int i = 0; i < 4; i++) y += sf[i] * _v[i]->y();
- double z = 0 ; for(int i = 0; i < 4; i++) z += sf[i] * _v[i]->z();
-
- int nbComplete = (ord + 1) * (ord + 2) * (ord + 3) / 6;
- int nbInterior = (ord - 3) * (ord - 2) * (ord - 1) / 6;
-
- const int N = nv ? nbComplete : nbComplete - nbInterior;
-
- for(int i = 4; i < N; i++) x += sf[i] * vs[i - 4]->x();
- for(int i = 4; i < N; i++) y += sf[i] * vs[i - 4]->y();
- for(int i = 4; i < N; i++) z += sf[i] * vs[i - 4]->z();
-
- p = SPoint3(x, y, z);
+ return NULL;
}
const int numSubEdges = 6;
@@ -1096,22 +931,22 @@ void MTriangleN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *
const double d = 1. / numSubEdges;
SPoint3 pnt1, pnt2, pnt3;
- double J1[2][3], J2[2][3], J3[2][3];
+ double J1[3][3], J2[3][3], J3[3][3];
if (ix % 2 == 0){
pnt(ix / 2 * d, iy * d, 0, pnt1);
pnt((ix / 2 + 1) * d, iy * d, 0, pnt2);
pnt(ix / 2 * d, (iy + 1) * d, 0, pnt3);
- jac(ix / 2 * d, iy * d, 0, J1);
- jac((ix / 2 + 1) * d, iy * d, 0, J2);
- jac(ix / 2 * d, (iy + 1) * d, 0, J3);
+ getJacobian(ix / 2 * d, iy * d, 0, J1);
+ getJacobian((ix / 2 + 1) * d, iy * d, 0, J2);
+ getJacobian(ix / 2 * d, (iy + 1) * d, 0, J3);
}
else{
pnt((ix / 2 + 1) * d, iy * d, 0, pnt1);
pnt((ix / 2 + 1) * d, (iy + 1) * d, 0, pnt2);
pnt(ix / 2 * d, (iy + 1) * d, 0, pnt3);
- jac((ix / 2 + 1) * d, iy * d, 0, J1);
- jac((ix / 2 + 1) * d, (iy + 1) * d, 0, J2);
- jac(ix / 2 * d, (iy + 1) * d, 0, J3);
+ getJacobian((ix / 2 + 1) * d, iy * d, 0, J1);
+ getJacobian((ix / 2 + 1) * d, (iy + 1) * d, 0, J2);
+ getJacobian(ix / 2 * d, (iy + 1) * d, 0, J3);
}
{
SVector3 d1(J1[0][0], J1[0][1], J1[0][2]);
@@ -1141,11 +976,11 @@ int MTriangleN::getNumEdgesRep(){ return 3 * numSubEdges; }
void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
- n[0] = n[1] = getFace(0).normal();
+ n[0] = n[1] = n[2] = getFace(0).normal();
int N = getNumEdgesRep() / 3;
if (num < N){
SPoint3 pnt1, pnt2;
- pnt((double)num / N, 0., 0,pnt1);
+ pnt((double)num / N, 0., 0.,pnt1);
pnt((double)(num + 1) / N, 0., 0, pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
@@ -1173,83 +1008,6 @@ void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
}
}
-void MTetrahedron10::getShapeFunction(int num, double uu, double vv, double ww, double &s)
-{
- double sf[256];
- gmshFunctionSpaces::find(MSH_TET_10).f(uu, vv, ww, sf);
- s = sf[num];
-}
-
-void MTetrahedron10::getGradShapeFunction(int num, double uu, double vv, double ww,
- double gs[3])
-{
- double gsf[256][3];
- gmshFunctionSpaces::find(MSH_TET_10).df(uu, vv, ww, gsf);
- for (int i = 0; i < 3; i++) gs[i] = gsf[num][i];
-}
-
-void MTetrahedronN::getShapeFunction(int num, double uu, double vv, double ww, double &s)
-{
- double sf[256];
-
- int nv = getNumVolumeVertices();
-
- if (!nv){
- switch(getPolynomialOrder()){
- case 1: gmshFunctionSpaces::find(MSH_TET_4).f(uu, vv, ww, sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TET_10).f(uu, vv, ww, sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TET_20).f(uu, vv, ww, sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TET_34).f(uu, vv, ww, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TET_52).f(uu, vv, ww, sf); break;
- default:
- Msg::Error("Order %d tetrahedron pnt not implemented", getPolynomialOrder());
- break;
- }
- }
- else{
- switch(getPolynomialOrder()){
- case 4: gmshFunctionSpaces::find(MSH_TET_35).f(uu, vv, ww, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TET_56).f(uu, vv, ww, sf); break;
- default:
- Msg::Error("Order %d tetrahedron pnt not implemented", getPolynomialOrder());
- break;
- }
- }
-
- s = sf[num];
-}
-
-void MTetrahedronN::getGradShapeFunction(int num, double uu, double vv, double ww,
- double gs[3])
-{
- double gsf[256][3];
-
- int nv = getNumVolumeVertices();
-
- if (!nv){
- switch(getPolynomialOrder()){
- case 1: gmshFunctionSpaces::find(MSH_TET_4) .df(uu, vv, ww, gsf); break;
- case 2: gmshFunctionSpaces::find(MSH_TET_10).df(uu, vv, ww, gsf); break;
- case 3: gmshFunctionSpaces::find(MSH_TET_20).df(uu, vv, ww, gsf); break;
- case 4: gmshFunctionSpaces::find(MSH_TET_34).df(uu, vv, ww, gsf); break;
- case 5: gmshFunctionSpaces::find(MSH_TET_52).df(uu, vv, ww, gsf); break;
- default:
- Msg::Error("Order %d tetrahedron pnt not implemented", getPolynomialOrder());
- break;
- }
- }
- else{
- switch(getPolynomialOrder()){
- case 4: gmshFunctionSpaces::find(MSH_TET_35).df(uu, vv, ww, gsf); break;
- case 5: gmshFunctionSpaces::find(MSH_TET_56).df(uu, vv, ww, gsf); break;
- default:
- Msg::Error("Order %d tetrahedron pnt not implemented", getPolynomialOrder());
- break;
- }
- }
- for (int i = 0; i < 3; i++) gs[i] = gsf[num][i];
-}
-
int MTetrahedronN::getNumEdgesRep(){ return 6 * numSubEdges; }
void MTetrahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -1300,7 +1058,7 @@ void MTetrahedronN::getFaceRep(int num, double *x, double *y, double *z, SVector
0 1 2 3
0 1 2 3 4
0 1 2 3 4 5
- */
+ */
// on the first layer, we have (numSubEdges-1) * 2 + 1 triangles
// on the second layer, we have (numSubEdges-2) * 2 + 1 triangles
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 1c5d2111c1d0bb1ec3960c92a6f36782835025b2..785e14ec06808a61c0f5e867d062df64021d54d1 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -14,6 +14,7 @@
#include "MEdge.h"
#include "MFace.h"
#include "Message.h"
+#include "FunctionSpace.h"
struct IntPt{
double pt[3];
@@ -66,7 +67,7 @@ class MElement
virtual int getDim() = 0;
// return the polynomial order the element
- virtual int getPolynomialOrder(){ return 1; }
+ virtual int getPolynomialOrder() const { return 1; }
// get/set the partition to which the element belongs
virtual int getPartition(){ return _partition; }
@@ -77,7 +78,7 @@ class MElement
virtual void setVisibility(char val){ _visible = val; }
// get the vertices
- virtual int getNumVertices() = 0;
+ virtual int getNumVertices() const = 0;
virtual MVertex *getVertex(int num) = 0;
// get the vertex using the I-deas UNV ordering
@@ -94,15 +95,17 @@ class MElement
// get the number of vertices associated with edges, faces and
// volumes (nonzero only for higher order elements)
- virtual int getNumEdgeVertices(){ return 0; }
- virtual int getNumFaceVertices(){ return 0; }
- virtual int getNumVolumeVertices(){ return 0; }
+ virtual int getNumEdgeVertices() const { return 0; }
+ virtual int getNumFaceVertices() const { return 0; }
+ virtual int getNumVolumeVertices() const { return 0; }
// get the number of primary vertices (first-order element)
int getNumPrimaryVertices()
{
- return getNumVertices() - getNumEdgeVertices() -
- getNumFaceVertices() - getNumVolumeVertices();
+ return getNumVertices() -
+ getNumEdgeVertices() -
+ getNumFaceVertices() -
+ getNumVolumeVertices();
}
// get the edges
@@ -138,6 +141,7 @@ class MElement
virtual double minEdge();
// get the quality measures
+
virtual double rhoShapeMeasure();
virtual double gammaShapeMeasure(){ return 0.; }
virtual double etaShapeMeasure(){ return 0.; }
@@ -154,26 +158,32 @@ class MElement
virtual double getVolume(){ return 0.; }
virtual int getVolumeSign(){ return 1; }
virtual void setVolumePositive(){ if(getVolumeSign() < 0) revert(); }
-
+
// return an information string for the element
virtual std::string getInfoString();
+ virtual const gmshFunctionSpace* getFunctionSpace(int ord = -1) const {
+ Msg::Fatal("Function space not implemented for element %s",getStringForPOS());
+ return NULL;
+ }
+
// return the interpolating nodal shape function associated with
// node num, evaluated at point (u,v,w) in parametric coordinates
- virtual void getShapeFunction(int num, double u, double v, double w,
- double &s) = 0;
+ virtual void getShapeFunction(int num, double u, double v, double w,double &s);
// return the gradient of of the nodal shape function associated
// with node num, evaluated at point (u,v,w) in parametric
// coordinates
- virtual void getGradShapeFunction(int num, double u, double v, double w,
- double s[3]) = 0;
+ virtual void getGradShapeFunction(int num, double u, double v, double w,double s[3]);
// return the Jacobian of the element evaluated at point (u,v,w) in
// parametric coordinates
+
double getJacobian(double u, double v, double w, double jac[3][3]);
-
- virtual void pnt(double u, double v, double w, SPoint3 &p);
+ double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
+
+ virtual void pnt (double u,double v,double w,SPoint3 &p);
+ virtual void primaryPnt(double u,double v,double w,SPoint3 &p);
// invert the parametrisation
virtual void xyz2uvw(double xyz[3], double uvw[3]);
@@ -209,11 +219,11 @@ class MElement
// info for specific IO formats (returning 0 means that the element
// is not implemented in that format)
- virtual int getTypeForMSH(){ return 0; }
- virtual int getTypeForUNV(){ return 0; }
- virtual int getTypeForVTK(){ return 0; }
- virtual const char *getStringForPOS(){ return 0; }
- virtual const char *getStringForBDF(){ return 0; }
+ virtual int getTypeForMSH() const { return 0; }
+ virtual int getTypeForUNV() const { return 0; }
+ virtual int getTypeForVTK() const { return 0; }
+ virtual const char *getStringForPOS() const { return 0; }
+ virtual const char *getStringForBDF() const { return 0; }
// return the number of vertices, as well as the element name if
// 'name' != 0
@@ -261,7 +271,7 @@ class MPoint : public MElement {
}
~MPoint(){}
virtual int getDim(){ return 0; }
- virtual int getNumVertices(){ return 1; }
+ virtual int getNumVertices() const { return 1; }
virtual MVertex *getVertex(int num){ return _v[0]; }
virtual int getNumEdges(){ return 0; }
virtual MEdge getEdge(int num){ return MEdge(); }
@@ -271,8 +281,8 @@ class MPoint : public MElement {
virtual MFace getFace(int num){ return MFace(); }
virtual int getNumFacesRep(){ return 0; }
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n){}
- virtual int getTypeForMSH(){ return MSH_PNT; }
- virtual const char *getStringForPOS(){ return "SP"; }
+ virtual int getTypeForMSH() const { return MSH_PNT; }
+ virtual const char *getStringForPOS() const { return "SP"; }
virtual void getShapeFunction(int num, double u, double v, double w, double &s)
{
s = 1.;
@@ -314,7 +324,7 @@ class MLine : public MElement {
}
~MLine(){}
virtual int getDim(){ return 1; }
- virtual int getNumVertices(){ return 2; }
+ virtual int getNumVertices() const { return 2; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual int getNumEdges(){ return 1; }
virtual MEdge getEdge(int num){ return MEdge(_v[0], _v[1]); }
@@ -332,38 +342,26 @@ class MLine : public MElement {
virtual MFace getFace(int num){ return MFace(); }
virtual int getNumFacesRep(){ return 0; }
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n){}
- virtual int getTypeForMSH(){ return MSH_LIN_2; }
- virtual int getTypeForUNV(){ return 21; } // linear beam
- virtual int getTypeForVTK(){ return 3; }
- virtual const char *getStringForPOS(){ return "SL"; }
- virtual const char *getStringForBDF(){ return "CBAR"; }
+ virtual int getTypeForMSH() const { return MSH_LIN_2; }
+ virtual int getTypeForUNV() const { return 21; } // linear beam
+ virtual int getTypeForVTK() const { return 3; }
+ virtual const char *getStringForPOS() const { return "SL"; }
+ virtual const char *getStringForBDF() const { return "CBAR"; }
virtual void revert()
{
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
}
- virtual void getShapeFunction(int num, double u, double v, double w, double &s)
- {
- switch(num) {
- case 0 : s = 0.5 * (1. - u); break;
- case 1 : s = 0.5 * (1. + u); break;
- default : s = 0.; break;
- }
- }
- virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
- {
- switch(num) {
- case 0 : s[0] = -0.5; s[1] = 0.; s[2] = 0.; break;
- case 1 : s[0] = 0.5; s[1] = 0.; s[2] = 0.; break;
- default : s[0] = s[1] = s[2] = 0.; break;
- }
- }
- virtual void pnt(double u, double v, double w, SPoint3 &p);
+
+
+ virtual const gmshFunctionSpace* getFunctionSpace(int=-1) const;
+
virtual bool isInside(double u, double v, double w, double tol=1.e-8)
{
if(u < -(1. + tol) || u > (1. + tol))
return false;
return true;
}
+
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
};
@@ -390,15 +388,15 @@ class MLine3 : public MLine {
_vs[0]->setPolynomialOrder(2);
}
~MLine3(){}
- virtual int getPolynomialOrder(){ return 2; }
- virtual int getNumVertices(){ return 3; }
+ virtual int getPolynomialOrder() const { return 2; }
+ virtual int getNumVertices() const { return 3; }
virtual MVertex *getVertex(int num){ return num < 2 ? _v[num] : _vs[num - 2]; }
virtual MVertex *getVertexUNV(int num)
{
static const int map[3] = {0, 2, 1};
return getVertex(map[num]);
}
- virtual int getNumEdgeVertices(){ return 1; }
+ virtual int getNumEdgeVertices() const { return 1; }
virtual int getNumEdgesRep(){ return 2; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -414,13 +412,11 @@ class MLine3 : public MLine {
MLine::_getEdgeVertices(v);
v[2] = _vs[0];
}
- virtual void getShapeFunction(int num, double u, double v, double w, double &s);
- virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]);
- virtual void pnt(double u, double v, double w, SPoint3 &p);
- virtual int getTypeForMSH(){ return MSH_LIN_3; }
- virtual int getTypeForUNV(){ return 24; } // parabolic beam
- virtual int getTypeForVTK(){ return 21; }
- virtual const char *getStringForPOS(){ return "SL2"; }
+
+ virtual int getTypeForMSH() const { return MSH_LIN_3; }
+ virtual int getTypeForUNV() const { return 24; } // parabolic beam
+ virtual int getTypeForVTK() const { return 21; }
+ virtual const char *getStringForPOS() const { return "SL2"; }
};
/*
@@ -448,10 +444,10 @@ class MLineN : public MLine {
_vs[i]->setPolynomialOrder(_vs.size() + 1);
}
~MLineN(){}
- virtual int getPolynomialOrder(){ return _vs.size() + 1; }
- virtual int getNumVertices(){ return _vs.size() + 2; }
+ virtual int getPolynomialOrder() const { return _vs.size() + 1; }
+ virtual int getNumVertices() const { return _vs.size() + 2; }
virtual MVertex *getVertex(int num){ return num < 2 ? _v[num] : _vs[num - 2]; }
- virtual int getNumEdgeVertices(){ return _vs.size(); }
+ virtual int getNumEdgeVertices() const { return _vs.size(); }
virtual int getNumEdgesRep(){ return _vs.size() + 1; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -465,10 +461,7 @@ class MLineN : public MLine {
MLine::_getEdgeVertices(v);
for(unsigned int i = 0; i != _vs.size(); ++i) v[i+2] = _vs[i];
}
- virtual void getShapeFunction(int num, double u, double v, double w, double &s);
- virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]);
- virtual void pnt(double u, double v, double w, SPoint3 &p);
- virtual int getTypeForMSH()
+ virtual int getTypeForMSH() const
{
if(_vs.size() == 2) return MSH_LIN_4;
if(_vs.size() == 3) return MSH_LIN_5;
@@ -521,7 +514,7 @@ class MTriangle : public MElement {
virtual int getDim(){ return 2; }
virtual double gammaShapeMeasure();
virtual double distoShapeMeasure();
- virtual int getNumVertices(){ return 3; }
+ virtual int getNumVertices() const { return 3; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual MVertex *getVertexMED(int num)
{
@@ -566,35 +559,26 @@ class MTriangle : public MElement {
v.resize(3);
_getFaceVertices(v);
}
- virtual int getTypeForMSH(){ return MSH_TRI_3; }
- virtual int getTypeForUNV(){ return 91; } // thin shell linear triangle
- virtual int getTypeForVTK(){ return 5; }
- virtual const char *getStringForPOS(){ return "ST"; }
- virtual const char *getStringForBDF(){ return "CTRIA3"; }
+ virtual int getTypeForMSH() const { return MSH_TRI_3; }
+ virtual int getTypeForUNV() const { return 91; } // thin shell linear triangle
+ virtual int getTypeForVTK() const { return 5; }
+ virtual const char *getStringForPOS() const { return "ST"; }
+ virtual const char *getStringForBDF() const { return "CTRIA3"; }
virtual void revert()
{
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
}
- virtual void getShapeFunction(int num, double u, double v, double w, double &s);
- virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]);
- virtual void jac(int order, MVertex *verts[], double u, double v, double w,
- double j[2][3]);
- virtual void jac(double u, double v, double w, double j[2][3])
- {
- jac(1, 0, u, v, w, j);
- }
- virtual void pnt(int order, MVertex *verts[], double u, double v, double w,
- SPoint3 &p);
- virtual void pnt(double u, double v, double w, SPoint3 &p)
- {
- pnt(1, 0, u, v, w, p);
- }
+
+
+ virtual const gmshFunctionSpace* getFunctionSpace(int=-1) const;
+
virtual bool isInside(double u, double v, double w, double tol=1.e-8)
{
if(u < (-tol) || v < (-tol) || u > ((1. + tol) - v))
return false;
return true;
}
+
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
virtual SPoint3 circumcenter();
private:
@@ -639,8 +623,8 @@ class MTriangle6 : public MTriangle {
for(int i = 0; i < 3; i++) _vs[i]->setPolynomialOrder(2);
}
~MTriangle6(){}
- virtual int getPolynomialOrder(){ return 2; }
- virtual int getNumVertices(){ return 6; }
+ virtual int getPolynomialOrder() const { return 2; }
+ virtual int getNumVertices() const { return 6; }
virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
virtual MVertex *getVertexUNV(int num)
{
@@ -652,7 +636,7 @@ class MTriangle6 : public MTriangle {
static const int map[6] = {0, 2, 1, 5, 4, 3};
return getVertex(map[num]);
}
- virtual int getNumEdgeVertices(){ return 3; }
+ virtual int getNumEdgeVertices() const { return 3; }
virtual int getNumEdgesRep(){ return 6; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -686,25 +670,17 @@ class MTriangle6 : public MTriangle {
v[4] = _vs[1];
v[5] = _vs[2];
}
- virtual int getTypeForMSH(){ return MSH_TRI_6; }
- virtual int getTypeForUNV(){ return 92; } // thin shell parabolic triangle
- virtual int getTypeForVTK(){ return 22; }
- virtual const char *getStringForPOS(){ return "ST2"; }
- virtual const char *getStringForBDF(){ return "CTRIA6"; }
+ virtual int getTypeForMSH() const { return MSH_TRI_6; }
+ virtual int getTypeForUNV() const { return 92; } // thin shell parabolic triangle
+ virtual int getTypeForVTK() const { return 22; }
+ virtual const char *getStringForPOS() const { return "ST2"; }
+ virtual const char *getStringForBDF() const { return "CTRIA6"; }
virtual void revert()
{
MVertex *tmp;
tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
tmp = _vs[0]; _vs[0] = _vs[2]; _vs[2] = tmp;
}
- virtual void jac(double u, double v, double w, double j[2][3])
- {
- MTriangle::jac(2, _vs, u, v, w, j);
- }
- virtual void pnt(double u, double v, double w, SPoint3 &p)
- {
- MTriangle::pnt(2, _vs, u, v, w, p);
- }
};
/*
@@ -750,10 +726,10 @@ class MTriangleN : public MTriangle {
for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
}
~MTriangleN(){}
- virtual int getPolynomialOrder(){ return _order; }
- virtual int getNumVertices(){ return 3 + _vs.size(); }
+ virtual int getPolynomialOrder() const { return _order; }
+ virtual int getNumVertices() const { return 3 + _vs.size(); }
virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
- virtual int getNumFaceVertices()
+ virtual int getNumFaceVertices() const
{
if(_order == 3 && _vs.size() == 6) return 0;
if(_order == 3 && _vs.size() == 7) return 1;
@@ -763,17 +739,17 @@ class MTriangleN : public MTriangle {
if(_order == 5 && _vs.size() == 18) return 6;
return 0;
}
- virtual int getNumEdgeVertices() const { return _order - 1; }
+ virtual int getNumEdgeVertices() const { return 3 * (_order - 1); }
virtual int getNumEdgesRep();
virtual int getNumFacesRep();
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
{
- v.resize(2 + getNumEdgeVertices());
+ v.resize(_order + 1);
MTriangle::_getEdgeVertices(num, v);
int j = 2;
- const int ie = (num + 1) * getNumEdgeVertices();
- for(int i = num*getNumEdgeVertices(); i != ie; ++i) v[j++] = _vs[i];
+ const int ie = (num + 1) * (_order - 1);
+ for(int i = num * (_order-1); i != ie; ++i) v[j++] = _vs[i];
}
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
@@ -782,7 +758,7 @@ class MTriangleN : public MTriangle {
MTriangle::_getFaceVertices(v);
for(unsigned int i = 0; i != _vs.size(); ++i) v[i + 3] = _vs[i];
}
- virtual int getTypeForMSH()
+ virtual int getTypeForMSH() const
{
if(_order == 3 && _vs.size() == 6) return MSH_TRI_9;
if(_order == 3 && _vs.size() == 7) return MSH_TRI_10;
@@ -800,14 +776,6 @@ class MTriangleN : public MTriangle {
inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
_vs = inv;
}
- virtual void jac(double u, double v, double w, double j[2][3])
- {
- MTriangle::jac(_order, &_vs[0], u, v, w, j);
- }
- virtual void pnt(double u, double v, double w, SPoint3 &p)
- {
- MTriangle::pnt(_order, &_vs[0], u, v, w, p);
- }
};
/*
@@ -853,7 +821,7 @@ class MQuadrangle : public MElement {
}
~MQuadrangle(){}
virtual int getDim(){ return 2; }
- virtual int getNumVertices(){ return 4; }
+ virtual int getNumVertices() const { return 4; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual MVertex *getVertexMED(int num)
{
@@ -892,35 +860,15 @@ class MQuadrangle : public MElement {
v.resize(4);
_getFaceVertices(v);
}
- virtual int getTypeForMSH(){ return MSH_QUA_4; }
- virtual int getTypeForUNV(){ return 94; } // thin shell linear quadrilateral
- virtual int getTypeForVTK(){ return 8; }
- virtual const char *getStringForPOS(){ return "SQ"; }
- virtual const char *getStringForBDF(){ return "CQUAD4"; }
+ virtual int getTypeForMSH() const { return MSH_QUA_4; }
+ virtual int getTypeForUNV() const { return 94; } // thin shell linear quadrilateral
+ virtual int getTypeForVTK() const { return 8; }
+ virtual const char *getStringForPOS() const { return "SQ"; }
+ virtual const char *getStringForBDF() const { return "CQUAD4"; }
virtual void revert()
{
MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
}
- virtual void getShapeFunction(int num, double u, double v, double w, double &s)
- {
- switch(num) {
- case 0 : s = 0.25 * (1. - u) * (1. - v); break;
- case 1 : s = 0.25 * (1. + u) * (1. - v); break;
- case 2 : s = 0.25 * (1. + u) * (1. + v); break;
- case 3 : s = 0.25 * (1. - u) * (1. + v); break;
- default : s = 0.; break;
- }
- }
- virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
- {
- switch(num) {
- case 0 : s[0] = -0.25 * (1. - v); s[1] = -0.25 * (1. - u); s[2] = 0.; break;
- case 1 : s[0] = 0.25 * (1. - v); s[1] = -0.25 * (1. + u); s[2] = 0.; break;
- case 2 : s[0] = 0.25 * (1. + v); s[1] = 0.25 * (1. + u); s[2] = 0.; break;
- case 3 : s[0] = -0.25 * (1. + v); s[1] = 0.25 * (1. - u); s[2] = 0.; break;
- default : s[0] = s[1] = s[2] = 0.; break;
- }
- }
virtual bool isInside(double u, double v, double w, double tol=1.e-8)
{
if(u < -(1. + tol) || v < -(1. + tol) || u > (1. + tol) || v > (1. + tol))
@@ -971,8 +919,8 @@ class MQuadrangle8 : public MQuadrangle {
for(int i = 0; i < 4; i++) _vs[i]->setPolynomialOrder(2);
}
~MQuadrangle8(){}
- virtual int getPolynomialOrder(){ return 2; }
- virtual int getNumVertices(){ return 8; }
+ virtual int getPolynomialOrder() const { return 2; }
+ virtual int getNumVertices() const { return 8; }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
virtual MVertex *getVertexUNV(int num)
{
@@ -984,7 +932,7 @@ class MQuadrangle8 : public MQuadrangle {
static const int map[8] = {0, 3, 2, 1, 7, 6, 5, 4};
return getVertex(map[num]);
}
- virtual int getNumEdgeVertices(){ return 4; }
+ virtual int getNumEdgeVertices() const { return 4; }
virtual int getNumEdgesRep(){ return 8; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -1020,10 +968,10 @@ class MQuadrangle8 : public MQuadrangle {
v[6] = _vs[2];
v[7] = _vs[3];
}
- virtual int getTypeForMSH(){ return MSH_QUA_8; }
- virtual int getTypeForUNV(){ return 95; } // shell parabolic quadrilateral
- virtual int getTypeForVTK(){ return 23; }
- virtual const char *getStringForBDF(){ return "CQUAD8"; }
+ virtual int getTypeForMSH() const { return MSH_QUA_8; }
+ virtual int getTypeForUNV() const { return 95; } // shell parabolic quadrilateral
+ virtual int getTypeForVTK() const { return 23; }
+ virtual const char *getStringForBDF() const { return "CQUAD8"; }
virtual void revert()
{
MVertex *tmp;
@@ -1063,11 +1011,11 @@ class MQuadrangle9 : public MQuadrangle {
for(int i = 0; i < 5; i++) _vs[i]->setPolynomialOrder(2);
}
~MQuadrangle9(){}
- virtual int getPolynomialOrder(){ return 2; }
- virtual int getNumVertices(){ return 9; }
+ virtual int getPolynomialOrder() const { return 2; }
+ virtual int getNumVertices() const { return 9; }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
- virtual int getNumEdgeVertices(){ return 4; }
- virtual int getNumFaceVertices(){ return 1; }
+ virtual int getNumEdgeVertices() const { return 4; }
+ virtual int getNumFaceVertices() const { return 6; }
virtual int getNumEdgesRep(){ return 8; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -1105,8 +1053,8 @@ class MQuadrangle9 : public MQuadrangle {
v[7] = _vs[3];
v[8] = _vs[4];
}
- virtual int getTypeForMSH(){ return MSH_QUA_9; }
- virtual const char *getStringForPOS(){ return "SQ2"; }
+ virtual int getTypeForMSH() const { return MSH_QUA_9; }
+ virtual const char *getStringForPOS() const { return "SQ2"; }
virtual void revert()
{
MVertex *tmp;
@@ -1166,7 +1114,7 @@ class MTetrahedron : public MElement {
}
~MTetrahedron(){}
virtual int getDim(){ return 3; }
- virtual int getNumVertices(){ return 4; }
+ virtual int getNumVertices() const { return 4; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual MVertex *getVertexMED(int num)
{
@@ -1208,11 +1156,11 @@ class MTetrahedron : public MElement {
v.resize(3);
_getFaceVertices(num, v);
}
- virtual int getTypeForMSH(){ return MSH_TET_4; }
- virtual int getTypeForUNV(){ return 111; } // solid linear tetrahedron
- virtual int getTypeForVTK(){ return 10; }
- virtual const char *getStringForPOS(){ return "SS"; }
- virtual const char *getStringForBDF(){ return "CTETRA"; }
+ virtual int getTypeForMSH() const { return MSH_TET_4; }
+ virtual int getTypeForUNV() const { return 111; } // solid linear tetrahedron
+ virtual int getTypeForVTK() const { return 10; }
+ virtual const char *getStringForPOS() const { return "SS"; }
+ virtual const char *getStringForBDF() const { return "CTETRA"; }
virtual void revert()
{
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
@@ -1235,54 +1183,16 @@ class MTetrahedron : public MElement {
virtual double distoShapeMeasure();
virtual double etaShapeMeasure();
void xyz2uvw(double xyz[3], double uvw[3]);
- virtual void getShapeFunction(int num, double u, double v, double w, double &s)
- {
- switch(num) {
- case 0 : s = 1. - u - v - w; break;
- case 1 : s = u ; break;
- case 2 : s = v ; break;
- case 3 : s = w; break;
- default : s = 0.; break;
- }
- }
- virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
- {
- switch(num) {
- case 0 : s[0] = -1.; s[1] = -1.; s[2] = -1.; break;
- case 1 : s[0] = 1.; s[1] = 0.; s[2] = 0.; break;
- case 2 : s[0] = 0.; s[1] = 1.; s[2] = 0.; break;
- case 3 : s[0] = 0.; s[1] = 0.; s[2] = 1.; break;
- default : s[0] = s[1] = s[2] = 0.; break;
- }
- }
+
+ virtual const gmshFunctionSpace* getFunctionSpace(int=-1) const;
+
virtual bool isInside(double u, double v, double w, double tol=1.e-8)
{
if(u < (-tol) || v < (-tol) || w < (-tol) || u > ((1. + tol) - v - w))
return false;
return true;
}
- virtual void jac(int order, MVertex *verts[], double u, double v, double w,
- double j[3][3]);
- virtual void jac(int order, std::vector<MVertex*>& verts, double u, double v, double w,
- double j[3][3])
- {
- jac(order, &verts[0], u, v, w, j);
- }
- virtual void jac(double u, double v, double w, double j[3][3])
- {
- jac(1, 0, u, v, w, j);
- }
- virtual void pnt(int order, MVertex *verts[], double u, double v, double w,
- SPoint3 &p);
- virtual void pnt(int order, std::vector<MVertex*> &verts, double u, double v, double w,
- SPoint3 &p)
- {
- pnt(order, &verts[0], u, v, w, p);
- }
- virtual void pnt(double u, double v, double w, SPoint3 &p)
- {
- pnt(1, 0, u, v, w, p);
- }
+
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
virtual SPoint3 circumcenter();
private:
@@ -1346,8 +1256,8 @@ class MTetrahedron10 : public MTetrahedron {
for(int i = 0; i < 6; i++) _vs[i]->setPolynomialOrder(2);
}
~MTetrahedron10(){}
- virtual int getPolynomialOrder(){ return 2; }
- virtual int getNumVertices(){ return 10; }
+ virtual int getPolynomialOrder() const { return 2; }
+ virtual int getNumVertices() const { return 10; }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
virtual MVertex *getVertexUNV(int num)
{
@@ -1364,7 +1274,7 @@ class MTetrahedron10 : public MTetrahedron {
static const int map[10] = {0, 2, 1, 3, 6, 5, 4, 7, 8, 9};
return getVertex(map[num]);
}
- virtual int getNumEdgeVertices(){ return 6; }
+ virtual int getNumEdgeVertices() const { return 6; }
virtual int getNumEdgesRep(){ return 12; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -1411,11 +1321,11 @@ class MTetrahedron10 : public MTetrahedron {
v[4] = _vs[f[num][1]];
v[5] = _vs[f[num][2]];
}
- virtual int getTypeForMSH(){ return MSH_TET_10; }
- virtual int getTypeForUNV(){ return 118; } // solid parabolic tetrahedron
- virtual int getTypeForVTK(){ return 24; }
- virtual const char *getStringForPOS(){ return "SS2"; }
- virtual const char *getStringForBDF(){ return "CTETRA"; }
+ virtual int getTypeForMSH() const { return MSH_TET_10; }
+ virtual int getTypeForUNV() const { return 118; } // solid parabolic tetrahedron
+ virtual int getTypeForVTK() const { return 24; }
+ virtual const char *getStringForPOS() const { return "SS2"; }
+ virtual const char *getStringForBDF() const { return "CTETRA"; }
virtual void revert()
{
MVertex *tmp;
@@ -1423,19 +1333,6 @@ 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 getShapeFunction(int num, double u, double v, double w, double &s);
- virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]);
-
- virtual void jac(double u,double v,double w, double jac[3][3])
- {
- MTetrahedron::jac(2,_vs,u,v,w,jac);
- }
- virtual void pnt(double u,double v,double w, SPoint3& p)
- {
- MTetrahedron::pnt(2, _vs, u, v, w, p);
- }
};
/*
@@ -1531,31 +1428,32 @@ class MTetrahedronN : public MTetrahedron {
for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
}
~MTetrahedronN(){}
- virtual int getPolynomialOrder(){ return _order; }
- virtual int getNumVertices(){ return 4 + _vs.size(); }
+ virtual int getPolynomialOrder() const { return _order; }
+ virtual int getNumVertices() const { return 4 + _vs.size(); }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
- virtual int getNumEdgeVertices() const { return _order - 1; }
+ virtual int getNumEdgeVertices() const { return 6*(_order - 1); }
virtual int getNumFaceVertices() const
{
- return ((_order - 1) * (_order - 2)) / 2;
+ return 4 * ((_order - 1) * (_order - 2)) / 2;
}
virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
{
- v.resize(2 + getNumEdgeVertices());
+ v.resize(_order + 1);
MTetrahedron::_getEdgeVertices(num, v);
int j = 2;
- const int ie = (num + 1) * getNumEdgeVertices();
- for(int i = num * getNumEdgeVertices(); i != ie; ++i) v[j++] = _vs[i];
+ const int ie = (num + 1) * (_order -1);
+ for(int i = num * (_order -1); i != ie; ++i) v[j++] = _vs[i];
}
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
{
- v.resize(3 + getNumFaceVertices());
+ v.resize(3 + 3 * (_order - 1) + (_order-1) * (_order - 2) /2);
MTetrahedron::_getFaceVertices(num, v);
int j = 3;
- const int ie = (num+1)*getNumFaceVertices();
- for(int i = num*getNumFaceVertices(); i != ie; ++i) v[j++] = _vs[i];
+ int nbV = (_order - 1) * (_order - 2) / 2;
+ const int ie = (num+1)*nbV;
+ for(int i = num*nbV; i != ie; ++i) v[j++] = _vs[i];
}
- virtual int getNumVolumeVertices()
+ virtual int getNumVolumeVertices() const
{
switch(getTypeForMSH()){
case MSH_TET_35 : return 1;
@@ -1563,8 +1461,7 @@ class MTetrahedronN : public MTetrahedron {
default : return 0;
}
}
- virtual int getNumEdgeVertices(){ return _order - 1; }
- virtual int getTypeForMSH()
+ virtual int getTypeForMSH() const
{
// (p+1)*(p+2)*(p+3)/6
if(_order == 3 && _vs.size() + 4 == 20) return MSH_TET_20;
@@ -1608,20 +1505,10 @@ class MTetrahedronN : public MTetrahedron {
}
}
}
- virtual void getShapeFunction(int num, double u, double v, double w, double &s);
- virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]);
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual int getNumEdgesRep();
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual int getNumFacesRep();
- virtual void jac(double u, double v, double w, double j[3][3])
- {
- MTetrahedron::jac(_order, _vs , u, v, w, j);
- }
- virtual void pnt(double u, double v, double w, SPoint3 &p)
- {
- MTetrahedron::pnt(_order, _vs, u, v, w, p);
- }
};
/*
@@ -1671,7 +1558,7 @@ class MHexahedron : public MElement {
}
~MHexahedron(){}
virtual int getDim(){ return 3; }
- virtual int getNumVertices(){ return 8; }
+ virtual int getNumVertices() const { return 8; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual MVertex *getVertexMED(int num)
{
@@ -1722,11 +1609,11 @@ class MHexahedron : public MElement {
v.resize(4);
_getFaceVertices(num, v);
}
- virtual int getTypeForMSH(){ return MSH_HEX_8; }
- virtual int getTypeForUNV(){ return 115; } // solid linear brick
- virtual int getTypeForVTK(){ return 12; }
- virtual const char *getStringForPOS(){ return "SH"; }
- virtual const char *getStringForBDF(){ return "CHEXA"; }
+ virtual int getTypeForMSH() const { return MSH_HEX_8; }
+ virtual int getTypeForUNV() const { return 115; } // solid linear brick
+ virtual int getTypeForVTK() const { return 12; }
+ virtual const char *getStringForPOS() const { return "SH"; }
+ virtual const char *getStringForBDF() const { return "CHEXA"; }
virtual void revert()
{
MVertex *tmp;
@@ -1858,8 +1745,8 @@ class MHexahedron20 : public MHexahedron {
for(int i = 0; i < 12; i++) _vs[i]->setPolynomialOrder(2);
}
~MHexahedron20(){}
- virtual int getPolynomialOrder(){ return 2; }
- virtual int getNumVertices(){ return 20; }
+ virtual int getPolynomialOrder() const { return 2; }
+ virtual int getNumVertices() const { return 20; }
virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
virtual MVertex *getVertexUNV(int num)
{
@@ -1879,7 +1766,7 @@ class MHexahedron20 : public MHexahedron {
8, 17, 19, 18, 16, 10, 15, 14, 12};
return getVertex(map[num]);
}
- virtual int getNumEdgeVertices(){ return 12; }
+ virtual int getNumEdgeVertices() const { return 12; }
virtual int getNumEdgesRep(){ return 24; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -1937,10 +1824,10 @@ class MHexahedron20 : public MHexahedron {
v[6] = _vs[f[num][2]];
v[7] = _vs[f[num][3]];
}
- virtual int getTypeForMSH(){ return MSH_HEX_20; }
- virtual int getTypeForUNV(){ return 116; } // solid parabolic brick
- virtual int getTypeForVTK(){ return 25; }
- virtual const char *getStringForBDF(){ return "CHEXA"; }
+ virtual int getTypeForMSH() const { return MSH_HEX_20; }
+ virtual int getTypeForUNV() const { return 116; } // solid parabolic brick
+ virtual int getTypeForVTK() const { return 25; }
+ virtual const char *getStringForBDF() const { return "CHEXA"; }
virtual void revert()
{
MVertex *tmp;
@@ -1997,12 +1884,12 @@ class MHexahedron27 : public MHexahedron {
for(int i = 0; i < 19; i++) _vs[i]->setPolynomialOrder(2);
}
~MHexahedron27(){}
- virtual int getPolynomialOrder(){ return 2; }
- virtual int getNumVertices(){ return 27; }
+ virtual int getPolynomialOrder() const { return 2; }
+ virtual int getNumVertices() const { return 27; }
virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
- virtual int getNumEdgeVertices(){ return 12; }
- virtual int getNumFaceVertices(){ return 6; }
- virtual int getNumVolumeVertices(){ return 1; }
+ virtual int getNumEdgeVertices() const { return 12; }
+ virtual int getNumFaceVertices() const { return 6; }
+ virtual int getNumVolumeVertices() const { return 1; }
virtual int getNumEdgesRep(){ return 24; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -2067,8 +1954,8 @@ class MHexahedron27 : public MHexahedron {
v[7] = _vs[f[num][3]];
v[8] = _vs[12+num];
}
- virtual int getTypeForMSH(){ return MSH_HEX_27; }
- virtual const char *getStringForPOS(){ return "SH2"; }
+ virtual int getTypeForMSH() const { return MSH_HEX_27; }
+ virtual const char *getStringForPOS() const { return "SH2"; }
virtual void revert()
{
MVertex *tmp;
@@ -2141,7 +2028,7 @@ class MPrism : public MElement {
}
~MPrism(){}
virtual int getDim(){ return 3; }
- virtual int getNumVertices(){ return 6; }
+ virtual int getNumVertices() const { return 6; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual MVertex *getVertexMED(int num)
{
@@ -2196,11 +2083,11 @@ class MPrism : public MElement {
v.resize((num < 2) ? 3 : 4);
_getFaceVertices(num, v);
}
- virtual int getTypeForMSH(){ return MSH_PRI_6; }
- virtual int getTypeForUNV(){ return 112; } // solid linear wedge
- virtual int getTypeForVTK(){ return 13; }
- virtual const char *getStringForPOS(){ return "SI"; }
- virtual const char *getStringForBDF(){ return "CPENTA"; }
+ virtual int getTypeForMSH() const { return MSH_PRI_6; }
+ virtual int getTypeForUNV() const { return 112; } // solid linear wedge
+ virtual int getTypeForVTK() const { return 13; }
+ virtual const char *getStringForPOS() const { return "SI"; }
+ virtual const char *getStringForBDF() const { return "CPENTA"; }
virtual void revert()
{
MVertex *tmp;
@@ -2321,8 +2208,8 @@ class MPrism15 : public MPrism {
for(int i = 0; i < 9; i++) _vs[i]->setPolynomialOrder(2);
}
~MPrism15(){}
- virtual int getPolynomialOrder(){ return 2; }
- virtual int getNumVertices(){ return 15; }
+ virtual int getPolynomialOrder() const { return 2; }
+ virtual int getNumVertices() const { return 15; }
virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num - 6]; }
virtual MVertex *getVertexUNV(int num)
{
@@ -2339,7 +2226,7 @@ class MPrism15 : public MPrism {
static const int map[15] = {0, 2, 1, 3, 5, 4, 7, 9, 6, 13, 14, 12, 8, 11, 10};
return getVertex(map[num]);
}
- virtual int getNumEdgeVertices(){ return 9; }
+ virtual int getNumEdgeVertices() const { return 9; }
virtual int getNumEdgesRep(){ return 18; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -2393,9 +2280,9 @@ class MPrism15 : public MPrism {
v[i+2] = _vs[f[num][2]];
if (num >= 2) v[7] = _vs[f[num][3]];
}
- virtual int getTypeForMSH(){ return MSH_PRI_15; }
- virtual int getTypeForUNV(){ return 113; } // solid parabolic wedge
- virtual const char *getStringForBDF(){ return "CPENTA"; }
+ virtual int getTypeForMSH() const { return MSH_PRI_15; }
+ virtual int getTypeForUNV() const { return 113; } // solid parabolic wedge
+ virtual const char *getStringForBDF() const { return "CPENTA"; }
virtual void revert()
{
MVertex *tmp;
@@ -2449,11 +2336,11 @@ class MPrism18 : public MPrism {
for(int i = 0; i < 12; i++) _vs[i]->setPolynomialOrder(2);
}
~MPrism18(){}
- virtual int getPolynomialOrder(){ return 2; }
- virtual int getNumVertices(){ return 18; }
+ virtual int getPolynomialOrder() const { return 2; }
+ virtual int getNumVertices() const { return 18; }
virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num - 6]; }
- virtual int getNumEdgeVertices(){ return 9; }
- virtual int getNumFaceVertices(){ return 3; }
+ virtual int getNumEdgeVertices() const { return 9; }
+ virtual int getNumFaceVertices() const { return 3; }
virtual int getNumEdgesRep(){ return 18; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -2513,8 +2400,8 @@ class MPrism18 : public MPrism {
v[8] = _vs[7+num];
}
}
- virtual int getTypeForMSH(){ return MSH_PRI_18; }
- virtual const char *getStringForPOS(){ return "SI2"; }
+ virtual int getTypeForMSH() const { return MSH_PRI_18; }
+ virtual const char *getStringForPOS() const { return "SI2"; }
virtual void revert()
{
MVertex *tmp;
@@ -2586,7 +2473,7 @@ class MPyramid : public MElement {
}
~MPyramid(){}
virtual int getDim(){ return 3; }
- virtual int getNumVertices(){ return 5; }
+ virtual int getNumVertices() const { return 5; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual MVertex *getVertexMED(int num)
{
@@ -2638,10 +2525,10 @@ class MPyramid : public MElement {
v.resize((num < 4) ? 3 : 4);
_getFaceVertices(num, v);
}
- virtual int getTypeForMSH(){ return MSH_PYR_5; }
- virtual int getTypeForVTK(){ return 14; }
- virtual const char *getStringForPOS(){ return "SY"; }
- virtual const char *getStringForBDF(){ return "CPYRAM"; }
+ virtual int getTypeForMSH() const { return MSH_PYR_5; }
+ virtual int getTypeForVTK() const { return 14; }
+ virtual const char *getStringForPOS() const { return "SY"; }
+ virtual const char *getStringForBDF() const { return "CPYRAM"; }
virtual void revert()
{
MVertex *tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
@@ -2762,15 +2649,15 @@ class MPyramid13 : public MPyramid {
for(int i = 0; i < 8; i++) _vs[i]->setPolynomialOrder(2);
}
~MPyramid13(){}
- virtual int getPolynomialOrder(){ return 2; }
- virtual int getNumVertices(){ return 13; }
+ virtual int getPolynomialOrder() const { return 2; }
+ virtual int getNumVertices() const { return 13; }
virtual MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
virtual MVertex *getVertexMED(int num)
{
static const int map[13] = {0, 3, 2, 1, 4, 6, 10, 8, 5, 7, 12, 11, 9};
return getVertex(map[num]);
}
- virtual int getNumEdgeVertices(){ return 8; }
+ virtual int getNumEdgeVertices() const { return 8; }
virtual int getNumEdgesRep(){ return 16; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -2828,7 +2715,7 @@ class MPyramid13 : public MPyramid {
v[7] = _vs[0];
}
}
- virtual int getTypeForMSH(){ return MSH_PYR_13; }
+ virtual int getTypeForMSH() const { return MSH_PYR_13; }
virtual void revert()
{
MVertex *tmp;
@@ -2879,11 +2766,11 @@ class MPyramid14 : public MPyramid {
for(int i = 0; i < 9; i++) _vs[i]->setPolynomialOrder(2);
}
~MPyramid14(){}
- virtual int getPolynomialOrder(){ return 2; }
- virtual int getNumVertices(){ return 14; }
+ virtual int getPolynomialOrder() const { return 2; }
+ virtual int getNumVertices() const { return 14; }
virtual MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
- virtual int getNumEdgeVertices(){ return 8; }
- virtual int getNumFaceVertices(){ return 1; }
+ virtual int getNumEdgeVertices() const { return 8; }
+ virtual int getNumFaceVertices() const { return 1; }
virtual int getNumEdgesRep(){ return 16; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -2943,8 +2830,8 @@ class MPyramid14 : public MPyramid {
v[8] = _vs[8];
}
}
- virtual int getTypeForMSH(){ return MSH_PYR_14; }
- virtual const char *getStringForPOS(){ return "SY2"; }
+ virtual int getTypeForMSH() const { return MSH_PYR_14; }
+ virtual const char *getStringForPOS() const { return "SY2"; }
virtual void revert()
{
MVertex *tmp;
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index da2a62726ff8c19e27b3dbcfa843bc3b0a3a265b..54793e1b86f187c73626f36a6772a20573820b85 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -28,6 +28,294 @@ typedef std::map<std::pair<MVertex*,MVertex*>, std::vector<MVertex*> > edgeConta
// typedef std::map<std::vector<MVertex*>, std::vector<MVertex*>> faceContainer;
typedef std::map<MFace, std::vector<MVertex*>,Less_Face> faceContainer;
+bool mappingIsInvertible (MTetrahedron *e)
+{
+ if (e->getPolynomialOrder() == 1)return 1.0;
+
+ double mat[3][3];
+ e->getPrimaryJacobian(0.,0.,0., mat);
+ double det0 = det3x3(mat);
+
+ IntPt *pts;
+ int npts;
+
+ e->getIntegrationPoints(e->getPolynomialOrder(),&npts, &pts);
+
+ for (int i=0;i<npts;i++){
+
+ const double u = pts[i].pt[0];
+ const double v = pts[i].pt[1];
+ const double w = pts[i].pt[2];
+
+ e->getJacobian(u, v, w, mat);
+
+ double detN = det3x3(mat);
+
+ if (det0 * detN <= 0.) return false;
+ }
+
+ const Double_Matrix& points = e->getFunctionSpace()->points;
+
+ for (int i=0;i<e->getNumPrimaryVertices();i++) {
+
+ const double u = points(i,0);
+ const double v = points(i,1);
+ const double w = points(i,2);
+
+ e->getJacobian(u,v,w,mat);
+
+ double detN = det3x3(mat);
+
+ if (det0 * detN <= 0.) return false;
+ }
+
+ return true;
+}
+
+bool deformElement(MElement *ele,GEntity* ge,const std::set<MVertex*>& blocked) {
+
+ int nbNodes = ele->getNumVertices();
+
+ double _E = 1.;
+ double _nu = 0.1;
+
+ double FACT = _E / (1. + _nu) / (1. - 2. * _nu);
+ double C11 = FACT * (1. - _nu);
+ double C12 = FACT * _nu;
+ double C44 = (C11 - C12) / 2;
+
+ Double_Matrix stiffness(3*nbNodes,3*nbNodes);
+ Double_Vector rhs(3*nbNodes);
+
+ int npts;
+ IntPt *pts;
+ ele->getIntegrationPoints(2*ele->getPolynomialOrder(),&npts, &pts);
+ double d0;
+
+ const gmshFunctionSpace* fs = ele->getFunctionSpace();
+
+ const Double_Matrix& points = fs->points;
+
+
+ double jac[3][3];
+ ele->getPrimaryJacobian(0,0,0,jac);
+ double invjac[3][3];
+ inv3x3(jac,invjac);
+
+ double gsf[256][3];
+
+
+ // cheaper version : quadrature only used for parametric stiffness
+ // paramGrads should be stored in functionspace.
+
+// for (int j=0;j<nbNodes;j++) {
+// for (int k=0;k<nbNodes;k++) {
+
+// double xx = 0;
+// double xy = 0;
+// double xz = 0;
+// double yx = 0;
+// double yy = 0;
+// double yz = 0;
+// double zx = 0;
+// double zy = 0;
+// double zz = 0;
+
+
+// for (int m=0;m<3;m++) {
+// for (int n=0;n<3;n++) {
+
+// double dd = paramGrads[m][n](j,k);
+
+// xx += dd * invJac[0][m] * invJac[0][n];
+// xy += dd * invJac[0][m] * invJac[1][n];
+// xz += dd * invJac[0][m] * invJac[2][n];
+// yx += dd * invJac[1][m] * invJac[0][n];
+// yy += dd * invJac[1][m] * invJac[1][n];
+// yz += dd * invJac[1][m] * invJac[2][n];
+// zx += dd * invJac[2][m] * invJac[0][n];
+// zy += dd * invJac[2][m] * invJac[1][n];
+// zz += dd * invJac[2][m] * invJac[2][n];
+// }
+// }
+
+// stiffness(3*j,3*k) += (C11 * xx + // dphidx . tauxx
+// C44 * yy + // dphidy . tauxy
+// C44 * zz); // dphidz . tauxz
+
+// stiffness(3*j,3*k+1) += (C12 * xy + // dphidx . tauxx
+// C44 * yx); // dphidy . tauxy
+
+// stiffness(3*j,3*k+2) += (C12 * xz + // dphidx . tauxx
+// C44 * zx); // dphidz . tauxz
+
+// stiffness(3*j+1,3*k) += (C44 * xy + // dphidx . tauxy
+// C12 * yx); // dphidy . tauyy
+
+// stiffness(3*j+1,3*k+1) += (C44 * xx + // dphidx . tauxy
+// C11 * yy + // dphidy . tauyy
+// C44 * zz ); // dphidz . tauzy
+
+// stiffness(3*j+1,3*k+2) += (C12 * yz + // dphidy . tauyy
+// C44 * zy); // dphidz . tauzy
+
+// stiffness(3*j+2,3*k) += (C44 * xz + // dphidx . tauxz
+// C12 * zx); // dphidz . tauzz
+
+// stiffness(3*j+2,3*k+1) += (C44 * yz + // dphidy . tauyz
+// C12 * zy); // dphidz . tauzz
+
+// stiffness(3*j+2,3*k+2) += (C44 * xx + // dphidx . tauxz
+// C44 * yy + // dphidy . tauyz
+// C11 * zz); // dphidz . tauzz
+// }
+// }
+
+
+ // more expensive version
+
+ double Grads[256][3];
+
+ stiffness.scale(0.);
+
+ for (int i=0;i<npts;i++){
+
+ const double u = pts[i].pt[0];
+ const double v = pts[i].pt[1];
+ const double w = pts[i].pt[2];
+
+ const double weight = pts[i].weight;
+
+ fs->df(u,v,w,gsf);
+
+ for (int j=0;j<nbNodes;j++) {
+ Grads[j][0] = invjac[0][0] * gsf[j][0] + invjac[0][1] * gsf[j][1] + invjac[0][2] * gsf[j][2];
+ Grads[j][1] = invjac[1][0] * gsf[j][0] + invjac[1][1] * gsf[j][1] + invjac[1][2] * gsf[j][2];
+ Grads[j][2] = invjac[2][0] * gsf[j][0] + invjac[2][1] * gsf[j][1] + invjac[2][2] * gsf[j][2];
+ }
+
+ for (int j=0;j<nbNodes;j++) {
+ for (int k=0;k<nbNodes;k++) {
+
+ // R_x = dphidx . tauxx + dphidy . tauxy + dphidz . tauxz
+ // tau_xx = C11 . e_xx + C12 . (e_yy + e_zz)
+ // tau_xy = C44 . e_xy
+ // e_xy = 0.5 * (dX/dx + dY/dy)
+
+ double xx = Grads[j][0] * Grads[k][0];
+ double xy = Grads[j][0] * Grads[k][1];
+ double xz = Grads[j][0] * Grads[k][2];
+
+ double yx = Grads[j][1] * Grads[k][0];
+ double yy = Grads[j][1] * Grads[k][1];
+ double yz = Grads[j][1] * Grads[k][2];
+
+ double zx = Grads[j][2] * Grads[k][0];
+ double zy = Grads[j][2] * Grads[k][1];
+ double zz = Grads[j][2] * Grads[k][2];
+
+
+ // Poisson
+
+ stiffness(3*j,3*k) += weight * (xx + yy + zz);
+ stiffness(3*j+1,3*k+1) += weight * (xx + yy + zz);
+ stiffness(3*j+2,3*k+2) += weight * (xx + yy + zz);
+
+
+ // Elasticity
+
+ // stiffness(3*j ,3*k ) += weight * (C11 * xx + C44 * yy + C44 * zz);
+// stiffness(3*j ,3*k+1) += weight * (C12 * xy + C44 * yx);
+// stiffness(3*j ,3*k+2) += weight * (C12 * xz + C44 * zx);
+
+// stiffness(3*j+1,3*k ) += weight * (C12 * yx + C44 * xy);
+// stiffness(3*j+1,3*k+1) += weight * (C44 * xx + C11 * yy + C44 * zz);
+// stiffness(3*j+1,3*k+2) += weight * (C12 * yz + C44 * zy);
+
+// stiffness(3*j+2,3*k ) += weight * (C12 * zx + C44 * xz);
+// stiffness(3*j+2,3*k+1) += weight * (C12 * zy + C44 * yz);
+// stiffness(3*j+2,3*k+2) += weight * (C44 * xx + C44 * yy + C11 * zz);
+
+ }
+ }
+ }
+
+ int nbDirichlet = 0;
+
+ Double_Vector original(nbNodes*3);
+
+ for (int i=0;i<nbNodes;i++) {
+
+ MVertex* v = ele->getVertex(i);
+
+ SPoint3 primary;
+ ele->primaryPnt(points(i,0),points(i,1),points(i,2),primary);
+ original(3*i ) = primary.x();
+ original(3*i+1) = primary.y();
+ original(3*i+2) = primary.z();
+
+ MFaceVertex* vf = dynamic_cast<MFaceVertex*> (v);
+ MFaceVertex* ve = dynamic_cast<MFaceVertex*> (v);
+
+ if (v->onWhat() != ge || blocked.find(v) != blocked.end()) {
+
+ nbDirichlet++;
+
+ double dx = v->x() - original(3*i );
+ double dy = v->y() - original(3*i+1);
+ double dz = v->z() - original(3*i+2);
+
+ for (int j=0;j<3*nbNodes;j++) {
+ stiffness(3*i ,j) = 0;
+ stiffness(3*i+1,j) = 0;
+ stiffness(3*i+2,j) = 0;
+ }
+
+ for (int k=0;k<3;k++) stiffness(i*3+k,i*3+k) = 1.;
+
+ rhs(i*3 ) = dx;
+ rhs(i*3+1) = dy;
+ rhs(i*3+2) = dz;
+ }
+ else {
+ rhs(i*3 ) = 0.;
+ rhs(i*3+1) = 0.;
+ rhs(i*3+2) = 0.;
+ }
+ }
+
+ if (nbDirichlet < 3) {
+ Msg::Warning("Could not deform element due to lack of constraints\n");
+ return false;
+ }
+
+ if (nbDirichlet == nbNodes) {
+ Msg::Warning("Could not deform element because fully constrained\n");
+ return false;
+ }
+
+
+ Msg::Warning("Deforming element - have %d constrained points of which %d from previous positioning",nbDirichlet,(int) blocked.size());
+
+
+ Double_Vector displacement(nbNodes*3);
+
+
+ stiffness.lu_solve(rhs,displacement);
+
+ for (int i=0;i<nbNodes;i++) {
+ MVertex* v = ele->getVertex(i);
+
+ v->x() = original(3*i) + displacement(3*i );
+ v->y() = original(3*i+1) + displacement(3*i+1);
+ v->z() = original(3*i+2) + displacement(3*i+2);
+
+ }
+
+ return true;
+}
+
+
bool reparamOnFace(MVertex *v, GFace *gf, SPoint2 ¶m)
{
if(v->onWhat()->dim() == 0){
@@ -364,8 +652,11 @@ void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
}
}
-void getEdgeVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &ve,
- edgeContainer &edgeVertices, bool linear, int nPts = 1)
+void getEdgeVertices(GRegion *gr, MElement *ele,
+ std::vector<MVertex*> &ve,
+ std::set<MVertex*>& blocked,
+ edgeContainer &edgeVertices, bool linear,
+ int nPts = 1)
{
for(int i = 0; i < ele->getNumEdges(); i++){
MEdge edge = ele->getEdge(i);
@@ -375,6 +666,9 @@ void getEdgeVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &ve,
ve.insert(ve.end(), edgeVertices[p].begin(), edgeVertices[p].end());
else
ve.insert(ve.end(), edgeVertices[p].rbegin(), edgeVertices[p].rend());
+ blocked.insert(edgeVertices[p].begin(),edgeVertices[p].end());
+ blocked.insert(edge.getMinVertex());
+ blocked.insert(edge.getMaxVertex());
}
else{
std::vector<MVertex*> temp;
@@ -689,7 +983,9 @@ void reorientTrianglePoints(std::vector<MVertex*>& vtcs,int orientation,bool swa
// KH: check face orientation wrt element ...
-void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
+void getFaceVertices(GRegion *gr, MElement *ele,
+ std::vector<MVertex*> &vf,
+ std::set<MVertex*>& blocked,
faceContainer &faceVertices,
edgeContainer& edgeVertices,
bool linear, int nPts = 1)
@@ -734,6 +1030,11 @@ void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
}
else Msg::Error("Error in face lookup for recuperation of high order face nodes");
vf.insert(vf.end(), vtcs.begin(), vtcs.end());
+
+ blocked.insert(vtcs.begin(),vtcs.end());
+ blocked.insert(face.getVertex(0));
+ blocked.insert(face.getVertex(1));
+ blocked.insert(face.getVertex(2));
}
else{
std::vector<MVertex*>& vtcs = faceVertices[face];
@@ -981,36 +1282,68 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
faceContainer &faceVertices, bool linear, bool incomplete,
int nPts = 1)
{
+ int nbCorr = 0;
+
std::vector<MTetrahedron*> tetrahedra2;
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
MTetrahedron *t = gr->tetrahedra[i];
+
+ std::set<MVertex*> blocked;
+
std::vector<MVertex*> ve,vf,vr;
- getEdgeVertices(gr, t, ve, edgeVertices, linear, nPts);
+ getEdgeVertices(gr, t, ve, blocked, edgeVertices, linear, nPts);
if (nPts == 1){
tetrahedra2.push_back
(new MTetrahedron10(t->getVertex(0), t->getVertex(1), t->getVertex(2),
t->getVertex(3), ve[0], ve[1], ve[2], ve[3], ve[4], ve[5]));
}
else{
- getFaceVertices(gr, t, vf, faceVertices, edgeVertices, linear, nPts);
+ getFaceVertices(gr, t, vf, blocked, faceVertices, edgeVertices, linear, nPts);
ve.insert(ve.end(), vf.begin(), vf.end());
MTetrahedronN incpl (t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
ve, nPts + 1);
getRegionVertices(gr, &incpl, t, vr, linear, nPts);
- ve.insert(ve.end(), vr.begin(), vr.end());
- tetrahedra2.push_back
- (new MTetrahedronN(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
- ve, nPts + 1));
+ ve.insert(ve.end(), vr.begin(), vr.end());
+
+ MTetrahedron* n = new MTetrahedronN(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),ve, nPts + 1);
+
+ if (!mappingIsInvertible(n)) {
+ Msg::Warning("Found invalid curved volume element (# %d in list) ",i);
+
+ // if (deformElement(n,gr,blocked)) {
+// if (mappingIsInvertible(n)) printf(" - corrected using Poisson smoothing\n");
+// else printf(" - could not correct the mapping\n");
+// nbCorr++;
+// }
+ }
+ tetrahedra2.push_back (n);
}
delete t;
}
gr->tetrahedra = tetrahedra2;
+ std::vector<int> invalid;
+
+ if (nbCorr != 0) {
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
+ if (!mappingIsInvertible(gr->tetrahedra[i])) invalid.push_back(i);
+ }
+
+ if (invalid.size()!=0) {
+ Msg::Warning("We have %d invalid elements remaining\n",(int) invalid.size());
+ std::vector<int>::iterator iIter = invalid.begin();
+ for (;iIter!=invalid.end();++iIter) printf("%d ",*iIter);
+ printf("\n");
+ }
+ }
+
+
std::vector<MHexahedron*> hexahedra2;
for(unsigned int i = 0; i < gr->hexahedra.size(); i++){
MHexahedron *h = gr->hexahedra[i];
std::vector<MVertex*> ve, vf;
- getEdgeVertices(gr, h, ve, edgeVertices, linear, nPts);
+ std::set<MVertex*> blocked;
+ getEdgeVertices(gr, h, ve, blocked, edgeVertices, linear, nPts);
if(incomplete){
hexahedra2.push_back
(new MHexahedron20(h->getVertex(0), h->getVertex(1), h->getVertex(2),
@@ -1020,7 +1353,7 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
ve[11]));
}
else{
- getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts);
+ getFaceVertices(gr, h, vf, blocked, faceVertices, edgeVertices, linear, nPts);
SPoint3 pc = h->barycenter();
MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
gr->mesh_vertices.push_back(v);
@@ -1039,7 +1372,8 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
for(unsigned int i = 0; i < gr->prisms.size(); i++){
MPrism *p = gr->prisms[i];
std::vector<MVertex*> ve, vf;
- getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
+ std::set<MVertex*> blocked;
+ getEdgeVertices(gr, p, ve, blocked,edgeVertices, linear, nPts);
if(incomplete){
prisms2.push_back
(new MPrism15(p->getVertex(0), p->getVertex(1), p->getVertex(2),
@@ -1047,7 +1381,7 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8]));
}
else{
- getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
+ getFaceVertices(gr, p, vf, blocked, faceVertices, edgeVertices, linear, nPts);
prisms2.push_back
(new MPrism18(p->getVertex(0), p->getVertex(1), p->getVertex(2),
p->getVertex(3), p->getVertex(4), p->getVertex(5),
@@ -1062,7 +1396,8 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
for(unsigned int i = 0; i < gr->pyramids.size(); i++){
MPyramid *p = gr->pyramids[i];
std::vector<MVertex*> ve, vf;
- getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
+ std::set<MVertex*> blocked;
+ getEdgeVertices(gr, p, ve, blocked, edgeVertices, linear, nPts);
if(incomplete){
pyramids2.push_back
(new MPyramid13(p->getVertex(0), p->getVertex(1), p->getVertex(2),
@@ -1070,7 +1405,7 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
ve[3], ve[4], ve[5], ve[6], ve[7]));
}
else{
- getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
+ getFaceVertices(gr, p, vf, blocked, faceVertices, edgeVertices, linear, nPts);
pyramids2.push_back
(new MPyramid14(p->getVertex(0), p->getVertex(1), p->getVertex(2),
p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2],
@@ -1164,7 +1499,9 @@ static double mesh_functional_distorsion(MTriangle *t, double u, double v)
{
// compute uncurved element jacobian d_u x and d_v x
double mat[2][3];
- t->jac(1, 0, 0, 0, 0, mat);
+ // t->getPrimaryJacobian(0, 0, 0, mat);
+ t->getJacobian(0, 0, 0, mat);
+
double v1[3] = {mat[0][0], mat[0][1], mat[0][2]};
double v2[3] = {mat[1][0], mat[1][1], mat[1][2]};
double normal1[3];
@@ -1172,7 +1509,7 @@ static double mesh_functional_distorsion(MTriangle *t, double u, double v)
double nn = sqrt(SQU(normal1[0]) + SQU(normal1[1]) + SQU(normal1[2]));
// compute uncurved element jacobian d_u x and d_v x
- t->jac(u, v, 0, mat);
+ t->getJacobian(u, v, 0, mat);
double v1b[3] = {mat[0][0], mat[0][1], mat[0][2]};
double v2b[3] = {mat[1][0], mat[1][1], mat[1][2]};
double normal[3];
@@ -1191,7 +1528,8 @@ void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
{
double mat[2][3];
int n = 3;
- t->jac(1, 0, 0, 0, 0, mat);
+ // t->getPrimaryJacobian(0, 0, 0, mat);
+ t->getJacobian(0, 0, 0, mat);
double v1[3] = {mat[0][0], mat[0][1], mat[0][2]};
double v2[3] = {mat[1][0], mat[1][1], mat[1][2]};
double normal1[3], normal[3];
@@ -1199,7 +1537,7 @@ void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
double nn = sqrt(SQU(normal1[0]) + SQU(normal1[1]) + SQU(normal1[2]));
for(int i = 0; i < n; i++){
for(int k = 0; k < n - i; k++){
- t->jac((double)i / (n - 1), (double)k / (n - 1), 0, mat);
+ t->getJacobian((double)i / (n - 1), (double)k / (n - 1), 0, mat);
double v1b[3] = {mat[0][0], mat[0][1], mat[0][2]};
double v2b[3] = {mat[1][0], mat[1][1], mat[1][2]};
prodve(v1b, v2b, normal);
diff --git a/Numeric/FunctionSpace.cpp b/Numeric/FunctionSpace.cpp
index 85647e93e8dd1d59ffddfed14385d9acbedb477f..419692981360053bc57d535e0aa4ecd07cc8753d 100644
--- a/Numeric/FunctionSpace.cpp
+++ b/Numeric/FunctionSpace.cpp
@@ -7,6 +7,7 @@
#include "GmshDefines.h"
#include "Message.h"
+
Double_Matrix generate1DMonomials(int order)
{
Double_Matrix monomials(order + 1, 1);
@@ -20,7 +21,7 @@ Double_Matrix generate1DPoints(int order)
line(0, 0) = -1.;
line(1, 0) = 1.;
double dd = 2. / order;
- for (int i = 2; i < order; i++) line(i, 0) = -1. + dd * (i - 1);
+ for (int i = 2; i < order+1; i++) line(i, 0) = -1. + dd * (i - 1);
return line;
}
@@ -585,6 +586,10 @@ const gmshFunctionSpace &gmshFunctionSpaces::find(int tag)
gmshFunctionSpace F;
switch (tag){
+ case MSH_PNT:
+ F.monomials = generate1DMonomials(0);
+ F.points = generate1DPoints(0);
+ break;
case MSH_LIN_2 :
F.monomials = generate1DMonomials(1);
F.points = generate1DPoints(1);
@@ -675,3 +680,31 @@ const gmshFunctionSpace &gmshFunctionSpaces::find(int tag)
fs.insert(std::make_pair(tag, F));
return fs[tag];
}
+
+std::map<std::pair<int,int>, Double_Matrix> gmshFunctionSpaces::injector;
+
+
+const Double_Matrix &gmshFunctionSpaces::findInjector(int tag1,int tag2) {
+
+ std::pair<int,int> key(tag1,tag2);
+ std::map<std::pair<int,int>,Double_Matrix>::const_iterator it = injector.find(key);
+ if (it != injector.end()) return it->second;
+
+
+ const gmshFunctionSpace& fs1 = find(tag1);
+ const gmshFunctionSpace& fs2 = find(tag2);
+
+ Double_Matrix inj(fs1.points.size1(),fs2.points.size1());
+
+ double sf[256];
+
+ for (int i=0;i<fs1.points.size1();i++) {
+ fs2.f(fs1.points(i,0),fs1.points(i,1),fs1.points(i,2),sf);
+ for (int j=0;j<fs2.points.size1();j++) inj(i,j) = sf[j];
+ }
+
+ injector.insert(std::make_pair(key,inj));
+ return injector[key];
+}
+
+
diff --git a/Numeric/FunctionSpace.h b/Numeric/FunctionSpace.h
index 9976604c06971346dcb800804e2137d3c0fbb428..3495e7a1ba573074d020ca5ebca858bab55421ba 100644
--- a/Numeric/FunctionSpace.h
+++ b/Numeric/FunctionSpace.h
@@ -10,85 +10,173 @@
#include <map>
#include "GmshMatrix.h"
+#define MAX_FUNCTIONSPACESIZE 56
+
struct gmshFunctionSpace
{
Double_Matrix points;
Double_Matrix monomials;
Double_Matrix coefficients;
- inline void computePows(double uu, double vv, double p[][2]) const
- {
- for(int j = 0; j < coefficients.size2(); j++){
- p[j][0] = pow(uu,monomials(j, 0));
- p[j][1] = pow(vv,monomials(j, 1));
+
+ inline void evaluateMonomials(double u,double v,double w,double p[]) const {
+ for (int j=0;j<monomials.size1();j++) {
+ p[j] = pow(u,(int) monomials(j,0));
+ if (monomials.size2() > 1) p[j] *= pow(v,(int) monomials(j,1));
+ if (monomials.size2() > 2) p[j] *= pow(w,(int) monomials(j,2));
}
}
+
+/* inline void evaluateMonomialDerivatives(double u,double v,double w,double p[][3]) const */
+/* { */
+
+
+
+/* for (int j=0;j<monomials.size1();j++) { */
+
+/* int expU = (int) monomials(j,0); */
+/* double uu = pow(u,expU); */
+
+/* p[j][0] = pow(u,expU - 1) * expU; */
+/* p[j][1] = 0.; */
+/* p[j][2] = 0.; */
+
+
+/* if (monomials.size2() > 1) { */
+
+/* int expV = (int) monomials(j,1); */
+/* double vv = pow(v,expV); */
+
+/* p[j][0] *= vv; */
+/* p[j][1] = uu * pow(v,expV-1) * expV; */
+
+/* if (monomials.size2() > 2) { */
+
+/* int expW = (int) monomials(j,2); */
+/* double ww = pow(w,expW); */
+
+/* p[j][0] *= ww; */
+/* p[j][1] *= ww; */
+/* p[j][2] = uu * vv * pow(w,expW-1) * expW; */
+/* } */
+/* } */
+/* } */
+/* } */
+
inline void f(double u, double v, double w, double *sf) const
{
- for(int i = 0; i < coefficients.size1(); i++){
- sf[i] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- sf[i] += coefficients(i, j) * pow(u, monomials(j, 0)) *
- pow(v, monomials(j, 1)) * pow(w, monomials(j, 2));
- }
- }
- }
- inline void f(double u, double v, double *sf) const
- {
- double p[256][2];
- computePows(u, v, p);
- for (int i = 0; i < coefficients.size1(); i++){
+
+ double p[256];
+ evaluateMonomials(u,v,w,p);
+
+ for (int i=0;i<coefficients.size1();i++) {
sf[i] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- sf[i] += coefficients(i, j) * p[j][0] * p[j][1];
+ for (int j=0;j<coefficients.size2();j++) {
+ sf[i] += coefficients(i,j) * p[j];
}
}
}
+
inline void df(double u, double v, double w, double grads[][3]) const
{
- for (int i = 0; i < coefficients.size1(); i++){
- grads[i][0] = 0;
- grads[i][1] = 0;
- grads[i][2] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- if ((monomials)(j, 0) > 0)
- grads[i][0] += (coefficients)(i, j) * pow(u, (monomials)(j, 0) - 1) *
- (monomials)(j, 0) * pow(v, (monomials)(j, 1)) *
- pow(w, (monomials)(j, 2));
- if ((monomials)(j, 1) > 0)
- grads[i][1] += (coefficients)(i, j) * pow(u,(monomials)(j, 0)) *
- pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1) *
- pow(w, (monomials)(j, 2));
- if ((monomials)(j, 2) > 0)
- grads[i][2] += (coefficients)(i, j) * pow(u, (monomials)(j, 0)) *
- pow(v, (monomials)(j, 1)) * pow(w, (monomials)(j, 2) - 1) *
- (monomials)(j, 2);
+
+ switch (monomials.size2()) {
+
+ case 1:
+ {
+ for (int i = 0; i < coefficients.size1(); i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ grads[i][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if ((monomials)(j, 0) > 0)
+ grads[i][0] += (coefficients)(i, j) * pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0);
+ }
+ }
+ break;
}
- }
- }
- inline void df(double u, double v, double w,double grads[][2]) const
- {
- double p[256][2];
- computePows(u, v, p);
- for (int i = 0; i < coefficients.size1(); i++){
- grads[i][0] = 0;
- grads[i][1] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- if ((monomials)(j, 0) > 0)
- grads[i][0] += (coefficients)(i, j) *
- pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0) * p[j][1];
- if ((monomials)(j, 1) > 0)
- grads[i][1] += (coefficients)(i, j) * p[j][0] *
- pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1);
+ case 2:
+ {
+
+ for (int i = 0; i < coefficients.size1(); i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ grads[i][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if ((monomials)(j, 0) > 0)
+ grads[i][0] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0) *
+ pow(v, (monomials)(j, 1));
+ if ((monomials)(j, 1) > 0)
+ grads[i][1] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1);
+ }
+ }
+
+ break;
+ }
+
+ case 3:
+ {
+
+ for (int i = 0; i < coefficients.size1(); i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ grads[i][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if ((monomials)(j, 0) > 0)
+ grads[i][0] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0) *
+ pow(v, (monomials)(j, 1)) *
+ pow(w, (monomials)(j, 2));
+ if ((monomials)(j, 1) > 0)
+ grads[i][1] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1) *
+ pow(w, (monomials)(j, 2));
+ if ((monomials)(j, 2) > 0)
+ grads[i][2] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1)) *
+ pow(w, (monomials)(j, 2) - 1) * (monomials)(j, 2);
+ }
+ }
+ break;
}
}
}
+
+
+/* inline void df(double u, double v, double w, double grads[][3]) const */
+/* { */
+
+/* double p[256][3]; */
+
+/* evaluateMonomialDerivatives(u,v,w,p); */
+
+/* for (int i=0;i<coefficients.size1();i++) { */
+
+/* grads[i][0] = 0.; */
+/* grads[i][1] = 0.; */
+/* grads[i][2] = 0.; */
+
+/* for (int j=0;j<coefficients.size2();j++) { */
+/* grads[i][0] += coefficients(i,j) * p[j][0]; */
+/* grads[i][1] += coefficients(i,j) * p[j][1]; */
+/* // grads[i][2] += coefficients(i,j) * p[j][2]; */
+/* } */
+/* } */
+/* } */
};
class gmshFunctionSpaces
{
static std::map<int, gmshFunctionSpace> fs;
+ static std::map<std::pair<int,int>,Double_Matrix> injector;
+
public :
- static const gmshFunctionSpace &find(int);
+ static const gmshFunctionSpace &find(int);
+ static const Double_Matrix &findInjector(int,int);
};
#endif