CirclePattern.cpp

Go to the documentation of this file.

#include "CirclePattern.h"
#include "util.h"
#include "Mesh.h"

using namespace std;

static double M2_PI = 2*M_PI;

double CirclePattern::ImLi2Sum (const double& theta, const double& x, 
                                    const double& tanHalfThetaStar, const double& clThetaStar)
{
   double tStar = M_PI - theta;
   double p = 2*atan(tanh(0.5*x)*tanHalfThetaStar);

   return  p*x + Cl2(p + tStar) + Cl2(-p + tStar) - clThetaStar;
}

double CirclePattern::evaluateEnergy (Mesh * mesh, double *xx) {
   int i, j;
   double energy = 0;
   for (i = 0; i < mesh->numFaces(); i++) {
      energy += xx[i];
   }

   energy = M2_PI*energy;   // Total contribution from faces

   Mesh::HalfEdgeIterator hIter = mesh->halfEdgeIterator();
   for (hIter.reset(); !hIter.end(); hIter++) {
      Edge * e = hIter.half_edge();

      if (!e->isBoundary()) {
          i = e->face->ID;
          j = e->pair->face->ID;
          if (i > j) {
             double theta = e->theta();
             energy += ImLi2Sum(theta, xx[i] - xx[j], e->tanHalfTStar(), e->Cl2TStar());
             energy -= (M_PI - theta)*(xx[i] + xx[j]);
          }
      }
      else {
          if (e->face) {
             i = e->face->ID;
             energy -= 2*(M_PI - e->theta())*(xx[i]);
          }
      }
   }
   return energy;
}

void CirclePattern::computeGradient (Mesh * mesh, double * xx, double * grad) {

   Mesh::FaceIterator fIter = mesh->faceIterator();

   /* Compute and store the gradient of the f. */
   int i, j;
   for (i = 0; i < mesh->numFaces(); i++)
      grad[i] = M2_PI;

   // sum over edges
   Mesh::HalfEdgeIterator hIter = mesh->halfEdgeIterator();
   for (hIter.reset(); !hIter.end(); hIter++) {
      Edge * e = hIter.half_edge();
      if (!e->isBoundary()) {
          i = e->face->ID;
          j = e->pair->face->ID;
          if (i > j) {
             grad[i] -= fTheta2(e->cosTheta(), e->sinTheta(), xx[j] - xx[i]);
             grad[j] -= fTheta2(e->cosTheta(), e->sinTheta(), xx[i] - xx[j]);
          }
      }
      else if (e->face) {
          i = e->face->ID;
          grad[i] -= 2*(M_PI - e->theta());
      }
   }

}

void CirclePattern::buildSparsityForHessian (Mesh * mesh, int * hessI, int * hessJ) {
   int i, k, j;

   for (k = 0; k < mesh->numFaces(); k++) {
      hessI[k] = k; 
      hessJ[k] = k;
   }

   k = mesh->numFaces();
   Mesh::HalfEdgeIterator hIter = mesh->halfEdgeIterator();
   for (hIter.reset(); !hIter.end(); hIter++) {
      if (!hIter.half_edge()->isBoundary()) {
          i = hIter.half_edge()->face->ID;
          j  = hIter.half_edge()->pair->face->ID;
          if (i > j) {
             hessI[k] = i; 
             hessJ[k] = j;
             k++;
          }       
      }
   }
}

void CirclePattern::computeHessian (Mesh * mesh, double * xx, double * hessval) {
   for (int k = 0; k < mesh->numFaces(); k++) 
      hessval[k] = 0;

   int count = mesh->numFaces();
   int i,j;
   Mesh::HalfEdgeIterator hIter = mesh->halfEdgeIterator();
   for (hIter.reset(); !hIter.end(); hIter++) {
      Edge * e = hIter.half_edge();
      if (!e->isBoundary()) {
          i = e->face->ID;
          j = e->pair->face->ID;
          if (i > j) {
             double ximj = xx[i] - xx[j];
             double fThPrime2 = fThetaPrime2(e->cosTheta(), e->sinTheta(), ximj);
             hessval[count] = -fThPrime2;
             hessval[i] += fThPrime2;
             hessval[j] += fThPrime2;
             count++;
          }
      }
   }
}

---