$search

mesh.cc

Go to the documentation of this file.

/*
  James R. Diebel
  Stanford University

  Started: 24 August 2004
  Last revised: 27 Sugust 2004

  mesh.cc - class implementation of Mesh class (mesh.hh)

  Depends on:
  - Mesh class (mesh.hh)
*/

#include "bmtk/triangulate.h"
#include "bmtk/mesh.hh"
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <list>
#include <vector>
#include <algorithm>

namespace bmtk {


// Constructor
Mesh::Mesh() {
  po = true;
  ut = true;
  ib = false;

  gs2 = 10;
  gsn2 = 0.001;

  segThresh = 0.02;

  //psiSpline = new USpline(50,0.0,2.0,2.0,2*(1-4/gsn2)*exp(-4/gsn2));
  psiSpline = new USpline(500,0.0,2.0,2.0,2.0);
  setNormalVar(gsn2);

  /*ofstream fout1("psi.dat");
  ofstream fout2("psi0.dat");
  ofstream fout3("dpsi.dat");
  ofstream fout4("dpsi0.dat");
  for (float y=0;y<2;y+=0.001) {
    fout1 << y << " " << psiSpline->evalf(y) << endl;
    fout2 << y << " " << 2*y*exp(-2*y/gsn2) << endl;
    fout3 << y << " " << psiSpline->evaldf(y) << endl;
    fout4 << y << " " << 2*(1-2*y/gsn2)*exp(-2*y/gsn2) << endl;
  }
  fout1.close();
  fout2.close();
  fout3.close();
  fout4.close();*/

  nv = ne = nf = nn = 0;
  nd = 0;
  psi = 0;
  normGradPsi = 0;
  l = 1;
  iter = 1;
  nStepsForReset = 100;
  wroteIter = 0;
  maxIter = 100;

  rho[0] = rho[1] = 0;
  ut = (po || ut);
  tCG = new Timer("","\n");
  tLS = new Timer("","\n");
  tMP = new Timer("","\n");
  tVG = new Timer("","\n");
  tFN = new Timer("","\n");
  callBackExists = 0;

  foutConv.open("convergence.dat");
}

// Destructor
Mesh::~Mesh() {
  cleanUp();
  delete psiSpline;
  foutConv.close();
  delete tCG;
  delete tLS;
  delete tMP;
  delete tVG;
  delete tFN;
  if (po) cout << "Done." << endl << flush;
}

void Mesh::cleanUp() {
  if (po) cout << "Deleting mesh data..." << flush;
  if (nv) delete [] v; if (po) cout << " vertices," << flush;
  if (nv) delete [] qv; if (po) cout << "vertex queue," << flush;
  if (ne) delete [] e; if (po) cout << "edges," << flush;
  if (nf) delete [] f; if (po) cout << "faces," << flush;
  if (nf) delete [] r; if (po) cout << "regions," << flush;
  if (nf) delete [] q; if (po) cout << "face queue, " << flush;
}

void Mesh::registerCB(void (*callBack_)(void)) {
  callBackExists = 1;
  callBack = callBack_;
}

// update calls callback function to refresh GUI display
void Mesh::update() {
  if (callBackExists)
    (*callBack)();
}

void Mesh::setVertVar(float gs2_) {
  gs2 = gs2_;
}

void Mesh::setNormalVar(float gsn2_) {
  gsn2 = gsn2_;
  for (int i=0;i<psiSpline->n;i++) {
    float p = 2.0*float(i)/float(psiSpline->n - 1);
    psiSpline->y[i] = sqrt(p+gsn2);
  }
  psiSpline->b0 = 0.5/sqrt(gsn2);
  psiSpline->bn = 0.5/sqrt(2.0f+gsn2);
  /*float a,b,c,d;
  c = 4; d = 0.5;
  b = (d*d-1)*c/(2.0f*(d-d*d));
  a = (-c-2*b)/3.0f;

  for (int i=0;i<psiSpline->n;i++) {
    float p = 2.0*float(i)/float(psiSpline->n - 1);
    psiSpline->y[i] = a*p*p*p + b*p*p + c*p;
  }
  psiSpline->b0 = c;
  psiSpline->bn = 12.0f*a + 4.0f*b + c;*/

  psiSpline->update();
}

// General Mesh Functions //
void Mesh::findBBox() {
  if (po) cout << "- Finding bounding box..." << flush;
  boxMin = boxMax = v[0].x;
  for (int i=1;i<nv;i++) {
    for (int j=0;j<3;j++) {
      if (boxMin[j] < v[i].x[j]) boxMin[j] = v[i].x[j];
      if (boxMax[j] > v[i].x[j]) boxMax[j] = v[i].x[j];
    }
  }
  box = ~(boxMax - boxMin);
  if (po) cout << "Done." << endl << flush;
}

void Mesh::normalize() {
  if (po) cout << "- Scaling mesh for unity median edge length..." << endl
         << flush;
  for (int i=0;i<nv;i++) {
    v[i].x /= l;
    v[i].x0 /= l;
    if (nd) for (int k=0;k<nd;k++) v[i].xs[k] /= l;
    v[i].l /= l;
  }
  for (int i=0;i<ne;i++) e[i].l /= l;
  for (int i=0;i<nf;i++) {
    f[i].d /= l;
    f[i].x /= l;
  }
  cout << "  ";
  findBBox();
  cout << "  ";
  findFaceAreas();
  if (po) cout << "  Done." << endl << flush;
}

void Mesh::findLength() {
  if (po) cout << "- Computing the median edge length..." << flush;
  vector<float> lengths;
  for (int i=0;i<ne;i++) lengths.push_back(e[i].l);
  nth_element(lengths.begin(),
        lengths.begin()+lengths.size()/2,
        lengths.end());
  l = lengths[lengths.size()/2];
  if (po) cout << "Done." << endl << flush;
}

Mesh* Mesh::pruneFaces(float factor) {
  if (po) cout << "- Flagging triangles with edges longer than "
         << factor << " times the median edge length" << flush << endl;

  resetFaceFlags();
  findLength();

  // flag faces with long edges, keeping count
  int count = 0;
  for (int i=0;i<nf;i++) {
    if (f[i].e[0]->l > l*factor ||
  f[i].e[1]->l > l*factor ||
  f[i].e[2]->l > l*factor)
      f[i].flag = 1;
    else count++;
  }

  if (po) cout << "  - Keeping " << count << " of " << nf << " faces"
         << flush << endl;

  // fill new vertex and face lists
  float* vertices = new float[3*nv];
  int* faces = new int[3*count];
  for (int i=0;i<nv;i++) for (int j=0;j<3;j++) vertices[3*i+j] = v[i].x[j];
  for (int i=0,k=0;i<nf;i++) if (f[i].flag != 1) {
    for (int j=0;j<3;j++) faces[3*k+j] = f[i].v[j]->i;
    k++;
  }

  // and build the new Mesh
  Mesh* nm = new Mesh();
  nm->buildFrom(vertices,nv,faces,count);

  delete [] vertices;
  delete [] faces;

  if (po) cout << "  Done." << flush << endl;
  return nm;
}

void Mesh::findMeshProps() {
  findBBox();
  findFaceProps();
  findEdgeProps();
  findVertProps();
  findLength();
}

void Mesh::reset() {
  if (po) cout << "- Reverting to original mesh..." << flush;
  bool temp = po;
  po = false;
  for (int i=0;i<nv;i++) v[i].x = v[i].x0;
  findFaceProps();
  findEdgeProps();
  findVertProps();
  findMeshPotential();
  findLocalEdgePotentials();
  findVertGradients();
  findSearchDirs(true);
  rho[0] = rho[1] = 0;
  po = temp;
  if (po) cout << "Done." << endl << flush;
}

void Mesh::findMeshPotential() {
  if (ut) tMP->start();
  if (po) cout << "- Computing mesh potential..." << endl << flush;
  if (po) cout << "  ";
  findEdgePotentials();
  if (po) cout << "  ";
  findVertPotentials();
  psi = psiVert = psiEdge = 0;
  for (int i=0;i<ne;i++) if (ib || e[i].nf > 1) psiEdge += e[i].psi;
  for (int i=0;i<nv;i++) if (ib || !v[i].bound) psiVert += v[i].psi;
  psi = psiVert + psiEdge;
  if (po) cout << "  Psi = " << psiVert << " (v) + " << psiEdge
         << " (e) = " << psi << " ...Done." << flush;
  if (po) tMP->printMark();
  if (ut) tMP->mark();
}

void Mesh::findPotentials() {
  findEdgePotentials();
  findVertPotentials();
}

void Mesh::doBisectionLineSearch() {
  if (ut) tLS->start();
  if (po) cout << "  - Performing Bisection line search..." << endl << flush;
  float psi0 = psi;
  float step = 0.001f;
  float total = 0;
  bool temp = po;
  po = false;
  while (fabs(step) > 1e-5) {
    if (psi>psi0) {
      step *= -0.5;
    }
    psi0 = psi;
    total += step;
    moveVerts(total);
    findFaceNormals();
    findMeshPotential();
    if (ut) {
      tCG->stop();
      tLS->stop();
    }
    update();
    if (ut) {
      tCG->start();
      tLS->start();
    }
  }
  po = temp;
  if (po) cout << "    Total Movement = " << total << endl << flush;
  if (po) cout << "    Done." << flush;
  if (po) tLS->printMark();
  if (ut) tLS->mark();
}

void Mesh::doNewtonLineSearch() {
  if (ut) tLS->start();
  if (po) cout << "  - Performing Newton line search..." << endl << flush;
  float normGradPsi0 = normGradPsi;
  //float step = nv*0.001f/normDir;
  float step = 0.001f;
  float istep = 1/step;
  float total = 0;
  bool temp = po;
  po = false;
  do {
    moveVerts(total+step);
    findFaceNormals();
    findEdgePotentials();
    findVertGradients();
    findNormGradPsi();
    float diff = istep*normGradPsi - istep*normGradPsi0;

    total -= normGradPsi0/diff;
    moveVerts(total);

    findFaceNormals();
    findEdgePotentials();
    findVertGradients();
    findNormGradPsi();
    normGradPsi0 = normGradPsi;

    if (ut) {
      tCG->stop();
      tLS->stop();
    }
    update();
    if (ut) {
      tCG->start();
      tLS->start();
    }
  } while (fabs(normGradPsi - normGradPsi0)/normGradPsi > 1e-5);

  po = temp;
  if (po) cout << "    Total Movement = " << total << endl << flush;
  if (po) cout << "    Done." << flush;
  if (po) tLS->printMark();
  if (ut) tLS->mark();
}

void Mesh::runCG() {
  float psi0;
  int iter0 = iter;
  writeConvDetails();
  do {
    psi0 = psi;
    if (ut) tCG->start();
    if (po) cout << endl << "- Performing CG iteration #" << iter;
    if (po) if ((iter-iter0)%nStepsForReset == 0)
      cout << " with restart...";
    if (po) cout << endl << flush;
    if (po) cout << "  ";
    findVertGradients();
    if (po) cout << "  ";
    if ((iter-iter0)%nStepsForReset == 0) findSearchDirs(true);
    else findSearchDirs();
    if (po) cout << "  ";
    findNormGradPsi();
    if (po) cout << "  ";
    saveRefVerts();
    doNewtonLineSearch();
    findMeshPotential();
    findLocalEdgePotentials();
    iter++;
    writeConvDetails();
    if (po) cout << "  Done. " << flush;
    if (po) tCG->printMark();
    if (ut) tCG->mark();
    cout << "Relative change: " << fabs(psi-psi0)/psi << endl << flush;
  } while ((fabs(psi-psi0)/psi > 1e-3 && psi < psi0) || iter < 20);

  if (po) {
    cout << endl
   << "Summary of time spent so far:" << endl
   << "-----------------------------" << endl;
    tCG->printStore("Cojugate Gradient Search");
    tLS->printStore("             Line Search");
    tMP->printStore("Computing Mesh Potential");
    tVG->printStore("Computing Vert Gradients");
    tFN->printStore("   Comuting Face Normals");
  }
}

void Mesh::anneal() {
  if (po) cout << endl << "- Performing Annealing iteration #" << iter
         << "..." << endl <<  flush;
  writeConvDetails();
  // perturb vertex field
  //blurFaceNormals(2);
  //findVertNormals();
  vertFaceConsistency(1);

  // fix everything up for next iteration
  findFaceNormals();
  findMeshPotential();
  findLocalEdgePotentials();
  findVertGradients();
  findSearchDirs(true);
  iter++;
  writeConvDetails();
  update();
}

void Mesh::writeConvDetails() {
  if (wroteIter < iter) {
    foutConv << wroteIter << " "
       << psi << " "
       << psiVert << " "
       << psiEdge << " "
       << endl;
    wroteIter++;
  }
}

// Equality tests
bool Mesh::operator == (const Mesh &m) const {
  if (nv != m.nv || nf != m.nf || ne != m.ne) return false;
  for (int i=0;i<nv;i++) {
    if ((v[i].x!=m.v[i].x) || (v[i].nf!=m.v[i].nf) || (v[i].ne!=m.v[i].ne))
      return false;
    for (int j=0;j<v[i].nf;j++)
      if (v[i].f[j]->i != m.v[i].f[j]->i) return false;
    for (int j=0;j<v[i].ne;j++)
      if ((v[i].e[j]->i!=m.v[i].e[j]->i) || (v[i].v[j]->i!=m.v[i].v[j]->i))
  return false;
  }
  for (int i=0;i<nf;i++) {
    for (int j=0;j<3;j++)
      if ((f[i].v[j]->i!=m.f[i].v[j]->i) || (f[i].e[j]->i!=m.f[i].e[j]->i))
  return false;
    for (int j=0;j<f[i].nf;j++)
      if (f[i].f[j]->i != m.f[i].f[j]->i) return false;
  }
  for (int i=0;i<ne;i++) {
    for (int j=0;j<2;j++) if (e[i].v[j]->i != m.e[i].v[j]->i) return false;
    for (int j=0;j<e[i].nf;j++)
      if (e[i].f[j]->i != m.e[i].f[j]->i) return false;
  }
  return true;
}

bool Mesh::operator != (const Mesh &m) const {
  return (!(m==*this));
}

} // namespace bmtk

---

bmtk
Author(s): Benjamin Pitzer
autogenerated on Mon Mar 4 11:05:51 2013