simulation.C

Go to the documentation of this file.

#ifndef _SIMULATION_C
#define _SIMULATION_C

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


#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cstddef>
#include <iostream>
#include <list>

#include "simulation.h"

#include "beziermesh.h"
#include "boundarymesh.h"
#include "cell.h"
#include "celltuple.h"
#include "cleaner.h"
#include "datastore.h"
#include "epswrite.h"
#include "globals.h"
#include "meshio.h"
#include "spline.h"
#include "timehelp.h"
#include "util.h"

using namespace std;

Simulation::Simulation()
{
    initialize_structures();
}

Simulation::Simulation( char *filename )
{
    initialize_structures();
    int len = strlen(filename);
    if(len>=4 && !strncmp(filename+(len-4),".geo",4))
        read_from_binary_file( filename );
    else
        read_from_text_file( filename );

}

void Simulation::initialize_structures()
{
    persistant_store = new PersistantStore();
    data_store = new DataStore(*persistant_store);
    bdry_mesh = new BoundaryMesh(*persistant_store, data_store);
    bezier_mesh = new BezierMesh(*persistant_store, data_store);
    bdry_mesh->set_bezier_mesh(bezier_mesh);
    bezier_mesh->set_bdry_mesh(bdry_mesh);
    positions_set = false;
}

BezierMesh* Simulation::get_bezier_mesh()
{
    return bezier_mesh;
}

BoundaryMesh* Simulation::get_boundary_mesh()
{
    return bdry_mesh;
}

DataStore* Simulation::get_data_store()
{
    return data_store;
}


int Simulation::read_from_text_file( char *filename )
{
    MeshInput meshin(bezier_mesh, bdry_mesh);
    if(! meshin.read( filename ) ){
        cout<<"Error reading in mesh from file: "<<filename<<endl;
        exit(-1);
    }
    return 0;
}


int Simulation::read_from_binary_file( char *filename )
{
    MeshBinaryInput meshin(this);
    if(! meshin.read( filename ) ){
        cout<<"Error reading in mesh from file: "<<filename<<endl;
        exit(-1);
    }
    return 0;
}

const LinearData& Simulation::get_data(const ControlPoint& cp) const
{
  return *data_store->get_data(cp);
}


static Point2D compute_displacement(const LinearData& data, int fixed,
                                    int v1x, int v1y, int v0x, int v0y, double dt)
{
  Point2D vel0(data.fields[v0x], data.fields[v0y] );
  Point2D vel1(data.fields[v1x], data.fields[v1y] );

  /* zero out fixed components */
  if (fixed == 1 || fixed == 3){
    vel0.coords[0] = 0.0;
    vel1.coords[0] = 0.0;
  }
  if (fixed == 2 || fixed == 3){
    vel0.coords[1] = 0.0;
    vel1.coords[1] = 0.0;
  }
  return ( vel0 * 0.5 + vel1 * 0.5 ) * dt;
}

int Simulation::to_text_file( char *filename )
{
    int ret;
    MeshOutput meshout(bezier_mesh, bdry_mesh);
    if( !meshout.write( filename ) ){
        cout<<"Error writing mesh to file: "<<filename<<endl;
        ret = -1;
    } else {
        ret = 0;
    }
    return ret;
}
int Simulation::to_binary_file( char *filename )
{
    int ret;
    MeshBinaryOutput meshout(this);
    if( !meshout.write( filename ) ){
        cout<<"Error writing binary mesh file: "<<filename<<endl;
        ret = -1;
    } else {
        ret = 0;
    }
    return ret;
}

int Simulation::to_eps( char *epsfile )
{
    EPSWrite eps(bezier_mesh, bdry_mesh);
    if( !eps.write( epsfile ) ) {
        cout<<"Error writing to eps file: "<<epsfile<<endl;
        return -1;
    }
    return 0;
}

void Simulation::print()
{
    cout<<data_store<<endl;
    cout<<bdry_mesh<<endl;
    cout<<bezier_mesh<<endl;
}

void Simulation::move_(double dt, unsigned vxCoord, unsigned vyCoord, unsigned xCoord, unsigned yCoord)
{
  move_(dt, vxCoord, vyCoord, vxCoord, vyCoord, xCoord, yCoord);
}

void Simulation::move_(double dt, unsigned vxCoord, unsigned xCoord)
{
  move_(dt, vxCoord, vxCoord+1, vxCoord, vxCoord+1, xCoord, xCoord+1);
}

void Simulation::move_(double dt, unsigned v1xCoord, unsigned v0xCoord, unsigned xCoord)
{
  move_(dt, v1xCoord, v1xCoord+1, v0xCoord, v0xCoord+1, xCoord, xCoord+1);
}



void Simulation::move_(double dt, unsigned vxCoord1, unsigned vyCoord1, unsigned vxCoord0, unsigned vyCoord0, unsigned xCoord, unsigned yCoord)
{
  //timeval start;
  //gettimeofday(&start,NULL);
  
  //check to make sure requested data values are here
  unsigned data_length = bezier_mesh->data_length();
  if((vxCoord1 >= data_length) || (vxCoord0 >= data_length)
     || (vyCoord1 >= data_length) || (vyCoord0 >= data_length) 
     || (xCoord >= data_length) || (yCoord >= data_length)) {
       cout<<"Simulation::move --- Requested velocity or position components do not exist in mesh data. Mesh Data Size:"<<data_length<<endl;
       cout<<"vxc1 = "<<vxCoord1<<"vyc1 = "<<vyCoord1<<"vxc0 = "<<vxCoord0<<"vyc0 = "<<vyCoord0<<"xc = "<<xCoord<<"yc = "<<yCoord<<endl;
       assert(0);
       exit(1);
     }
     
     /* Step 1: move BoundaryEdge's.*/
     /* Note spline->move will NOT move d0 vertices, we will do this seperately */
     for(BoundaryMesh::Edge_Hash_T::iterator be = bdry_mesh->get_edges_begin();
     be != bdry_mesh->get_edges_end(); ++be) {
       
       BoundaryEdge *e = (*be);
       QBSpline *spline = e->get_spline();
       int size_dp = spline->k.size();
       int size_dq = size_dp - 1;
       int fixed = e->get_fixed();
       Point2D dp[size_dp];
       Point2D dq[size_dq];
       
       //compute dp and dq
       for (int i = 0; i < size_dp; i++) {
     dp[i] = compute_displacement(get_data(spline->b[2*i]), fixed,
                      vxCoord1, vyCoord1,
                      vxCoord0, vyCoord0, dt);
       }
       for (int i = 0; i < size_dq; i++){
     LinearData data;
     data  = get_data(spline->b[2*i  ]) * 0.25;
     data += get_data(spline->b[2*i+1]) * 0.50;
     data += get_data(spline->b[2*i+2]) * 0.25;
     dq[i] = compute_displacement(data, fixed, vxCoord1, vyCoord1,
                      vxCoord0, vyCoord0, dt);
       }

       //       cout << "sim:move --> Calling spline-> move" <<endl;
       spline->move(dp,dq, xCoord, yCoord);
     }

     //     cout << "sim::move --> finished step 1" << endl;
     
     Point2D displacement;

     /* Step 2: Move d0 vertices*/
     for (BoundaryMesh::Vertex_Hash_T::iterator bv = bdry_mesh->get_vertices_begin();
      bv != bdry_mesh->get_vertices_end(); ++bv) {
       BoundaryVertex *v = *bv;
       displacement = v->get_bezier()->get_cp() + 
     compute_displacement(v->get_bezier()->get_dp(),v->get_fixed(), 
                  vxCoord1, vyCoord1,
                  vxCoord0, vyCoord0, dt);

       //cout << "v0 is cp:"<<v->get_bezier()->get_cp()<<" ld: "<<v->get_bezier()->get_dp()<<endl;

       v -> get_bezier() -> set_data(xCoord, displacement.x());
       v -> get_bezier() -> set_data(yCoord, displacement.y());

       //cout << "now v0 has ld: "<<v->get_bezier()->get_dp()<<endl;
     }


     //     cout << "sim::move --> finished step 2" << endl;
     
     /* Step 3: Iterate over all non-boundary vertex and non-boundary edge control points */
     for(BezierMesh::Vertex_Hash_T::iterator vi = bezier_mesh->get_vertices_begin();
     vi!=bezier_mesh->get_vertices_end(); ++vi) {
       
       BezierVertex *v = *vi;
       if(!v->is_boundary()) {
     displacement = v->get_cp() + 
       compute_displacement(v->get_dp(), 0, vxCoord1, vyCoord1,
                vxCoord0, vyCoord0, dt);

     // cout << "v is cp:"<<v->get_cp()<<" ld: "<<v->get_dp()<<endl;

     v -> set_data(xCoord, displacement.x());
     v -> set_data(yCoord, displacement.y());

     // cout << "now v has ld: "<<v->get_dp()<<endl;
       }
     }

    /* Displace the edge midpoints */
    for(BezierMesh::Edge_Hash_T::iterator ei = bezier_mesh->get_edges_begin();
        ei!=bezier_mesh->get_edges_end(); ++ei){
        BezierEdge *e = *ei;
        if(!e->is_boundary()){
      
      LinearData ld = (*(e->get_data_point(0)) + *(e->get_data_point(2))) * 0.5;

      //displacement = Point2D(ld.get(xCoord), ld.get(yCoord));
          displacement = e->get_cp(1) +  compute_displacement(e->get_dp(1), 0,
                              vxCoord1, vyCoord1,
                              vxCoord0, vyCoord0, dt);

      // cout << "e has cp:"<<e->get_cp(1)<<" ld: "<<e->get_dp(1)<<endl;

      e -> set_data(1, xCoord, displacement.x());
      e -> set_data(1, yCoord, displacement.y());

      // cout << "now e has ld: "<<e->get_dp(1)<<endl;
        }
    }

    //    cout << "sim::move --> finished step 3" << endl;

    assert(bezier_mesh->check_consistency());

    // cout<<"Quadratic Moving complete: Time step = "<<dt<<"sec    Elapsed Time = "<<elapsed_time_ms(start)<<"ms"<<endl;

}

void Simulation::set_positions(unsigned xPosI, unsigned yPosI)
{
  XposIndex = xPosI;
  YposIndex = yPosI;
  positions_set = true;
  return;
}

// 3 arguments, t, vx, vy
void Simulation::move(double time_step, unsigned vx_pos, unsigned vy_pos)
{
  move(time_step, vx_pos, vy_pos, vx_pos, vy_pos);
}
// Shorthand, assume y index is always x+1
// 2 arguments, t, vxy
void Simulation::move(double time_step, unsigned vx_pos)
{
  move(time_step,vx_pos, vx_pos+1, vx_pos, vx_pos+1);
}
// 5 arguments, t, v0x, v0y, v1x, v1y
void Simulation::move(double time_step, unsigned vx1_pos, unsigned vy1_pos, unsigned vx0_pos, unsigned vy0_pos)
{
  if (!positions_set)
    {
      cout << "Positions to store the x and y data have not been set nor initialized." << endl;
      assert(0);
    }

  // Put new position in the data
  //cout << "Putting new position into the dangling data vector" << endl;

  move_(time_step, vx1_pos, vy1_pos, vx0_pos, vy0_pos, XposIndex, YposIndex);

  // Swap it into the geometry
  //cout << "Swapping new position into the dangling geometry" << endl;

  data_store->swap_geo_with ( XposIndex, YposIndex);

  return;
}

// 3 arguments, t, vx, vy
double Simulation::move_safe(double time_step, unsigned vx_pos, unsigned vy_pos)
{
  return move_safe(time_step, vx_pos, vy_pos, vx_pos, vy_pos);
}
// Shorthand, assume y index is always x+1
// 2 arguments, t, vxy
double Simulation::move_safe(double time_step, unsigned vx_pos)
{
  return move_safe(time_step,vx_pos, vx_pos+1, vx_pos, vx_pos+1);
}

double Simulation::move_safe(double time_step, unsigned vx1_pos, unsigned vy1_pos, unsigned vx0_pos, unsigned vy0_pos)
{
  move(time_step, vx1_pos, vy1_pos, vx0_pos, vy0_pos);
  
  vector<BezierTriangle*> inverted;
  bezier_mesh->find_inverted_triangles(inverted);
  if (inverted.size() != 0) {
    cout << "[MoveSafe] Found at least one inverting triangle:" << endl;
    (inverted[0]) -> print_mathematica();
    data_store->swap_geo_with ( XposIndex, YposIndex);
    return -1.0; 
  }

  return 1.0;
}

void Simulation::store_state(unsigned vIndexFrom, unsigned vIndexTo, unsigned xIndexTo)
{
  data_store->copy_data_2(vIndexFrom,  vIndexTo);
  data_store->write_geo_into(xIndexTo);
}

void Simulation::restore_state(unsigned vIndexFrom, unsigned vIndexTo, unsigned xIndexFrom)
{
  data_store->copy_data_2(vIndexFrom, vIndexTo);
  data_store->read_geo_from(xIndexFrom);
}

void Simulation::print_statistics()
{
    cout<<"************ Simulation ****************"<<endl;
    cout<<"BoundaryMesh:"<<endl;
    bdry_mesh->print_statistics();
    cout<<"BezierMesh:"<<endl;
    bezier_mesh->print_statistics();
    data_store->print();
}




#endif /* _SIMULATION_C */

---