cleanup

CMakeLists.txt

@@ -1419,8 +1419,7 @@ endif()
 # don't issue warnings for contributed libraries
 check_cxx_compiler_flag("-w" NOWARN)
 if(NOWARN)
-  file(GLOB_RECURSE NOWARN_SRC contrib/*.cpp contrib/*.cc contrib/*.cxx contrib/*.c
-      ${CMAKE_CURRENT_BINARY_DIR}/contrib/taucs/*.c)
+  file(GLOB_RECURSE NOWARN_SRC contrib/*.cpp contrib/*.cc contrib/*.cxx contrib/*.c)
   set_compile_flags(NOWARN_SRC "-w")
 endif()
 

Solver/elasticityTerm.cpp

@@ -37,8 +37,6 @@ void elasticityTerm::createData(MElement *e) const
     d.w.push_back(w);
     d.weight.push_back(weight);
   }
-  //  printf("coucou creation of a data for %d with %d
-  //  points\n",e->getTypeForMSH(),npts);
   _data[e->getTypeForMSH()] = d;
 }
 
@@ -87,7 +85,6 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
     BT.setAll(0.);
 
     if(se->getShapeEnrichement() == se->getTestEnrichement()) {
-      //      printf("coucou\n");
       for(int j = 0; j < nbSF; j++) {
         BT(j, 0) = B(0, j) = Grads[j][0];
         BT(j, 3) = B(3, j) = Grads[j][1];
@@ -125,7 +122,6 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
     BTH.gemm(BT, H);
     m.gemm(BTH, B, weight * detJ, 1.);
   }
-  return;
 }
 
 void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const

Solver/linearSystemPETSc.hpp

@@ -124,7 +124,7 @@ template <class scalar> void linearSystemPETSc<scalar>::preAllocateEntries()
   int blockSize = _getBlockSizeFromParameters();
   std::vector<int> nByRowDiag(_localSize), nByRowOffDiag(_localSize);
   if(_sparsity.getNbRows() == 0) {
-    PetscInt prealloc = 500;
+    PetscInt prealloc = 100;
     PetscBool set;
     PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set);
     prealloc = std::min(prealloc, _localSize);

contrib/HighOrderMeshOptimizer/HighOrderMeshElasticAnalogy.cpp

@@ -100,7 +100,7 @@ void HighOrderMeshElasticAnalogy(GModel *m, bool onlyVisible)
                  t2 - t1);
 }
 
-void highOrderTools::moveToStraightSidedLocation(MElement *e) const
+void highOrderTools::_moveToStraightSidedLocation(MElement *e) const
 {
   for(int i = 0; i < e->getNumVertices(); i++) {
     MVertex *v = e->getVertex(i);
@@ -239,21 +239,16 @@ void highOrderTools::ensureMinimumDistorsion(double threshold)
   ensureMinimumDistorsion(v, threshold);
 }
 
-void highOrderTools::computeMetricInfo(GFace *gf, MElement *e,
-                                       fullMatrix<double> &J,
-                                       fullMatrix<double> &JT,
-                                       fullVector<double> &D)
+void highOrderTools::_computeMetricInfo(GFace *gf, MElement *e,
+                                        fullMatrix<double> &J,
+                                        fullMatrix<double> &JT,
+                                        fullVector<double> &D)
 {
   int nbNodes = e->getNumVertices();
-  //  printf("ELEMENT --\n");
   for(int j = 0; j < nbNodes; j++) {
     SPoint2 param;
     reparamMeshVertexOnFace(e->getVertex(j), gf, param);
-    //    printf("%g %g vs %g %g %g\n",param.x(),param.y(),
-    //	   e->getVertex(j)->x(),e->getVertex(j)->y(),e->getVertex(j)->z());
-
     Pair<SVector3, SVector3> der = gf->firstDer(param);
-
     int XJ = j;
     int YJ = j + nbNodes;
     int ZJ = j + 2 * nbNodes;
@@ -283,11 +278,12 @@ void highOrderTools::computeMetricInfo(GFace *gf, MElement *e,
 
 void highOrderTools::applySmoothingTo(std::vector<MElement *> &all, GFace *gf)
 {
-#ifdef HAVE_TAUCS
-  linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+#if !defined(HAVE_PETSC)
+  Msg::Error("Elastic analogy smoother on surface %d requires PETSc",
+             gf->tag());
 #else
   linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
-#endif
+
   // compute the straight sided positions of high order nodes that are
   // on the edges of the face in the UV plane
   dofManager<double> myAssembler(lsys);
@@ -303,6 +299,7 @@ void highOrderTools::applySmoothingTo(std::vector<MElement *> &all, GFace *gf)
   }
 
   if(!v.size()) return;
+
   Msg::Info(
     "Smoothing surface %d (%d elements considered in elastic analogy, "
     "worst mapping %12.5E, %3d bad elements)",
@@ -323,7 +320,7 @@ void highOrderTools::applySmoothingTo(std::vector<MElement *> &all, GFace *gf)
 
   // fix all vertices that cannot move
   for(std::size_t i = 0; i < v.size(); i++) {
-    moveToStraightSidedLocation(v[i]);
+    _moveToStraightSidedLocation(v[i]);
     for(int j = 0; j < v[i]->getNumVertices(); j++) {
       MVertex *vert = v[i]->getVertex(j);
       if(vert->onWhat()->dim() < 2) {
@@ -344,12 +341,12 @@ void highOrderTools::applySmoothingTo(std::vector<MElement *> &all, GFace *gf)
     }
   }
 
-  double dx0 = smooth_metric_(v, gf, myAssembler, verticesToMove, El);
+  double dx0 = _smoothMetric(v, gf, myAssembler, verticesToMove, El);
   double dx = dx0;
   Msg::Debug(" dx0 = %12.5E", dx0);
   int iter = 0;
   while(0) {
-    double dx2 = smooth_metric_(v, gf, myAssembler, verticesToMove, El);
+    double dx2 = _smoothMetric(v, gf, myAssembler, verticesToMove, El);
     Msg::Debug(" dx2  = %12.5E", dx2);
     if(fabs(dx2 - dx) < 1.e-4 * dx0) break;
     if(iter++ > 2) break;
@@ -374,12 +371,13 @@ void highOrderTools::applySmoothingTo(std::vector<MElement *> &all, GFace *gf)
     }
   }
   delete lsys;
+#endif
 }
 
-double highOrderTools::smooth_metric_(std::vector<MElement *> &v, GFace *gf,
-                                      dofManager<double> &myAssembler,
-                                      std::set<MVertex *> &verticesToMove,
-                                      elasticityTerm &El)
+double highOrderTools::_smoothMetric(std::vector<MElement *> &v, GFace *gf,
+                                     dofManager<double> &myAssembler,
+                                     std::set<MVertex *> &verticesToMove,
+                                     elasticityTerm &El)
 {
   std::set<MVertex *>::iterator it;
   double dx = 0.0;
@@ -402,7 +400,7 @@ double highOrderTools::smooth_metric_(std::vector<MElement *> &v, GFace *gf,
       K33.setAll(0.0);
       SElement se(e);
       El.elementMatrix(&se, K33);
-      computeMetricInfo(gf, e, J32, J23, D3);
+      _computeMetricInfo(gf, e, J32, J23, D3);
       J23K33.gemm(J23, K33, 1, 0);
       K22.gemm(J23K33, J32, 1, 0);
       J23K33.mult(D3, R2);
@@ -440,7 +438,7 @@ double highOrderTools::smooth_metric_(std::vector<MElement *> &v, GFace *gf,
 
 highOrderTools::highOrderTools(GModel *gm) : _gm(gm), _tag(111), _dim(2)
 {
-  computeStraightSidedPositions();
+  _computeStraightSidedPositions();
 }
 
 highOrderTools::highOrderTools(GModel *gm, GModel *mesh, int order)
@@ -449,10 +447,10 @@ highOrderTools::highOrderTools(GModel *gm, GModel *mesh, int order)
   GeomMeshMatcher::instance()->forceTomatch(gm, mesh, 1.e-5);
   GeomMeshMatcher::instance()->destroy();
   SetOrderN(gm, order, false, false);
-  computeStraightSidedPositions();
+  _computeStraightSidedPositions();
 }
 
-void highOrderTools::computeStraightSidedPositions()
+void highOrderTools::_computeStraightSidedPositions()
 {
   _clean();
   // compute straigh sided positions that are actually X,Y,Z positions
@@ -605,26 +603,31 @@ void highOrderTools::computeStraightSidedPositions()
 
 // apply a displacement that does not create elements that are distorted over a
 // value "thres"
-double highOrderTools::apply_incremental_displacement(
+double highOrderTools::_applyIncrementalDisplacement(
   double max_incr, std::vector<MElement *> &v, bool mixed, double thres,
   char *meshName, std::vector<MElement *> &disto)
 {
-#ifdef HAVE_PETSC
+#if !defined(HAVE_PETSC)
+  Msg::Error("Elastic analogy smoother requires PETSc");
+#else
+  if(v.empty()) return 1.;
+
   // assume that the mesh is OK, yet already curved
-  // linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
   linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
-  lsys->setParameter("petscOptions", "-pc_type ilu");
-  lsys->setParameter("petscOptions", "-ksp_monitor");
-
+  char opt[256];
+  sprintf(opt, "-pc_type ilu -ksp_monitor -petsc_prealloc %d",
+          100 * (v[0]->getPolynomialOrder() + 2));
+  lsys->setParameter("petscOptions", opt);
+  Msg::Info("Assembling linear system...");
   dofManager<double> myAssembler(lsys);
   elasticityMixedTerm El_mixed(0, 1.0, .333, _tag);
   elasticityTerm El(0, 1.0, .333, _tag);
 
   std::set<MVertex *> _vertices;
 
-  //+++++++++ Boundary Conditions & Numbering +++++++++++++++++++++++++++++++
-  // fix all dof that correspond to vertices on the boundary
-  // the value is equal
+  // Boundary Conditions & Numbering
+
+  // fix all dof that correspond to vertices on the boundary the value is equal
   for(std::size_t i = 0; i < v.size(); i++) {
     for(int j = 0; j < v[i]->getNumVertices(); j++) {
       MVertex *vert = v[i]->getVertex(j);
@@ -632,7 +635,7 @@ double highOrderTools::apply_incremental_displacement(
     }
   }
 
-  //+++++++++ Fix d tr(eps) = 0 +++++++++++++++++++++++++++++++
+  // Fix d tr(eps) = 0
   if(mixed) {
     for(std::size_t i = 0; i < disto.size(); i++) {
       for(int j = 0; j < disto[i]->getNumVertices(); j++) {
@@ -658,7 +661,6 @@ double highOrderTools::apply_incremental_displacement(
       myAssembler.fixVertex(vert, 1, _tag, 0);
       myAssembler.fixVertex(vert, 2, _tag, 0);
     }
-    //    }
     if(_dim == 2) myAssembler.fixVertex(vert, 2, _tag, 0);
     // number vertices
     myAssembler.numberVertex(vert, 0, _tag);
@@ -679,6 +681,7 @@ double highOrderTools::apply_incremental_displacement(
       else
         El.addToMatrix(myAssembler, &se);
     }
+    Msg::Info("Solving linear system (%d unknowns)...", myAssembler.sizeOfR());
     // solve the system
     lsys->systemSolve();
   }
@@ -757,17 +760,21 @@ double highOrderTools::applySmoothingTo(std::vector<MElement *> &all,
   int ITER = 0;
   double minD;
   std::vector<MElement *> disto;
+
   // move the points to their straight sided locations
   for(std::size_t i = 0; i < all.size(); i++)
-    moveToStraightSidedLocation(all[i]);
+    _moveToStraightSidedLocation(all[i]);
+
   // apply the displacement
 
   // _gm->writeMSH("straightSided.msh");
 
   char sm[] = "sm.msh";
-  double percentage_of_what_is_left =
-    apply_incremental_displacement(1., all, mixed, -100000000, sm, all);
-  //  ensureMinimumDistorsion (all, threshold);
+  double percentage_of_what_is_left = _applyIncrementalDisplacement
+    (1., all, mixed, -100000000, sm, all);
+
+  // ensureMinimumDistorsion (all, threshold);
+
   return 1.;
 
   double percentage = 0.0;
@@ -775,7 +782,7 @@ double highOrderTools::applySmoothingTo(std::vector<MElement *> &all,
     char NN[256];
     sprintf(NN, "smoothing-%d.msh", ITER++);
     percentage_of_what_is_left =
-      apply_incremental_displacement(1., all, mixed, threshold, NN, all);
+      _applyIncrementalDisplacement(1., all, mixed, threshold, NN, all);
     percentage += (1. - percentage) * percentage_of_what_is_left / 100.;
     Msg::Info("Smoother was able to do %3d percent of the motion",
               (int)(percentage * 100.));

contrib/HighOrderMeshOptimizer/HighOrderMeshElasticAnalogy.h

@@ -44,12 +44,12 @@ class GFace;
 class GRegion;
 
 class highOrderTools {
+private:
   GModel *_gm;
   const int _tag;
   // contains the position of vertices in a straight sided mesh
   std::map<MVertex *, SVector3> _straightSidedLocation;
-  // contains the position of vertices in the best curvilinear mesh
-  // available
+  // contains the position of vertices in the best curvilinear mesh available
   std::map<MVertex *, SVector3> _targetLocation;
   int _dim;
   void _clean()
@@ -57,18 +57,18 @@ class highOrderTools {
     _straightSidedLocation.clear();
     _targetLocation.clear();
   }
-  double smooth_metric_(std::vector<MElement *> &v, GFace *gf,
-                        dofManager<double> &myAssembler,
-                        std::set<MVertex *> &verticesToMove,
-                        elasticityTerm &El);
-  void computeMetricInfo(GFace *gf, MElement *e, fullMatrix<double> &J,
-                         fullMatrix<double> &JT, fullVector<double> &D);
-  double apply_incremental_displacement(double max_incr,
-                                        std::vector<MElement *> &v, bool mixed,
-                                        double thres, char *meshName,
-                                        std::vector<MElement *> &disto);
-  void moveToStraightSidedLocation(MElement *) const;
-  void computeStraightSidedPositions();
+  double _smoothMetric(std::vector<MElement *> &v, GFace *gf,
+                       dofManager<double> &myAssembler,
+                       std::set<MVertex *> &verticesToMove,
+                       elasticityTerm &El);
+  void _computeMetricInfo(GFace *gf, MElement *e, fullMatrix<double> &J,
+                          fullMatrix<double> &JT, fullVector<double> &D);
+  double _applyIncrementalDisplacement(double max_incr,
+                                       std::vector<MElement *> &v, bool mixed,
+                                       double thres, char *meshName,
+                                       std::vector<MElement *> &disto);
+  void _moveToStraightSidedLocation(MElement *) const;
+  void _computeStraightSidedPositions();
 
 public:
   highOrderTools(GModel *gm);

contrib/HighOrderMeshOptimizer/HighOrderMeshOptimizer.cpp

@@ -307,7 +307,8 @@ void HighOrderMeshOptimizer(std::vector<GEntity *> &entities,
   bool order1 = false;
   for(std::size_t i = 0; i < entities.size(); i++){
     for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++){
-      if(entities[i]->getMeshElement(j)->getPolynomialOrder() < 2){
+      if(entities[i]->dim() > 0 &&
+         entities[i]->getMeshElement(j)->getPolynomialOrder() < 2){
         order1 = true;
         break;
       }

contrib/MeshOptimizer/MeshOptimizer.cpp

@@ -467,6 +467,10 @@ namespace {
                            bndElts[iPatch], par);
     }
     if(par.nCurses) displayResultTable(nbPatchSuccess, toOptimize.size());
+
+#if defined(_OPENMP)
+//#pragma omp parallel for schedule(dynamic)
+#endif
     for(int iPatch = 0; iPatch < toOptimize.size(); ++iPatch) {
       // Initialize optimization and output if asked
       if(par.nCurses) {
@@ -502,25 +506,37 @@ namespace {
 
       // Evaluate mesh and update it if (partial) success
       opt.updateResults();
-      if(newObjFunctionRange.size() == 0) {
-        newObjFunctionRange = opt.objFunction()->minMax();
-        objFunctionNames = opt.objFunction()->names();
-      }
-      else {
-        for(int i = 0; i < newObjFunctionRange.size(); i++) {
-          newObjFunctionRange[i].first = std::min(
-            newObjFunctionRange[i].first, opt.objFunction()->minMax()[i].first);
-          newObjFunctionRange[i].second =
-            std::max(newObjFunctionRange[i].second,
-                     opt.objFunction()->minMax()[i].second);
+
+#if defined(_OPENMP)
+//#pragma omp critical
+#endif
+      {
+        if(newObjFunctionRange.size() == 0) {
+          newObjFunctionRange = opt.objFunction()->minMax();
+          objFunctionNames = opt.objFunction()->names();
+        }
+        else {
+          for(int i = 0; i < newObjFunctionRange.size(); i++) {
+            newObjFunctionRange[i].first =
+              std::min(newObjFunctionRange[i].first,
+                       opt.objFunction()->minMax()[i].first);
+            newObjFunctionRange[i].second =
+              std::max(newObjFunctionRange[i].second,
+                       opt.objFunction()->minMax()[i].second);
+          }
         }
       }
+
       if(success >= 0) opt.patch.updateGEntityPositions();
 
-      //#pragma omp critical
-      par.success = std::min(par.success, success);
+#if defined(_OPENMP)
+//#pragma omp critical
+#endif
+      {
+        par.success = std::min(par.success, success);
+        nbPatchSuccess[success + 1]++;
+      }
 
-      nbPatchSuccess[success + 1]++;
       if(par.nCurses) {
         displayMinMaxVal(nbPatchSuccess, objFunctionNames, newObjFunctionRange);
         displayResultTable(nbPatchSuccess, toOptimize.size());
@@ -528,6 +544,7 @@ namespace {
                                   iPatch, -1);
       }
     }
+
     while(_patchHistory.size() > 0) {
       delete[] _patchHistory.back();
       _patchHistory.pop_back();