AnglesOptimization.cpp

Go to the documentation of this file.

#include "AnglesOptimization.h"
#include <iostream>
#include "util.h"
#include "Mesh.h"
#include <fstream>
#include <iomanip>

using namespace std;

const double AnglesOptimization::epsRange = 1.0e-3;

AnglesOptimization::AnglesOptimization(Mesh *_mesh, const MSKenv_t & env) {
   mesh = _mesh;

   NUMVAR = mesh->numFaces() * 3;
   NUMCON = mesh->numFaces() + mesh->numVertices() + mesh->numEdges() - mesh->numBoundary();
   NUMANZ = 3*NUMVAR - mesh->numBoundary();
   NUMQNZ = NUMVAR;

   int r = MSK_RES_OK;
   r = MSK_maketask(env, NUMCON, NUMVAR, &task);
   check_error (r != MSK_RES_OK, "Cannot make MOSEK task for angle optimization.");
   r = MSK_linkfiletotaskstream(task, MSK_STREAM_LOG, "_AngleOptMosekOut.txt", 0);
   check_error (r != MSK_RES_OK, "Cannot link MOSEK \"_AngleOptMosekOut.txt\" output file.");
}

void AnglesOptimization::doSolve() {
   allocateMosekArrays();
   setDataStructres();
   optimize();
}

void AnglesOptimization::setDataStructres() {

   int k = 0;
   int kBefore = 0;

   Mesh::FaceIterator fIter = mesh->faceIterator();
   for ( ; !fIter.end(); fIter++) {
      k = fIter.face()->ID;
      bkc[k] = MSK_BK_FX;
      blc[k] = M_PI;
      buc[k] = M_PI;
   }

   kBefore = mesh->numFaces();

   Mesh::VertexIterator vIter = mesh->vertexIterator();
   for (; !vIter.end(); vIter++) {
      Vertex * v = vIter.vertex();
      k = kBefore + v->ID;

      if (!v->constrained) {
          if (!v->isBoundary()) {
             bkc[k] = MSK_BK_FX;
             blc[k] = 2*M_PI;
             buc[k] = 2*M_PI;
          }
          else {
             bkc[k] = MSK_BK_RA;
             blc[k] = 0;
             buc[k] = 2*M_PI;  
          }
      }
      else {
          if (v->max_cone_angle == v->min_cone_angle) 
             bkc[k] = MSK_BK_FX;
          else 
             bkc[k] = MSK_BK_RA;
          blc[k] = v->min_cone_angle*M_PI;
          buc[k] = v->max_cone_angle*M_PI;
      }
   }

   kBefore += mesh->numVertices();
   k = kBefore;

   Mesh::EdgeIterator eIter = mesh->edgeIterator();
   for ( ; !eIter.end(); eIter++) {
      Edge * e = eIter.edge();
      if (!e->isBoundary()) {
          e->ID = k;
          e->pair->ID = e->ID;

          bkc[k] = MSK_BK_RA;
          blc[k] = epsRange; 
          buc[k] = M_PI - epsRange; 
          k++;
      }
   }

   /* Because we solve for difference between initial angles and valid angles ,   */
   /* we need to subtruct initial angles from sum conditions.                     */
   Mesh::HalfEdgeIterator hIter = mesh->halfEdgeIterator();
   for (; !hIter.end(); hIter++) {
      Edge * e = hIter.half_edge();
      if (e->face) {
          int vId = e->oppositeVertex()->ID;
          int fId = e->face->ID;
          int eId = e->ID;

          vId += mesh->numFaces();

          blc[fId] -= e->alphaOposite;
          buc[fId] -= e->alphaOposite;
          blc[vId] -= e->alphaOposite;
          buc[vId] -= e->alphaOposite;

          if (!e->isBoundary()) {
             blc[eId] -= e->alphaOposite;
             buc[eId] -= e->alphaOposite;
          }
      }
   }


   int currPtr = 0;
   int count = 0;

   for (hIter.reset(); !hIter.end(); hIter++) {
      Edge * e = hIter.half_edge();
      if (e->face) {
          Edge * e = hIter.half_edge();

          /* Linear minim part*/
          c[count]    = 0.0;

          /* Constraint matrix (our A). */  // assigning over the columns
          ptrb[count] = currPtr;

          if (e->isBoundary())
             ptre[count] = currPtr + 2;
          else
             ptre[count] = currPtr + 3;

          int fId = e->face->ID;
          sub[currPtr]  = fId;
          val[currPtr]  = 1.0;  
          currPtr++;

          int vId = e->oppositeVertex()->ID;
          vId += mesh->numFaces();
          sub[currPtr]  = vId;
          val[currPtr]  = 1.0;  
          currPtr++;

          if (!e->isBoundary()) {
             int eId = e->ID;
             sub[currPtr]  = eId;
             val[currPtr]  = 1.0;  
             currPtr++;
          }

          /* Bounds. */
          bkx[count]  = MSK_BK_RA;
          blx[count]  = epsRange - e->alphaOposite;
          bux[count]  = (M_PI - epsRange) - e->alphaOposite;
          count++;
      }
   }
}

void AnglesOptimization::optimize() {
   int r = MSK_RES_OK;
   r = MSK_inputdata(task, NUMCON,NUMVAR, NUMCON,NUMVAR, c, 
      0.0, ptrb, ptre, sub, val, bkc, blc, buc, bkx, blx, bux);    

   check_error ( r != MSK_RES_OK, 
      "Could not input constraints for angle optimization.\n Check  \"AngleOptMosekOut.txt\" for output." );

   int count = 0;
   for (count = 0; count < NUMVAR; count++) {
      qsubi[count] = count;   
      qsubj[count] = count;  
      qval[count] = 1;
   }

   r = MSK_putqobj(task, NUMQNZ, qsubi, qsubj, qval);
   check_error ( r != MSK_RES_OK, 
      "Could not input objective for angle optimization.\n Check  \"AngleOptMosekOut.txt\" for output." );


   r = MSK_optimize(task);
   check_error ( r != MSK_RES_OK, 
      "Could not solve angle optimization. Possibly the problem is infisible.\n Check  \"AngleOptMosekOut.txt\" for output." );

   r = MSK_solutionsummary(task, MSK_STREAM_ERR);
   check_error ( r != MSK_RES_OK, 
      "Could not output solution summary.\n Check  \"AngleOptMosekOut.txt\" for output." );

   //MSK_writedata(task,"problem.mbt");
   int ps, ss;
   r = MSK_getsolution (task, MSK_SOL_ITR, &ps, &ss, NULL, NULL, NULL, NULL, xx, NULL, NULL, NULL, NULL, NULL, NULL);

   check_error ( r != MSK_RES_OK, "Could not get solution of angle optimization.\n Check  \"AngleOptMosekOut.txt\" for output." );
   check_error ( ps != MSK_PRO_STA_PRIM_AND_DUAL_FEAS, 
      "Either problem is infisible or could not find optimal solution for some other reason.\n Check  \"AngleOptMosekOut.txt\" for output." );

   Mesh::HalfEdgeIterator hIter = mesh->halfEdgeIterator();
   count = 0;
   for (hIter.reset(); !hIter.end(); hIter++) {
      Edge * e = hIter.half_edge();
      if (e->face) 
          e->alphaOposite = e->alphaOposite + xx[count++];
      else
          e->alphaOposite = 0;
   }


   Mesh::EdgeIterator eIter = mesh->edgeIterator();
   for (eIter.reset(); !eIter.end(); eIter++){ 
      Edge * e = eIter.edge();
      e->setTheta (M_PI - e->alphaOposite - e->pair->alphaOposite);
      e->pair->setTheta (e->theta());
   }

   Mesh::VertexIterator vIter = mesh->vertexIterator();
   for (vIter.reset(); !vIter.end(); vIter++) {
      vIter.vertex()->cone_angle = 0;
   }

   for (eIter.reset(); !eIter.end(); eIter++) {
      Edge * e = eIter.edge();
      e->vertex->cone_angle += e->theta();
      e->pair->vertex->cone_angle += e->theta();
   }

   for (vIter.reset(); !vIter.end(); vIter++) {
      vIter.vertex()->cone_angle /= M_PI;
   }
}


void AnglesOptimization::allocateMosekArrays() {
   bkc = new int[NUMCON];  bkx = new int[NUMVAR];
   ptrb = new int[NUMVAR]; ptre = new int[NUMVAR];
   sub = new int[NUMANZ];
   qsubi = new int[NUMQNZ]; qsubj = new int[NUMQNZ];

   blc = new double[NUMCON];  buc = new double[NUMCON];
   c = new double[NUMVAR];    blx = new double[NUMVAR];
   bux = new double[NUMVAR];  val = new double[NUMANZ];
   xx = new double[NUMVAR];   qval = new double[NUMQNZ];
}

void AnglesOptimization::clearMosekArrays() {

   if (task) {
      delete [] bkc;  delete [] bkx;  delete [] ptrb; 
      delete [] ptre;  delete [] sub;  delete [] qsubi; delete [] qsubj;
      delete [] buc; delete [] c;   delete [] blx; delete [] blc;
      delete [] bux;  delete [] val; delete [] xx;  delete [] qval;
      MSK_deletetask(&task);
      task = NULL;
   }
}

---