VertexMeshWriter.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
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 "VertexMeshWriter.hpp"
00030 #include "Version.hpp"
00031 
00035 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00036 struct MeshWriterIterators
00037 {
00039     typename AbstractMesh<ELEMENT_DIM,SPACE_DIM>::NodeIterator* pNodeIter;
00041     typename VertexMesh<ELEMENT_DIM,SPACE_DIM>::VertexElementIterator* pElemIter;
00042 };
00043 
00045 // Implementation
00047 
00048 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00049 VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::VertexMeshWriter(const std::string& rDirectory,
00050                                                            const std::string& rBaseName,
00051                                                            const bool clearOutputDir)
00052     : AbstractMeshWriter<ELEMENT_DIM, SPACE_DIM>(rDirectory, rBaseName, clearOutputDir),
00053       mpMesh(NULL),
00054       mpIters(new MeshWriterIterators<ELEMENT_DIM, SPACE_DIM>),
00055       mpNodeMap(NULL),
00056       mNodeMapCurrentIndex(0)
00057 {
00058     mpIters->pNodeIter = NULL;
00059     mpIters->pElemIter = NULL;
00060 
00061 #ifdef CHASTE_VTK
00062      // Dubious, since we shouldn't yet know what any details of the mesh are.
00063      mpVtkUnstructedMesh = vtkUnstructuredGrid::New();
00064 #endif //CHASTE_VTK
00065 }
00066 
00067 
00068 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00069 VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::~VertexMeshWriter()
00070 {
00071     if (mpIters->pNodeIter)
00072     {
00073         delete mpIters->pNodeIter;
00074         delete mpIters->pElemIter;
00075     }
00076 
00077     delete mpIters;
00078 
00079     if (mpNodeMap)
00080     {
00081         delete mpNodeMap;
00082     }
00083 
00084 #ifdef CHASTE_VTK
00085      // Dubious, since we shouldn't yet know what any details of the mesh are.
00086      mpVtkUnstructedMesh->Delete(); // Reference counted
00087 #endif //CHASTE_VTK
00088 }
00089 
00090 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00091 std::vector<double> VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::GetNextNode()
00092 {
00093     if (mpMesh)
00094     {
00095         // Sanity check
00096         assert(this->mNumNodes == mpMesh->GetNumNodes());
00097 
00098         std::vector<double> coordinates(SPACE_DIM+1);
00099 
00100         // Get the node coordinates using the node iterator (thus skipping deleted nodes)
00101         for (unsigned j=0; j<SPACE_DIM; j++)
00102         {
00103             coordinates[j] = (*(mpIters->pNodeIter))->GetPoint()[j];
00104         }
00105         coordinates[SPACE_DIM] = (*(mpIters->pNodeIter))->IsBoundaryNode();
00106 
00107         ++(*(mpIters->pNodeIter));
00108 
00109         return coordinates;
00110     }
00111     else
00112     {
00113         return AbstractMeshWriter<ELEMENT_DIM,SPACE_DIM>::GetNextNode();
00114     }
00115 }
00116 
00117 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00118 VertexElementData VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::GetNextElementWithFaces()
00119 {
00120     // This method should only be called in 3D
00121     assert(SPACE_DIM == 3);
00122     assert(mpMesh);
00123     assert(this->mNumElements == mpMesh->GetNumElements());
00124 
00125     // Create data structure for this element
00126     VertexElementData elem_data;
00127 
00128     // Store node indices owned by this element
00129     elem_data.NodeIndices.resize((*(mpIters->pElemIter))->GetNumNodes());
00130     for (unsigned j=0; j<elem_data.NodeIndices.size(); j++)
00131     {
00132         unsigned old_index = (*(mpIters->pElemIter))->GetNodeGlobalIndex(j);
00133         elem_data.NodeIndices[j] = mpMesh->IsMeshChanging() ? mpNodeMap->GetNewIndex(old_index) : old_index;
00134     }
00135 
00136     // Store faces owned by this element
00137     elem_data.Faces.resize((*(mpIters->pElemIter))->GetNumFaces());
00138     for (unsigned i=0; i<elem_data.Faces.size(); i++)
00139     {
00140         // Get pointer to this face
00141         VertexElement<ELEMENT_DIM-1, SPACE_DIM>* p_face = (*(mpIters->pElemIter))->GetFace(i);
00142 
00143         // Create data structure for this face
00144         ElementData face_data;
00145 
00146         // Store this face's index
00147         face_data.AttributeValue = p_face->GetIndex();
00148 
00149         // Store node indices owned by this face
00150         face_data.NodeIndices.resize(p_face->GetNumNodes());
00151         for (unsigned j=0; j<face_data.NodeIndices.size(); j++)
00152         {
00153             unsigned old_index = p_face->GetNodeGlobalIndex(j);
00154             face_data.NodeIndices[j] = mpMesh->IsMeshChanging() ? mpNodeMap->GetNewIndex(old_index) : old_index;
00155         }
00156 
00157         // Store this face
00158         elem_data.Faces[i] = face_data;
00159 
00161     }
00162 
00164     ++(*(mpIters->pElemIter));
00165     return elem_data;
00166 }
00167 
00168 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00169 ElementData VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::GetNextElement()
00170 {
00172 
00173     if (mpMesh)
00174     {
00175         assert(this->mNumElements == mpMesh->GetNumElements());
00176 
00177         ElementData elem_data;
00178         elem_data.NodeIndices.resize((*(mpIters->pElemIter))->GetNumNodes());
00179         for (unsigned j=0; j<elem_data.NodeIndices.size(); j++)
00180         {
00181             unsigned old_index = (*(mpIters->pElemIter))->GetNodeGlobalIndex(j);
00182             elem_data.NodeIndices[j] = mpMesh->IsMeshChanging() ? mpNodeMap->GetNewIndex(old_index) : old_index;
00183         }
00185         ++(*(mpIters->pElemIter));
00186         return elem_data;
00187     }
00188     else
00189     {
00190         return AbstractMeshWriter<ELEMENT_DIM, SPACE_DIM>::GetNextElement();
00191     }
00192 }
00193 
00194 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00195 void VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::WriteVtkUsingMesh(VertexMesh<ELEMENT_DIM, SPACE_DIM>& rMesh, std::string stamp)
00196 {
00197 #ifdef CHASTE_VTK
00198     //Make the Vtk mesh
00199     assert(SPACE_DIM==3 || SPACE_DIM == 2);
00200     vtkPoints* p_pts = vtkPoints::New(VTK_DOUBLE);
00201     p_pts->GetData()->SetName("Vertex positions");
00202     for (unsigned node_num=0; node_num<rMesh.GetNumNodes(); node_num++)
00203     {
00204         c_vector<double, SPACE_DIM> position = rMesh.GetNode(node_num)->rGetLocation();
00205         if (SPACE_DIM==2)
00206         {
00207             p_pts->InsertPoint(node_num, position[0], position[1], 0.0);
00208         }
00209         else
00210         {
00211             p_pts->InsertPoint(node_num, position[0], position[1], position[2]);
00212         }
00213     }
00214 
00215     mpVtkUnstructedMesh->SetPoints(p_pts);
00216     p_pts->Delete(); // Reference counted
00217     for (typename VertexMesh<ELEMENT_DIM,SPACE_DIM>::VertexElementIterator iter = rMesh.GetElementIteratorBegin();
00218          iter != rMesh.GetElementIteratorEnd();
00219          ++iter)
00220     {
00221         vtkCell* p_cell;
00222         if (SPACE_DIM == 2)
00223         {
00224             p_cell = vtkPolygon::New();
00225         }
00226         else
00227         {
00228             p_cell = vtkConvexPointSet::New();
00229         }
00230         vtkIdList* p_cell_id_list = p_cell->GetPointIds();
00231         p_cell_id_list->SetNumberOfIds(iter->GetNumNodes());
00232         for (unsigned j = 0; j < iter->GetNumNodes(); ++j)
00233         {
00234             p_cell_id_list->SetId(j, iter->GetNodeGlobalIndex(j));
00235         }
00236         mpVtkUnstructedMesh->InsertNextCell(p_cell->GetCellType(), p_cell_id_list);
00237         p_cell->Delete(); //Reference counted
00238     }
00239 
00240     //Vtk mesh is now made
00241     assert(mpVtkUnstructedMesh->CheckAttributes() == 0);
00242     vtkXMLUnstructuredGridWriter* p_writer = vtkXMLUnstructuredGridWriter::New();
00243     p_writer->SetInput(mpVtkUnstructedMesh);
00244     //Uninitialised stuff arises (see #1079), but you can remove
00245     //valgrind problems by removing compression:
00246     // **** REMOVE WITH CAUTION *****
00247     p_writer->SetCompressor(NULL);
00248     // **** REMOVE WITH CAUTION *****
00249 
00250     std::string vtk_file_name = this->mpOutputFileHandler->GetOutputDirectoryFullPath() + this->mBaseName;
00251     if (stamp != "")
00252     {
00253         vtk_file_name += "_" + stamp;
00254     }
00255     vtk_file_name += ".vtu";
00256 
00257     p_writer->SetFileName(vtk_file_name.c_str());
00258     //p_writer->PrintSelf(std::cout, vtkIndent());
00259     p_writer->Write();
00260     p_writer->Delete(); //Reference counted
00261 #endif //CHASTE_VTK
00262 }
00263 
00264 
00265 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00266 void VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::AddCellData(std::string dataName, std::vector<double> dataPayload)
00267 {
00268 #ifdef CHASTE_VTK
00269     vtkDoubleArray* p_scalars = vtkDoubleArray::New();
00270     p_scalars->SetName(dataName.c_str());
00271     for (unsigned i=0; i<dataPayload.size(); i++)
00272     {
00273         p_scalars->InsertNextValue(dataPayload[i]);
00274     }
00275 
00276     vtkCellData* p_cell_data = mpVtkUnstructedMesh->GetCellData();
00277     p_cell_data->AddArray(p_scalars);
00278     p_scalars->Delete(); //Reference counted
00279 #endif //CHASTE_VTK
00280 }
00281 
00282 
00283 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00284 void VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::AddPointData(std::string dataName, std::vector<double> dataPayload)
00285 {
00286 #ifdef CHASTE_VTK
00287     vtkDoubleArray* p_scalars = vtkDoubleArray::New();
00288     p_scalars->SetName(dataName.c_str());
00289     for (unsigned i=0; i<dataPayload.size(); i++)
00290     {
00291         p_scalars->InsertNextValue(dataPayload[i]);
00292     }
00293 
00294     vtkPointData* p_point_data = mpVtkUnstructedMesh->GetPointData();
00295     p_point_data->AddArray(p_scalars);
00296     p_scalars->Delete(); //Reference counted
00297 #endif //CHASTE_VTK
00298 }
00300 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00301 void VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::WriteFilesUsingMesh(VertexMesh<ELEMENT_DIM,SPACE_DIM>& rMesh)
00302 {
00303 
00304     this->mpMeshReader = NULL;
00305     mpMesh = &rMesh;
00306 
00307     this->mNumNodes = mpMesh->GetNumNodes();
00308     this->mNumElements = mpMesh->GetNumElements();
00309 
00310     typedef typename AbstractMesh<ELEMENT_DIM,SPACE_DIM>::NodeIterator NodeIterType;
00311     mpIters->pNodeIter = new NodeIterType(mpMesh->GetNodeIteratorBegin());
00312 
00313     typedef typename VertexMesh<ELEMENT_DIM,SPACE_DIM>::VertexElementIterator ElemIterType;
00314     mpIters->pElemIter = new ElemIterType(mpMesh->GetElementIteratorBegin());
00315 
00316     // Set up node map if we might have deleted nodes
00317     mNodeMapCurrentIndex = 0;
00318     if (mpMesh->IsMeshChanging())
00319     {
00320         mpNodeMap = new NodeMap(mpMesh->GetNumAllNodes());
00321         for (NodeIterType it = mpMesh->GetNodeIteratorBegin(); it != mpMesh->GetNodeIteratorEnd(); ++it)
00322         {
00323             mpNodeMap->SetNewIndex(it->GetIndex(), mNodeMapCurrentIndex++);
00324         }
00325     }
00326     WriteFiles();
00327 }
00328 
00329 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00330 void VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::WriteFiles()
00331 {
00332     std::string comment = "# " + ChasteBuildInfo::GetProvenanceString();
00333 
00334     // Write node file
00335     std::string node_file_name = this->mBaseName + ".node";
00336     out_stream p_node_file = this->mpOutputFileHandler->OpenOutputFile(node_file_name);
00337 
00338     // Write the node header
00339     unsigned num_attr = 0;
00340     unsigned max_bdy_marker = 1; // as we include boundary node information in the node file
00341     unsigned num_nodes = this->GetNumNodes();
00342 
00343     *p_node_file << num_nodes << "\t";
00344     *p_node_file << SPACE_DIM << "\t";
00345     *p_node_file << num_attr << "\t";
00346     *p_node_file << max_bdy_marker << "\n";
00347     *p_node_file << std::setprecision(6);
00348 
00349     // Write each node's data
00350     for (unsigned item_num=0; item_num<num_nodes; item_num++)
00351     {
00352         std::vector<double> current_item = this->GetNextNode();
00353         *p_node_file << item_num;
00354         for (unsigned i=0; i<SPACE_DIM+1; i++)
00355         {
00356             *p_node_file << "\t" << current_item[i];
00357         }
00358         *p_node_file << "\n";
00359     }
00360     *p_node_file << comment << "\n";
00361     p_node_file->close();
00362 
00363     // Write element file
00364     std::string element_file_name = this->mBaseName + ".cell";
00365     out_stream p_element_file = this->mpOutputFileHandler->OpenOutputFile(element_file_name);
00366 
00367     // Write the element header
00368     unsigned num_elements = this->GetNumElements();
00369 
00370     *p_element_file << num_elements << "\t";
00371     *p_element_file << num_attr << "\n";
00372 
00373     // Write each element's data
00374     for (unsigned item_num=0; item_num<num_elements; item_num++)
00375     {
00376         if (SPACE_DIM == 2) // In 2D, write the node indices owned by this element
00377         {
00378             // Get data for this element
00379             ElementData elem_data = this->GetNextElement();
00380 
00381             // Get the node indices owned by this element
00382             std::vector<unsigned> node_indices = elem_data.NodeIndices;
00383 
00384             // Write this element's index and the number of nodes owned by it to file
00385             *p_element_file << item_num <<  "\t" << node_indices.size();
00386 
00387             // Write the node indices owned by this element to file
00388             for (unsigned i=0; i<node_indices.size(); i++)
00389             {
00390                 *p_element_file << "\t" << node_indices[i];
00391             }
00392 
00394 
00395             // New line
00396             *p_element_file << "\n";
00397         }
00398         else // 3D
00399         {
00400             assert(SPACE_DIM == 3);
00401 
00402             // Get data for this element
00403             VertexElementData elem_data_with_faces = this->GetNextElementWithFaces();
00404 
00405             // Get the node indices owned by this element
00406             std::vector<unsigned> node_indices = elem_data_with_faces.NodeIndices;
00407 
00408             // Write this element's index and the number of nodes owned by it to file
00409             *p_element_file << item_num <<  "\t" << node_indices.size();
00410 
00411             // Write the node indices owned by this element to file
00412             for (unsigned i=0; i<node_indices.size(); i++)
00413             {
00414                 *p_element_file << "\t" << node_indices[i];
00415             }
00416 
00417             // Get the faces owned by this element (if any)
00418             std::vector<ElementData> faces = elem_data_with_faces.Faces;
00419 
00420             // Write the number of faces owned by this element to file
00421             *p_element_file << "\t" << faces.size();
00422 
00423             for (unsigned j=0; j<faces.size(); j++)
00424             {
00425                 // Get the node indices owned by this face
00426                 std::vector<unsigned> face_node_indices = faces[j].NodeIndices;
00427 
00428                 // Write this face's index and the number of nodes owned by it to file
00429                 *p_element_file << "\t" << faces[j].AttributeValue <<  "\t" << face_node_indices.size();
00430 
00431                 // Write the node indices owned by this element to file
00432                 for (unsigned i=0; i<face_node_indices.size(); i++)
00433                 {
00434                     *p_element_file << "\t" << face_node_indices[i];
00435                 }
00436 
00438             }
00439 
00441 
00442             // New line
00443             *p_element_file << "\n";
00444         }
00445     }
00446 
00447     *p_element_file << comment << "\n";
00448     p_element_file->close();
00449 }
00450 
00452 // Explicit instantiation
00454 
00455 template class VertexMeshWriter<1,1>;
00456 template class VertexMeshWriter<1,2>;
00457 template class VertexMeshWriter<1,3>;
00458 template class VertexMeshWriter<2,2>;
00459 template class VertexMeshWriter<2,3>;
00460 template class VertexMeshWriter<3,3>;

Generated on Mon Nov 1 12:35:24 2010 for Chaste by  doxygen 1.5.5