Hdf5ToCmguiConverter.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include "Hdf5ToCmguiConverter.hpp"
00030 #include "CmguiMeshWriter.hpp"
00031 #include "UblasCustomFunctions.hpp"
00032 #include "HeartConfig.hpp"
00033 #include "PetscTools.hpp"
00034 #include "Exception.hpp"
00035 #include "ReplicatableVector.hpp"
00036 #include "DistributedVector.hpp"
00037 #include "DistributedVectorFactory.hpp"
00038 #include "Version.hpp"
00039 #include "GenericMeshReader.hpp"    
00040 
00041 
00042 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00043 void Hdf5ToCmguiConverter<ELEMENT_DIM,SPACE_DIM>::Write(std::string type)
00044 {
00045     out_stream p_file=out_stream(NULL);
00046 
00047     unsigned num_nodes = this->mpReader->GetNumberOfRows();
00048     unsigned num_timesteps = this->mpReader->GetUnlimitedDimensionValues().size();
00049 
00050     DistributedVectorFactory factory(num_nodes);
00051 
00052     Vec data = factory.CreateVec();//for V
00053     Vec data_phie = factory.CreateVec();//for phi_e
00054     Vec data_second_cell = factory.CreateVec();//for the V of the second cell, used in extended bidomain problems.
00055 
00056     for (unsigned time_step=0; time_step<num_timesteps; time_step++)
00057     {
00058         //create the file for this time step
00059         std::stringstream time_step_string;
00060         //unsigned to string
00061         time_step_string << time_step;
00062         if (PetscTools::AmMaster())
00063         {
00064             p_file = this->mpOutputFileHandler->OpenOutputFile(this->mFileBaseName + "_" + time_step_string.str() + ".exnode");
00065         }
00066 
00067         std::vector<ReplicatableVector*> all_data;
00068         unsigned num_vars = this->mpReader->GetVariableNames().size();
00069         for (unsigned var=0; var<num_vars; var++)
00070         {
00071             //read the data for this time step
00072             this->mpReader->GetVariableOverNodes(data, this->mpReader->GetVariableNames()[var], time_step);
00073             ReplicatableVector* p_repl_data = new ReplicatableVector(data);
00074             assert(p_repl_data->GetSize()==num_nodes);
00075             all_data.push_back(p_repl_data);
00076         }
00077 
00078         if(PetscTools::AmMaster())
00079         {
00080             //write provenance info
00081             std::string comment = "! " + ChasteBuildInfo::GetProvenanceString();
00082             *p_file << comment;
00083             //The header first
00084             *p_file << "Group name: " << this->mFileBaseName << "\n";
00085             *p_file << "#Fields=" << num_vars << "\n";
00086             for (unsigned var=0; var<num_vars; var++)
00087             {
00088                 *p_file << " " << var+1 << ") " <<this->mpReader->GetVariableNames()[var]<< " , field, rectangular cartesian, #Components=1" << "\n" << "x.  Value index=1, #Derivatives=0, #Versions=1"<<"\n";
00089                 if (var!=num_vars-1)
00090                 {
00091                     *p_file << "\n";
00092                 }
00093             }
00094 
00095             //write the data
00096             for(unsigned i=0; i<num_nodes; i++)
00097             {
00098                 //cmgui counts nodes from 1
00099                 *p_file << "Node: "<< i+1 << "\n";
00100                 for (unsigned var=0; var<num_vars; var++)
00101                 {
00102                     *p_file  << (*(all_data[var]))[i] << "\n";
00103                 }
00104             }
00105         }
00106         
00107         for (unsigned var=0; var<num_vars; var++)
00108         {
00109            delete all_data[var];
00110         }
00111     }
00112     VecDestroy(data);
00113     VecDestroy(data_phie);
00114     VecDestroy(data_second_cell);
00115 
00116     if(PetscTools::AmMaster())
00117     {
00118         p_file->close();
00119     }
00120 }
00121 
00122 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00123 Hdf5ToCmguiConverter<ELEMENT_DIM,SPACE_DIM>::Hdf5ToCmguiConverter(std::string inputDirectory,
00124                           std::string fileBaseName,
00125                           AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM> *pMesh,
00126                           bool hasBath) :
00127                     AbstractHdf5Converter<ELEMENT_DIM,SPACE_DIM>(inputDirectory, fileBaseName, pMesh, "cmgui_output")
00128 {
00129     // Write the node data out
00130     Write("");
00131 
00132     //Write mesh in a suitable form for cmgui
00133     std::string output_directory =  HeartConfig::Instance()->GetOutputDirectory() + "/cmgui_output";
00134     
00135     CmguiMeshWriter<ELEMENT_DIM,SPACE_DIM> cmgui_mesh_writer(output_directory, HeartConfig::Instance()->GetOutputFilenamePrefix(), false);
00136     //Used to inform the mesh of the data names
00137     std::vector<std::string> field_names = this->mpReader->GetVariableNames();
00138     cmgui_mesh_writer.SetAdditionalFieldNames(field_names);
00139     if (hasBath)
00140     {
00141         std::vector<std::string> names;
00142         names.push_back("tissue");
00143         names.push_back("bath");
00144         cmgui_mesh_writer.SetRegionNames(names);
00145     }
00146     
00147    
00148     // Normally the in-memory mesh is converted:
00149     if (HeartConfig::Instance()->GetOutputUsingOriginalNodeOrdering() == false )
00150     {
00151         cmgui_mesh_writer.WriteFilesUsingMesh(*(this->mpMesh));
00152     }
00153     else
00154     {
00155         //In this case we expect the mesh to have been read in from file
00157         //Note that the next line will throw if the mesh has not been read from file
00158         std::string original_file=this->mpMesh->GetMeshFileBaseName();
00159         GenericMeshReader<ELEMENT_DIM, SPACE_DIM> original_mesh_reader(original_file);
00160         cmgui_mesh_writer.WriteFilesUsingMeshReader(original_mesh_reader);
00161     }
00162     
00163     WriteCmguiScript();
00164     PetscTools::Barrier("Hdf5ToCmguiConverter");
00165 }
00166 
00167 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00168 void Hdf5ToCmguiConverter<ELEMENT_DIM,SPACE_DIM>::WriteCmguiScript()
00169 {
00170     unsigned num_timesteps = this->mpReader->GetUnlimitedDimensionValues().size();
00171     assert(this->mpReader->GetVariableNames().size() > 0);//seg fault guard
00172     std::string variable_name = this->mpReader->GetVariableNames()[0];
00173 
00174     if(PetscTools::AmMaster())
00175     {
00176         out_stream p_script_file = this->mpOutputFileHandler->OpenOutputFile("script.com");
00177         //write provenance info, note the # instead of ! because this is - essentially - a PERL script that Cmgui interpretes
00178         std::string comment = "# " + ChasteBuildInfo::GetProvenanceString();
00179         *p_script_file << comment;
00180 
00181         *p_script_file << "# Read the mesh \n"
00182                        << "gfx read node "<<HeartConfig::Instance()->GetOutputFilenamePrefix()<<".exnode \n"
00183                        << "gfx read elem "<<HeartConfig::Instance()->GetOutputFilenamePrefix()<<".exelem generate_faces_and_lines \n" //note the mesh file name is taken from HeartConfig
00184                        << "# Create a window \n"
00185                        << "gfx cre win 1 \n"
00186                        << "# Modify the scene (obtained by gfx list g_element XXXX commands) to visualize first var on lines and nodes \n"
00187                        << "gfx modify g_element "<< HeartConfig::Instance()->GetOutputFilenamePrefix()<<" general clear circle_discretization 6 default_coordinate coordinates element_discretization \"4*4*4\" native_discretization none; \n"
00188                        << "gfx modify g_element "<< HeartConfig::Instance()->GetOutputFilenamePrefix()<<" lines select_on material default data "<<variable_name<<" spectrum default selected_material default_selected; \n"
00189                        << "gfx modify g_element "<< HeartConfig::Instance()->GetOutputFilenamePrefix()<<" node_points glyph point general size \"1*1*1\" centre 0,0,0 font default select_on material default data "<<variable_name<<" spectrum default selected_material default_selected; \n"
00190                        << "# Load the data \n"
00191                        << "for ($i=0; $i<" << num_timesteps << "; $i++) { \n"
00192                        << "    gfx read node " << this->mFileBaseName << "_$i.exnode time $i\n" //..while the data file from mFileBaseName...
00193                        << "}\n";
00194         p_script_file->close();
00195     }
00196 }
00198 // Explicit instantiation
00200 
00201 template class Hdf5ToCmguiConverter<1,1>;
00202 template class Hdf5ToCmguiConverter<1,2>;
00203 template class Hdf5ToCmguiConverter<2,2>;
00204 template class Hdf5ToCmguiConverter<1,3>;
00205 template class Hdf5ToCmguiConverter<2,3>;
00206 template class Hdf5ToCmguiConverter<3,3>;

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5