Skip to content
Snippets Groups Projects
Select Git revision
  • 83a1a852b7405612e563efb02171b0ceb3bc2f5e
  • master default protected
  • overlaps_tags_and_distributed_export
  • overlaps_tags_and_distributed_export_rebased
  • relaying
  • alphashapes
  • patches-4.14
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

HierarchicalBasisH1Quad.cpp

Blame
  • OptHOM.cpp 12.54 KiB
    // Copyright (C) 2013 ULg-UCL
    //
    // Permission is hereby granted, free of charge, to any person
    // obtaining a copy of this software and associated documentation
    // files (the "Software"), to deal in the Software without
    // restriction, including without limitation the rights to use, copy,
    // modify, merge, publish, distribute, and/or sell copies of the
    // Software, and to permit persons to whom the Software is furnished
    // to do so, provided that the above copyright notice(s) and this
    // permission notice appear in all copies of the Software and that
    // both the above copyright notice(s) and this permission notice
    // appear in supporting documentation.
    //
    // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
    // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
    // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
    // NONINFRINGEMENT OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE
    // COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR
    // ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY
    // DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
    // WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
    // ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
    // OF THIS SOFTWARE.
    //
    // Please report all bugs and problems to the public mailing list
    // <gmsh@geuz.org>.
    //
    // Contributors: Thomas Toulorge, Jonathan Lambrechts
    
    #include <iostream>
    #include <algorithm>
    #include "OptHomMesh.h"
    #include "OptHOM.h"
    #include "OptHomRun.h"
    #include "GmshMessage.h"
    #include "GmshConfig.h"
    
    #if defined(HAVE_BFGS)
    
    #include "ap.h"
    #include "alglibinternal.h"
    #include "alglibmisc.h"
    #include "linalg.h"
    #include "optimization.h"
    
    static inline double compute_f(double v, double barrier)
    {
      if ((v > barrier && barrier < 1) || (v < barrier && barrier > 1)) {
        const double l = log((v - barrier) / (1 - barrier));
        const double m = (v - 1);
        return l * l + m * m;
      }
      else return 1.e300;
      // if (v < 1.) return pow(1.-v,powM);
      // if (v < 1.) return exp((long double)pow(1.-v,3));
      // else return pow(v-1.,powP);
    }
    
    static inline double compute_f1(double v, double barrier)
    {
      if ((v > barrier && barrier < 1) || (v < barrier && barrier > 1)) {
        return 2 * (v - 1) + 2 * log((v - barrier) / (1 - barrier)) / (v - barrier);
      }
      else return barrier < 1 ? -1.e300 : 1e300;
      // if (v < 1.) return -powM*pow(1.-v,powM-1.);
      // if (v < 1.) return -3.*pow(1.-v,2)*exp((long double)pow(1.-v,3));
      // else return powP*pow(v-1.,powP-1.);
    }
    
    OptHOM::OptHOM(const std::map<MElement*,GEntity*> &element2entity,
                   const std::set<MElement*> &els, std::set<MVertex*> & toFix,
                   bool fixBndNodes, bool fastJacEval) :
      mesh(element2entity, els, toFix, fixBndNodes, fastJacEval)
    {
      _optimizeMetricMin = false;
    }
    
    // Contribution of the element Jacobians to the objective function value and
    // gradients (2D version)
    bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
    {
      minJac = 1.e300;
      maxJac = -1.e300;
    
      for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
        // Scaled Jacobians
        std::vector<double> sJ(mesh.nBezEl(iEl));
        // Gradients of scaled Jacobians
        std::vector<double> gSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl));
        mesh.scaledJacAndGradients (iEl,sJ,gSJ);
    
        for (int l = 0; l < mesh.nBezEl(iEl); l++) {
          double f1 = compute_f1(sJ[l], jacBar);
          Obj += compute_f(sJ[l], jacBar);
          if (_optimizeBarrierMax) {
            Obj += compute_f(sJ[l], barrier_min);
            f1 += compute_f1(sJ[l], barrier_min);
          }
          for (int iPC = 0; iPC < mesh.nPCEl(iEl); iPC++)
            gradObj[mesh.indPCEl(iEl,iPC)] += f1*gSJ[mesh.indGSJ(iEl,l,iPC)];
          minJac = std::min(minJac, sJ[l]);
          maxJac = std::max(maxJac, sJ[l]);
        }
      }
    
      return true;
    }
    
    bool OptHOM::addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj)
    {
    
      minJac = 1.e300;
      maxJac = -1.e300;
    
      for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
        // Scaled Jacobians
        std::vector<double> sJ(mesh.nBezEl(iEl));
        // Gradients of scaled Jacobians
        std::vector<double> gSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl));
        mesh.metricMinAndGradients (iEl,sJ,gSJ);
    
        for (int l = 0; l < mesh.nBezEl(iEl); l++) {
          Obj += compute_f(sJ[l], jacBar);
          const double f1 = compute_f1(sJ[l], jacBar);
          for (int iPC = 0; iPC < mesh.nPCEl(iEl); iPC++)
            gradObj[mesh.indPCEl(iEl,iPC)] += f1*gSJ[mesh.indGSJ(iEl,l,iPC)];
          minJac = std::min(minJac, sJ[l]);
          maxJac = std::max(maxJac, sJ[l]);
        }
      }
      return true;
    }
    
    // Contribution of the vertex distance to the objective function value and
    // gradients (2D version)
    bool OptHOM::addDistObjGrad(double Fact, double Fact2, double &Obj,
                                alglib::real_1d_array &gradObj)
    {
      maxDist = 0;
      avgDist = 0;
      int nbBnd = 0;
    
      for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
        const double Factor = invLengthScaleSq*(mesh.forced(iFV) ? Fact : Fact2);
        const double dSq = mesh.distSq(iFV), dist = sqrt(dSq);
        Obj += Factor * dSq;
        std::vector<double> gDSq(mesh.nPCFV(iFV));
        mesh.gradDistSq(iFV,gDSq);
        for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++)
          gradObj[mesh.indPCFV(iFV,iPC)] += Factor*gDSq[iPC];
        maxDist = std::max(maxDist, dist);
        avgDist += dist;
        nbBnd++;
      }
      if (nbBnd != 0) avgDist /= nbBnd;
    
      return true;
    }
    
    void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj,
                             alglib::real_1d_array &gradObj)
    {
      mesh.updateMesh(x.getcontent());
    
      Obj = 0.;
      for (int i = 0; i < gradObj.length(); i++) gradObj[i] = 0.;
    
      addJacObjGrad(Obj, gradObj);
      addDistObjGrad(lambda, lambda2, Obj, gradObj);
      if(_optimizeMetricMin)
        addMetricMinObjGrad(Obj, gradObj);
      if ((minJac > barrier_min) && (maxJac < barrier_max || !_optimizeBarrierMax)) {
        Msg::Info("Reached %s (%g %g) requirements, setting null gradient",
                  _optimizeMetricMin ? "svd" : "jacobian", minJac, maxJac);
        Obj = 0.;
        for (int i = 0; i < gradObj.length(); i++) gradObj[i] = 0.;
      }
    }
    
    void evalObjGradFunc(const alglib::real_1d_array &x, double &Obj,
                         alglib::real_1d_array &gradObj, void *HOInst)
    {
      (static_cast<OptHOM*>(HOInst))->evalObjGrad(x, Obj, gradObj);
    }
    
    void OptHOM::recalcJacDist()
    {
      maxDist = 0;
      avgDist = 0;
      int nbBnd = 0;
    
      for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
        if (mesh.forced(iFV)) {
          double dSq = mesh.distSq(iFV);
          maxDist = std::max(maxDist, sqrt(dSq));
          avgDist += sqrt(dSq);
          nbBnd++;
        }
      }
      if (nbBnd != 0) avgDist /= nbBnd;
    
      minJac = 1.e300;
      maxJac = -1.e300;
      for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
        // Scaled Jacobians
        std::vector<double> sJ(mesh.nBezEl(iEl));
        // (Dummy) gradients of scaled Jacobians
        std::vector<double> dumGSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl));
        mesh.scaledJacAndGradients (iEl,sJ,dumGSJ);
        if(_optimizeMetricMin)
          mesh.metricMinAndGradients (iEl,sJ,dumGSJ);
        for (int l = 0; l < mesh.nBezEl(iEl); l++) {
          minJac = std::min(minJac, sJ[l]);
          maxJac = std::max(maxJac, sJ[l]);
        }
      }
    }
    
    void OptHOM::printProgress(const alglib::real_1d_array &x, double Obj)
    {
      iter++;
    
      if (iter % progressInterv == 0) {
        Msg::Info("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E) "
                  "-- minJ = %12.5E  maxJ = %12.5E Max D = %12.5E Avg D = %12.5E",
                  iter, Obj, Obj/initObj, minJac, maxJac, maxDist, avgDist);
      }
    }
    
    void printProgressFunc(const alglib::real_1d_array &x, double Obj, void *HOInst)
    {
      ((OptHOM*)HOInst)->printProgress(x,Obj);
    }
    
    void OptHOM::calcScale(alglib::real_1d_array &scale)
    {
      scale.setlength(mesh.nPC());
    
      // Calculate scale
      for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
        std::vector<double> scaleFV(mesh.nPCFV(iFV),1.);
        mesh.pcScale(iFV,scaleFV);
        for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++)
          scale[mesh.indPCFV(iFV,iPC)] = scaleFV[iPC];
      }
    }
    
    void OptHOM::OptimPass(alglib::real_1d_array &x,
                           const alglib::real_1d_array &initGradObj, int itMax)
    {
      static const double EPSG = 0.;
      static const double EPSF = 0.;
      static const double EPSX = 0.;
      static int OPTMETHOD = 1;
    
      Msg::Info("--- Optimization pass with initial jac. range (%g, %g), jacBar = %g",
                minJac, maxJac, jacBar);
    
      iter = 0;
    
      alglib::real_1d_array scale;
      calcScale(scale);
    
      int iterationscount = 0, nfev = 0, terminationtype = -1;
      if (OPTMETHOD == 1) {
        alglib::mincgstate state;
        alglib::mincgreport rep;
        mincgcreate(x, state);
        mincgsetscale(state,scale);
        mincgsetprecscale(state);
        mincgsetcond(state, EPSG, EPSF, EPSX, itMax);
        mincgsetxrep(state, true);
        alglib::mincgoptimize(state, evalObjGradFunc, printProgressFunc, this);
        mincgresults(state, x, rep);
        iterationscount = rep.iterationscount;
        nfev = rep.nfev;
        terminationtype = rep.terminationtype;
      }
      else {
        alglib::minlbfgsstate state;
        alglib::minlbfgsreport rep;
        minlbfgscreate(3, x, state);
        minlbfgssetscale(state,scale);
        minlbfgssetprecscale(state);
        minlbfgssetcond(state, EPSG, EPSF, EPSX, itMax);
        minlbfgssetxrep(state, true);
        alglib::minlbfgsoptimize(state, evalObjGradFunc, printProgressFunc, this);
        minlbfgsresults(state, x, rep);
        iterationscount = rep.iterationscount;
        nfev = rep.nfev;
        terminationtype = rep.terminationtype;
      }
    
      Msg::Info("Optimization finalized after %d iterations (%d function evaluations),",
                iterationscount, nfev);
      switch(int(terminationtype)) {
      case 1: Msg::Info("because relative function improvement is no more than EpsF"); break;
      case 2: Msg::Info("because relative step is no more than EpsX"); break;
      case 4: Msg::Info("because gradient norm is no more than EpsG"); break;
      case 5: Msg::Info("because the maximum number of steps was taken"); break;
      default: Msg::Info("with code %d", int(terminationtype)); break;
      }
    }
    
    int OptHOM::optimize(double weightFixed, double weightFree, double b_min,
                         double b_max, bool optimizeMetricMin, int pInt,
                         int itMax, int optPassMax)
    {
      barrier_min = b_min;
      barrier_max = b_max;
      progressInterv = pInt;
    //  powM = 4;
    //  powP = 3;
    
      _optimizeMetricMin = optimizeMetricMin;
      // Set weights & length scale for non-dimensionalization
      lambda = weightFixed;
      lambda2 = weightFree;
      std::vector<double> dSq(mesh.nEl());
      mesh.distSqToStraight(dSq);
      const double maxDSq = *max_element(dSq.begin(),dSq.end());
      if (maxDSq < 1.e-10) {                                        // Length scale for non-dim. distance
        std::vector<double> sSq(mesh.nEl());
        mesh.elSizeSq(sSq);
        const double maxSSq = *max_element(sSq.begin(),sSq.end());
        invLengthScaleSq = 1./maxSSq;
      }
      else invLengthScaleSq = 1./maxDSq;
    
      // Set initial guess
      alglib::real_1d_array x;
      x.setlength(mesh.nPC());
      mesh.getUvw(x.getcontent());
    
      // Calculate initial performance
      recalcJacDist();
      initMaxDist = maxDist;
      initAvgDist = avgDist;
    
      const double jacBarStart = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
      jacBar = jacBarStart;
    
      _optimizeBarrierMax = false;
      // Calculate initial objective function value and gradient
      initObj = 0.;
      alglib::real_1d_array gradObj;
      gradObj.setlength(mesh.nPC());
      for (int i = 0; i < mesh.nPC(); i++) gradObj[i] = 0.;
      evalObjGrad(x, initObj, gradObj);
    
      Msg::Info("Start optimizing %d elements (%d vertices, %d free vertices, %d variables) "
                "with min barrier %g and max. barrier %g", mesh.nEl(), mesh.nVert(),
                mesh.nFV(), mesh.nPC(), barrier_min, barrier_max);
    
      int ITER = 0;
      bool minJacOK = true;
      while (minJac < barrier_min) {
        const double startMinJac = minJac;
        OptimPass(x, gradObj, itMax);
        recalcJacDist();
        jacBar = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
        if (ITER ++ > optPassMax) {
          minJacOK = (minJac > barrier_min);
          break;
        }
        if (fabs((minJac-startMinJac)/startMinJac) < 0.01) {
          Msg::Info("Stagnation in minJac detected, stopping optimization");
          minJacOK = false;
          break;
        }
      }
    
      ITER = 0;
      if (minJacOK && (!_optimizeMetricMin)) {
        _optimizeBarrierMax = true;
        jacBar =  1.1 * maxJac;
        while (maxJac > barrier_max ) {
          const double startMaxJac = maxJac;
          OptimPass(x, gradObj, itMax);
          recalcJacDist();
          jacBar =  1.1 * maxJac;
          if (ITER ++ > optPassMax) break;
          if (fabs((maxJac-startMaxJac)/startMaxJac) < 0.01) {
            Msg::Info("Stagnation in maxJac detected, stopping optimization");
            break;
          }
        }
      }
    
      Msg::Info("Optimization done Range (%g,%g)", minJac, maxJac);
    
      if (minJac > barrier_min && maxJac < barrier_max) return 1;
      if (minJac > 0.0) return 0;
      return -1;
    }
    
    #endif