Basso: /Users/jack/Code/basso_dev/inc/WriteEnsight.h Source File

00001 
00012 #ifndef _WRITE_ENSIGHT_BASSO_H_
00013 #define _WRITE_ENSIGHT_BASSO_H_
00014 
00015 // std includes
00016 #include <stdio.h>
00017 #include <map>
00018 #include <list>
00019 
00020 // Basso includes
00021 #include "basso.h"
00022 #include "Point.h"
00023 #include "Element.h"
00024 
00025 namespace basso {
00026 
00031 class WriteEnsight {
00032 
00033 public:
00034         
00035         // PUBLIC DEFS
00036         
00037         // CONSTRUCTORS 
00038         
00042         WriteEnsight( const String &jobname ) { fJobname=jobname; fGeomFilename=jobname + ".geom"; }
00043         
00044         ~WriteEnsight() { }
00045         
00046         // ACCESSORS
00047         
00048         // MEMBER FUNCTIONS
00055         template < class NoDe, class ElEmLiSt >
00056         void WriteGeometry( const Array<NoDe> &node, const ElEmLiSt &element,
00057                 int pid=1, const String &description="mesh generated by Basso" ) const;
00058         
00065         void WriteScalarField2( const nArray &d, const String &varname, int step ) const;
00066                 
00073         template < class ElEmLiSt >
00074         void WriteElementScalarField2( const nArray &d, const ElEmLiSt &elements, 
00075                 const String &varname, int step=0 );
00076 
00084         void WriteVectorField2( const nArray &d, const String &varname, int sdim, int step ) const;
00085 
00095         template < class ElEmLiSt >
00096         void WriteElementVectorField2( const nArray &d, const ElEmLiSt &elements, 
00097                 const String &varname, int sdim, int step ) const;
00098         
00104         void WriteCase( const list<String> &scalar_vars, 
00105                         const list<String> &vect_vars, const list<Numeric> &times) const;
00106                         
00108         void AddNodalScalarField( const String &fieldname ) { sclr_pernode.insert(fieldname); }
00109         
00111         void AddElementScalarField( const String &fieldname ) { sclr_perelem.insert(fieldname); }
00112         
00114         void AddNodalVectorField( const String &fieldname ) { vctr_pernode.insert(fieldname); }
00115         
00117         void AddElementVectorField( const String &fieldname ) { vctr_perelem.insert(fieldname); }
00118         
00124         void WriteNodalScalarField( const nArray &d, const String &fieldname, Numeric timestep=0.0 );
00125 
00131         template < class ElEmLiSt >
00132         void WriteElementScalarField( const nArray &d, const ElEmLiSt &element,
00133                 const String &fieldname, Numeric timestep=0.0 );
00134 
00141         void WriteNodalVectorField( const nArray &d, int sdim, const String &fieldname, 
00142                 Numeric timestep=0.0 );
00143 
00150         template < class ElEmLiSt >
00151         void WriteElementVectorField( const nArray &d, int sdim, const ElEmLiSt &element,
00152                 const String &fieldname, Numeric timestep=0.0 );
00153                         
00154         void WriteEnsight::WriteCase() const;
00155         
00156 protected:
00157         typedef map< String, list<Element> > sorted_element_type; 
00158         
00160         String ElementName( const Element &e ) const;
00161         
00163         template < class ElEmLiSt >
00164         void SortElements( const ElEmLiSt &element, sorted_element_type &sortedElem ) const;
00165         
00167         template < class ElEmLiSt >
00168         void SortElementId( const ElEmLiSt &element );
00169         
00171         template < class ElEmLiSt >
00172         void WriteEnsight::SortElementData( const ElEmLiSt &element, 
00173                 const nArray &val, map< String, list<Numeric> > &sortedVal, int sdim=1 ) const;
00174                 
00176         void WriteElement( FILE *outFile, const Element &elem ) const;  
00177 
00178 protected:
00179         String fJobname;
00180         String fGeomFilename;
00181         
00182         set<String> sclr_pernode;
00183         set<String> vctr_pernode;
00184         set<String> sclr_perelem;
00185         set<String> vctr_perelem;
00186         
00187         list<Numeric> time_steps;
00188         
00189         map< String, list<int> > sorted_elemid;
00190         
00191         list< Element > *elementDefPtr;
00192         
00193 };
00194 
00195 template < class NoDe, class ElEmLiSt >
00196 void WriteEnsight::WriteGeometry ( const Array<NoDe> &node, 
00197         const ElEmLiSt &element, int pid, const String &description ) const
00198 {
00199         FILE * outFile;
00200         outFile = fopen( fGeomFilename.c_str(), "w" );
00201         if ( outFile == NULL ) {
00202                 WarningMessage("void WriteEnsight::WriteGeometry","Could not open file");
00203                 return;
00204         }
00205         fprintf (outFile, "Ensight Gold Geometry File\n");
00206         fprintf (outFile, "%-79s\n",fGeomFilename.c_str());
00207         fprintf (outFile, "node id assign\n");
00208         fprintf (outFile, "element id assign\n");
00209         fprintf (outFile, "part\n");
00210         fprintf (outFile, "%10d\n",pid);
00211         fprintf (outFile, "%-79s\n",description.c_str());
00212         fprintf (outFile, "coordinates\n");
00213         fprintf (outFile, "%10d\n",node.size());        
00214                 
00215         for ( int i=0; i<3; ++i )  // loop over nodes
00216                 for ( int I=0; I<node.size(); ++I )
00217                         fprintf ( outFile, "%12.5e\n", node[I].x(i) );
00218                 
00219         sorted_element_type sortedElem;
00220         SortElements( element, sortedElem );
00221         sorted_element_type::const_iterator sortItr;
00222         for ( sortItr=sortedElem.begin(); sortItr!=sortedElem.end(); ++sortItr ) { // loop over element lists
00223                 list<Element>::const_iterator eItr=sortItr->second.begin();
00224                 fprintf (outFile, "%-79s\n",(sortItr->first).c_str());
00225                 fprintf (outFile, "%10d\n", sortItr->second.size() );
00226                 for ( ; eItr!=sortItr->second.end(); ++eItr)   // loop over elements in list
00227                         WriteElement( outFile, *eItr );
00228         }       
00229                 
00230         fclose(outFile);        
00231 }
00232 
00233 void WriteEnsight::WriteElement( FILE *outFile, const Element &elem ) const 
00234 {
00235         iArray econn;
00236         econn=elem.Connectivity();
00237 
00238         switch ( elem.Type() ) {
00239         
00240                 case kQUAD9:  // write as a quad8
00241                 for ( int i=0; i<8; ++i ) 
00242                         fprintf ( outFile, "%10d", econn[i]+1 );
00243                 fprintf ( outFile, "\n" );
00244                 break;
00245 
00246                 default:
00247                 for ( int i=0; i<econn.size(); ++i )
00248                         fprintf ( outFile, "%10d", econn[i]+1 );
00249                 fprintf ( outFile, "\n" );
00250                 break;
00251         }       
00252 }
00253 
00254 
00255 String WriteEnsight::ElementName( const Element &e ) const
00256 {
00257         String ename;
00258         switch ( e.Type() ) {
00259                 
00260                 case kPOINT1:
00261                 ename="point";
00262                 break;
00263                 
00264                 case kLINE2:
00265                 ename="bar2";
00266                 break;
00267                 
00268                 case kTRIA3:
00269                 ename="tria3";
00270                 break;
00271                 
00272                 case kQUAD4:
00273                 ename="quad4";
00274                 break;
00275                 
00276                 case kTETRA4:
00277                 ename="tetra4";
00278                 break;
00279                 
00280                 case kHEXA8:
00281                 ename="hexa8";
00282                 break;
00283                 
00284                 case kPRISM6:
00285                 ename="penta6";
00286                 break;
00287                 
00288                 case kLINE3:
00289                 ename="bar3";
00290                 break;
00291                 
00292                 case kTRIA6:
00293                 ename="tria6";
00294                 break;
00295                 
00296                 case kQUAD8:
00297                 ename="quad8";
00298                 break;
00299                 
00300                 case kQUAD9:
00301                 ename="quad8";
00302                 break;
00303                 
00304                 case kTETRA10:
00305                 ename="tetra10";
00306                 break;
00307                 
00308                 case kHEXA20:
00309                 ename="hexa20";
00310                 break;
00311                 
00312                 default:
00313                 ename="unknown";
00314                 WarningMessage("WriteEnsight::WriteElement","unknown or unsupported element type");
00315                 
00316         }
00317         return ename;
00318 }
00319         
00320 template < class ElEmLiSt >
00321 void WriteEnsight::SortElements( const ElEmLiSt &element, 
00322          sorted_element_type &sortedElem ) const
00323 {
00324         typename ElEmLiSt::const_iterator eItr;
00325         for ( eItr=element.begin(); eItr!=element.end(); ++eItr ) {
00326                 String ename=ElementName(*eItr);
00327                 sortedElem[ename].push_back(*eItr);
00328         }
00329 }
00330 
00331 template < class ElEmLiSt >
00332 void WriteEnsight::SortElementData( const ElEmLiSt &element,  
00333                 const nArray &val, map< String, list<Numeric> > &sortedVal, int sdim ) const
00334 {
00335         int ev=0;
00336         typename ElEmLiSt::const_iterator eItr;
00337         for ( eItr=element.begin(); eItr!=element.end(); ++eItr ) {
00338                 String ename=ElementName(*eItr);
00339                 for ( int s=0; s<sdim; ++s )
00340                         sortedVal[ename].push_back(val[ev++]);
00341         }
00342 }
00343 
00344 template < class ElEmLiSt >
00345 void WriteEnsight::SortElementId( const ElEmLiSt &element )  
00346 {
00347         typename ElEmLiSt::const_iterator eItr;
00348         for ( eItr=element.begin(); eItr!=element.end(); ++eItr ) {
00349                 String ename=ElementName(*eItr);
00350                 sorted_elemid[ename].push_back(e);
00351         }
00352 }
00353 
00354 void WriteEnsight::WriteScalarField2( const nArray &d, const String &varname, int step ) const
00355 {
00356         int pid=1;
00357         char stepStr[5];
00358         sprintf(stepStr,"%04d",step); 
00359         String fileName = fJobname + ".ps" + stepStr + "." + varname;
00360         FILE * outFile;
00361         outFile = fopen( fileName.c_str(), "w" );
00362         if ( outFile == NULL ) {
00363                 WarningMessage("void WriteEnsight::WriteScalarField","Could not open file");
00364                 return;
00365         }
00366         fprintf (outFile, "Ensight Gold scalar variable per node file\n");
00367         fprintf (outFile, "part\n");
00368         fprintf (outFile, "%10d\n",pid);
00369         fprintf (outFile, "coordinates\n");
00370         for ( int i=0; i<d.size(); ++i )
00371                 fprintf ( outFile, "%12.5e\n", d[i] ); 
00372         fclose(outFile);
00373 }
00374 
00375 template < class ElEmLiSt >
00376 void WriteEnsight::WriteElementScalarField2( const nArray &d, const ElEmLiSt &elements, 
00377         const String &varname, int step )  
00378 {
00379         int pid=1;
00380         char stepStr[5];
00381         sprintf(stepStr,"%04d",step); 
00382         String fileName = fJobname + ".ps" + stepStr + "." + varname;
00383         FILE * outFile;
00384         outFile = fopen( fileName.c_str(), "w" );
00385         if ( outFile == NULL ) {
00386                 WarningMessage("void WriteEnsight::WriteElementScalarField","Could not open file");
00387                 return;
00388         }
00389         
00390         if ( sorted_elemid.empty() )  // sort element ids to speed thngs up
00391                 SortElementId( elements );
00392         
00393         fprintf (outFile, "Ensight Gold scalar variable per element file\n");
00394         fprintf (outFile, "part\n");
00395         fprintf (outFile, "%10d\n",pid);
00396         
00397         map< String, list<int> >::const_iterator sortItr;
00398         for ( sortItr=sorted_elemid.begin(); sortItr!=sorted_elemid.end(); ++sortItr ) { // loop over element lists
00399                 list<int>::const_iterator eidItr=sortItr->second.begin();
00400                 fprintf (outFile, "%-79s\n",(sortItr->first).c_str());
00401                 for ( ; eidItr!=sortItr->second.end(); ++eidItr)   // loop over elements in list
00402                         fprintf ( outFile, "%12.5e\n", d[*eidItr] );
00403         }
00404         
00405         fclose(outFile);
00406 }
00407 
00408 
00409 void WriteEnsight::WriteVectorField2( const nArray &d, const String &varname, int sdim, int step ) const
00410 {
00411         int pid=1;
00412         char stepStr[5];
00413         sprintf(stepStr,"%04d",step);  
00414         String fileName = fJobname + ".ps" + stepStr + "." + varname;
00415         FILE * outFile;
00416         outFile = fopen( fileName.c_str(), "w" );
00417         if ( outFile == NULL ) {
00418                 WarningMessage("void WriteEnsight::WriteField","Could not open file");
00419                 return;
00420         }
00421         fprintf ( outFile, "Ensight Gold Vector Variable per Node File\n" );
00422         fprintf ( outFile, "part\n" );
00423         fprintf ( outFile, "%10d\n", pid );
00424         fprintf ( outFile, "coordinates\n" );
00425         if ( sdim==1 ) {
00426                 for ( int i=0; i<d.size(); ++i ) 
00427                         fprintf ( outFile, "%12.5e\n", d[i] );
00428                 for ( int i=0; i<2*d.size(); ++i ) 
00429                         fprintf ( outFile, "%12.5e\n", 0.0 );
00430         }
00431         else if ( sdim==2 ) {
00432                 for ( int i=0; i<d.size(); i+=2 ) 
00433                         fprintf ( outFile, "%12.5e\n", d[i] );
00434                 for ( int i=1; i<d.size(); i+=2 ) 
00435                         fprintf ( outFile, "%12.5e\n", d[i] );
00436                 for ( int i=0; i<d.size()/2; ++i ) 
00437                         fprintf ( outFile, "%12.5e\n", 0.0 );
00438         }
00439         else  { // sdim==3
00440                 for ( int i=0; i<d.size(); i+=3 )
00441                         fprintf ( outFile, "%12.5e\n", d[i] );
00442                 for ( int i=1; i<d.size(); i+=3 )
00443                         fprintf ( outFile, "%12.5e\n", d[i] );
00444                 for ( int i=2; i<d.size(); i+=3 )
00445                         fprintf ( outFile, "%12.5e\n", d[i] );
00446         }
00447         fclose(outFile);
00448 }
00449 
00450 void WriteEnsight::WriteCase( const list<String> &scalar_vars, 
00451                 const list<String> &vect_vars, const list<Numeric> &times ) const
00452 {
00453         String fileName = fJobname + ".case";
00454         FILE * outFile;
00455         outFile = fopen( fileName.c_str(), "w" );
00456         if ( outFile == NULL ) {
00457                 WarningMessage("void WriteEnsight::WriteCase","Could not open file");
00458                 return;
00459         }
00460         
00461         fprintf ( outFile, "FORMAT\n" );
00462         fprintf ( outFile, "type: ensight gold\n" );
00463         
00464         fprintf ( outFile, "\nGEOMETRY\n" );
00465         fprintf ( outFile, "model: %-72s\n", fGeomFilename.c_str() );
00466         
00467         fprintf ( outFile, "\nVARIABLE\n" );
00468         list<String>::const_iterator itr;
00469         for ( itr=scalar_vars.begin(); itr!=scalar_vars.end(); ++itr ) {
00470                 String varFile = fJobname + ".ps****." + *itr; 
00471                 fprintf ( outFile, "scalar per node: %s %s\n", itr->c_str(), varFile.c_str() );
00472         }
00473         for ( itr=vect_vars.begin(); itr!=vect_vars.end(); ++itr ) {
00474                 String varFile = fJobname + ".ps****." + *itr; 
00475                 fprintf ( outFile, "vector per node: %s %s\n", itr->c_str(), varFile.c_str() );
00476         }
00477         
00478         int ts=1, start_val=0, inc=1;
00479         fprintf ( outFile, "\nTIME\n" );
00480         fprintf ( outFile, "time set:              %10d\n", ts );
00481         fprintf ( outFile, "number of steps:       %10d\n", times.size() );
00482         fprintf ( outFile, "filename start number: %10d\n", start_val );
00483         fprintf ( outFile, "filename increment:    %10d\n", inc );
00484         fprintf ( outFile, "time values: " );
00485         int npl=1;
00486         list<Numeric>::const_iterator tItr;
00487         for ( tItr=times.begin(); tItr!=times.end(); ++tItr, ++npl ) {
00488                 if ( npl==6 ) {
00489                         npl=0;
00490                         fprintf ( outFile, "\n" );
00491                 }
00492                 fprintf ( outFile, "%12.5e", *tItr );
00493         }
00494                 
00495         fclose(outFile);
00496 }
00497 
00498 void WriteEnsight::WriteNodalScalarField
00499         ( const nArray &d, const String &fieldname, Numeric timestep )  
00500 {
00501         if ( timestep != time_steps.back() )
00502                 time_steps.push_back(timestep);
00503 
00504         if ( sclr_pernode.find(fieldname) == sclr_pernode.end() )
00505                 sclr_pernode.insert(fieldname);
00506                 
00507         WriteScalarField2( d, fieldname, time_steps.size()-1 ); 
00508 }
00509 
00510 template< class ElEmLiSt >
00511 void WriteEnsight::WriteElementScalarField
00512         ( const nArray &d, const ElEmLiSt &element, const String &fieldname, Numeric timestep )  
00513 {
00514                 
00515         if ( timestep != time_steps.back() )
00516                 time_steps.push_back(timestep);
00517 
00518         if ( sclr_perelem.find(fieldname) == sclr_perelem.end() )
00519                 sclr_perelem.insert(fieldname);
00520         
00521         int ts=time_steps.size()-1;
00522         WriteElementScalarField2( d, element, fieldname, ts );          
00523                 
00524 }
00525 
00526 void WriteEnsight::WriteNodalVectorField
00527         ( const nArray &d, int sdim, const String &fieldname, Numeric timestep ) 
00528 {
00529         if ( timestep != time_steps.back() )
00530                 time_steps.push_back(timestep);
00531 
00532         if ( vctr_pernode.find(fieldname) == vctr_pernode.end() )
00533                 vctr_pernode.insert(fieldname);
00534                 
00535         WriteVectorField2( d, fieldname, sdim, time_steps.size()-1 );
00536 }
00537 
00538 template < class ElEmLiSt >
00539 void WriteEnsight::WriteElementVectorField      ( const nArray &d, int sdim, 
00540         const ElEmLiSt &element, const String &fieldname, Numeric timestep )  
00541 {
00542 /*      int pid=1;
00543         char stepStr[5];
00544         sprintf(stepStr,"%04d",step); 
00545         String fileName = fJobname + ".ps" + stepStr + "." + varname;
00546         FILE * outFile;
00547         outFile = fopen( fileName.c_str(), "w" );
00548         if ( outFile == NULL ) {
00549                 WarningMessage("void WriteEnsight::WriteElementVectorField","Could not open file");
00550                 return;
00551         }
00552         
00553         if ( sorted_elemid.empty() )  // sort element ids to speed thngs up
00554                 SortElementId( elements );
00555         
00556         fprintf (outFile, "Ensight Gold vector variable per element file\n");
00557         fprintf (outFile, "part\n");
00558         fprintf (outFile, "%10d\n",pid);
00559         
00560         map< String, list<int> >::const_iterator sortItr;
00561         for ( sortItr=sorted_elemid.begin(); sortItr!=sorted_elemid.end(); ++sortItr ) { // loop over element lists
00562                 list<int>::const_iterator eidItr=sortItr->second.begin();
00563                 fprintf (outFile, "%-79s\n",(sortItr->first).c_str());
00564                 for ( ; eidItr!=sortItr->second.end(); ++eidItr)   // loop over elements in list
00565                         fprintf ( outFile, "%12.5e\n", d[*eidItr] );
00566         }
00567         
00568         fclose(outFile); */
00569         WarningMessage("WriteEnsight::WriteElementVectorField","Not yet implemented");
00570 }
00571 
00572 void WriteEnsight::WriteCase() const
00573 {
00574         String fileName = fJobname + ".case";
00575         FILE * outFile;
00576         outFile = fopen( fileName.c_str(), "w" );
00577         if ( outFile == NULL ) {
00578                 WarningMessage("void WriteEnsight::WriteCase","Could not open file");
00579                 return;
00580         }
00581         
00582         fprintf ( outFile, "FORMAT\n" );
00583         fprintf ( outFile, "type: ensight gold\n" );
00584         
00585         fprintf ( outFile, "\nGEOMETRY\n" );
00586         fprintf ( outFile, "model: %-72s\n", fGeomFilename.c_str() );
00587         
00588         fprintf ( outFile, "\nVARIABLE\n" );
00589         set<String>::const_iterator itr;
00590         String varFile;
00591         for ( itr=sclr_pernode.begin(); itr!=sclr_pernode.end(); ++itr ) {
00592                 varFile = fJobname + ".ps****." + *itr; 
00593                 fprintf ( outFile, "scalar per node: %s %s\n", itr->c_str(), varFile.c_str() );
00594         }
00595         for ( itr=vctr_pernode.begin(); itr!=vctr_pernode.end(); ++itr ) {
00596                 varFile = fJobname + ".ps****." + *itr; 
00597                 fprintf ( outFile, "vector per node: %s %s\n", itr->c_str(), varFile.c_str() );
00598         }
00599         for ( itr=sclr_perelem.begin(); itr!=sclr_perelem.end(); ++itr ) {
00600                 varFile = fJobname + ".ps****." + *itr; 
00601                 fprintf ( outFile, "scalar per element: %s %s\n", itr->c_str(), varFile.c_str() );
00602         }
00603         for ( itr=vctr_perelem.begin(); itr!=vctr_perelem.end(); ++itr ) {
00604                 varFile = fJobname + ".ps****." + *itr; 
00605                 fprintf ( outFile, "vector per element: %s %s\n", itr->c_str(), varFile.c_str() );
00606         }
00607         
00608         int ts=1, start_val=0, inc=1;
00609         fprintf ( outFile, "\nTIME\n" );
00610         fprintf ( outFile, "time set:              %10d\n", ts );
00611         fprintf ( outFile, "number of steps:       %10d\n", time_steps.size() );
00612         fprintf ( outFile, "filename start number: %10d\n", start_val );
00613         fprintf ( outFile, "filename increment:    %10d\n", inc );
00614         fprintf ( outFile, "time values: " );
00615         int npl=1;
00616         list<Numeric>::const_iterator tItr;
00617         for ( tItr=time_steps.begin(); tItr!=time_steps.end(); ++tItr, ++npl ) {
00618                 if ( npl==6 ) {
00619                         npl=0;
00620                         fprintf ( outFile, "\n" );
00621                 }
00622                 fprintf ( outFile, "%12.5e", *tItr );
00623         }
00624                 
00625         fclose(outFile);
00626 }
00627 
00628 } // end namespace
00629 
00630 #endif
00631 
00632