diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index d58c857965fd0a9986be6bad026e7143738b4b11..eda15b5d7c5c3b34bcc4227bf2d84bf3010e5a56 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -337,36 +337,31 @@ std::string MElement::getInfoString()
return std::string(tmp);
}
-double MElement::getJacobian(double u, double v, double w, double jac[3][3])
+static double _computeDeterminantAndRegularize(MElement *ele, 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 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()) {
+ switch (ele->getDim()) {
case 1:
{
- dJ = sqrt(jac[0][0] * jac[0][0] + jac[1][1] * jac[0][0] + jac[2][2] * jac[2][2]);
+ dJ = sqrt(SQU(jac[0][0]) + SQU(jac[0][1]) + SQU(jac[0][2]));
+
+ // regularize matrix
+ double a[3], b[3], c[3];
+ a[0] = ele->getVertex(1)->x() - ele->getVertex(0)->x();
+ a[1] = ele->getVertex(1)->y() - ele->getVertex(0)->y();
+ a[2] = ele->getVertex(1)->z() - ele->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];
break;
}
case 2:
@@ -374,6 +369,18 @@ double MElement::getJacobian(double u, double v, double w, double jac[3][3])
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]));
+
+ // regularize matrix
+ 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];
break;
}
case 3:
@@ -387,112 +394,88 @@ double MElement::getJacobian(double u, double v, double w, double jac[3][3])
return dJ;
}
-
-
-double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+double MElement::getJacobian(double u, double v, double w, double jac[3][3])
{
-
- const gmshFunctionSpace* fs = getFunctionSpace(1);
+ 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 gsf[256][3];
- fs->df(u,v,w,gsf);
-
- for (int i=0;i<getNumPrimaryVertices();i++) {
-
+ 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++) {
+ 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;
+ return _computeDeterminantAndRegularize(this, jac);
+}
+
+double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+{
+ const gmshFunctionSpace* fs = getFunctionSpace(1);
- switch (fs->monomials.size2()) {
+ 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.;
- case 1:
- {
- dJ = sqrt(jac[0][0] * jac[0][0] + jac[0][0] * jac[0][0] + jac[0][0] * jac[0][0]);
- break;
- }
- 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]);
- break;
+ 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];
}
}
- return dJ;
-}
-
-void MElement::pnt(double uu, double vv, double ww, SPoint3 &p) {
+ return _computeDeterminantAndRegularize(this, jac);
+}
- double x=0.;
- double y=0.;
- double z=0.;
-
+void MElement::pnt(double uu, double vv, double ww, SPoint3 &p)
+{
+ double x = 0., y = 0., z = 0.;
const gmshFunctionSpace* fs = getFunctionSpace();
-
- if (fs) {
-
+ if(fs){
double sf[256];
- fs->f(uu,vv,ww,sf);
-
- for (int j=0;j<getNumVertices();j++) {
+ 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);
-
+ else
+ Msg::Error("Could not find function space");
+ 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.;
-
+void MElement::primaryPnt(double uu, double vv, double ww, SPoint3 &p)
+{
+ double x = 0., y = 0., z = 0.;
const gmshFunctionSpace* fs = getFunctionSpace(1);
-
- if (fs) {
-
+ 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++) {
+ fs->f(uu, vv, ww, sf);
+ if (getNumPrimaryVertices() != 4)
+ Msg::Error("Incorrect number of vertices %d", 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])
@@ -902,7 +885,10 @@ const gmshFunctionSpace* MLine::getFunctionSpace(int o) const {
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;
+ default:
+ Msg::Error("Order %d line point interpolation not implemented",
+ getPolynomialOrder());
+ break;
}
return NULL;
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 785e14ec06808a61c0f5e867d062df64021d54d1..3ca34e483c465ead312025c9edec9d7d48ebacd5 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -102,9 +102,7 @@ class MElement
// get the number of primary vertices (first-order element)
int getNumPrimaryVertices()
{
- return getNumVertices() -
- getNumEdgeVertices() -
- getNumFaceVertices() -
+ return getNumVertices() - getNumEdgeVertices() - getNumFaceVertices() -
getNumVolumeVertices();
}
@@ -141,11 +139,10 @@ class MElement
virtual double minEdge();
// get the quality measures
-
virtual double rhoShapeMeasure();
virtual double gammaShapeMeasure(){ return 0.; }
virtual double etaShapeMeasure(){ return 0.; }
- virtual double distoShapeMeasure() {return 1.0;}
+ virtual double distoShapeMeasure(){ return 1.0; }
// compute the barycenter
virtual SPoint3 barycenter();
@@ -162,28 +159,31 @@ class MElement
// 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;
+ // get the function space for the element
+ virtual const gmshFunctionSpace* getFunctionSpace(int ord = -1) const
+ {
+ Msg::Fatal("Function space not implemented for element %s", getStringForPOS());
+ return 0;
}
// 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);
+ 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]);
+ 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]);
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);
+ // get the point in cartesian coordinates corresponding to the point
+ // (u,v,w) in parametric coordinates
+ 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]);
@@ -351,17 +351,13 @@ class MLine : public MElement {
{
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
}
-
-
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;
};
@@ -405,14 +401,12 @@ class MLine3 : public MLine {
};
_getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n);
}
-
virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(3);
MLine::_getEdgeVertices(v);
v[2] = _vs[0];
}
-
virtual int getTypeForMSH() const { return MSH_LIN_3; }
virtual int getTypeForUNV() const { return 24; } // parabolic beam
virtual int getTypeForVTK() const { return 21; }
@@ -568,17 +562,13 @@ class MTriangle : public MElement {
{
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
}
-
-
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:
@@ -1192,7 +1182,6 @@ class MTetrahedron : public MElement {
return false;
return true;
}
-
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
virtual SPoint3 circumcenter();
private: