Changed handling of contributions in mesh optimizer (to make future parallelization easier)

contrib/HighOrderMeshOptimizer/OptHomObjContribMetricMin.h

@@ -13,6 +13,7 @@ class ObjContribMetricMin : public ObjContrib, public FuncType
 public:
   ObjContribMetricMin(double weight);
   virtual ~ObjContribMetricMin() {}
+  virtual ObjContrib *copy() const;
   virtual void initialize(Patch *mesh);
   virtual bool fail() { return false; }
   virtual bool addContrib(double &Obj, alglib::real_1d_array &gradObj);
@@ -36,6 +37,13 @@ ObjContribMetricMin<FuncType>::ObjContribMetricMin(double weight) :
 }
 
 
+template<class FuncType>
+ObjContrib *ObjContribMetricMin<FuncType>::copy() const
+{
+  return new ObjContribMetricMin<FuncType>(*this);
+}
+
+
 template<class FuncType>
 void ObjContribMetricMin<FuncType>::initialize(Patch *mesh)
 {

contrib/HighOrderMeshOptimizer/OptHomObjContribScaledJac.h

@@ -13,6 +13,7 @@ class ObjContribScaledJac : public ObjContrib, public FuncType
 public:
   ObjContribScaledJac(double weight);
   virtual ~ObjContribScaledJac() {}
+  virtual ObjContrib *copy() const;
   virtual void initialize(Patch *mesh);
   virtual bool fail() { return _min <= 0.; }
   virtual bool addContrib(double &Obj, alglib::real_1d_array &gradObj);
@@ -36,6 +37,13 @@ ObjContribScaledJac<FuncType>::ObjContribScaledJac(double weight) :
 }
 
 
+template<class FuncType>
+ObjContrib *ObjContribScaledJac<FuncType>::copy() const
+{
+  return new ObjContribScaledJac<FuncType>(*this);
+}
+
+
 template<class FuncType>
 void ObjContribScaledJac<FuncType>::initialize(Patch *mesh)
 {

contrib/HighOrderMeshOptimizer/OptHomRun.cpp

@@ -745,7 +745,6 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
 //#include "OptHomRun.h"
 #include "MeshOptCommon.h"
 #include "MeshOptObjContribFunc.h"
-#include "MeshOptObjectiveFunction.h"
 #include "MeshOptObjContrib.h"
 #include "MeshOptObjContribScaledNodeDispSq.h"
 #include "OptHomObjContribScaledJac.h"
@@ -808,27 +807,35 @@ void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p)
   par.patchDef = &patchDef;
   par.optDisplay = 30;
   par.verbose = 4;
+
   ObjContribScaledNodeDispSq<ObjContribFuncSimple> nodeDistFunc(p.weightFixed, p.weightFree);
   ObjContribScaledJac<ObjContribFuncBarrierMin> minJacBarFunc(1.);
+  minJacBarFunc.setTarget(p.BARRIER_MIN, 1.);
   ObjContribScaledJac<ObjContribFuncBarrierMinMax> minMaxJacBarFunc(1.);
+  minMaxJacBarFunc.setTarget(p.BARRIER_MAX, 1.);
+  par.allContrib.resize(3);
+  par.allContrib[0] = &nodeDistFunc;
+  par.allContrib[1] = &minJacBarFunc;
+  par.allContrib[2] = &minMaxJacBarFunc;
+
   MeshOptParameters::PassParameters minJacPass;
   minJacPass.barrierIterMax = p.optPassMax;
   minJacPass.optIterMax = p.itMax;
-  minJacPass.objFunc = ObjectiveFunction();
-  minJacPass.objFunc.push_back(&nodeDistFunc);
-  minJacBarFunc.setTarget(p.BARRIER_MIN, 1.);
-  minJacPass.objFunc.push_back(&minJacBarFunc);
+  minJacPass.contribInd.resize(2);
+  minJacPass.contribInd[0] = 0;
+  minJacPass.contribInd[1] = 1;
   par.pass.push_back(minJacPass);
+
   if (p.BARRIER_MAX > 0.) {
     MeshOptParameters::PassParameters minMaxJacPass;
     minMaxJacPass.barrierIterMax = p.optPassMax;
     minMaxJacPass.optIterMax = p.itMax;
-    minMaxJacPass.objFunc = ObjectiveFunction();
-    minMaxJacPass.objFunc.push_back(&nodeDistFunc);
-    minMaxJacBarFunc.setTarget(p.BARRIER_MAX, 1.);
-    minMaxJacPass.objFunc.push_back(&minMaxJacBarFunc);
+    minMaxJacPass.contribInd.resize(2);
+    minMaxJacPass.contribInd[0] = 0;
+    minMaxJacPass.contribInd[1] = 2;
     par.pass.push_back(minMaxJacPass);
   }
+
   MeshOptResults res;
 
   meshOptimizer(gm, par, res);

contrib/MeshOptimizer/MeshOpt.cpp

@@ -66,10 +66,18 @@ void printProgressFunc(const alglib::real_1d_array &x, double Obj, void *MOInst)
 
 MeshOpt::MeshOpt(const std::map<MElement*,GEntity*> &element2entity,
                  const std::set<MElement*> &els, std::set<MVertex*> &toFix,
-                 bool fixBndNodes) :
+                 bool fixBndNodes, const std::vector<ObjContrib*> &allContrib) :
   patch(element2entity, els, toFix, fixBndNodes), _verbose(0),
-  _iPass(0), _objFunc(0), _iter(0), _intervDisplay(0), _initObj(0)
+  _iPass(0), _objFunc(), _iter(0), _intervDisplay(0), _initObj(0)
 {
+  _allContrib.resize(allContrib.size());
+  for (int i = 0; i < allContrib.size(); i++) _allContrib[i] = allContrib[i]->copy();
+}
+
+
+MeshOpt::~MeshOpt()
+{
+  for (int i = 0; i < _allContrib.size(); i++) delete _allContrib[i];
 }
 
 
@@ -77,8 +85,8 @@ void MeshOpt::evalObjGrad(const alglib::real_1d_array &x, double &obj,
                                          alglib::real_1d_array &gradObj)
 {
   patch.updateMesh(x.getcontent());
-  _objFunc->compute(obj, gradObj);
-  if (_objFunc->targetReached()) {
+  _objFunc.compute(obj, gradObj);
+  if (_objFunc.targetReached()) {
     if (_verbose > 2) Msg::Info("Reached target values, setting null gradient");
     obj = 0.;
     for (int i = 0; i < gradObj.length(); i++) gradObj[i] = 0.;
@@ -92,7 +100,7 @@ void MeshOpt::printProgress(const alglib::real_1d_array &x, double Obj)
 
   if ((_verbose > 2) && (_iter % _intervDisplay == 0))
     Msg::Info(("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E)" +
-              _objFunc->minMaxStr()).c_str(), _iter, Obj, Obj/_initObj);
+              _objFunc.minMaxStr()).c_str(), _iter, Obj, Obj/_initObj);
 }
 
 
@@ -112,7 +120,7 @@ void MeshOpt::calcScale(alglib::real_1d_array &scale)
 
 void MeshOpt::updateResults(MeshOptResults &res)
 {
-  _objFunc->updateResults(res);
+  _objFunc.updateResults(res);
 }
 
 
@@ -180,15 +188,15 @@ int MeshOpt::optimize(MeshOptParameters &par)
   for (_iPass=0; _iPass<par.pass.size(); _iPass++) {
 
     // Set objective function Output if required
-    _objFunc = &par.pass[_iPass].objFunc;
+    _objFunc = ObjectiveFunction(_allContrib, par.pass[_iPass].contribInd);
     if (_verbose > 2) {
-      std::string msgStr = "* Pass %d with contributions: " + _objFunc->contribNames();
+      std::string msgStr = "* Pass %d with contributions: " + _objFunc.contribNames();
       Msg::Info(msgStr.c_str(), _iPass);
     }
 
     // Initialize contributions
-    _objFunc->initialize(&patch);
-    _objFunc->updateParameters();
+    _objFunc.initialize(&patch);
+    _objFunc.updateParameters();
 
     // Calculate initial objective function value and gradient
    _initObj = 0.;
@@ -198,14 +206,14 @@ int MeshOpt::optimize(MeshOptParameters &par)
     evalObjGrad(x, _initObj, gradObj);
 
     // Loop for update of objective function parameters (barrier movement)
-    bool targetReached = _objFunc->targetReached();
+    bool targetReached = _objFunc.targetReached();
     for (int iBar=0; (iBar<par.pass[_iPass].barrierIterMax) && (!targetReached); iBar++) {
       if (_verbose > 2) Msg::Info("--- Optimization run %d", iBar);
-      _objFunc->updateParameters();
+      _objFunc.updateParameters();
       runOptim(x, gradObj, par.pass[_iPass].optIterMax);
-      _objFunc->updateMinMax();
-      targetReached = _objFunc->targetReached();
-      if (_objFunc->stagnated()) {
+      _objFunc.updateMinMax();
+      targetReached = _objFunc.targetReached();
+      if (_objFunc.stagnated()) {
         if (_verbose > 2) Msg::Info("Stagnation detected, stopping optimization");
         break;
       }
@@ -213,7 +221,7 @@ int MeshOpt::optimize(MeshOptParameters &par)
 
     // Check results of pass and output if required
     if (!targetReached) result = 0;
-    std::string failMeasures = _objFunc->failMeasures();
+    std::string failMeasures = _objFunc.failMeasures();
     if (!failMeasures.empty()) {
       result = -1;
       if (_verbose > 2)
@@ -224,7 +232,7 @@ int MeshOpt::optimize(MeshOptParameters &par)
       if (targetReached)
         Msg::Info("Target reached for pass %d", _iPass);
       else {
-        std::string failedTargets = _objFunc->targetsNotReached();
+        std::string failedTargets = _objFunc.targetsNotReached();
         Msg::Warning("Failed to reach target in pass %d for "
                       "contributions %s", _iPass, failedTargets.c_str());
       }

contrib/MeshOptimizer/MeshOpt.h

@@ -48,7 +48,9 @@ class MeshOpt
 public:
   Patch patch;
   MeshOpt(const std::map<MElement*,GEntity*> &element2entity,
-          const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBndNodes);
+          const std::set<MElement*> &els, std::set<MVertex*> &toFix,
+          bool fixBndNodes, const std::vector<ObjContrib*> &allContrib);
+  ~MeshOpt();
   int optimize(MeshOptParameters &par);
   void updateMesh(const alglib::real_1d_array &x);
   void updateResults(MeshOptResults &res);
@@ -59,7 +61,8 @@ public:
  private:
   int _verbose;
   int _iPass;
-  ObjectiveFunction *_objFunc;                                                      // Current contributions to objective function
+  std::vector<ObjContrib*> _allContrib;                                             // All contributions to objective function for all passes
+  ObjectiveFunction _objFunc;                                                       // Contributions to objective function for current pass
   int _iter, _intervDisplay;                                                        // Current iteration, interval of iterations for reporting
   double _initObj;                                                                  // Values for reporting
   void calcScale(alglib::real_1d_array &scale);

contrib/MeshOptimizer/MeshOptCommon.h

@@ -31,10 +31,10 @@
 #define _MESHOPTCOMMON_H_
 
 #include <vector>
-#include "MeshOptObjectiveFunction.h"
 
 
 class MElement;
+class ObjContrib;
 
 
 struct MeshOptResults {                             // Output of mesh optimization
@@ -52,7 +52,7 @@ private:
 struct MeshOptParameters {                             // Parameters controlling the strategy
   enum { STRAT_CONNECTED, STRAT_ONEBYONE };
   struct PassParameters {                              // Parameters controlling the optimization procedure in each pass
-    ObjectiveFunction objFunc;                         // Contributions to objective function
+    std::vector<int> contribInd;                       // Indices of contributions to objective function
     int optIterMax;                                    // Max. number of opt. iterations each time the barrier is moved
     int barrierIterMax;                                // Max. number of times the barrier is moved
   };
@@ -74,6 +74,7 @@ struct MeshOptParameters {                             // Parameters controlling
   bool onlyVisible ;                                    // Apply optimization to visible entities ONLY
   bool fixBndNodes;                                     // If points can move on boundaries
   PatchDefParameters *patchDef;
+  std::vector<ObjContrib*> allContrib;                  // All contributions to objective functions for all passes
   std::vector<PassParameters> pass;
   int optDisplay;                                       // Sampling rate in opt. iterations for display
   int verbose;                                          // Level of information displayed and written to disk

contrib/MeshOptimizer/MeshOptObjContrib.h

@@ -16,6 +16,7 @@ class ObjContrib
 public:
   ObjContrib(std::string mesName, std::string name);
   virtual ~ObjContrib() {}
+  virtual ObjContrib *copy() const = 0;
   const double getMin() { return _min; }
   const double getMax() { return _max; }
   const std::string &getName() const { return _name; }

contrib/MeshOptimizer/MeshOptObjContribScaledNodeDispSq.h

@@ -13,6 +13,7 @@ class ObjContribScaledNodeDispSq : public ObjContrib, public FuncType
 public:
   ObjContribScaledNodeDispSq(double weightFixed, double weightFree);
   virtual ~ObjContribScaledNodeDispSq() {}
+  virtual ObjContrib *copy() const;
   virtual void initialize(Patch *mesh);
   virtual bool fail() { return false; }
   virtual bool addContrib(double &Obj, alglib::real_1d_array &gradObj);
@@ -38,6 +39,13 @@ ObjContribScaledNodeDispSq<FuncType>::ObjContribScaledNodeDispSq(double weightFi
 }
 
 
+template<class FuncType>
+ObjContrib *ObjContribScaledNodeDispSq<FuncType>::copy() const
+{
+  return new ObjContribScaledNodeDispSq<FuncType>(*this);
+}
+
+
 template<class FuncType>
 void ObjContribScaledNodeDispSq<FuncType>::initialize(Patch *mesh)
 {

contrib/MeshOptimizer/MeshOptObjectiveFunction.cpp

@@ -5,6 +5,14 @@
 #include "MeshOptObjectiveFunction.h"
 
 
+ObjectiveFunction::ObjectiveFunction(const std::vector<ObjContrib*> &allContrib,
+                                              const std::vector<int> &contribInd)
+{
+  resize(contribInd.size());
+  for (int i=0; i<contribInd.size(); i++) operator[](i) = allContrib[contribInd[i]];
+}
+
+
 void ObjectiveFunction::initialize(Patch *mesh)
 {
    for (std::vector<ObjContrib*>::iterator it=begin(); it!=end(); it++)

contrib/MeshOptimizer/MeshOptObjectiveFunction.h

@@ -12,9 +12,12 @@ class MeshOptResults;
 class Patch;
 
 
-class ObjectiveFunction : public  std::vector<ObjContrib*>                         // Contributions to objective function in each pass
+class ObjectiveFunction : protected std::vector<ObjContrib*>    // Contributions to objective function in each pass
 {
 public:
+  ObjectiveFunction() {}
+  ObjectiveFunction(const std::vector<ObjContrib*> &model,
+                        const std::vector<int> &contribInd);
   void initialize(Patch *mesh);
   std::string contribNames();
   std::string minMaxStr();

contrib/MeshOptimizer/MeshOptimizer.cpp

@@ -267,7 +267,7 @@ static void optimizeConnectedPatches
       Msg::Info("Optimizing patch %i/%i composed of %i elements", iPatch,
                       toOptimize.size()-1, toOptimize[iPatch].first.size());
     MeshOpt opt(element2entity, toOptimize[iPatch].first, toOptimize[iPatch].second,
-                                                                      par.fixBndNodes);
+                                                      par.fixBndNodes, par.allContrib);
     if (par.verbose > 3) {
       std::ostringstream ossI1;
       ossI1 << "initial_patch-" << iPatch << ".msh";
@@ -356,7 +356,7 @@ static void optimizeOneByOne
       if (par.verbose > 1)
         Msg::Info("Optimizing patch %i (max. %i remaining) composed of %4d elements",
                                             iBadEl, badElts.size(), toOptimize.size());
-      MeshOpt opt(element2entity, toOptimize, toFix, par.fixBndNodes);
+      MeshOpt opt(element2entity, toOptimize, toFix, par.fixBndNodes, par.allContrib);
       if (par.verbose > 3) {
         std::ostringstream ossI1;
         ossI1 << "initial_patch-" << iBadEl << ".msh";