amroc/amr/MultiStatCPTLevelSet.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2013 Ralf Deiterding
00004 // German Aerospace Center (DLR)
00005 
00006 #ifndef AMROC_MULTISTATCPTLEVELSET_H
00007 #define AMROC_MULTISTATCPTLEVELSET_H
00008 
00016 #include "CPTLevelSet.h"
00017 
00018 #define MAXFILES 40
00019 
00026 template <class DataType, int dim>
00027 class MultiStatCPTLevelSet : public CPTLevelSet<DataType,dim> {
00028   typedef CPTLevelSet<DataType,dim> base;
00029 
00030 public:
00031   typedef typename base::grid_fct_type grid_fct_type;   
00032   typedef typename base::grid_data_type grid_data_type;   
00033   typedef typename cpt::State<dim,DataType> cpt_type;
00034   typedef ads::FixedArray<dim,DataType> point_type;
00035   typedef ads::FixedArray<dim,int> multi_index_type;
00036 
00037   MultiStatCPTLevelSet() : base(), nbreps(0), num_vertices(0), num_connections(0) {
00038     for (register int i=0; i<MAXFILES; i++) {
00039       brep_filetype[i] = 0;
00040       origin[i] = 0.;
00041       morigin[i] = 0.;
00042       rotate[i] = 0.;
00043       translate[i] = 0.;
00044       scale[i] = 1.;
00045       vertices[i] = 0;
00046       connections[i] = 0;
00047     }
00048     base::_Stationary = 1;
00049     pi = 4.0*std::atan(1.0);
00050   }
00051 
00052   virtual ~MultiStatCPTLevelSet() {    
00053     if (vertices[0]!=0) delete[] vertices[0];
00054     if (connections[0]!=0) delete[] connections[0];
00055   }
00056 
00057   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00058     base::LocCtrl = Ctrl.getSubDevice(prefix+"MultiStatCPT");
00059     char VariableName[32];
00060     for (int i=0; i<MAXFILES; i++) {
00061       std::sprintf(VariableName,"Brep(%d)",i+1);    
00062       RegisterAt(base::LocCtrl,VariableName,brep_filename[i]);
00063       std::sprintf(VariableName,"FileType(%d)",i+1);    
00064       RegisterAt(base::LocCtrl,VariableName,brep_filetype[i]);
00065       for (int d=0; d<dim; d++) {
00066         std::sprintf(VariableName,"Origin(%d,%d)",i+1,d+1);
00067         RegisterAt(base::LocCtrl,VariableName,origin[i](d));
00068         std::sprintf(VariableName,"MOrigin(%d,%d)",i+1,d+1);
00069         RegisterAt(base::LocCtrl,VariableName,morigin[i](d));
00070         if (dim==3) {
00071           std::sprintf(VariableName,"Rotate(%d,%d)",i+1,d+1);
00072           RegisterAt(base::LocCtrl,VariableName,rotate[i](d));
00073         }
00074         else if (dim==2) {
00075           std::sprintf(VariableName,"Rotate(%d)",i+1);
00076           RegisterAt(base::LocCtrl,VariableName,rotate[i](0));
00077         }         
00078         std::sprintf(VariableName,"Translate(%d,%d)",i+1,d+1);
00079         RegisterAt(base::LocCtrl,VariableName,translate[i](d));
00080       }
00081       std::sprintf(VariableName,"Scale(%d)",i+1);    
00082       RegisterAt(base::LocCtrl,VariableName,scale[i]);
00083     }
00084     RegisterAt(base::LocCtrl,"NBreps",nbreps);
00085     RegisterAt(base::LocCtrl,"PlotPhi",base::_PlotPhi);
00086     RegisterAt(base::LocCtrl,"Inverse",base::inverse);
00087     RegisterAt(base::LocCtrl,"Unsigned",base::unsign);
00088     RegisterAt(base::LocCtrl,"FillWidth",base::_FillWidth);
00089   } 
00090   virtual void register_at(ControlDevice& Ctrl) {
00091     register_at(Ctrl, "");
00092   }
00093 
00094   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00095     base::SetupData(gh, ghosts);
00096 
00097     if (nbreps > MAXFILES) nbreps=MAXFILES;
00098 
00099     num_vertices = 0; num_connections = 0;
00100     for (int i=0; i<nbreps; i++) {
00101       std::ifstream brep_file( brep_filename[i].c_str() );
00102       if ( ! brep_file )
00103         std::cerr << "Could not read " << brep_filename[i] << std::endl;
00104  
00105       if (brep_filetype[i]==1) {
00106         int fdim, mdim;
00107         brep_file >> fdim >> mdim;
00108         if ( fdim!=dim || mdim!=dim-1 ) {
00109           std::cerr << "Mesh in b-rep file has wrong dimension.\n";
00110           continue;
00111         }
00112         brep_file >> nvertices[i];
00113         vertices[i] = new point_type[ nvertices[i] ];
00114         const point_type* point_iter_end = vertices[i] + nvertices[i];
00115         for ( point_type* piter = vertices[i]; piter != point_iter_end; ++piter )
00116           brep_file >> *piter;
00117         brep_file >> nconnections[i];
00118         connections[i] = new multi_index_type[ nconnections[i] ];
00119       }
00120       else {
00121         brep_file >> nvertices[i] >> nconnections[i];
00122       vertices[i] = new point_type[ nvertices[i] ];
00123       connections[i] = new multi_index_type[ nconnections[i] ];
00124       const point_type* point_iter_end = vertices[i] + nvertices[i];
00125       for ( point_type* piter = vertices[i]; piter != point_iter_end; ++piter )
00126         brep_file >> *piter;
00127       }
00128       
00129       const multi_index_type* connections_iter_end = connections[i] + nconnections[i];
00130       for ( multi_index_type* eiter = connections[i]; eiter != connections_iter_end; ++eiter ) 
00131         brep_file >> *eiter;
00132       
00133       assert( brep_file.good() );
00134       brep_file.close();
00135 
00136       num_vertices += nvertices[i];
00137       num_connections += nconnections[i];
00138     }
00139 
00140     point_type* ver = new point_type[num_vertices];
00141     multi_index_type* con = new multi_index_type[num_connections];
00142     
00143     int iver=0, icon=0, idbase=0;
00144     for (int i=0; i<nbreps; i++) {
00145       if (nvertices[i]<=0) continue;
00146       point_type* current_vertices=vertices[i];
00147       vertices[i] = ver;
00148       for (register int j=0; j<nvertices[i]; j++, iver++)
00149         ver[iver] = ForwardTransformN(i,current_vertices[j]);
00150       delete [] current_vertices;
00151       
00152       multi_index_type* current_connections=connections[i];
00153       connections[i] = con;
00154       for (register int j=0; j<nconnections[i]; j++, icon++) {
00155         con[icon] = current_connections[j];
00156         con[icon] += idbase;
00157       } 
00158       delete [] current_connections;
00159       
00160       idbase += nvertices[i];
00161 
00162       morigin[i] = ForwardTransformN(i,morigin);
00163     }
00164     
00165     base::SetBrep(num_vertices, reinterpret_cast<const DataType*>(vertices[0]), 
00166                   num_connections, reinterpret_cast<const int*>(connections[0]));
00167     
00168   }
00169 
00170   point_type ForwardTransformN(const int i, const point_type& p) {
00171     point_type r=p;
00172     r -= origin[i];
00173     point_type a=rotate[i]/180.*pi;
00174     if (dim==3) {
00175       DataType x,y,z,sx,sy,sz,cx,cy,cz;
00176       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; } 
00177       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; } 
00178       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; }
00179     }
00180     else if (dim==2) { 
00181       DataType x,y,sz,cz;
00182       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; }
00183     }
00184     r *= scale[i];
00185     r += translate[i];
00186     return r;
00187   }
00188 
00189   point_type InverseTransformN(const int i, const point_type& p) {
00190     point_type r=p;
00191     r -= translate[i];
00192     r /= scale[i];
00193     point_type a=-rotate[i]/180.*pi;
00194     if (dim==3) {
00195       DataType x,y,z,sx,sy,sz,cx,cy,cz;
00196       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; } 
00197       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; } 
00198       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; }
00199     }
00200     else if (dim==2) { 
00201       DataType x,y,sz,cz;
00202       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; }
00203     }
00204     r += origin[i];
00205     return r;
00206   }
00207 
00208   inline void SetNBreps(const int _nbreps) { nbreps = _nbreps; }
00209   inline int NBreps() const { return nbreps; }
00210   inline void SetFilename(const int i, const std::string filename) { brep_filename[i] = filename; }
00211   inline std::string Filename(const int i) const { return brep_filename[i]; }
00212 
00213   bool GetNthBrep(const int i, int& nver, const DataType*& ver, int& nconn, const int*& conn) {
00214     nver = nvertices[i]; ver = reinterpret_cast<const DataType*>(vertices[i]);
00215     nconn = nconnections[i]; conn = reinterpret_cast<const int*>(connections[i]);
00216     return ((nver>0) && (nconn>0) && (ver!=0) && (conn!=0));
00217   }
00218 
00219 private:
00220   DataType pi;
00221 
00222 protected:
00223   int nbreps, num_vertices, num_connections;
00224   int nvertices[MAXFILES], nconnections[MAXFILES];
00225   point_type* vertices[MAXFILES];
00226   multi_index_type* connections[MAXFILES];
00227 
00228 public:
00229   std::string brep_filename[MAXFILES];
00230   int brep_filetype[MAXFILES];
00231   point_type origin[MAXFILES], morigin[MAXFILES], rotate[MAXFILES], translate[MAXFILES];
00232   DataType scale[MAXFILES];
00233 };
00234 
00235 
00236 #endif
Generated on Thu Jun 30 02:19:28 2016 for AMROC Fluid-solver Framework by 1.6.3