meshconstructor.C

Go to the documentation of this file.

#ifdef HAVE_CONFIG_H
#include <tumble-conf.h>
#endif /* HAVE_CONFIG_H */

#include "meshconstructor.h"

#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <sys/time.h>

#include "boundarymesh.h"
#include "beziermesh.h"
#include "conformalmesher.h"
#include "globals.h"
#include "cleaner.h"
#include "spline.h"

using namespace std;
using namespace hashers;

MeshConstructor::MeshConstructor(BezierMesh* _bezier_mesh, BoundaryMesh* _bdry_mesh)
{
  bezier = _bezier_mesh;
  bdry = _bdry_mesh;
  mesher = new ConformalMesher(bezier, bdry);
}

MeshConstructor::~MeshConstructor()
{
  delete mesher;
}

BoundaryVertex* MeshConstructor::add_vertex(const Point2D &p)
{
    return bdry->add_boundary_vertex(p);
}

BoundaryVertex* MeshConstructor::add_vertex()
{
    Point2D p;
    return bdry->add_boundary_vertex(p);
}

BoundaryFace* MeshConstructor::add_face(const vector<BoundaryEdge*> &edges,
                                        const Point2D &point)
{
  /* refine angle is set to zero by default */
    BoundaryFace *face;
    face=bdry->add_boundary_face( edges, 0.0, 0);
    mesher->add_face_point(face, point);
    return face;
}

BoundaryEdge* MeshConstructor::create_edge(BoundaryVertex *v0, BoundaryVertex *v1,
                                           const vector<Point2D> &gps, const Point2D &initial)
{
    if(!v0 || !v1 || gps.size() < 2) return NULL;

    vector<Point2D> d;
    vector<double> k;
    unsigned j;

    /* compute deboor points */
    d.push_back( gps.front() );
    d.push_back( initial );
    for(j=1; j < gps.size() - 1; j++ ) d.push_back( gps[j]*2.0 - d[j]);
    d.push_back( gps.back() );

    /* compute knot vector */
    double temp = 0.0;
    double inc = 1.0/(gps.size()-1);
    for(j=0; j < gps.size(); j++){
        k.push_back( temp );
        temp += inc;
    }

    /* create new edge */
    BoundaryEdge *e = bdry->add_boundary_edge(v0, v1, d, k);
    if( !e ) cout<<"MeshConstructor::create_edge -- Error creating edge!"<<endl;
    return e;
}


BoundaryEdge*
MeshConstructor::add_line(BoundaryVertex* v0, BoundaryVertex* v1, int segments)
{
    return add_line(v0, v1, v0->get_cp(), v1->get_cp(), segments);
}

BoundaryEdge*
MeshConstructor::add_line(BoundaryVertex* v0, BoundaryVertex* v1, const Point2D &p1, int segments)
{
    return add_line(v0, v1, v0->get_cp(), p1, segments);
}

BoundaryEdge*
MeshConstructor::add_line(BoundaryVertex* v0, BoundaryVertex* v1, const Point2D &p0, const Point2D &p1, int segments)
{
    vector<Point2D>  gps;
    Point2D inc = (p1 - p0)/segments;
    unsigned ngps = segments + 1;

    v0->set_cp(p0);
    v1->set_cp(p1);

    Point2D temp = p0;
    for(unsigned i=0; i<ngps; i++){
        gps.push_back( temp );
        temp += inc;
    }

    Point2D initial = (gps[0] + gps[1]) * 0.5;

    return create_edge(v0, v1, gps, initial);
}

BoundaryEdge* MeshConstructor::add_circle(BoundaryVertex *v0, const Point2D &center, double radius, Radians start, int segments)
{
    if(!v0 || radius<=0.0 || segments < 4) return NULL;
    vector<Point2D>  gps;
    Radians inc = TWOPI / (double) segments;
    unsigned ngps = segments + 1;


    double theta = start;
    Point2D diff;
    for(unsigned i=0; i<ngps; i++){
        diff.assign(cos(theta), sin(theta));
        gps.push_back( center + diff*radius );
        theta += inc;
    }

    /* override v0 control point */
    v0->set_cp(gps[0]);

    /* compute initial deboor */
    /* where tangents to circle at gps[0] and gps[1] intersect */
    Point2D initial;
    double r;
    theta = start + inc/2.0;
    initial.assign(cos(theta), sin(theta));
    r = radius * sqrt(1.0 + sin(inc/2.0)*sin(inc/2.0) );
    initial = center + initial * r;

    return create_edge(v0, v0, gps, initial);
}

BoundaryEdge* MeshConstructor::add_arc(BoundaryVertex *v0,
                                       BoundaryVertex *v1,
                                       const Point2D &center,
                                       double radius,
                                       Radians start,
                                       Radians stop,
                                       int segments)
{
    Radians span = fabs(stop - start);
    if(!v0 || !v1 || radius<=0.0 || segments < (int)(2.0*span/PI) ) return NULL;
    vector<Point2D>  gps;
    Radians inc = span / (double) segments;
    unsigned ngps = segments + 1;


    double theta = start;
    Point2D diff;
    for(unsigned i=0; i<ngps; i++){
        diff.assign(cos(theta), sin(theta));
        gps.push_back( center + diff*radius );
        theta += inc;
    }

    /* override v0 and v1 control point */
    v0->set_cp(gps[0]);
    v1->set_cp(gps[segments]);

    /* compute initial deboor */
    /* where tangents to circle at gps[0] and gps[1] intersect */
    Point2D initial;
    double r;
    theta = start + inc/2.0;
    initial.assign(cos(theta), sin(theta));
    r = radius * sqrt(1.0 + sin(inc/2.0)*sin(inc/2.0) );
    initial = center + initial * r;

    return create_edge(v0, v1, gps, initial);
}

BoundaryEdge* MeshConstructor::add_spiral( BoundaryVertex *v0,
                                           BoundaryVertex *v1,
                                           const Point2D &center,
                                           double divergence,
                                           double radius,
                                           Radians start,
                                           Radians stop,
                                           int segments)
{
    Radians span = fabs(stop - start);
    if(!v0 || !v1 || stop<start || divergence<=0.0 ||
       radius<=0.0 || segments < (int)(2.0*span/PI) ) return NULL;

    vector<Point2D>  gps;
    Radians inc = span / (double) segments;
    unsigned ngps = segments + 1;


    double theta = start;
    double r = radius;
    Point2D diff;
    for(unsigned i=0; i<ngps; i++){
        diff.assign(cos(theta), sin(theta));
        gps.push_back( center + diff * r );
        theta += inc;
        r += divergence * (inc/TWOPI);
    }

    /* override v0 and v1 control point */
    v0->set_cp(gps[0]);
    v1->set_cp(gps[segments]);

    /* compute initial deboor */
    /* use parametric representation of sprial where
     * x(phi) = (radius + divergence*phi) cos(phi)
     * y(phi) = (radius + divergence*phi) sin(phi)
     * m0 = dx/dy(phi0)
     * m1 = dx/dy(phi1)
     */
    Point2D initial;
    double phi0 = start;
    double phi1 = start+inc;

    double m0 = spiral_derivative(phi0, radius, divergence, start);
    double x0 = (radius + divergence * (phi0 - start) / TWOPI)*cos(phi0);
    double y0 = (radius + divergence * (phi0 - start) / TWOPI)*sin(phi0);
    //cout<<"phi0: "<<phi0<<" x0: "<<x0<<" y0: "<<y0<<" m0: "<<m0<<endl;
    double m1 = spiral_derivative(phi1, radius, divergence, start);
    double x1 = (radius + divergence * (phi1 - start) / TWOPI)*cos(phi1);
    double y1 = (radius + divergence * (phi1 - start) / TWOPI)*sin(phi1);
    //cout<<"phi1: "<<phi1<<" x1: "<<x1<<" y1: "<<y1<<" m1: "<<m1<<endl;
    initial.coords[0] = -(y0 - y1 + m1*x1 - m0*x0) / (m0 - m1);
    initial.coords[1] = y0 + m0*(initial.x() - x0);
    //cout<<"initial: "<<initial<<endl;
    return create_edge(v0, v1, gps, initial);
}

/* phi is angular parameter to parametric equations */
/* r0 is initial radius (at phi=0.0) */
/* delta is divergence of spiral in distance/radian */
/* start is starting radians */
double MeshConstructor::spiral_derivative(Radians phi, double r0, Radians delta, Radians start)
{
    double cosphi = cos(phi);
    double sinphi = sin(phi);
    delta /= TWOPI;
    double alpha = r0 + (phi-start) * delta;
    return (alpha*cosphi + delta * sinphi)/ (delta*cosphi - alpha * sinphi);
}


BoundaryEdge* MeshConstructor::add_bezier(BoundaryVertex *v0,
                                          BoundaryVertex *v1,
                                          const Point2D &p0,
                                          const Point2D &p1,
                                          const Point2D &p2 )
{
    vector<Point2D>  gps;
    v0->set_cp(p0);
    v1->set_cp(p2);
    gps.push_back( p0 );
    gps.push_back( p2 );
    return create_edge(v0, v1, gps, p1);
}

void
MeshConstructor::add_box_void(const Point2D &p0, const Point2D &p1, const Point2D &p2, const Point2D &p3, int segments)
{
  add_box(p0,p1,p2,p3,segments);
  return;
}

vector<BoundaryEdge*>
MeshConstructor::add_box(const Point2D &p0, const Point2D &p1, const Point2D &p2, const Point2D &p3, int segments)
{
    BoundaryVertex *v0, *v1, *v2, *v3;
    vector<BoundaryEdge*> es;

    v0 = add_vertex(p0);
    v1 = add_vertex(p1);
    v2 = add_vertex(p2);
    v3 = add_vertex(p3);

    es.push_back( add_line(v0, v1, segments) );
    es.push_back( add_line(v1, v2, segments) );
    es.push_back( add_line(v2, v3, segments) );
    es.push_back( add_line(v3, v0, segments) );

    return es;
}

void MeshConstructor::add_remove_point(const Point2D &p)
{
  mesher->add_hole(p);
}


int MeshConstructor::mesh(bool refine)
{
  /* create a conformal mesh */
  mesher->mesh();

  /* angle should already have bet set in degrees on the boundary facets*/
  if (refine) {
    remove_small_angles();
  }

  return 0;
}

void MeshConstructor::remove_small_angles()
{
    cout << "Begin small angle removal!  Initial size: "
      << bezier->get_num_faces() << endl;
    bezier->set_linear_flips(true);
    bezier->make_delaunay();
    bezier->remove_small_angles();
    bezier->set_linear_flips(false);
    cout << "End small angle removal!  Final size: "
      << bezier->get_num_faces() << endl;
}


void MeshConstructor::add_red_blood_cells_box(int num, const Point2D &box_top, const Point2D &box_bot,
                                              int cell_type, double cell_size, vector<BoundaryVertex*> &verts,
                                              vector<BoundaryEdge*> &cells, vector<Point2D> &centers)
{
    vector<Point2D> initial_cells;
    Point2D p;

    double x_begin = box_bot.x();
    double x_span = box_top.x() - box_bot.x();
    double y_begin = box_bot.y();
    double y_span = box_top.y() - box_bot.y();

    /* generate 30*n initial_cell locations */
    timeval tv;
    gettimeofday(&tv, NULL);
    srand48(tv.tv_sec);
    for(int i=0; i < 30*num; ++i){
        p.assign(x_begin + x_span * drand48(), y_begin + y_span * drand48());
        initial_cells.push_back(p);
    }

    insert_cells(num, initial_cells, cell_type, cell_size, verts, cells, centers);
}

void MeshConstructor::add_red_blood_cells_triangle(int num,
                                                   const Point2D &p1, const Point2D &p2, const Point2D &p3,
                                                   int cell_type, double cell_size,
                                                   vector<BoundaryVertex*> &verts,
                                                   vector<BoundaryEdge*> &cells,
                                                   vector<Point2D> &centers)
{
    vector<Point2D> initial_cells;
    Point2D p;
    double x_begin = min(p1.x(), min(p2.x(), p3.x()));
    double x_end   = max(p1.x(), max(p2.x(), p3.x()));
    double x_span  = x_end - x_begin;
    double y_begin = min(p1.y(), min(p2.y(), p3.y()));
    double y_end   = max(p1.y(), max(p2.y(), p3.y()));
    double y_span  = y_end - y_begin;

    /* generate 30*n initial_cell locations */
    timeval tv;
    gettimeofday(&tv, NULL);
    srand48(tv.tv_sec);
    for(int i=0;;){
        p.assign(x_begin + x_span * drand48(), y_begin + y_span * drand48());
        if(p.line_side_test(p1, p2) > 0.0 && p.line_side_test(p2, p3) > 0.0 && p.line_side_test(p3, p1) > 0.0){
            i++;
            initial_cells.push_back(p);
        }
        if(i >= 30*num -1) break;
    }

    insert_cells(num, initial_cells, cell_type, cell_size, verts, cells, centers);
}

void MeshConstructor::insert_cells(int num, vector<Point2D> initial_cells,
                                   int cell_type, double cell_size, vector<BoundaryVertex*> &verts,
                                   vector<BoundaryEdge*> &cells, vector<Point2D> &centers)
{
    double safe_dist = cell_size*1.1;
    hash_set<int> keep;
    int i;
    hash_set<int>::iterator j;
    bool should_keep;
    Point2D p;
    if(cell_type < 0) return;
    /* find a maximal independant set of centers
    * This is O(n^2) for number of cells, but it is
    * an 'offline' calculation
    */
    for(i=0; i < (int) initial_cells.size(); ++i){
        if(keep.find(i) != keep.end()) continue;
        should_keep = true;
        for(j=keep.begin(); j != keep.end(); ++j){
            if((i != *j) && ( (initial_cells[i] - initial_cells[*j]).mag() < safe_dist) ){
                should_keep=false;
                break;
            }
        }
        if(should_keep) keep.insert(i);
        if((int)keep.size() >= num) break;
        if(i%1000==0 && i!=0) cout<<"created "<<i<<" cells..."<<endl;
    }


    /* add no more than n cells, record Edges and centers */
    BoundaryVertex *v;
    for(i=0, j = keep.begin(); (i < num) && (j != keep.end()); ++i, ++j){
        v=add_vertex(p);
        verts.push_back(v);
        centers.push_back( initial_cells[ *j ] );
        cells.push_back( add_circle(v, centers.back(), cell_size*0.5, 0.0, 12) );
    }
}

---