amroc/amr/MultiCPTLevelSet.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2016 Ralf Deiterding
// University of Southampton

#ifndef AMROC_MULTICPTLEVELSET_H
#define AMROC_MULTICPTLEVELSET_H

#include "MultiStatCPTLevelSet.h"

#define MAXFILES 40

template <class DataType, int dim>
class MultiCPTLevelSet : public MultiStatCPTLevelSet<DataType,dim> {
  typedef MultiStatCPTLevelSet<DataType,dim> base;

public:
  typedef typename base::grid_fct_type grid_fct_type;   
  typedef typename base::grid_data_type grid_data_type;   
  typedef typename cpt::State<dim,DataType> cpt_type;
  typedef ads::FixedArray<dim,DataType> point_type;
  typedef ads::FixedArray<dim,int> multi_index_type;

  MultiCPTLevelSet() : base() {
    base::_Stationary = 0;
  }

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl,prefix);
  } 
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }

  virtual void SetPhi(grid_fct_type& phi, const int Time, const int Level, double t) {
    point_type* ver = new point_type[num_vertices];
    multi_index_type* con = new multi_index_type[num_connections];
    
    int iver=0, icon=0, idbase=0;
    for (int i=0; i<nbreps; i++) {
      if (nvertices[i]<=0) continue;
      point_type* current_vertices=vertices[i];
      vertices[i] = ver;
      for (register int j=0; j<nvertices[i]; j++, iver++)
        ver[iver] = ForwardTransformN(i,current_vertices[j]);
      delete [] current_vertices;
      
      multi_index_type* current_connections=connections[i];
      connections[i] = con;
      for (register int j=0; j<nconnections[i]; j++, icon++) {
        con[icon] = current_connections[j];
        con[icon] += idbase;
      } 
      delete [] current_connections;
      
      idbase += nvertices[i];

      morigin[i] = ForwardTransformN(i,morigin);
    }
    
    base::SetBrep(num_vertices, reinterpret_cast<const DataType*>(vertices[0]), 
                  num_connections, reinterpret_cast<const int*>(connections[0]));
    
    base::SetPhi(phi, Time, Level, t);
  }

  virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
    base::SetupData(gh, ghosts);

    if (nbreps > MAXFILES) nbreps=MAXFILES;

    num_vertices = 0; num_connections = 0;
    for (int i=0; i<nbreps; i++) {
      std::ifstream brep_file( brep_filename[i].c_str() );
      if ( ! brep_file )
        std::cerr << "Could not read " << brep_filename[i] << std::endl;
 
      if (brep_filetype[i]==1) {
        int fdim, mdim;
        brep_file >> fdim >> mdim;
        if ( fdim!=dim || mdim!=dim-1 ) {
          std::cerr << "Mesh in b-rep file has wrong dimension.\n";
          continue;
        }
        brep_file >> nvertices[i];
        vertices[i] = new point_type[ nvertices[i] ];
        const point_type* point_iter_end = vertices[i] + nvertices[i];
        for ( point_type* piter = vertices[i]; piter != point_iter_end; ++piter )
          brep_file >> *piter;
        brep_file >> nconnections[i];
        connections[i] = new multi_index_type[ nconnections[i] ];
      }
      else {
        brep_file >> nvertices[i] >> nconnections[i];
      vertices[i] = new point_type[ nvertices[i] ];
      connections[i] = new multi_index_type[ nconnections[i] ];
      const point_type* point_iter_end = vertices[i] + nvertices[i];
      for ( point_type* piter = vertices[i]; piter != point_iter_end; ++piter )
        brep_file >> *piter;
      }
      
      const multi_index_type* connections_iter_end = connections[i] + nconnections[i];
      for ( multi_index_type* eiter = connections[i]; eiter != connections_iter_end; ++eiter ) 
        brep_file >> *eiter;
      
      assert( brep_file.good() );
      brep_file.close();

      num_vertices += nvertices[i];
      num_connections += nconnections[i];
    }

    point_type* ver = new point_type[num_vertices];
    multi_index_type* con = new multi_index_type[num_connections];
    
    int iver=0, icon=0, idbase=0;
    for (int i=0; i<nbreps; i++) {
      if (nvertices[i]<=0) continue;
      point_type* current_vertices=vertices[i];
      vertices[i] = ver;
      for (register int j=0; j<nvertices[i]; j++, iver++)
        ver[iver] = ForwardTransformN(i,current_vertices[j]);
      delete [] current_vertices;
      
      multi_index_type* current_connections=connections[i];
      connections[i] = con;
      for (register int j=0; j<nconnections[i]; j++, icon++) {
        con[icon] = current_connections[j];
        con[icon] += idbase;
      } 
      delete [] current_connections;
      
      idbase += nvertices[i];

      morigin[i] = ForwardTransformN(i,morigin);
    }
    
    base::SetBrep(num_vertices, reinterpret_cast<const DataType*>(vertices[0]), 
                  num_connections, reinterpret_cast<const int*>(connections[0]));
    
  }

  point_type ForwardTransformN(const int i, const point_type& p) {
    point_type r=p;
    r -= origin[i];
    point_type a=rotate[i]/180.*pi;
    if (dim==3) {
      DataType x,y,z,sx,sy,sz,cx,cy,cz;
      if (a(0)!=0.) { cx=std::cos(a(0)); sx=std::sin(a(0)); y=r(1); z=r(2); r(1) = cx*y-sx*z; r(2) = sx*y+cx*z; } 
      if (a(1)!=0.) { cy=std::cos(a(1)); sy=std::sin(a(1)); x=r(0); z=r(2); r(0) = cy*x-sy*z; r(2) = sy*x+cy*z; } 
      if (a(2)!=0.) { cz=std::cos(a(2)); sz=std::sin(a(2)); x=r(0); y=r(1); r(0) = cz*x-sz*y; r(1) = sz*x+cz*y; }
    }
    else if (dim==2) { 
      DataType x,y,sz,cz;
      if (a(0)!=0.) { cz=std::cos(a(0)); sz=std::sin(a(0)); x=r(0); y=r(1); r(0) = cz*x-sz*y; r(1) = sz*x+cz*y; }
    }
    r *= scale[i];
    r += translate[i];
    return r;
  }

  point_type InverseTransformN(const int i, const point_type& p) {
    point_type r=p;
    r -= translate[i];
    r /= scale[i];
    point_type a=-rotate[i]/180.*pi;
    if (dim==3) {
      DataType x,y,z,sx,sy,sz,cx,cy,cz;
      if (a(0)!=0.) { cx=std::cos(a(0)); sx=std::sin(a(0)); y=r(1); z=r(2); r(1) = cx*y-sx*z; r(2) = sx*y+cx*z; } 
      if (a(1)!=0.) { cy=std::cos(a(1)); sy=std::sin(a(1)); x=r(0); z=r(2); r(0) = cy*x-sy*z; r(2) = sy*x+cy*z; } 
      if (a(2)!=0.) { cz=std::cos(a(2)); sz=std::sin(a(2)); x=r(0); y=r(1); r(0) = cz*x-sz*y; r(1) = sz*x+cz*y; }
    }
    else if (dim==2) { 
      DataType x,y,sz,cz;
      if (a(0)!=0.) { cz=std::cos(a(0)); sz=std::sin(a(0)); x=r(0); y=r(1); r(0) = cz*x-sz*y; r(1) = sz*x+cz*y; }
    }
    r += origin[i];
    return r;
  }

  inline void SetNBreps(const int _nbreps) { nbreps = _nbreps; }
  inline int NBreps() const { return nbreps; }
  inline void SetFilename(const int i, const std::string filename) { brep_filename[i] = filename; }
  inline std::string Filename(const int i) const { return brep_filename[i]; }

  bool GetNthBrep(const int i, int& nver, const DataType*& ver, int& nconn, const int*& conn) {
    nver = nvertices[i]; ver = reinterpret_cast<const DataType*>(vertices[i]);
    nconn = nconnections[i]; conn = reinterpret_cast<const int*>(connections[i]);
    return ((nver>0) && (nconn>0) && (ver!=0) && (conn!=0));
  }

private:
  DataType pi;

protected:
  int nbreps, num_vertices, num_connections;
  int nvertices[MAXFILES], nconnections[MAXFILES];
  point_type* vertices[MAXFILES];
  multi_index_type* connections[MAXFILES];

public:
  std::string brep_filename[MAXFILES];
  int brep_filetype[MAXFILES];
  point_type origin[MAXFILES], morigin[MAXFILES], rotate[MAXFILES], translate[MAXFILES];
  DataType scale[MAXFILES];
};


#endif

---