AdaptiveTetrahedralMesh.cpp

00001 /*
00002 
00003 Copyright (C) Fujitsu Laboratories of Europe, 2009-2010
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_ADAPTIVITY
00038 
00039 #include "AdaptiveTetrahedralMesh.hpp"
00040 #include "HeartConfig.hpp"
00041 #include "OutputFileHandler.hpp"
00042 
00043 AdaptiveTetrahedralMesh::AdaptiveTetrahedralMesh() :
00044     mNumNodes(0),
00045     mNumElements(0),
00046     mNumLocalNodes(0),
00047     mAdaptSuccess(false),
00048     mpDiscreteGeometryConstraints(NULL),
00049     mpErrorMeasure(NULL),
00050     mpAdapt(NULL),
00051     mGoodEdgeRange(0.0),
00052     mBadEdgeCriterion(0.0),
00053     mVerbose(false)
00054 
00055 {
00056         mpVtkUnstructuredGrid = vtkUnstructuredGrid::New();
00057 }
00058 
00059 AdaptiveTetrahedralMesh::~AdaptiveTetrahedralMesh()
00060 {
00061     Reset();//Delete multiple-use pointers
00062     mpVtkUnstructuredGrid->Delete();
00063 }
00064 
00065 void AdaptiveTetrahedralMesh::ConstructFromVtuFile(std::string fileName)
00066 {
00067     vtkXMLUnstructuredGridReader *p_vtk_reader = vtkXMLUnstructuredGridReader::New();
00068     p_vtk_reader->SetFileName(fileName.c_str());
00069     p_vtk_reader->Update();
00070 
00071     mpVtkUnstructuredGrid->DeepCopy(p_vtk_reader->GetOutput());
00072     mpVtkUnstructuredGrid->Update();
00073 
00074     p_vtk_reader->Delete();
00075 
00076     mNumNodes = mpVtkUnstructuredGrid->GetNumberOfPoints();
00077     mNumElements = mpVtkUnstructuredGrid->GetNumberOfCells();
00078     mNumLocalNodes = mNumNodes;       // Since each process has a complete copy of the vtkUnstructuredGrid
00079 }
00080 
00081 void AdaptiveTetrahedralMesh::ConstructFromMesh(AbstractTetrahedralMesh<3,3>* rMesh)
00082 {
00083     vtkPoints *p_pts = vtkPoints::New();
00084     p_pts->SetDataTypeToDouble();
00085     for (unsigned i=0; i<rMesh->GetNumNodes(); i++)
00086     {
00087         p_pts->InsertPoint(i, rMesh->GetNode(i)->rGetLocation().data());
00088     }
00089 
00090     mpVtkUnstructuredGrid->Allocate(rMesh->GetNumNodes(), rMesh->GetNumNodes());
00091     mpVtkUnstructuredGrid->SetPoints(p_pts);
00092     mpVtkUnstructuredGrid->Update();
00093     p_pts->Delete();    // Reference counted
00094 
00095     for (unsigned i=0; i<rMesh->GetNumElements(); i++)
00096     {
00097         vtkTetra *p_tetra = vtkTetra::New();
00098         for (int j = 0; j < 4; ++j)
00099         {
00100             p_tetra->GetPointIds()->SetId(j, rMesh->GetElement(i)->GetNodeGlobalIndex(j));
00101         }
00102         mpVtkUnstructuredGrid->InsertNextCell(p_tetra->GetCellType(), p_tetra->GetPointIds());
00103         mpVtkUnstructuredGrid->Update();
00104         p_tetra->Delete();    // Reference counted
00105     }
00106 
00107     mNumNodes = mpVtkUnstructuredGrid->GetNumberOfPoints();
00108     mNumElements = mpVtkUnstructuredGrid->GetNumberOfCells();
00109     mNumLocalNodes = mNumNodes;       // Since each process has a complete copy of the vtkUnstructuredGrid
00110 }
00111 
00112 void AdaptiveTetrahedralMesh::ConstructFromDistributedMesh(DistributedTetrahedralMesh<3,3>* rMesh)
00113 {
00114     vtkPoints *p_pts = vtkPoints::New();
00115     p_pts->SetDataTypeToDouble();
00116 
00117     std::vector<unsigned> global_node_numbers, halo_nodes;
00118     std::map<unsigned, unsigned> global_to_local_index_map;
00119 
00120 //    vtkUnsignedIntArray *p_scalars = vtkUnsignedIntArray::New();
00121 //    p_scalars->SetName("GlobalNodeNumbers");
00122 
00123     unsigned index = 0;
00124     for (DistributedTetrahedralMesh<3,3>::NodeIterator it=rMesh->GetNodeIteratorBegin();
00125          it != rMesh->GetNodeIteratorEnd();
00126          ++it)
00127     {
00128         p_pts->InsertPoint(index, it->rGetLocation().data());
00129         global_node_numbers.push_back(it->GetIndex());
00130         global_to_local_index_map[it->GetIndex()] = index;
00131         index++;
00132     }
00133 
00134     rMesh->GetHaloNodeIndices(halo_nodes);
00135     for(unsigned i=0; i<halo_nodes.size(); i++)
00136     {
00137         global_node_numbers.push_back(halo_nodes[i]);
00138         p_pts->InsertPoint(index, rMesh->GetNodeOrHaloNode(halo_nodes[i])->rGetLocation().data());
00139         global_to_local_index_map[halo_nodes[i]] = index;
00140         index++;
00141     }
00142 
00143     AddPointData("GlobalNodeNumbers", global_node_numbers);
00144 
00145     mpVtkUnstructuredGrid->Allocate(global_node_numbers.size(), global_node_numbers.size());
00146     mpVtkUnstructuredGrid->SetPoints(p_pts);
00147     mpVtkUnstructuredGrid->Update();
00148     p_pts->Delete();    // Reference counted
00149 
00150     for (DistributedTetrahedralMesh<3,3>::ElementIterator it=rMesh->GetElementIteratorBegin();
00151          it != rMesh->GetElementIteratorEnd();
00152          ++it)
00153     {
00154         vtkTetra *p_tetra = vtkTetra::New();
00155         for (int j = 0; j < 4; ++j)
00156         {
00157             p_tetra->GetPointIds()->SetId(j, global_to_local_index_map[it->GetNodeGlobalIndex(j)]);
00158         }
00159         mpVtkUnstructuredGrid->InsertNextCell(p_tetra->GetCellType(), p_tetra->GetPointIds());
00160         mpVtkUnstructuredGrid->Update();
00161         p_tetra->Delete();    // Reference counted
00162     }
00163 
00164     mNumNodes = rMesh->GetNumNodes();
00165     mNumElements = rMesh->GetNumElements();
00166     mNumLocalNodes = mpVtkUnstructuredGrid->GetNumberOfPoints() - halo_nodes.size();
00167 }
00168 
00169 void AdaptiveTetrahedralMesh::AddPointData(std::string dataName, std::vector<double> dataPayload)
00170 {
00171     vtkDoubleArray *p_scalars = vtkDoubleArray::New();
00172     p_scalars->SetName(dataName.c_str());
00173     for (unsigned i=0; i<dataPayload.size(); i++)
00174     {
00175         p_scalars->InsertNextValue(dataPayload[i]);
00176     }
00177 
00178     mpVtkUnstructuredGrid->GetPointData()->AddArray(p_scalars);
00179     mpVtkUnstructuredGrid->Update();
00180     p_scalars->Delete(); // Reference counted
00181 }
00182 
00183 void AdaptiveTetrahedralMesh::AddPointData(std::string dataName, std::vector<unsigned> dataPayload)
00184 {
00185     vtkUnsignedIntArray *p_scalars = vtkUnsignedIntArray::New();
00186     p_scalars->SetName(dataName.c_str());
00187     for (unsigned i=0; i<dataPayload.size(); i++)
00188     {
00189         p_scalars->InsertNextValue(dataPayload[i]);
00190     }
00191 
00192     mpVtkUnstructuredGrid->GetPointData()->AddArray(p_scalars);
00193     mpVtkUnstructuredGrid->Update();
00194     p_scalars->Delete(); // Reference counted
00195 }
00196 
00197 void AdaptiveTetrahedralMesh::RemoveArray(std::string dataName)
00198 {
00199     mpVtkUnstructuredGrid->GetPointData()->RemoveArray(dataName.c_str());
00200 }
00201 
00202 void AdaptiveTetrahedralMesh::WriteMeshToFile(std::string directory, std::string fileName)
00203 
00204 {
00205     std::string vtk_file_name = directory + fileName;
00206 
00207     vtkXMLUnstructuredGridWriter *vtk_writer = vtkXMLUnstructuredGridWriter::New();
00208     //Uninitialised stuff arises (see #1079), but you can remove
00209     //valgrind problems by removing compression:
00210     // **** REMOVE WITH CAUTION *****
00211     vtk_writer->SetCompressor(NULL);
00212     // **** REMOVE WITH CAUTION *****
00213     vtk_writer->SetFileName(vtk_file_name.c_str());
00214     vtk_writer->SetInput(mpVtkUnstructuredGrid);
00215     vtk_writer->Write();
00216     vtk_writer->Delete();
00217 }
00218 
00219 void AdaptiveTetrahedralMesh::WriteMeshToDistributedFile(std::string directory, std::string fileName)
00220 {
00221     std::string vtk_file_name = directory + fileName;
00222 
00223     vtkXMLPUnstructuredGridWriter *vtk_writer = vtkXMLPUnstructuredGridWriter::New();
00224     //Uninitialised stuff arises (see #1079), but you can remove
00225     //valgrind problems by removing compression:
00226     // **** REMOVE WITH CAUTION *****
00227     vtk_writer->SetCompressor(NULL);
00228     // **** REMOVE WITH CAUTION *****
00229     vtk_writer->SetFileName(vtk_file_name.c_str());
00230     vtk_writer->SetDataModeToBinary();
00231 
00232     vtk_writer->SetNumberOfPieces(PetscTools::GetNumProcs());
00233     vtk_writer->SetGhostLevel(1);
00234     vtk_writer->SetStartPiece(PetscTools::GetMyRank());
00235     vtk_writer->SetEndPiece(PetscTools::GetMyRank());
00236 
00237     vtk_writer->SetInput(mpVtkUnstructuredGrid);
00238     vtk_writer->Write();
00239     vtk_writer->Delete();
00240 }
00241 
00242 vtkUnstructuredGrid* AdaptiveTetrahedralMesh::GetVtkUnstructuredGrid()
00243 {
00244     return mpVtkUnstructuredGrid;
00245 }
00246 
00247 void AdaptiveTetrahedralMesh::SetAdaptCriterion(double range, double criterion)
00248 {
00249     mGoodEdgeRange = range;
00250     mBadEdgeCriterion = criterion;
00251 }
00252 
00253 unsigned AdaptiveTetrahedralMesh::GetNumNodes()
00254 {
00255     return mNumNodes;
00256 }
00257 
00258 unsigned AdaptiveTetrahedralMesh::GetNumLocalNodes()
00259 {
00260     return mNumLocalNodes;
00261 }
00262 
00263 unsigned AdaptiveTetrahedralMesh::GetNumLocalAndHaloNodes()
00264 {
00265     return mpVtkUnstructuredGrid->GetNumberOfPoints();
00266 }
00267 
00268 unsigned AdaptiveTetrahedralMesh::GetNumElements()
00269 {
00270     return mNumElements;
00271 }
00272 
00273 unsigned AdaptiveTetrahedralMesh::GetNumLocalElements()
00274 {
00275     return mpVtkUnstructuredGrid->GetNumberOfCells();
00276 }
00277 
00278 unsigned AdaptiveTetrahedralMesh::GetNumSurfaceElements()
00279 {
00280     return sids.size();
00281 }
00282 
00283 void AdaptiveTetrahedralMesh::CalculateSENListAndSids(double coplanarTolerance)
00284 {
00285     DiscreteGeometryConstraints constraints;
00286 
00287     if (mVerbose) constraints.verbose_on();
00288     constraints.set_coplanar_tolerance( coplanarTolerance );
00289     constraints.set_volume_input(mpVtkUnstructuredGrid);
00290 
00291     constraints.get_coplanar_ids(sids);
00292     constraints.get_surface(SENList);
00293     assert(sids.size()*3==SENList.size());
00294 }
00295 
00296 void AdaptiveTetrahedralMesh::GetGeometryConstraints()
00297 {
00298     assert(mpDiscreteGeometryConstraints==NULL);
00299     mpDiscreteGeometryConstraints = new DiscreteGeometryConstraints;
00300 
00301     if (mVerbose) mpDiscreteGeometryConstraints->verbose_on();
00302     mpDiscreteGeometryConstraints->set_surface_input(mpVtkUnstructuredGrid, SENList, sids);
00303     mpDiscreteGeometryConstraints->get_constraints(max_len);
00304 //    mpDiscreteGeometryConstraints->write_vtk(std::string("sids.vtu"));
00305 }
00306 
00307 void AdaptiveTetrahedralMesh::CalculateErrorMetric()
00308 {
00309     double error           = HeartConfig::Instance()->GetTargetErrorForAdaptivity();
00310     double sigma           = HeartConfig::Instance()->GetSigmaForAdaptivity();
00311     double max_edge_length = HeartConfig::Instance()->GetMaxEdgeLengthForAdaptivity();
00312     double min_edge_length = HeartConfig::Instance()->GetMinEdgeLengthForAdaptivity();
00313     double gradation       = HeartConfig::Instance()->GetGradationForAdaptivity();
00314     unsigned max_nodes     = HeartConfig::Instance()->GetMaxNodesForAdaptivity();
00315     assert(mpErrorMeasure == NULL);
00316     mpErrorMeasure = new ErrorMeasure;
00317 
00318     if (mVerbose) mpErrorMeasure->verbose_on();
00319     mpErrorMeasure->set_input(mpVtkUnstructuredGrid);
00320     mpErrorMeasure->add_field("Vm", error, false, sigma);
00321     mpErrorMeasure->set_max_length(max_edge_length);
00322     mpErrorMeasure->set_max_length(&(max_len[0]), mpVtkUnstructuredGrid->GetNumberOfPoints());
00323     mpErrorMeasure->set_min_length(min_edge_length);
00324     mpErrorMeasure->apply_gradation(gradation);
00325     mpErrorMeasure->set_max_nodes(max_nodes);
00326 
00327     mpErrorMeasure->diagnostics();
00328 }
00329 
00330 double AdaptiveTetrahedralMesh::GetEdgeLengthDistribution(double range)
00331 {
00332     //This is a temporary adaptivity class
00333     Adaptivity an_adapter;
00334     an_adapter.set_from_vtk(mpVtkUnstructuredGrid, true);
00335      return an_adapter.edgeLengthDistribution(range);
00336 }
00337 
00338 void AdaptiveTetrahedralMesh::Adapt()
00339 {
00340     mpAdapt = new Adaptivity;
00341     mAdaptSuccess = false;
00342 
00343     if (mVerbose) mpAdapt->verbose_on();
00344     mpAdapt->set_from_vtk(mpVtkUnstructuredGrid, true);
00345     mpAdapt->set_adapt_sweeps(HeartConfig::Instance()->GetNumberOfAdaptiveSweeps());
00346     mpAdapt->set_surface_mesh(SENList);
00347     mpAdapt->set_surface_ids(sids);
00348 
00349     if ( mpAdapt->edgeLengthDistribution( mGoodEdgeRange ) > mBadEdgeCriterion )
00350     {
00351         mpAdapt->adapt();
00352 
00353         mpAdapt->get_surface_ids(sids);
00354         mpAdapt->get_surface_mesh(SENList);
00355 
00356         vtkUnstructuredGrid *p_new_vtk_unstructured_grid = mpAdapt->get_adapted_vtu();
00357         mpVtkUnstructuredGrid->DeepCopy( p_new_vtk_unstructured_grid );
00358         mpVtkUnstructuredGrid->Update();
00359         p_new_vtk_unstructured_grid->Delete();
00360 
00361         mAdaptSuccess = true;
00362     }
00363     else
00364     {
00366         //mpAdapt->get_adapted_vtu()->Initialize();
00367     }
00368 
00369     mpVtkUnstructuredGrid->GetPointData()->RemoveArray("metric");
00370     mpVtkUnstructuredGrid->Squeeze();
00371 
00372     // Need to free these pointers as fresh instances are required for the next adapt (else odd things happen)
00373     Reset();
00374 
00375     // Update private member variables - note: these assume that the adapt happens sequentially
00376     mNumNodes = mpVtkUnstructuredGrid->GetNumberOfPoints();
00377     mNumElements = mpVtkUnstructuredGrid->GetNumberOfCells();
00378     mNumLocalNodes = mNumNodes;       // Since each process has a complete copy of the vtkUnstructuredGrid
00379 }
00380 
00381 void AdaptiveTetrahedralMesh::AdaptMesh()
00382 {
00383     GetGeometryConstraints();
00384 
00385     CalculateErrorMetric();
00386 
00387     RemoveArray("mean_desired_lengths");
00388     RemoveArray("desired_lengths");
00389 //    SetAdaptCriterion( mGoodEdgeRange, mBadEdgeCriterion );
00390     Adapt();
00391 }
00392 
00393 void AdaptiveTetrahedralMesh::Reset()
00394 {
00395     delete mpDiscreteGeometryConstraints;
00396     delete mpErrorMeasure;
00397     delete mpAdapt;
00398     mpDiscreteGeometryConstraints=NULL;
00399     mpErrorMeasure=NULL;
00400     mpAdapt=NULL;
00401 }
00402 
00403 bool AdaptiveTetrahedralMesh::GetAdaptSuccess()
00404 {
00405     return mAdaptSuccess;
00406 }
00407 
00408 void AdaptiveTetrahedralMesh::MakeVerbose(bool verbose)
00409 {
00410     mVerbose = verbose;
00411 }
00412 
00413 #endif // CHASTE_ADAPTIVITY

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