Fixed bug (wrong basis type) in OptHOM. Fixed bug in computation of Bergot...

Fixed bug (wrong basis type) in OptHOM. Fixed bug in computation of Bergot basis gradient, still wrong.

Numeric/BergotBasis.cpp

@@ -27,7 +27,7 @@ BergotBasis::BergotBasis(int p): order(p)
   for (int i=0;i<=order;i++) {
     for (int j=0;j<=order;j++) {
       int mIJ = std::max(i,j);
-      for (int k=0;k<= (int) (order - mIJ);k++,index++) {
+      for (int k=0; k<=order-mIJ; k++,index++) {
 
         iOrder[index] = i;
         jOrder[index] = j;
@@ -56,13 +56,13 @@ BergotBasis::~BergotBasis()
 
 void BergotBasis::f(double u, double v, double w, double* val) const
 {
-  double uhat = (w == 1.) ? 1 : u/(1.-w);
+  double uhat = (w == 1.) ? 1. : u/(1.-w);
 
   std::vector<double> uFcts(order+1);
   //double uFcts[order+1];
   legendre.f(uhat,&(uFcts[0]));
 
-  double vhat = (w == 1.) ? 1 : v/(1.-w);
+  double vhat = (w == 1.) ? 1. : v/(1.-w);
   std::vector<double> vFcts(order+1);
   legendre.f(vhat,&(vFcts[0]));
 
@@ -83,9 +83,9 @@ void BergotBasis::f(double u, double v, double w, double* val) const
   for (int i=0;i<=order;i++) {
     for (int j=0;j<=order;j++) {
       int mIJ = std::max(i,j);
-      double fact = pow_int(1-w,(int) mIJ);
+      double fact = pow(1-w, mIJ);
       double* wf = (wFcts.find(mIJ))->second;
-      for (int k=0;k<=(int) (order - mIJ);k++,index++) {
+      for (int k=0; k<=order-mIJ; k++,index++) {
         val[index] = uFcts[i] * vFcts[j] * wf[k] * fact;
       }
     }
@@ -140,11 +140,11 @@ void BergotBasis::df(double u, double v, double w, double grads[][3]) const
       double* wf = (wFcts .find(mIJ))->second;
       double* wg = (wGrads.find(mIJ))->second;
 
-      double wPowM2 = pow_int(1.-w,((int) mIJ) - 2);
+      double wPowM2 = pow(1.-w, mIJ-2);
       double wPowM1 = wPowM2*(1.-w);
       double wPowM0 = wPowM1*(1.-w);
 
-      for (int k=0;k<= (int) (order - mIJ);k++,index++) {
+      for (int k=0; k<=order-mIJ; k++,index++) {
 
         grads[index][0] = uGrads[i] * vFcts[j]  * wf[k] * wPowM1;
 

Numeric/polynomialBasis.cpp

@@ -422,7 +422,6 @@ void polynomialBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) con
 {
   double dfv[1256][3];
   dfm.resize (coefficients.size1(), coord.size1() * 3, false);
-  int ii = 0;
   int dimCoord = coord.size2();
   for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
     df(coord(iPoint,0), dimCoord > 1 ? coord(iPoint, 1) : 0., dimCoord > 2 ? coord(iPoint, 2) : 0., dfv);
@@ -430,7 +429,6 @@ void polynomialBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) con
       dfm(i, iPoint * 3 + 0) = dfv[i][0];
       dfm(i, iPoint * 3 + 1) = dfv[i][1];
       dfm(i, iPoint * 3 + 2) = dfv[i][2];
-      ++ii;
     }
   }
 }

Numeric/pyramidalBasis.cpp

@@ -43,7 +43,7 @@ pyramidalBasis::~pyramidalBasis()
 
 
 
-inline void pyramidalBasis::f(double u, double v, double w, double *val) const
+void pyramidalBasis::f(double u, double v, double w, double *val) const
 {
 
   const int N = bergot->size();
@@ -60,7 +60,7 @@ inline void pyramidalBasis::f(double u, double v, double w, double *val) const
 
 
 
-inline void pyramidalBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf)
+void pyramidalBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf)
 {
 
   const int N = bergot->size(), NPts = coord.size1();
@@ -81,7 +81,7 @@ inline void pyramidalBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf)
 
 
 
-inline void pyramidalBasis::df(double u, double v, double w, double grads[][3]) const
+void pyramidalBasis::df(double u, double v, double w, double grads[][3]) const
 {
 
   const int N = bergot->size();
@@ -102,18 +102,18 @@ inline void pyramidalBasis::df(double u, double v, double w, double grads[][3])
 
 
 
-inline void pyramidalBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const
+void pyramidalBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const
 {
 
   const int N = bergot->size(), NPts = coord.size1();
 
   double dfv[N][3];
-  dfm.resize (NPts, 3*N, false);
+  dfm.resize (N, 3*NPts, false);
 
   for (int iPt=0; iPt<NPts; iPt++) {
     df(coord(iPt,0), coord(iPt,1), coord(iPt,2), dfv);
     for (int i = 0; i < N; i++) {
-      dfm(i, 3*iPt+0) = dfv[i][0];
+      dfm(i, 3*iPt) = dfv[i][0];
       dfm(i, 3*iPt+1) = dfv[i][1];
       dfm(i, 3*iPt+2) = dfv[i][2];
     }

contrib/HighOrderMeshOptimizer/OptHomMesh.cpp

@@ -13,7 +13,7 @@ std::map<int, fullMatrix<double> > Mesh::_lag2Bez;
 
 
 
-fullMatrix<double> Mesh::computeGSF(const polynomialBasis *lagrange, const bezierBasis *bezier)
+fullMatrix<double> Mesh::computeGSF(const nodalBasis *lagrange, const bezierBasis *bezier)
 {
   // bezier points are defined in the [0,1] x [0,1] quad
   fullMatrix<double> bezierPoints = bezier->points;
@@ -70,7 +70,7 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
   for(std::set<MElement*>::const_iterator it = els.begin(); it != els.end(); ++it, ++iEl) {
     MElement *el = *it;
     _el[iEl] = el;
-    const polynomialBasis *lagrange = (polynomialBasis*) el->getFunctionSpace();
+    const nodalBasis *lagrange = el->getFunctionSpace();
     const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier;
     if (_lag2Bez.find(lagrange->type) == _lag2Bez.end()) {
       _gradShapeFunctions[lagrange->type] = computeGSF(lagrange, bezier);

contrib/HighOrderMeshOptimizer/OptHomMesh.h

@@ -85,7 +85,7 @@ private:
   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 polynomialBasis *lagrange, const bezierBasis *bezier);
+  static fullMatrix<double> computeGSF(const nodalBasis *lagrange, const bezierBasis *bezier);
   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; }