diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index a948b4b9b169f1069caa2e8c02b4acf38f9a99eb..cac071cbb16b5119d8924d684601f0019fe34c3e 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -317,17 +317,6 @@ double MQuadrangle::gammaShapeMeasure(){
}
-double MQuadrangle::distoShapeMeasure()
-{
-#if defined(HAVE_MESH)
- // return qmTriangleAngles(this);
- return qmDistorsionOfMapping(this);
-#else
- return 0.;
-#endif
-}
-
-
double MQuadrangle::angleShapeMeasure()
{
#if defined(HAVE_MESH)
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 5eb78ee5049b14cc7a8fade4d920bbcdddc37a98..1036738ad38fc01da64f9a33a058608d1dc3e6c9 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -148,7 +148,6 @@ class MQuadrangle : public MElement {
}
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
virtual double angleShapeMeasure();
- virtual double distoShapeMeasure();
// Computes the minimum inradius of the all the circles tangents to
// 3 of the 4 edges of the quad. If the 4 points of the quad are not
// planar, we compute the mean plane due to the least-square
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 23cf8e8e6c3cfcbe52995260c202bc2e609ff9cd..614fd85bf17b161b2c3a93f17654291d713f7310 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -44,82 +44,6 @@ double MTetrahedron::getCircumRadius()
#endif
}
-void MTetrahedron10::scaledJacRange(double &jmin, double &jmax){
-
- // the fastest available algo for computing scaled jacobians for
- // quadratic tetrahedron
- const int nbNodT = 10;
- const int nbBezT = 20;
- static int done = 0;
- static fullMatrix<double> gsf[nbBezT];
- const JacobianBasis *jac_basis = getJacobianFuncSpace();
- if (!done){
- double gs[nbNodT][3];
- for (int i=0;i<jac_basis->getNumJacNodes();i++){
- const double u = jac_basis->getPoints()(i,0);
- const double v = jac_basis->getPoints()(i,1);
- const double w = jac_basis->getPoints()(i,2);
- getGradShapeFunctions(u,v,w,gs);
- fullMatrix<double> a (nbNodT,3);
- for (int j=0;j < nbNodT;j++){
- a(j,0) = gs[j][0];
- a(j,1) = gs[j][1];
- a(j,2) = gs[j][2];
- }
- gsf[i]= a;
- }
- done=1;
- }
-
- double x[nbNodT];for (int i=0;i<nbNodT;i++)x[i] = getVertex(i)->x();
- double y[nbNodT];for (int i=0;i<nbNodT;i++)y[i] = getVertex(i)->y();
- double z[nbNodT];for (int i=0;i<nbNodT;i++)z[i] = getVertex(i)->z();
-
- SVector3 v1 (x[1] - x[0],y[1] - y[0],z[1] - z[0]);
- SVector3 v2 (x[2] - x[0],y[2] - y[0],z[2] - z[0]);
- SVector3 v3 (x[3] - x[0],y[3] - y[0],z[3] - z[0]);
-
- // take absolute value for quite complex reasons
- // straight sided element may be WRONG while curved
- // one is OK
- const double ss = fabs(1./dot(crossprod(v1,v2),v3));
-
- double jac[3][3];
- static fullVector<double>Ji(nbBezT);
- for (int i=0;i < nbBezT;i++){
- 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.;
- fullMatrix<double> &g = gsf[i];
- for (int j = 0; j < nbNodT; j++) {
- for (int k = 0; k < 3; k++) {
- const double gg = g(j, k);
- jac[k][0] += x[j] * gg;
- jac[k][1] += y[j] * gg;
- jac[k][2] += z[j] * gg;
- }
- }
- const double dJ = 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];
- Ji(i) = dJ * ss;
- }
- static fullVector<double> Bi( nbBezT );
- jac_basis->lag2Bez(Ji,Bi);
- // printf("%22.15E\n",*std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size()));
- jmin= *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
- jmax= *std::max_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
-}
-
-double MTetrahedron::distoShapeMeasure()
-{
-#if defined(HAVE_MESH)
- return qmDistorsionOfMapping(this);
-#else
- return 0.;
-#endif
-}
-
double MTetrahedron::getInnerRadius()
{
// radius of inscribed sphere = 3 * Volume / sum(Area_i)
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index d0e5e9cae82f934bf64b3baa9be8ccbfd10d7066..5eb4f6b89e796b4f4239c7f5900fda6f82f64629 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -123,7 +123,6 @@ class MTetrahedron : public MElement {
virtual double getVolume();
virtual int getVolumeSign(){ return (getVolume() >= 0) ? 1 : -1; }
virtual double gammaShapeMeasure();
- virtual double distoShapeMeasure();
virtual double getInnerRadius();
virtual double getCircumRadius();
virtual double etaShapeMeasure();
@@ -213,7 +212,6 @@ class MTetrahedron10 : public MTetrahedron {
for(int i = 0; i < 6; i++) _vs[i]->setPolynomialOrder(2);
}
~MTetrahedron10(){}
- virtual void scaledJacRange(double &jmin, double &jmax);
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]; }
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 518675aaf416896387c3ba13bebbb1807598c76e..6b7929109870528515c841b672b6757b0f72344e 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -36,15 +36,6 @@ double MTriangle::getVolume()
return norm(crossprod(v1, v2)) / 2.;
}
-double MTriangle::distoShapeMeasure()
-{
-#if defined(HAVE_MESH)
- return qmDistorsionOfMapping(this);
-#else
- return 0.;
-#endif
-}
-
double MTriangle::getInnerRadius()
{
// radius of inscribed circle = 2 * Area / sum(Line_i)
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 6a2168e4ef0b495f548ef0106e56f7f0fcd4a9ac..49cf7ecfeda33547d0fcc8a50742a98aa67ccf92 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -52,7 +52,6 @@ class MTriangle : public MElement {
virtual int getDim() const { return 2; }
virtual double etaShapeMeasure();
virtual double gammaShapeMeasure();
- virtual double distoShapeMeasure();
virtual double getInnerRadius();
virtual double getOuterRadius();
virtual double angleShapeMeasure();
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 46a8a908c283299c3f9e209a537a0e36da716138..3a6e4a5f9a77f64f779728bb8c1c6089954d7312 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -208,242 +208,6 @@ double conditionNumberAndDerivativeOfTet(const double &x1, const double &y1, con
}
*/
-double mesh_functional_distorsion_2D(MElement *t, double u, double v)
-{
- // compute uncurved element jacobian d_u x and d_v x
- double mat[3][3];
- t->getPrimaryJacobian(u, v, 0, mat);
- // t->getJacobian(u,v,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];
- prodve(v1, v2, normal1);
- double nn = norm3(normal1);
-
- // compute uncurved element jacobian d_u x and d_v x
-
- 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];
- prodve(v1b, v2b, normal);
-
- // printf("%g %g %g -- %g %g %g - %g\n",mat[0][0], mat[0][1], mat[0][2],mat[1][0], mat[1][1], mat[1][2],nn);
-
- double sign = 1.0;
- prosca(normal1, normal, &sign);
- // double det = (norm3(normal) / nn) * (sign > 0 ? 1. : -1.);
-
- // printf("%g %g : %g : %g n1 (%g,%g,%g)\n",u,v,sign,det, normal1[0], normal1[1], normal1[2]);
- // for (int i=0;i<t->getNumVertices();i++){
- // printf("COORD (%d) = %g %g %g\n",i,t->getVertex(i)->x(),t->getVertex(i)->y(),t->getVertex(i)->z());
- // }
- // printf("n (%g,%g,%g)\n", normal1[0], normal1[1], normal1[2]);
-
- return sign/(nn*nn);
-}
-
-static double MINQ (double a, double b, double c){
- if (a == 0) return std::min(a+b+c,c);
- double xmin = -b/(2*a);
- if (xmin < 0 || xmin > 1)return std::min(c,a+b+c);
- return std::min(c,std::min(a+b+c,a * xmin * xmin + b * xmin + c));
-}
-
-
-double mesh_functional_distorsion_p2_bezier_refined(MTriangle *t)
-{
- double J1 =mesh_functional_distorsion_2D(t,0.0,0.0);
- double J2 =mesh_functional_distorsion_2D(t,1.0,0.0);
- double J3 =mesh_functional_distorsion_2D(t,0.0,1.0);
- double J4 =mesh_functional_distorsion_2D(t,0.5,0.0);
- double J5 =mesh_functional_distorsion_2D(t,0.5,0.5);
- double J6 =mesh_functional_distorsion_2D(t,0.0,0.5);
-
- double J36 =mesh_functional_distorsion_2D(t,0.0,.75);
- double J35 =mesh_functional_distorsion_2D(t,0.25,.75);
- double J56 =mesh_functional_distorsion_2D(t,0.25,.5);
-
- double J16 =mesh_functional_distorsion_2D(t,0.0,.25);
- double J14 =mesh_functional_distorsion_2D(t,0.25,.0);
- double J46 =mesh_functional_distorsion_2D(t,0.25,.25);
-
-
- double J45 =mesh_functional_distorsion_2D(t,0.5,.25);
- double J52 =mesh_functional_distorsion_2D(t,0.75,.25);
- double J24 =mesh_functional_distorsion_2D(t,0.75,.0);
-
- double d[15] = {
- J1,J6,J4,2*J16-0.5*(J1+J6),2*J14-0.5*(J1+J4),2*J46-0.5*(J6+J4),
- J3,J5,2*J36-0.5*(J3+J6),2*J35-0.5*(J3+J5),2*J56-0.5*(J5+J6),
- J2,2*J45-0.5*(J4+J5),2*J52-0.5*(J5+J2),2*J24-0.5*(J2+J4)};
- return *std::min_element(d,d+15);
-}
-
-double mesh_functional_distorsion_p2_exact(MTriangle *t)
-{
- double J1 =mesh_functional_distorsion_2D(t,0.0,0.0);
- double J2 =mesh_functional_distorsion_2D(t,1.0,0.0);
- double J3 =mesh_functional_distorsion_2D(t,0.0,1.0);
- double J4 =mesh_functional_distorsion_2D(t,0.5,0.0);
- double J5 =mesh_functional_distorsion_2D(t,0.5,0.5);
- double J6 =mesh_functional_distorsion_2D(t,0.0,0.5);
-
- const double a = J1;
- const double b = -3*J1-J2+4*J4;
- const double c = -3*J1-J3+4*J6;
- const double d = 4*(J1-J4+J5-J6);
- const double e = 2*(J1+J2-2*J4);
- const double f = 2*(J1+J3-2*J6);
-
- double js[3] = {
- MINQ (2*(J1+J2-2*J4), -3*J1-J2+4*J4, J1),
- MINQ (2*(J1+J3-2*J6), -3*J1-J3+4*J6, J1),
- MINQ (2*(J3+J2-2*J5), -3*J2-J3+4*J5, J2)
- };
- double min_interm = *std::min_element(js,js+3);
-
- double mat[2][2] = {{2*e,d},{d,2*f}};
- double x[2], rhs[2] = {-b,-c};
-
- if (!sys2x2(mat,rhs,x))return min_interm;
-
- const double ximin = x[0];
- const double etamin = x[1];
-
- if (ximin> 0 && etamin > 0 && 1-ximin-etamin>0){
- const double m4 = a+b*ximin+c*etamin+d*ximin*etamin+
- e*ximin*ximin + f*etamin*etamin;
- /*
- if (m4 < min_interm && (m4 < .9 || m4 > 1.1)){
- printf("m4 = %g xi = %g eta = %g min_interm = %g min_edges = %g %g %g\n",m4,ximin,etamin,min_interm, MINQ (e,b,a), MINQ (f,c,a), MINQ (-d+e+f,b-c+d-2*f,a+c+f));
- FILE *f = fopen ("t.pos","w");
- fprintf(f,"ST2(%g,%g,0,%g,%g,0,%g,%g,0,%g,%g,0,%g,%g,0,%g,%g,0){%g,%g,%g,%g,%g,%g}\n",
- t->getVertex(0)->x(),t->getVertex(0)->y(),
- t->getVertex(1)->x(),t->getVertex(1)->y(),
- t->getVertex(2)->x(),t->getVertex(2)->y(),
- t->getVertex(3)->x(),t->getVertex(3)->y(),
- t->getVertex(4)->x(),t->getVertex(4)->y(),
- t->getVertex(5)->x(),t->getVertex(5)->y(),
- J1,J2,J3,J4,J5,J6);
- fclose(f);
- getchar();
- }
- */
- return std::min(m4, min_interm);
- }
- return min_interm;
-}
-
-double mesh_functional_distorsion_pN(MElement *t)
-{
- const JacobianBasis *jac = t->getJacobianFuncSpace();
- fullVector<double>Ji(jac->getNumJacNodes());
- for (int i=0;i<jac->getNumJacNodes();i++){
- double u = jac->getPoints()(i,0);
- double v = jac->getPoints()(i,1);
- if (t->getType() == TYPE_QUA){
- u = -1 + 2*u;
- v = -1 + 2*v;
- }
-
- Ji(i) = mesh_functional_distorsion_2D(t,u,v);
- }
-
- fullVector<double> Bi( jac->getNumJacNodes() );
- jac->lag2Bez(Ji,Bi);
- return *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
-}
-
-double qmDistorsionOfMapping (MTriangle *e)
-{
- // return 1.0;
- if (e->getPolynomialOrder() == 1) return 1.0;
- else if (e->getPolynomialOrder() == 2) {
- // const double exact = mesh_functional_distorsion_p2_exact(e);
- const double bezier= mesh_functional_distorsion_pN(e);
- // if (bezier < .1){
- // const double bezier_refined= mesh_functional_distorsion_p2_bezier_refined(e);
- // return bezier_refined;
- // }
- /*
- if (exact < .99 || exact > 1.01){
- FILE *f = fopen ("statistics.dat","a");
- fprintf(f,"%12.5E %12.5E %12.5E\n",exact,bezier,bezier_refined);
- fclose(f);
-
- if (exact > 0 && bezier < 0){
- f = fopen ("t.pos","w");
- double J1 =mesh_functional_distorsion(e,0.0,0.0);
- double J2 =mesh_functional_distorsion(e,1.0,0.0);
- double J3 =mesh_functional_distorsion(e,0.0,1.0);
- double J4 =mesh_functional_distorsion(e,0.5,0.0);
- double J5 =mesh_functional_distorsion(e,0.5,0.5);
- double J6 =mesh_functional_distorsion(e,0.0,0.5);
- fprintf(f,"ST2(%g,%g,0,%g,%g,0,%g,%g,0,%g,%g,0,%g,%g,0,%g,%g,0){%g,%g,%g,%g,%g,%g}\n",
- e->getVertex(0)->x(),e->getVertex(0)->y(),
- e->getVertex(1)->x(),e->getVertex(1)->y(),
- e->getVertex(2)->x(),e->getVertex(2)->y(),
- e->getVertex(3)->x(),e->getVertex(3)->y(),
- e->getVertex(4)->x(),e->getVertex(4)->y(),
- e->getVertex(5)->x(),e->getVertex(5)->y(),
- J1,J2,J3,J4,J5,J6);
- fclose(f);
- getchar();
- }
- }
- */
- return bezier;
- }
- else return mesh_functional_distorsion_pN(e);
-}
-
-double qmDistorsionOfMapping (MQuadrangle *e)
-{
- // return 1.0;
- if (e->getPolynomialOrder() == 1) return 1.0;
- else return mesh_functional_distorsion_pN(e);
-}
-
-
-double mesh_functional_distorsion_3D(MElement *t, double u, double v, double w)
-{
- // compute uncurved element jacobian d_u x and d_v x
- double mat[3][3];
- const double det1 = t->getPrimaryJacobian(u, v, w, mat);
-// const double det1 = det3x3(mat);
- const double detN = t->getJacobian(u, v, w, mat);
-// const double detN = det3x3(mat);
- if (det1 == 0 || detN == 0) return 0;
- double dist = detN / det1;
- return dist;
-}
-
-double qmDistorsionOfMapping(MTetrahedron *t)
-{
- const JacobianBasis *jac = t->getJacobianFuncSpace();
- fullVector<double>Ji(jac->getNumJacNodes());
- for (int i=0;i<jac->getNumJacNodes();i++){
- const double u = jac->getPoints()(i,0);
- const double v = jac->getPoints()(i,1);
- const double w = jac->getPoints()(i,2);
- Ji(i) = mesh_functional_distorsion_3D(t,u,v,w);
- }
-
- fullVector<double> Bi( jac->getNumJacNodes() );
- jac->lag2Bez(Ji,Bi);
- /*
- jac->matrixLag2Bez.print("Lag2Bez");
-
- jac->points.print("Points");
- t->getFunctionSpace(t->getPolynomialOrder())->points.print("lagrangianNodes");
- t->getFunctionSpace(t->getPolynomialOrder())->monomials.print("Monomials");
- t->getFunctionSpace(t->getPolynomialOrder())->coefficients.print("Coefficients");
- */
-
- return *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
-}
-
double qmTriangleAngles (MTriangle *e) {
double a = 500;
double worst_quality = std::numeric_limits<double>::max();
diff --git a/Mesh/qualityMeasures.h b/Mesh/qualityMeasures.h
index dfc8984e0f8780b078d86b9358e0ce8613535de0..ce740497283661a9ce16c850ec89dca16b6e1fd5 100644
--- a/Mesh/qualityMeasures.h
+++ b/Mesh/qualityMeasures.h
@@ -17,10 +17,6 @@ class MElement;
enum qualityMeasure4Triangle {QMTRI_RHO, QMTRI_COND};
enum qualityMeasure4Tet{QMTET_1, QMTET_2, QMTET_3, QMTET_ONE, QMTET_COND};
-double qmDistorsionOfMapping (MQuadrangle *e);
-double qmDistorsionOfMapping (MTriangle *e);
-double qmDistorsionOfMapping (MTetrahedron *e);
-
double qmTriangleAngles (MTriangle *e);
double qmQuadrangleAngles (MQuadrangle *e);
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 0ebc8979fc3733e48bcd152d325fc7c025d186cb..95f49ec23a58b35e57fefa3aac89ea48109c0434 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -11,7 +11,7 @@
#include "pyramidalBasis.h"
#include "pointsGenerators.h"
#include "BasisFactory.h"
-#include "MElement.h"
+#include "Numeric.h"
// Bezier Exponents
static fullMatrix<double> generate1DExponents(int order)
@@ -935,6 +935,25 @@ static fullMatrix<double> generateSubDivisor
+// Convert Bezier points to Lagrange points
+static fullMatrix<double> bez2LagPoints(bool dim1, bool dim2, bool dim3, const fullMatrix<double> &bezierPoints)
+{
+
+ const int numPt = bezierPoints.size1();
+ fullMatrix<double> lagPoints(numPt,3);
+
+ for (int i=0; i<numPt; i++) {
+ lagPoints(i,0) = dim1 ? -1. + 2*bezierPoints(i,0) : bezierPoints(i,0);
+ lagPoints(i,1) = dim2 ? -1. + 2*bezierPoints(i,1) : bezierPoints(i,1);
+ lagPoints(i,2) = dim3 ? -1. + 2*bezierPoints(i,2) : bezierPoints(i,2);
+ }
+
+ return lagPoints;
+
+}
+
+
+
std::map<int, bezierBasis> bezierBasis::_bbs;
const bezierBasis *bezierBasis::find(int tag)
{
@@ -948,9 +967,10 @@ const bezierBasis *bezierBasis::find(int tag)
if (MElement::ParentTypeFromTag(tag) == TYPE_PYR) {
B.dim = 3;
B.numLagPts = 5;
- B.points = gmshGeneratePointsPyramid(B.order,false);
- B.matrixLag2Bez.resize(B.points.size1(),B.points.size1(),0.);
- B.matrixBez2Lag.resize(B.points.size1(),B.points.size1(),0.);
+ B.bezierPoints = gmshGeneratePointsPyramid(B.order,false);
+ B.lagPoints = B.bezierPoints;
+ B.matrixLag2Bez.resize(B.bezierPoints.size1(),B.bezierPoints.size1(),0.);
+ B.matrixBez2Lag.resize(B.bezierPoints.size1(),B.bezierPoints.size1(),0.);
for (int i=0; i<B.matrixBez2Lag.size1(); ++i) {
B.matrixLag2Bez(i,i) = 1.;
B.matrixBez2Lag(i,i) = 1.;
@@ -965,14 +985,16 @@ const bezierBasis *bezierBasis::find(int tag)
B.dim = 0;
B.numLagPts = 1;
B.exponents = generate1DExponents(0);
- B.points = generate1DPoints(0);
+ B.bezierPoints = generate1DPoints(0);
+ B.lagPoints = B.bezierPoints;
dimSimplex = 0;
break;
case TYPE_LIN : {
B.dim = 1;
B.numLagPts = 2;
B.exponents = generate1DExponents(B.order);
- B.points = generate1DPoints(B.order);
+ B.bezierPoints = generate1DPoints(B.order);
+ B.lagPoints = bez2LagPoints(true,false,false,B.bezierPoints);
dimSimplex = 0;
subPoints = generateSubPointsLine(0);
break;
@@ -981,7 +1003,8 @@ const bezierBasis *bezierBasis::find(int tag)
B.dim = 2;
B.numLagPts = 3;
B.exponents = generateExponentsTriangle(B.order);
- B.points = gmshGeneratePointsTriangle(B.order,false);
+ B.bezierPoints = gmshGeneratePointsTriangle(B.order,false);
+ B.lagPoints = B.bezierPoints;
dimSimplex = 2;
subPoints = generateSubPointsTriangle(B.order);
break;
@@ -990,7 +1013,8 @@ const bezierBasis *bezierBasis::find(int tag)
B.dim = 2;
B.numLagPts = 4;
B.exponents = generateExponentsQuad(B.order);
- B.points = generatePointsQuad(B.order,false);
+ B.bezierPoints = generatePointsQuad(B.order,false);
+ B.lagPoints = bez2LagPoints(true,true,false,B.bezierPoints);
dimSimplex = 0;
subPoints = generateSubPointsQuad(B.order);
// B.points.print("points");
@@ -1000,7 +1024,8 @@ const bezierBasis *bezierBasis::find(int tag)
B.dim = 3;
B.numLagPts = 4;
B.exponents = generateExponentsTetrahedron(B.order);
- B.points = gmshGeneratePointsTetrahedron(B.order,false);
+ B.bezierPoints = gmshGeneratePointsTetrahedron(B.order,false);
+ B.lagPoints = B.bezierPoints;
dimSimplex = 3;
subPoints = generateSubPointsTetrahedron(B.order);
break;
@@ -1009,7 +1034,8 @@ const bezierBasis *bezierBasis::find(int tag)
B.dim = 3;
B.numLagPts = 6;
B.exponents = generateExponentsPrism(B.order);
- B.points = generatePointsPrism(B.order, false);
+ B.bezierPoints = generatePointsPrism(B.order, false);
+ B.lagPoints = bez2LagPoints(false,false,true,B.bezierPoints);
dimSimplex = 2;
subPoints = generateSubPointsPrism(B.order);
break;
@@ -1018,7 +1044,8 @@ const bezierBasis *bezierBasis::find(int tag)
B.dim = 3;
B.numLagPts = 8;
B.exponents = generateExponentsHex(B.order);
- B.points = generatePointsHex(B.order, false);
+ B.bezierPoints = generatePointsHex(B.order, false);
+ B.lagPoints = bez2LagPoints(true,true,true,B.bezierPoints);
dimSimplex = 0;
subPoints = generateSubPointsHex(B.order);
break;
@@ -1028,16 +1055,17 @@ const bezierBasis *bezierBasis::find(int tag)
B.dim = 3;
B.numLagPts = 4;
B.exponents = generateExponentsTetrahedron(0);
- B.points = gmshGeneratePointsTetrahedron(0, false);
+ B.bezierPoints = gmshGeneratePointsTetrahedron(0, false);
+ B.lagPoints = B.bezierPoints;
dimSimplex = 3;
subPoints = generateSubPointsTetrahedron(0);
break;
}
}
- B.matrixBez2Lag = generateBez2LagMatrix(B.exponents, B.points, B.order, dimSimplex);
+ B.matrixBez2Lag = generateBez2LagMatrix(B.exponents, B.bezierPoints, B.order, dimSimplex);
B.matrixBez2Lag.invert(B.matrixLag2Bez);
B.subDivisor = generateSubDivisor(B.exponents, subPoints, B.matrixLag2Bez, B.order, dimSimplex);
- B.numDivisions = (int) pow(2., (int) B.points.size2());
+ B.numDivisions = (int) pow(2., (int) B.bezierPoints.size2());
}
return &B;
@@ -1120,7 +1148,7 @@ JacobianBasis::JacobianBasis(int tag)
xBar /= numPrimMapNodes;
yBar /= numPrimMapNodes;
zBar /= numPrimMapNodes;
- double barDPsi[numPrimMapNodes][3];
+ double (*barDPsi)[3] = new double[numPrimMapNodes][3];
primMapBasis->df(xBar, yBar, zBar, barDPsi);
primGradShapeBarX.resize(numPrimMapNodes);
primGradShapeBarY.resize(numPrimMapNodes);
@@ -1130,6 +1158,7 @@ JacobianBasis::JacobianBasis(int tag)
primGradShapeBarY(j) = barDPsi[j][1];
primGradShapeBarZ(j) = barDPsi[j][2];
}
+ delete[] barDPsi;
}
@@ -1192,7 +1221,7 @@ double JacobianBasis::getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMa
dxyzdXbar(2) += primGradShapeBarX(j)*nodesXYZ(j,2);
dxyzdYbar(0) += primGradShapeBarY(j)*nodesXYZ(j,0);
dxyzdYbar(1) += primGradShapeBarY(j)*nodesXYZ(j,1);
- dxyzdYbar(2) += primGradShapeBarX(j)*nodesXYZ(j,2);
+ dxyzdYbar(2) += primGradShapeBarY(j)*nodesXYZ(j,2);
}
result(0,2) = dxyzdXbar(0) * dxyzdYbar(1) - dxyzdXbar(1) * dxyzdYbar(0);
@@ -1309,6 +1338,165 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ,
+// Calculate (signed) Jacobian at mapping's nodes for one element, given vectors for regularization
+// of line and surface Jacobians in 3D
+void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &JDJ) const
+{
+
+ switch (bezier->dim) {
+
+ case 0 : {
+ for (int i = 0; i < numJacNodes; i++) {
+ for (int j = 0; j < numMapNodes; j++) {
+ JDJ (i,j) = 0.;
+ JDJ (i,j+1*numMapNodes) = 0.;
+ JDJ (i,j+2*numMapNodes) = 0.;
+ }
+ JDJ(i,3*numMapNodes) = 1.;
+ }
+ break;
+ }
+
+ case 1 : {
+ fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3);
+ gradShapeMatX.mult(nodesXYZ, dxyzdX);
+ for (int i = 0; i < numJacNodes; i++) {
+ const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
+ const double &dxdY = normals(0,0), &dydY = normals(0,1), &dzdY = normals(0,2);
+ const double &dxdZ = normals(1,0), &dydZ = normals(1,1), &dzdZ = normals(1,2);
+ for (int j = 0; j < numMapNodes; j++) {
+ const double &dPhidX = gradShapeMatX(i,j);
+ JDJ (i,j) = dPhidX * dydY * dzdZ + dPhidX * dzdY * dydZ;
+ JDJ (i,j+1*numMapNodes) = dPhidX * dzdY * dxdZ - dPhidX * dxdY * dzdZ;
+ JDJ (i,j+2*numMapNodes) = dPhidX * dxdY * dydZ - dPhidX * dydY * dxdZ;
+ }
+ JDJ(i,3*numMapNodes) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ);
+ }
+ break;
+ }
+
+ case 2 : {
+ fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3);
+ gradShapeMatX.mult(nodesXYZ, dxyzdX);
+ gradShapeMatY.mult(nodesXYZ, dxyzdY);
+ for (int i = 0; i < numJacNodes; i++) {
+ const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
+ const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
+ const double &dxdZ = normals(0,0), &dydZ = normals(0,1), &dzdZ = normals(0,2);
+ for (int j = 0; j < numMapNodes; j++) {
+ const double &dPhidX = gradShapeMatX(i,j);
+ const double &dPhidY = gradShapeMatY(i,j);
+ JDJ (i,j) =
+ dPhidX * dydY * dzdZ + dzdX * dPhidY * dydZ +
+ dPhidX * dzdY * dydZ - dydX * dPhidY * dzdZ;
+ JDJ (i,j+1*numMapNodes) =
+ dxdX * dPhidY * dzdZ +
+ dPhidX * dzdY * dxdZ - dzdX * dPhidY * dxdZ
+ - dPhidX * dxdY * dzdZ;
+ JDJ (i,j+2*numMapNodes) =
+ dPhidX * dxdY * dydZ +
+ dydX * dPhidY * dxdZ - dPhidX * dydY * dxdZ -
+ dxdX * dPhidY * dydZ;
+ }
+ JDJ(i,3*numMapNodes) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ);
+ }
+ break;
+ }
+
+ case 3 : {
+ fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3), dxyzdZ(numJacNodes,3);
+ gradShapeMatX.mult(nodesXYZ, dxyzdX);
+ gradShapeMatY.mult(nodesXYZ, dxyzdY);
+ gradShapeMatZ.mult(nodesXYZ, dxyzdZ);
+ for (int i = 0; i < numJacNodes; i++) {
+ const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
+ const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
+ const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2);
+ for (int j = 0; j < numMapNodes; j++) {
+ const double &dPhidX = gradShapeMatX(i,j);
+ const double &dPhidY = gradShapeMatY(i,j);
+ const double &dPhidZ = gradShapeMatZ(i,j);
+ JDJ (i,j) =
+ dPhidX * dydY * dzdZ + dzdX * dPhidY * dydZ +
+ dydX * dzdY * dPhidZ - dzdX * dydY * dPhidZ -
+ dPhidX * dzdY * dydZ - dydX * dPhidY * dzdZ;
+ JDJ (i,j+1*numMapNodes) =
+ dxdX * dPhidY * dzdZ + dzdX * dxdY * dPhidZ +
+ dPhidX * dzdY * dxdZ - dzdX * dPhidY * dxdZ -
+ dxdX * dzdY * dPhidZ - dPhidX * dxdY * dzdZ;
+ JDJ (i,j+2*numMapNodes) =
+ dxdX * dydY * dPhidZ + dPhidX * dxdY * dydZ +
+ dydX * dPhidY * dxdZ - dPhidX * dydY * dxdZ -
+ dxdX * dPhidY * dydZ - dydX * dxdY * dPhidZ;
+ }
+ JDJ(i,3*numMapNodes) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ);
+ }
+ break;
+ }
+
+ }
+
+}
+
+
+
+void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &nodesXYZStraight,
+ fullVector<double> &lambdaJ , fullMatrix<double> &gradLambdaJ) const
+{
+
+ // jacobian of the straight elements (only triangles for now)
+ SPoint3 v0(nodesXYZ(0,0),nodesXYZ(0,1),nodesXYZ(0,2));
+ SPoint3 v1(nodesXYZ(1,0),nodesXYZ(1,1),nodesXYZ(1,2));
+ SPoint3 v2(nodesXYZ(2,0),nodesXYZ(2,1),nodesXYZ(2,2));
+ SPoint3 *IXYZ[3] = {&v0, &v1, &v2};
+ double jaci[2][2] = {
+ {IXYZ[1]->x() - IXYZ[0]->x(), IXYZ[2]->x() - IXYZ[0]->x()},
+ {IXYZ[1]->y() - IXYZ[0]->y(), IXYZ[2]->y() - IXYZ[0]->y()}
+ };
+ double invJaci[2][2];
+ inv2x2(jaci, invJaci);
+
+ for (int l = 0; l < numJacNodes; l++) {
+ double jac[2][2] = {{0., 0.}, {0., 0.}};
+ for (int i = 0; i < numMapNodes; i++) {
+ const double &dPhidX = gradShapeMatX(l,i);
+ const double &dPhidY = gradShapeMatY(l,i);
+ const double dpsidx = dPhidX * invJaci[0][0] + dPhidY * invJaci[1][0];
+ const double dpsidy = dPhidX * invJaci[0][1] + dPhidY * invJaci[1][1];
+ jac[0][0] += nodesXYZ(i,0) * dpsidx;
+ jac[0][1] += nodesXYZ(i,0) * dpsidy;
+ jac[1][0] += nodesXYZ(i,1) * dpsidx;
+ jac[1][1] += nodesXYZ(i,1) * dpsidy;
+ }
+ const double dxdx = jac[0][0] * jac[0][0] + jac[0][1] * jac[0][1];
+ const double dydy = jac[1][0] * jac[1][0] + jac[1][1] * jac[1][1];
+ const double dxdy = jac[0][0] * jac[1][0] + jac[0][1] * jac[1][1];
+ const double sqr = sqrt((dxdx - dydy) * (dxdx - dydy) + 4 * dxdy * dxdy);
+ const double osqr = sqr > 1e-8 ? 1/sqr : 0;
+ lambdaJ(l) = 0.5 * (dxdx + dydy - sqr);
+ const double axx = (1 - (dxdx - dydy) * osqr) * jac[0][0] - 2 * dxdy * osqr * jac[1][0];
+ const double axy = (1 - (dxdx - dydy) * osqr) * jac[0][1] - 2 * dxdy * osqr * jac[1][1];
+ const double ayx = -2 * dxdy * osqr * jac[0][0] + (1 - (dydy - dxdx) * osqr) * jac[1][0];
+ const double ayy = -2 * dxdy * osqr * jac[0][1] + (1 - (dydy - dxdx) * osqr) * jac[1][1];
+ const double axixi = axx * invJaci[0][0] + axy * invJaci[0][1];
+ const double aetaeta = ayx * invJaci[1][0] + ayy * invJaci[1][1];
+ const double aetaxi = ayx * invJaci[0][0] + ayy * invJaci[0][1];
+ const double axieta = axx * invJaci[1][0] + axy * invJaci[1][1];
+ for (int i = 0; i < numMapNodes; i++) {
+ const double &dPhidX = gradShapeMatX(l,i);
+ const double &dPhidY = gradShapeMatY(l,i);
+ gradLambdaJ(l, i + 0 * numMapNodes) = axixi * dPhidX + axieta * dPhidY;
+ gradLambdaJ(l, i + 1 * numMapNodes) = aetaxi * dPhidX + aetaeta * dPhidY;
+ }
+ }
+
+}
+
+
+
// Calculate (signed) Jacobian at mapping's nodes for several elements
// TODO: Optimize and test 1D & 2D
void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 801ccb63c908e842669fb6332369b332d82299d0..c24fc3e5c9b207d8e39bcc284a5b1aa713b04dc1 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -21,7 +21,7 @@ class bezierBasis {
int numDivisions;
// the 'numLagPts' first exponents are related to 'Lagrangian' values
fullMatrix<double> exponents; // exponents of Bezier FF
- fullMatrix<double> points; // sampling point
+ fullMatrix<double> bezierPoints, lagPoints; // sampling point
fullMatrix<double> matrixLag2Bez, matrixBez2Lag;
fullMatrix<double> subDivisor;
static const bezierBasis *find(int);
@@ -35,12 +35,14 @@ class JacobianBasis {
fullVector<double> primGradShapeBarX, primGradShapeBarY, primGradShapeBarZ;
fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
int numJacNodes, numMapNodes, numPrimJacNodes, numPrimMapNodes;
+ inline const fullMatrix<double> &getPoints() const { return bezier->lagPoints; }
public :
static const JacobianBasis *find(int);
JacobianBasis(int tag);
- inline int getNumJacNodes() const { return gradShapeMatX.size1(); }
- inline int getNumMapNodes() const { return gradShapeMatX.size2(); }
- inline const fullMatrix<double> &getPoints() const { return bezier->points; }
+ inline int getNumJacNodes() const { return numJacNodes; }
+ inline int getNumMapNodes() const { return numMapNodes; }
+ inline int getNumPrimJacNodes() const { return numPrimJacNodes; }
+ inline int getNumPrimMapNodes() const { return numPrimMapNodes; }
inline int getNumDivisions() const { return bezier->numDivisions; }
inline int getNumSubNodes() const { return bezier->subDivisor.size1(); }
inline int getNumLagPts() const { return bezier->numLagPts; }
@@ -48,6 +50,11 @@ class JacobianBasis {
double getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const;
void getSignedJacobian(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &normals, fullVector<double> &jacobian) const;
+ void getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const;
+ void getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &nodesXYZStraight,
+ fullVector<double> &lambdaJ , fullMatrix<double> &gradLambdaJ) const;
void getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
void getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 12c2c6dedfcacdf4ee4e2644021c47319bf8c1fb..08398951b89659e78805abe5d62e3a90f43b78e6 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -8,28 +8,6 @@
-std::map<int, fullMatrix<double> > Mesh::_gradShapeFunctions;
-
-
-
-fullMatrix<double> Mesh::computeGSF(const nodalBasis *lagrange, const JacobianBasis *jac)
-{
- // bezier points are defined in the [0,1] x [0,1] quad
- fullMatrix<double> bezierPoints = jac->getPoints();
- if (lagrange->parentType == TYPE_QUA) {
- for (int i = 0; i < bezierPoints.size1(); ++i) {
- bezierPoints(i, 0) = -1 + 2 * bezierPoints(i, 0);
- bezierPoints(i, 1) = -1 + 2 * bezierPoints(i, 1);
- }
- }
-
- fullMatrix<double> allDPsi;
- lagrange->df(bezierPoints, allDPsi);
- return allDPsi;
-}
-
-
-
Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFix, int method) :
_ge(ge)
{
@@ -67,16 +45,12 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
_indPCEl.resize(nElements);
int iEl = 0;
for(std::set<MElement*>::const_iterator it = els.begin(); it != els.end(); ++it, ++iEl) {
- MElement *el = *it;
- _el[iEl] = el;
- const nodalBasis *lagrange = el->getFunctionSpace();
- const JacobianBasis *jac = JacobianBasis::find(lagrange->type);
- if (_gradShapeFunctions.find(lagrange->type) == _gradShapeFunctions.end())
- _gradShapeFunctions[lagrange->type] = computeGSF(lagrange, jac);
+ _el[iEl] = *it;
+ const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace();
_nBezEl[iEl] = jac->getNumJacNodes();
- _nNodEl[iEl] = lagrange->points.size1();
- for (int iVEl = 0; iVEl < lagrange->points.size1(); iVEl++) {
- MVertex *vert = el->getVertex(iVEl);
+ _nNodEl[iEl] = jac->getNumMapNodes();
+ for (int iVEl = 0; iVEl < jac->getNumMapNodes(); iVEl++) {
+ MVertex *vert = _el[iEl]->getVertex(iVEl);
int iV = addVert(vert);
_el2V[iEl].push_back(iV);
const int nPCV = _pc->nCoord(vert);
@@ -105,8 +79,8 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
// Set normals to 2D elements (with magnitude of inverse Jacobian) or initial Jacobians of 3D elements
if (_dim == 2) {
- _normEl.resize(nEl());
- for (int iEl = 0; iEl < nEl(); iEl++) _normEl[iEl] = getNormalEl(iEl);
+ _scaledNormEl.resize(nEl());
+ for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(iEl);
}
else {
_invStraightJac.resize(nEl(),1.);
@@ -125,34 +99,25 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
-SVector3 Mesh::getNormalEl(int iEl)
+void Mesh::calcScaledNormalEl2D(int iEl)
{
- switch (_el[iEl]->getType()) {
- case TYPE_TRI: {
- const int iV0 = _el2V[iEl][0], iV1 = _el2V[iEl][1], iV2 = _el2V[iEl][2];
- SVector3 v10 = _xyz[iV1]-_xyz[iV0], v20 = _xyz[iV2]-_xyz[iV0];
- SVector3 n = crossprod(v10, v20);
- double xxx = n.norm();
- n *= 1./(xxx*xxx);
- return n;
- break;
- }
- case TYPE_QUA: {
- const int iV0 = _el2V[iEl][0], iV1 = _el2V[iEl][1], iV3 = _el2V[iEl][3];
- SVector3 v10 = _xyz[iV1]-_xyz[iV0], v30 = _xyz[iV3]-_xyz[iV0];
- SVector3 n = crossprod(v10, v30);
- double xxx = n.norm();
- n *= 4./(xxx*xxx);
- return n;
- break;
- }
- default:
- std::cout << "ERROR: getNormalEl: Unknown element type" << std::endl;
- break;
+ const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace();
+
+ fullMatrix<double> primNodesXYZ(jac->getNumPrimMapNodes(),3);
+ for (int i=0; i<jac->getNumPrimMapNodes(); i++) {
+ const int &iV = _el2V[iEl][i];
+ primNodesXYZ(i,0) = _xyz[iV].x();
+ primNodesXYZ(i,1) = _xyz[iV].y();
+ primNodesXYZ(i,2) = _xyz[iV].z();
}
- return SVector3(0.); // Just to avoid compilation warnings...
+ _scaledNormEl[iEl].resize(1,3);
+ const double norm = jac->getPrimNormal2D(primNodesXYZ,_scaledNormEl[iEl]);
+
+ _scaledNormEl[iEl](0,0) /= norm; // Re-scaling normal here is faster than an
+ _scaledNormEl[iEl](0,1) /= norm; // extra scaling operation on the Jacobian
+ _scaledNormEl[iEl](0,2) /= norm;
}
@@ -262,56 +227,32 @@ void Mesh::updateGEntityPositions()
}
+
void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda , std::vector<double> &gradLambda)
{
- fullMatrix<double> &gsf = _gradShapeFunctions[_el[iEl]->getTypeForMSH()];
- //const fullMatrix<double> &l2b = _lag2Bez[_el[iEl]->getTypeForMSH()];
- const int nbBez = _nBezEl[iEl];
- const int nbNod = _nNodEl[iEl];
- fullVector<double> lambdaJ(nbBez), lambdaB(nbBez);
- fullMatrix<double> gradLambdaJ(nbBez, 2 * nbNod);
- fullMatrix<double> gradLambdaB(nbBez, 2 * nbNod);
-
- // jacobian of the straight elements (only triangles for now)
- SPoint3 *IXYZ[3] = {&_ixyz[_el2V[iEl][0]], &_ixyz[_el2V[iEl][1]], &_ixyz[_el2V[iEl][2]]};
- double jaci[2][2] = {
- {IXYZ[1]->x() - IXYZ[0]->x(), IXYZ[2]->x() - IXYZ[0]->x()},
- {IXYZ[1]->y() - IXYZ[0]->y(), IXYZ[2]->y() - IXYZ[0]->y()}
- };
- double invJaci[2][2];
- inv2x2(jaci, invJaci);
+
+ const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace();
+ const int &numJacNodes = jacBasis->getNumJacNodes();
+ const int &numMapNodes = jacBasis->getNumMapNodes(), &numPrimMapNodes = jacBasis->getNumPrimMapNodes();
+ fullVector<double> lambdaJ(numJacNodes), lambdaB(numJacNodes);
+ fullMatrix<double> gradLambdaJ(numJacNodes, 2 * numMapNodes);
+ fullMatrix<double> gradLambdaB(numJacNodes, 2 * numMapNodes);
- for (int l = 0; l < nbBez; l++) {
- fullMatrix<double> dPsi(gsf, l * 3, 3);
- double jac[2][2] = {{0., 0.}, {0., 0.}};
- for (int i = 0; i < nbNod; i++) {
- int &iVi = _el2V[iEl][i];
- const double dpsidx = dPsi(i, 0) * invJaci[0][0] + dPsi(i, 1) * invJaci[1][0];
- const double dpsidy = dPsi(i, 0) * invJaci[0][1] + dPsi(i, 1) * invJaci[1][1];
- jac[0][0] += _xyz[iVi].x() * dpsidx;
- jac[0][1] += _xyz[iVi].x() * dpsidy;
- jac[1][0] += _xyz[iVi].y() * dpsidx;
- jac[1][1] += _xyz[iVi].y() * dpsidy;
- }
- const double dxdx = jac[0][0] * jac[0][0] + jac[0][1] * jac[0][1];
- const double dydy = jac[1][0] * jac[1][0] + jac[1][1] * jac[1][1];
- const double dxdy = jac[0][0] * jac[1][0] + jac[0][1] * jac[1][1];
- const double sqr = sqrt((dxdx - dydy) * (dxdx - dydy) + 4 * dxdy * dxdy);
- const double osqr = sqr > 1e-8 ? 1/sqr : 0;
- lambdaJ(l) = 0.5 * (dxdx + dydy - sqr);
- const double axx = (1 - (dxdx - dydy) * osqr) * jac[0][0] - 2 * dxdy * osqr * jac[1][0];
- const double axy = (1 - (dxdx - dydy) * osqr) * jac[0][1] - 2 * dxdy * osqr * jac[1][1];
- const double ayx = -2 * dxdy * osqr * jac[0][0] + (1 - (dydy - dxdx) * osqr) * jac[1][0];
- const double ayy = -2 * dxdy * osqr * jac[0][1] + (1 - (dydy - dxdx) * osqr) * jac[1][1];
- const double axixi = axx * invJaci[0][0] + axy * invJaci[0][1];
- const double aetaeta = ayx * invJaci[1][0] + ayy * invJaci[1][1];
- const double aetaxi = ayx * invJaci[0][0] + ayy * invJaci[0][1];
- const double axieta = axx * invJaci[1][0] + axy * invJaci[1][1];
- for (int i = 0; i < nbNod; i++) {
- gradLambdaJ(l, i + 0 * nbNod) = axixi * dPsi(i, 0) + axieta * dPsi(i, 1);
- gradLambdaJ(l, i + 1 * nbNod) = aetaxi * dPsi(i, 0) + aetaeta * dPsi(i, 1);
+ // Coordinates of nodes
+ fullMatrix<double> nodesXYZ(numMapNodes,3), nodesXYZStraight(numPrimMapNodes,3);
+ for (int i = 0; i < numMapNodes; i++) {
+ int &iVi = _el2V[iEl][i];
+ nodesXYZ(i,0) = _xyz[iVi].x();
+ nodesXYZ(i,1) = _xyz[iVi].y();
+ nodesXYZ(i,2) = _xyz[iVi].z();
+ if (i < numPrimMapNodes) {
+ nodesXYZStraight(i,0) = _ixyz[iVi].x();
+ nodesXYZStraight(i,1) = _ixyz[iVi].y();
+ nodesXYZStraight(i,2) = _ixyz[iVi].z();
}
}
+
+ jacBasis->getMetricMinAndGradients(nodesXYZ,nodesXYZStraight,lambdaJ,gradLambdaJ);
//l2b.mult(lambdaJ, lambdaB);
//l2b.mult(gradLambdaJ, gradLambdaB);
@@ -319,19 +260,19 @@ void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda , std::vec
gradLambdaB = gradLambdaJ;
int iPC = 0;
- std::vector<SPoint3> gXyzV(nbBez);
- std::vector<SPoint3> gUvwV(nbBez);
- for (int l = 0; l < nbBez; l++) {
+ std::vector<SPoint3> gXyzV(numJacNodes);
+ std::vector<SPoint3> gUvwV(numJacNodes);
+ for (int l = 0; l < numJacNodes; l++) {
lambda[l] = lambdaB(l);
}
- for (int i = 0; i < nbNod; i++) {
+ for (int i = 0; i < numMapNodes; i++) {
int &iFVi = _el2FV[iEl][i];
if (iFVi >= 0) {
- for (int l = 0; l < nbBez; l++) {
- gXyzV [l] = SPoint3(gradLambdaB(l,i+0*nbNod),gradLambdaB(l,i+1*nbNod),/*BDB(l,i+2*nbNod)*/ 0.);
+ for (int l = 0; l < numJacNodes; l++) {
+ gXyzV [l] = SPoint3(gradLambdaB(l,i+0*numMapNodes),gradLambdaB(l,i+1*numMapNodes),/*BDB(l,i+2*nbNod)*/ 0.);
}
_pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
- for (int l = 0; l < nbBez; l++) {
+ for (int l = 0; l < numJacNodes; l++) {
gradLambda[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
if (_nPCFV[iFVi] >= 2) gradLambda[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
if (_nPCFV[iFVi] == 3) gradLambda[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
@@ -341,120 +282,54 @@ void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda , std::vec
}
}
-/*
- A faster version that computes jacobians and their gradients
-
- Terms of jacobian are of the form J_{11} = dX / dU
-
- X = \sum \phi_j X_{j}
-
- d J_{11} / dX_k = d \phi_k / du
-
- */
-
void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ , std::vector<double> &gSJ) {
+
+ const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace();
+ const int &numJacNodes = jacBasis->getNumJacNodes();
+ const int &numMapNodes = jacBasis->getNumMapNodes();
+ fullMatrix<double> JDJ (numJacNodes,3*numMapNodes+1);
+ fullMatrix<double> BDB (numJacNodes,3*numMapNodes+1);
- // if (_dim == 2 && !projJac){
- // // printf("coucou\n");
- // gradScaledJac(iEl,gSJ);
- // scaledJac(iEl,sJ);
- // return;
- // }
-
- SVector3 n ;
- if (_dim == 2)
- n = _normEl[iEl];
-
- // std::vector<double> OLD = sJ;
- // std::vector<double> OLD = gSJ;
- // gradScaledJac(iEl,OLD);
- // scaledJac(iEl,OLD);
-
- const JacobianBasis *jacBasis = JacobianBasis::find(_el[iEl]->getTypeForMSH());
- fullMatrix<double> &gsf = _gradShapeFunctions[_el[iEl]->getTypeForMSH()];
- const int nbBez = _nBezEl[iEl];
- const int nbNod = _nNodEl[iEl];
- fullMatrix<double> JDJ (nbBez,3*nbNod+1);
- fullMatrix<double> BDB (nbBez,3*nbNod+1);
-
- double jac[3][3];
- for (int l = 0; l < nbBez; l++) {
- fullMatrix<double> dPsi(gsf, l * 3, 3);
-
- jac[0][0] = jac[0][1] = jac[0][2] = 0.;
- jac[1][0] = jac[1][1] = jac[1][2] = 0.;
- jac[2][0] = (_dim == 2) ? n.x() : 0.0; jac[2][1] = (_dim == 2) ? n.y() : 0.0; jac[2][2] = (_dim == 2) ? n.z() : 0.0;
-
- const double INVJ = (_dim == 3) ? _invStraightJac[iEl] : 1.0;
- for (int i = 0; i < nbNod; i++) {
- int &iVi = _el2V[iEl][i];
- const double x = _xyz[iVi].x();
- const double y = _xyz[iVi].y();
- const double z = _xyz[iVi].z();
- for (int k = 0; k < _dim; k++) {
- const double gg = dPsi(i, k);
- jac[k][0] += x * gg;
- jac[k][1] += y * gg;
- jac[k][2] += z * gg;
- }
- }
- for (int i = 0; i < nbNod; i++) {
- JDJ (l,i+0*nbNod) =
- ( dPsi(i, 0) * jac[1][1] * jac[2][2] + jac[0][2] * dPsi(i, 1) * jac[2][1] +
- jac[0][1] * jac[1][2] * dPsi(i, 2) - jac[0][2] * jac[1][1] * dPsi(i, 2) -
- dPsi(i, 0) * jac[1][2] * jac[2][1] - jac[0][1] * dPsi(i, 1) * jac[2][2])
- * INVJ;
- JDJ (l,i+1*nbNod) =
- (jac[0][0] * dPsi(i, 1) * jac[2][2] + jac[0][2] * jac[1][0] * dPsi(i, 2) +
- dPsi(i, 0) * jac[1][2] * jac[2][0] - jac[0][2] * dPsi(i, 1) * jac[2][0] -
- jac[0][0] * jac[1][2] * dPsi(i, 2) - dPsi(i, 0) * jac[1][0] * jac[2][2])
- * INVJ;
- JDJ (l,i+2*nbNod) =
- (jac[0][0] * jac[1][1] * dPsi(i, 2) + dPsi(i, 0) * jac[1][0] * jac[2][1] +
- jac[0][1] * dPsi(i, 1) * jac[2][0] - dPsi(i, 0) * jac[1][1] * jac[2][0] -
- jac[0][0] * dPsi(i, 1) * jac[2][1] - jac[0][1] * jac[1][0] * dPsi(i, 2))
- * INVJ;
- }
- const double dJ =
- 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];
- JDJ(l,3*nbNod) = dJ * INVJ;
+ // Coordinates of nodes
+ fullMatrix<double> nodesXYZ(numMapNodes,3), normals(_dim,3);
+ for (int i = 0; i < numMapNodes; i++) {
+ int &iVi = _el2V[iEl][i];
+ nodesXYZ(i,0) = _xyz[iVi].x();
+ nodesXYZ(i,1) = _xyz[iVi].y();
+ nodesXYZ(i,2) = _xyz[iVi].z();
}
+
+ // Calculate Jacobian and gradients, scale if 3D (already scaled by regularization normals in 2D)
+ jacBasis->getSignedJacAndGradients(nodesXYZ,_scaledNormEl[iEl],JDJ);
+ if (_dim == 3) JDJ.scale(_invStraightJac[iEl]);
- // (N_b x N_b) x (N_b x 3*N_n + 1)
+ // Transform Jacobian and gradients from Lagrangian to Bezier basis
jacBasis->lag2Bez(JDJ,BDB);
- // the scaled jacobian
- // printf("ELEMENT %d\n",iEl);
- for (int l = 0; l < nbBez; l++) {
- sJ [l] = BDB (l,3*nbNod);
- // printf("OLD %12.5E NEW %12.5E\n",OLD[l],sJ[l]);
- }
+ // Scaled jacobian
+ for (int l = 0; l < numJacNodes; l++) sJ [l] = BDB (l,3*numMapNodes);
- // gradients of the scaled jacobian
+ // Gradients of the scaled jacobian
int iPC = 0;
- std::vector<SPoint3> gXyzV(nbBez);
- std::vector<SPoint3> gUvwV(nbBez);
- for (int i = 0; i < nbNod; i++) {
+ std::vector<SPoint3> gXyzV(numJacNodes);
+ std::vector<SPoint3> gUvwV(numJacNodes);
+ for (int i = 0; i < numMapNodes; i++) {
int &iFVi = _el2FV[iEl][i];
if (iFVi >= 0) {
- for (int l = 0; l < nbBez; l++) {
- gXyzV [l] = SPoint3(BDB(l,i+0*nbNod),BDB(l,i+1*nbNod),BDB(l,i+2*nbNod));
- }
+ for (int l = 0; l < numJacNodes; l++)
+ gXyzV [l] = SPoint3(BDB(l,i+0*numMapNodes),BDB(l,i+1*numMapNodes),BDB(l,i+2*numMapNodes));
_pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
- for (int l = 0; l < nbBez; l++) {
- gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
- if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
- if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
- // if (fabs(OLD[indGSJ(iEl,l,iPC+2)]-gSJ[indGSJ(iEl,l,iPC+2)]) > 1.e-6)
- // printf("OLD = %12.5E new = %12.5E\n",OLD[indGSJ(iEl,l,iPC+2)],gSJ[indGSJ(iEl,l,iPC+2)]);
+ for (int l = 0; l < numJacNodes; l++) {
+ gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
+ if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
+ if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
}
iPC += _nPCFV[iFVi];
}
}
+
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index a7fa37a76848e2fb834036ea05b8d957423f3c67..0a716ae5a2760fc45126aaf2e3b4a21289d7e26a 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -60,31 +60,28 @@ private:
GEntity *_ge;
int _dim;
- int _nPC; // Total nb. of parametric coordinates
-
- std::vector<MElement*> _el; // List of elements
- std::vector<SVector3> _normEl; // Normals to 2D elements
- std::vector<double> _invStraightJac; // Initial Jacobians for 3D elements
- std::vector<MVertex*> _vert, _freeVert; // List of vert., free vert.
- std::vector<int> _fv2V; // Index of free vert. -> index of vert.
- std::vector<bool> _forced; // Is vertex forced?
- std::vector<SPoint3> _xyz, _ixyz; // Physical coord. of ALL vertices (current, straight, init.)
- std::vector<SPoint3> _uvw, _iuvw; // Parametric coord. of FREE vertices (current, straight, init.)
- std::vector<int> _startPCFV; // Start index of parametric coordinates for a free vertex
- std::vector<int> _nPCFV; // Number of parametric coordinates for a free vertex
- std::vector<std::vector<int> > _el2FV, _el2V; // Free vertices, vertices in element
- std::vector<int> _nBezEl, _nNodEl; // Number of Bezier poly. and nodes for an el.
- std::vector<std::vector<int> > _indPCEl; // Index of parametric coord. for an el.
+ int _nPC; // Total nb. of parametric coordinates
+
+ std::vector<MElement*> _el; // List of elements
+ std::vector<fullMatrix<double> > _scaledNormEl; // Normals to 2D elements for Jacobian regularization and scaling
+ std::vector<double> _invStraightJac; // Initial Jacobians for 3D elements
+ std::vector<MVertex*> _vert, _freeVert; // List of vert., free vert.
+ std::vector<int> _fv2V; // Index of free vert. -> index of vert.
+ std::vector<bool> _forced; // Is vertex forced?
+ std::vector<SPoint3> _xyz, _ixyz; // Physical coord. of ALL vertices (current, straight, init.)
+ std::vector<SPoint3> _uvw, _iuvw; // Parametric coord. of FREE vertices (current, straight, init.)
+ std::vector<int> _startPCFV; // Start index of parametric coordinates for a free vertex
+ std::vector<int> _nPCFV; // Number of parametric coordinates for a free vertex
+ std::vector<std::vector<int> > _el2FV, _el2V; // Free vertices, vertices in element
+ std::vector<int> _nBezEl, _nNodEl; // Number of Bezier poly. and nodes for an el.
+ std::vector<std::vector<int> > _indPCEl; // Index of parametric coord. for an el.
ParamCoord *_pc;
- bool projJac; // Using "projected" Jacobians or not
-
- static std::map<int, fullMatrix<double> > _gradShapeFunctions; // gradients of shape functions at Bezier points
+ bool projJac; // Using "projected" Jacobians or not
int addVert(MVertex* vert);
int addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVertex*> &toFix);
- SVector3 getNormalEl(int iEl);
- static fullMatrix<double> computeGSF(const nodalBasis *lagrange, const JacobianBasis *jac);
+ void calcScaledNormalEl2D(int iEl);
static inline int indJB2DBase(int nNod, int l, int i, int j) { return (l*nNod+i)*nNod+j; }
inline int indJB2D(int iEl, int l, int i, int j) { return indJB2DBase(_nNodEl[iEl],l,i,j); }
static inline int indJB3DBase(int nNod, int l, int i, int j, int m) { return ((l*nNod+i)*nNod+j)*nNod+m; }