VtkMeshReader.cpp

00001 /*
00002 
00003 Copyright (C) Fujitsu Laboratories of Europe, 2009
00004 
00005 */
00006 
00007 /*
00008 
00009 Copyright (C) University of Oxford, 2005-2011
00010 
00011 University of Oxford means the Chancellor, Masters and Scholars of the
00012 University of Oxford, having an administrative office at Wellington
00013 Square, Oxford OX1 2JD, UK.
00014 
00015 This file is part of Chaste.
00016 
00017 Chaste is free software: you can redistribute it and/or modify it
00018 under the terms of the GNU Lesser General Public License as published
00019 by the Free Software Foundation, either version 2.1 of the License, or
00020 (at your option) any later version.
00021 
00022 Chaste is distributed in the hope that it will be useful, but WITHOUT
00023 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00024 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00025 License for more details. The offer of Chaste under the terms of the
00026 License is subject to the License being interpreted in accordance with
00027 English Law and subject to any action against the University of Oxford
00028 being under the jurisdiction of the English Courts.
00029 
00030 You should have received a copy of the GNU Lesser General Public License
00031 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00032 
00033 */
00034 
00035 
00036 
00037 #ifdef CHASTE_VTK
00038 
00039 #include "VtkMeshReader.hpp"
00040 #include "Exception.hpp"
00041 
00042 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00043 VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::VtkMeshReader(std::string pathBaseName) :
00044     mIndexFromZero(true),
00045     mNumNodes(0),
00046     mNumElements(0),
00047     mNumFaces(0),
00048     mNodesRead(0),
00049     mElementsRead(0),
00050     mFacesRead(0),
00051     mBoundaryFacesRead(0),
00052     mNumElementAttributes(0),
00053     mNumFaceAttributes(0),
00054     mOrderOfElements(1),
00055     mNodesPerElement(4)
00056 {
00057     vtkXMLUnstructuredGridReader* vtk_xml_unstructured_grid_reader;
00058 
00059     // Check file exists
00060     mVtuFile.open(pathBaseName.c_str());
00061     if ( !mVtuFile.is_open() )
00062     {
00063         EXCEPTION("Could not open VTU file: " + pathBaseName);
00064     }
00065     mVtuFile.close();
00066 
00067     // Load the mesh geometry and data from a file
00068     vtk_xml_unstructured_grid_reader = vtkXMLUnstructuredGridReader::New();
00069     vtk_xml_unstructured_grid_reader->SetFileName( pathBaseName.c_str() );
00070     vtk_xml_unstructured_grid_reader->Update();
00071     mpVtkUnstructuredGrid = vtk_xml_unstructured_grid_reader->GetOutput();
00072 
00073     mNumNodes = vtk_xml_unstructured_grid_reader->GetNumberOfPoints();
00074     mNumElements = vtk_xml_unstructured_grid_reader->GetNumberOfCells();
00075 
00076     // Extract the surface faces
00077     mpVtkGeometryFilter = vtkGeometryFilter::New();
00078     mpVtkGeometryFilter->SetInput(mpVtkUnstructuredGrid);
00079     mpVtkGeometryFilter->Update();
00080 
00081     mNumFaces = mpVtkGeometryFilter->GetOutput()->GetNumberOfCells();
00082 
00083     vtk_xml_unstructured_grid_reader->Delete();
00084 }
00085 
00086 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00087 VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::VtkMeshReader(vtkUnstructuredGrid* p_vtkUnstructuredGrid) :
00088     mIndexFromZero(true),
00089     mNumNodes(0),
00090     mNumElements(0),
00091     mNumFaces(0),
00092     mNodesRead(0),
00093     mElementsRead(0),
00094     mFacesRead(0),
00095     mBoundaryFacesRead(0),
00096     mNumElementAttributes(0),
00097     mNumFaceAttributes(0),
00098     mOrderOfElements(1),
00099     mNodesPerElement(4)
00100 {
00101     mpVtkUnstructuredGrid = p_vtkUnstructuredGrid;
00102 
00103     mNumNodes = mpVtkUnstructuredGrid->GetNumberOfPoints();
00104     mNumElements = mpVtkUnstructuredGrid->GetNumberOfCells();
00105 
00106     // Extract the surface faces
00107     mpVtkGeometryFilter = vtkGeometryFilter::New();
00108     mpVtkGeometryFilter->SetInput(mpVtkUnstructuredGrid);
00109     mpVtkGeometryFilter->Update();
00110 
00111     mNumFaces = mpVtkGeometryFilter->GetOutput()->GetNumberOfCells();
00112 }
00113 
00114 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00115 VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::~VtkMeshReader()
00116 {
00117     mpVtkGeometryFilter->Delete();
00118 }
00119 
00120 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00121 unsigned VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetNumElements() const
00122 {
00123     return mNumElements;
00124 }
00125 
00126 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00127 unsigned VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetNumNodes() const
00128 {
00129     return mNumNodes;
00130 }
00131 
00132 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00133 unsigned VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetNumFaces() const
00134 {
00135     return mNumFaces;
00136 }
00137 
00138 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00139 unsigned VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetNumEdges() const
00140 {
00141     return mNumFaces;
00142 }
00143 
00144 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00145 unsigned VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetNumElementAttributes() const
00146 {
00147     return mNumElementAttributes;
00148 }
00149 
00150 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00151 unsigned VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetNumFaceAttributes() const
00152 {
00153     return mNumFaceAttributes;
00154 }
00155 
00156 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00157 void VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::Reset()
00158 {
00159 //    CloseFiles();
00160 //    OpenFiles();
00161 //    ReadHeaders();
00162 
00163     mNodesRead=0;
00164     mElementsRead=0;
00165     mFacesRead=0;
00166     mBoundaryFacesRead=0;
00167 }
00168 
00169 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00170 void VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::Initialize()
00171 {
00172     mpVtkUnstructuredGrid->Initialize();
00173 }
00174 
00175 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00176 std::vector<double> VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetNextNode()
00177 {
00178     if ( mNodesRead >= mNumNodes )
00179     {
00180         EXCEPTION( "Trying to read data for a node that doesn't exist" );
00181     }
00182 
00183     std::vector<double> next_node;
00184 
00185     for (unsigned i = 0; i < 3; i++)
00186     {
00187         next_node.push_back( mpVtkUnstructuredGrid->GetPoint(mNodesRead)[i] );
00188     }
00189 
00190     mNodesRead++;
00191     return next_node;
00192 }
00193 
00194 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00195 ElementData VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetNextElementData()
00196 {
00197     if ( mElementsRead >= mNumElements )
00198     {
00199         EXCEPTION( "Trying to read data for an element that doesn't exist" );
00200     }
00201 
00202     if( !mpVtkUnstructuredGrid->GetCell(mElementsRead)->IsA("vtkTetra") )
00203     {
00204         EXCEPTION("Element is not a vtkTetra");
00205     }
00206 
00207     ElementData next_element_data;
00208 
00209     for (unsigned i = 0; i < mNodesPerElement; i++)
00210     {
00211         next_element_data.NodeIndices.push_back(mpVtkUnstructuredGrid->GetCell(mElementsRead)->GetPointId(i));
00212     }
00213 
00214     // \todo implement method to read element data properly (currently returns zero always...)
00215     next_element_data.AttributeValue = 0;
00216 
00217     mElementsRead++;
00218     return next_element_data;
00219 }
00220 
00221 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00222 ElementData VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetNextFaceData()
00223 {
00224     if ( mBoundaryFacesRead >= mNumFaces)
00225     {
00226         EXCEPTION( "Trying to read data for a boundary element that doesn't exist");
00227     }
00228 
00229     ElementData next_face_data;
00230 
00231     for (unsigned i = 0; i < 3; i++)
00232     {
00233         next_face_data.NodeIndices.push_back(mpVtkGeometryFilter->GetOutput()->GetCell(mBoundaryFacesRead)->GetPointId(i));
00234     }
00235 
00236     mBoundaryFacesRead++;
00237     return next_face_data;
00238 }
00239 
00240 
00241 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00242 void VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetCellData(std::string dataName, std::vector<double>& dataPayload)
00243 {
00244     vtkCellData *p_cell_data = mpVtkUnstructuredGrid->GetCellData();
00245 
00246     if ( !p_cell_data->HasArray(dataName.c_str()) )
00247     {
00248         EXCEPTION("No cell data '" + dataName + "'");
00249     }
00250 
00251     vtkDataArray *p_scalars = p_cell_data->GetArray( dataName.c_str() );
00252     if (p_scalars->GetNumberOfComponents() != 1)
00253     {
00254         EXCEPTION("The cell data '" + dataName + "' is not scalar data.");
00255     }
00256 
00257     dataPayload.clear();
00258     for (unsigned i = 0; i < mNumElements; i++)
00259     {
00260         dataPayload.push_back( p_scalars->GetTuple(i)[0] );
00261     }
00262 }
00263 
00264 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00265 void VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetCellData(std::string dataName, std::vector<c_vector<double,SPACE_DIM> >& dataPayload)
00266 {
00267     vtkCellData *p_cell_data = mpVtkUnstructuredGrid->GetCellData();
00268 
00269     if ( !p_cell_data->HasArray(dataName.c_str()) )
00270     {
00271         EXCEPTION("No cell data '" + dataName + "'");
00272     }
00273 
00274     vtkDataArray *p_scalars = p_cell_data->GetArray( dataName.c_str() );
00275     if (p_scalars->GetNumberOfComponents() != 3)
00276     {
00277         EXCEPTION("The cell data '" + dataName + "' is not 3-vector data.");
00278     }
00279 
00280     dataPayload.clear();
00281     for (unsigned i = 0; i < mNumElements; i++)
00282     {
00283         c_vector <double, SPACE_DIM> data;
00284         for (unsigned j = 0; j < SPACE_DIM; j++)
00285         {
00286             data[j]=  p_scalars->GetTuple(i)[j];
00287         }
00288         dataPayload.push_back(data);
00289     }
00290 }
00291 
00292 
00293 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00294 void VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetPointData(std::string dataName, std::vector<double>& dataPayload)
00295 {
00296     vtkPointData *p_point_data = mpVtkUnstructuredGrid->GetPointData();
00297 
00298     if ( !p_point_data->HasArray(dataName.c_str()) )
00299     {
00300         EXCEPTION("No point data '" + dataName + "'");
00301     }
00302 
00303     vtkDataArray *p_scalars = p_point_data->GetArray( dataName.c_str() );
00304     if (p_scalars->GetNumberOfComponents() != 1)
00305     {
00306         EXCEPTION("The point data '" + dataName + "' is not scalar data.");
00307     }
00308 
00309     dataPayload.clear();
00310 
00311     for (unsigned i = 0; i < mNumNodes; i++)
00312     {
00313         dataPayload.push_back( p_scalars->GetTuple(i)[0] );
00314     }
00315 }
00316 
00317 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00318 void VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetPointData(std::string dataName, std::vector<c_vector<double,SPACE_DIM> >& dataPayload)
00319 {
00320    vtkPointData *p_point_data = mpVtkUnstructuredGrid->GetPointData();
00321 
00322     if ( !p_point_data->HasArray(dataName.c_str()) )
00323     {
00324         EXCEPTION("No point data '" + dataName + "'");
00325     }
00326 
00327     vtkDataArray *p_scalars = p_point_data->GetArray( dataName.c_str() );
00328 
00329     if (p_scalars->GetNumberOfComponents() != 3)
00330     {
00331         EXCEPTION("The point data '" + dataName + "' is not 3-vector data.");
00332     }
00333     dataPayload.clear();
00334     for (unsigned i = 0; i < mNumNodes; i++)
00335     {
00336         c_vector<double, SPACE_DIM> data;
00337         for (unsigned j = 0; j< SPACE_DIM; j++)
00338         {
00339             data[j] = p_scalars->GetTuple(i)[j];
00340         }
00341         dataPayload.push_back( data );
00342     }
00343 }
00344 
00345 
00346 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00347 vtkUnstructuredGrid* VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::OutputMeshAsVtkUnstructuredGrid()
00348 {
00349     return mpVtkUnstructuredGrid;
00350 }
00351 
00353 // Explicit instantiation
00355 
00356 template class VtkMeshReader<1,1>;
00357 template class VtkMeshReader<1,2>;
00358 template class VtkMeshReader<1,3>;
00359 template class VtkMeshReader<2,2>;
00360 template class VtkMeshReader<2,3>;
00361 template class VtkMeshReader<3,3>;
00362 #endif // CHASTE_VTK

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