coupling/Trash/MotionSolverDemo.h

Go to the documentation of this file.

// -*- C++ -*-

#ifndef MOTION_SOLVER_H
#define MOTION_SOLVER_H

#include "Solver.h"

template <class DataType, int dim>
class MotionSolver : public Solver {
public:
  typedef Solver base;
  typedef ads::FixedArray<dim,DataType> point_type;
  typedef ads::FixedArray<dim,int> multi_index_type;

  MotionSolver() : base(), Tr(0.), scale(1.), dtret(0.), brep_filetype(0), num_vertices(0), num_connections(0), steps(0), 
                   initial_vertices(0), vertices(0), velocities(0), connections(0), nodeids(0), pressures(0) {
    OutputName="-";
    CheckpointName="motion";
    for (int i=0; i<dim; i++) xs[i]=0.; 
    pi=4.0*std::atan(1.);
  }  

  ~MotionSolver() {
    if (initial_vertices) delete[] initial_vertices;
    if (vertices) delete[] vertices;
    if (velocities) delete [] velocities;
    if (connections) delete[] connections;
    if (nodeids) delete [] nodeids;
    if (pressures) delete [] pressures;
  }
  
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    LocCtrl = Ctrl.getSubDevice(prefix+"MotionSolver");
    RegisterAt(LocCtrl,"FileType",brep_filetype);
    RegisterAt(LocCtrl,"InputBrep",brep_filename);
    RegisterAt(LocCtrl,"OutputName",OutputName);
    RegisterAt(LocCtrl,"CheckpointName",CheckpointName);

    RegisterAt(LocCtrl,"Scale",scale);
    RegisterAt(LocCtrl,"Center(1)",xs[0]);
    RegisterAt(LocCtrl,"Center(2)",xs[1]);
    RegisterAt(LocCtrl,"Center(3)",xs[2]);
    RegisterAt(LocCtrl,"RotationDuration",Tr);
    RegisterAt(LocCtrl,"dt",dtret);    
  }
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }

  virtual bool setup() {
    std::ifstream brep_file( brep_filename.c_str() );
    if ( ! brep_file ) {
      std::cerr << "Could not read the b-rep file.  Exiting...\n";
      std::exit( 1 );
    }

    // ISS file
    if (brep_filetype==1) {
      int fdim, mdim;
      brep_file >> fdim >> mdim;
      if ( fdim!=dim || mdim!=dim-1 ) {
        std::cerr << "Mesh in b-rep file has wrong dimension.  Exiting...\n";
        std::exit( 1 );
      }
      brep_file >> num_vertices;
      initial_vertices = new point_type[ num_vertices ];
      const point_type* point_iter_end = initial_vertices + num_vertices;
      for ( point_type* piter = initial_vertices; piter != point_iter_end; ++piter ) {
        brep_file >> *piter; 
        *piter *= scale;
      }
      brep_file >> num_connections;
    }
    else {
      brep_file >> num_vertices >> num_connections;
      initial_vertices = new point_type[ num_vertices ];
      const point_type* point_iter_end = initial_vertices + num_vertices;
      for ( point_type* piter = initial_vertices; piter != point_iter_end; ++piter ) {
        brep_file >> *piter;
        *piter *= scale;
      }
    }
      
    connections = new multi_index_type[ num_connections ];
    const multi_index_type* connections_iter_end = connections + num_connections;
    for ( multi_index_type* eiter = connections; eiter != connections_iter_end; ++eiter ) 
      brep_file >> *eiter;
    brep_file.close();

    velocities = new point_type[ num_vertices ];
    vertices = new point_type[ num_vertices ];
    pressures = new DataType[ num_vertices ];  
    nodeids = new int[ num_vertices ];
    for (register int j=0; j<num_vertices; j++)
      nodeids[j] = j;
  
    return true;
  }    
  
  virtual void finish() {} 

  virtual void SetVertices(double& t) {
    for (int i=0; i<num_vertices; i++) {
      DataType x = initial_vertices[i](0)-xs[0];
      DataType y = initial_vertices[i](1)-xs[1];
      DataType rs = std::sqrt(x*x+y*y);
      DataType alpha = std::atan2(y,x);

      vertices[i](0) = xs[0]+rs*std::cos(2.*pi/Tr*t+alpha);
      vertices[i](1) = xs[1]+rs*std::sin(2.*pi/Tr*t+alpha);
      vertices[i](2) = initial_vertices[i](2);
    }
  }

  virtual void SetVelocities(double& t) {
    for (int i=0; i<num_vertices; i++) {
      DataType x = vertices[i](0)-xs[0];
      DataType y = vertices[i](1)-xs[1];
      DataType rs = std::sqrt(x*x+y*y);
      DataType alpha = std::atan2(y,x);
      if (Tr!=0.) {
        DataType vr = 2.*pi/Tr*rs;
        velocities[i](0) = -vr*std::sin(alpha);
        velocities[i](1) =  vr*std::cos(alpha);
        velocities[i](2) = 0.;
      }
      else 
        velocities[i] = 0.;
    }
  }

  virtual void Initialize(double& t, double& dt) {
    t=0.; time=t; dt=dtret;
    SetVertices(t);
    SetVelocities(t);
  }

  virtual void Advance(double& t, double& dt) {
    t+=dt;
    SetVertices(t);
    SetVelocities(t);
    time=t;
    steps++;
    dt=dtret;
  }

  virtual void Output() {
    if (OutputName.c_str()[0] == '-') 
      return;
    char ioname[256];
    sprintf(ioname,"%s_%d.tec",OutputName.c_str(),NSteps());
    std::ofstream outfile(ioname, std::ios::out);

    std::ostream* out;
    outfile.open(ioname, std::ios::out);
    out = new std::ostream(outfile.rdbuf());

    *out << "title = mesh" << std::endl;
    *out << "variables = x y z Pressures Velocity_u Velocity_v Velocity_w VelocityVec" << std::endl;
    *out << "zone t=\"zone 1\", I=" << num_vertices << ", J="
         << num_connections << ", F=FEPOINT, ET=TRIANGLE" << std::endl;
    for (register int j=0;j<num_vertices; j++)
      *out << vertices[j](0) << " " << vertices[j](1) << " " << vertices[j](2) << " " 
           << pressures[j] << " " << velocities[j](0) << " " << velocities[j](1) << " " << velocities[j](2) << " "
           << std::sqrt(velocities[j](0)*velocities[j](0)+velocities[j](1)*velocities[j](1)+velocities[j](2)*velocities[j](2)) 
           << std::endl;
    for (register int j=0;j<num_connections; j++)
      *out << connections[j](0)+1 << " " << connections[j](1)+1 << " " << connections[j](2)+1 << std::endl; 

    outfile.close();
    delete out; 
  }    

  virtual void Restart(double& t, double& dt) {
    char fname[256];
    sprintf(fname,"%s.cp",CheckpointName.c_str());
    std::ifstream ifs(fname, std::ios::in);
    ifs.read((char*)&time,sizeof(double));
    ifs.read((char*)&steps,sizeof(int));
    ifs.close();
    dt=dtret;
  }

  virtual void Checkpointing() {
    char fname[256];
    sprintf(fname,"%s.cp",CheckpointName.c_str());
    std::ofstream ofs(fname, std::ios::out);
    ofs.write((char*)&time,sizeof(double));
    ofs.write((char*)&steps,sizeof(int));
    ofs.close();
  }

  inline double CurrentTime() const { return time; }
  virtual  int NSteps() { return steps; }
 
protected:
  DataType xs[dim], Tr, pi, scale;
  double dtret, time;
  std::string brep_filename;
  std::string OutputName;
  std::string CheckpointName;
  int brep_filetype, num_vertices, num_connections, steps;
  point_type *initial_vertices, *vertices, *velocities;
  multi_index_type* connections;
  int *nodeids;
  DataType *pressures;
  ControlDevice LocCtrl;
};


#endif

---