FineCoarseMeshPair.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 
00030 #include "FineCoarseMeshPair.hpp"
00031 
00032 
00033 template<unsigned DIM>
00034 FineCoarseMeshPair<DIM>::FineCoarseMeshPair(TetrahedralMesh<DIM,DIM>& rFineMesh, QuadraticMesh<DIM>& rCoarseMesh)
00035     : mrFineMesh(rFineMesh),
00036       mrCoarseMesh(rCoarseMesh)
00037 {
00038     mpFineMeshBoxCollection = NULL;
00039     mpCoarseMeshBoxCollection = NULL;
00040     ResetStatisticsVariables();
00041 }
00042 
00043 
00044 
00045 template<unsigned DIM>
00046 FineCoarseMeshPair<DIM>::~FineCoarseMeshPair()
00047 {
00048     DeleteFineBoxCollection();
00049     DeleteCoarseBoxCollection();
00050 }
00051 
00052 template<unsigned DIM>
00053 void FineCoarseMeshPair<DIM>::DeleteFineBoxCollection()
00054 {
00055     if(mpFineMeshBoxCollection != NULL)
00056     {
00057         delete mpFineMeshBoxCollection;
00058         mpFineMeshBoxCollection = NULL;
00059     }
00060 }
00061 
00062 
00063 template<unsigned DIM>
00064 void FineCoarseMeshPair<DIM>::DeleteCoarseBoxCollection()
00065 {
00066     if(mpCoarseMeshBoxCollection != NULL)
00067     {
00068         delete mpCoarseMeshBoxCollection;
00069         mpCoarseMeshBoxCollection = NULL;
00070     }
00071 }
00072 
00074 //
00075 //
00076 //   Setting up boxes methods
00077 //
00078 //
00080 
00081 template<unsigned DIM>
00082 void FineCoarseMeshPair<DIM>::SetUpBoxesOnFineMesh(double boxWidth)
00083 {
00084     SetUpBoxes(mrFineMesh, boxWidth, mpFineMeshBoxCollection);
00085 }
00086 
00087 template<unsigned DIM>
00088 void FineCoarseMeshPair<DIM>::SetUpBoxesOnCoarseMesh(double boxWidth)
00089 {
00090     SetUpBoxes(mrCoarseMesh, boxWidth, mpCoarseMeshBoxCollection);
00091 }
00092 
00093 
00094 template<unsigned DIM>
00095 void FineCoarseMeshPair<DIM>::SetUpBoxes(TetrahedralMesh<DIM,DIM>& rMesh,
00096                                          double boxWidth,
00097                                          BoxCollection<DIM>*& rpBoxCollection)
00098 {
00099     if(rpBoxCollection)
00100     {
00101         delete rpBoxCollection;
00102         rpBoxCollection = NULL;
00103     }
00104     
00105     // compute min and max values for the fine mesh nodes
00106     ChasteCuboid<DIM> bounding_box = rMesh.CalculateBoundingBox();
00107     
00108     // set up the boxes. Use a domain which is a touch larger than the fine mesh
00109     c_vector<double,2*DIM> extended_min_and_max;
00110     for(unsigned i=0; i<DIM; i++)
00111     {
00112         double width = bounding_box.GetWidth(i);
00113  
00114         // subtract from the minima
00115         extended_min_and_max(2*i) = bounding_box.rGetLowerCorner()[i] - 0.05*width;
00116         // add to the maxima
00117         extended_min_and_max(2*i+1) = bounding_box.rGetUpperCorner()[i] + 0.05*width;
00118     }
00119 
00120 
00121 
00122     if(boxWidth < 0)
00123     {
00124         // use default value = max(max_edge_length, w20),  where w20 is the width corresponding to 
00125         // 20 boxes in the x-direction 
00126 
00127         // BoxCollection creates an extra box so divide by 19 not 20.  Add a little bit on to ensure
00128         // minor numerical fluctuations don't change the answer.
00129         boxWidth = (extended_min_and_max(1) - extended_min_and_max(0))/19.000000001;
00130 
00131         // determine the maximum edge length
00132         double max_edge_length = -1;
00133 
00134         for (typename TetrahedralMesh<DIM,DIM>::EdgeIterator edge_iterator = rMesh.EdgesBegin();
00135              edge_iterator!=rMesh.EdgesEnd();
00136              ++edge_iterator)
00137         {
00138             c_vector<double, 3> location1 = edge_iterator.GetNodeA()->rGetLocation();
00139             c_vector<double, 3> location2 = edge_iterator.GetNodeB()->rGetLocation();
00140             double edge_length = norm_2(location1-location2);
00141 
00142             if(edge_length>max_edge_length)
00143             {
00144                 max_edge_length = edge_length;
00145             }
00146         }
00147         
00148         if(boxWidth < max_edge_length)
00149         {
00150             boxWidth = 1.1*max_edge_length;
00151         }
00152     }
00153  
00154     rpBoxCollection = new BoxCollection<DIM>(boxWidth, extended_min_and_max);
00155     rpBoxCollection->SetupAllLocalBoxes();
00156 
00157     // for each element, if ANY of its nodes are physically in a box, put that element
00158     // in that box
00159     for(unsigned i=0; i<rMesh.GetNumElements(); i++)
00160     {
00161         Element<DIM,DIM>* p_element = rMesh.GetElement(i);
00162 
00163         std::set<unsigned> box_indices_each_node_this_elem;
00164         for(unsigned j=0; j<DIM+1; j++) // num vertices per element
00165         {
00166             Node<DIM>* p_node = p_element->GetNode(j);
00167             unsigned box_index = rpBoxCollection->CalculateContainingBox(p_node);
00168             box_indices_each_node_this_elem.insert(box_index);
00169         }
00170 
00171         for(std::set<unsigned>::iterator iter = box_indices_each_node_this_elem.begin();
00172             iter != box_indices_each_node_this_elem.end();
00173             ++iter)
00174         {
00175             rpBoxCollection->rGetBox( *iter ).AddElement(p_element);
00176         }
00177     }
00178 }
00179 
00180 
00182 // 
00183 // 
00184 // ComputeFineElementsAndWeightsForCoarseQuadPoints() 
00185 // and 
00186 // ComputeFineElementsAndWeightsForCoarseNodes()
00187 // and 
00188 // common method
00189 //
00190 //
00192 
00193 template<unsigned DIM>
00194 void FineCoarseMeshPair<DIM>::ComputeFineElementsAndWeightsForCoarseQuadPoints(GaussianQuadratureRule<DIM>& rQuadRule,
00195                                                                                bool safeMode)
00196 {
00197     if(mpFineMeshBoxCollection==NULL)
00198     {
00199         EXCEPTION("Call SetUpBoxesOnFineMesh() before ComputeFineElementsAndWeightsForCoarseQuadPoints()");
00200     }
00201 
00202     // get the quad point (physical) positions
00203     QuadraturePointsGroup<DIM> quad_point_posns(mrCoarseMesh, rQuadRule);
00204 
00205     // resize the elements and weights vector.
00206     mFineMeshElementsAndWeights.resize(quad_point_posns.Size());
00207 
00208     #ifdef FINECOARSEMESHPAIR_VERBOSE
00209     std::cout << "\nComputing fine elements and weights for coarse quad points\n";
00210     #endif
00211 
00212     ResetStatisticsVariables();
00213     for(unsigned i=0; i<quad_point_posns.Size(); i++)
00214     {
00215         #ifdef FINECOARSEMESHPAIR_VERBOSE
00216         std::cout << "\r " << i << " of " << quad_point_posns.Size() << std::flush;
00217         #endif
00218         
00219         // get the box this point is in
00220         unsigned box_for_this_point = mpFineMeshBoxCollection->CalculateContainingBox( quad_point_posns.Get(i) );
00221 
00222         // a chaste point version of the c-vector is needed for the GetContainingElement call.
00223         ChastePoint<DIM> point(quad_point_posns.Get(i));
00224 
00225         ComputeFineElementAndWeightForGivenPoint(point, safeMode, box_for_this_point, i);
00226     }
00227 
00228     if(mStatisticsCounters[1]>0)
00229     {
00230         std::stringstream stream;
00231         stream << mStatisticsCounters[1] << " of " << quad_point_posns.Size() << " coarse-mesh quadrature points were outside the fine mesh"; 
00232         WARNING(stream.str());
00233     }
00234 }
00235 
00236 
00237 template<unsigned DIM>
00238 void FineCoarseMeshPair<DIM>::ComputeFineElementsAndWeightsForCoarseNodes(bool safeMode)
00239 {
00240     if(mpFineMeshBoxCollection==NULL)
00241     {
00242         EXCEPTION("Call SetUpBoxesOnFineMesh() before ComputeFineElementsAndWeightsForCoarseNodes()");
00243     }
00244 
00245     // resize the elements and weights vector.
00246     mFineMeshElementsAndWeights.resize(mrCoarseMesh.GetNumNodes());
00247 
00248     #ifdef FINECOARSEMESHPAIR_VERBOSE
00249     std::cout << "\nComputing fine elements and weights for coarse nodes\n";
00250     #endif
00251 
00252     ResetStatisticsVariables();
00253     for(unsigned i=0; i<mrCoarseMesh.GetNumNodes(); i++)
00254     {
00255         #ifdef FINECOARSEMESHPAIR_VERBOSE
00256         std::cout << "\r " << i << " of " << mrCoarseMesh.GetNumNodes() << std::flush;
00257         #endif
00258         
00259         Node<DIM>* p_node = mrCoarseMesh.GetNode(i);
00260         
00261         // get the box this point is in
00262         unsigned box_for_this_point = mpFineMeshBoxCollection->CalculateContainingBox( p_node->rGetModifiableLocation() );
00263 
00264         // a chaste point version of the c-vector is needed for the GetContainingElement call.
00265         ChastePoint<DIM> point(p_node->rGetLocation());
00266 
00267         ComputeFineElementAndWeightForGivenPoint(point, safeMode, box_for_this_point, i);
00268     }
00269 }
00270 
00271 
00276 template<unsigned DIM>
00277 void FineCoarseMeshPair<DIM>::ComputeFineElementAndWeightForGivenPoint(ChastePoint<DIM>& rPoint, 
00278                                                                        bool safeMode,
00279                                                                        unsigned boxForThisPoint,
00280                                                                        unsigned index)
00281 {
00282     std::set<unsigned> test_element_indices;
00283 
00284     // the elements to try (initially) are those contained in the box the point is in
00285     // NOTE: it is possible the point to be in an element inot 'in' this box, as it is possible
00286     // for all element nodes to be in different boxes.
00287     CollectElementsInContainingBox(mpFineMeshBoxCollection, boxForThisPoint, test_element_indices);
00288 
00289     unsigned elem_index;
00290     c_vector<double,DIM+1> weight;
00291 
00292     try
00293     {
00294         //std::cout << "\n" << "# test elements initially " << test_element_indices.size() << "\n";
00295         // try these elements only, initially
00296         elem_index = mrFineMesh.GetContainingElementIndex(rPoint,
00297                                                           false,
00298                                                           test_element_indices,
00299                                                           true /* quit if not in test_elements */);
00300         weight = mrFineMesh.GetElement(elem_index)->CalculateInterpolationWeights(rPoint);
00301 
00302         mStatisticsCounters[0]++;
00303     }
00304     catch(Exception& e)
00305     {
00306         // now try all the elements, but trying the elements contained in the boxes local to this
00307         // element first
00308         std::set<unsigned> test_element_indices;
00309 
00310         CollectElementsInLocalBoxes(mpFineMeshBoxCollection, boxForThisPoint, test_element_indices);
00311 
00312         try
00313         {
00314             elem_index = mrFineMesh.GetContainingElementIndex(rPoint,
00315                                                               false,
00316                                                               test_element_indices,
00317                                                               true);
00318             weight = mrFineMesh.GetElement(elem_index)->CalculateInterpolationWeights(rPoint);
00319 
00320             mStatisticsCounters[0]++;
00321 
00322         }
00323         catch(Exception& e)
00324         {
00325             if(safeMode)
00326             {
00327                 // try the remaining elements
00328                 try
00329                 {
00330                     elem_index = mrFineMesh.GetContainingElementIndex(rPoint,
00331                                                                       false);
00332                     weight = mrFineMesh.GetElement(elem_index)->CalculateInterpolationWeights(rPoint);
00333                     mStatisticsCounters[0]++;
00334 
00335                 }
00336                 catch (Exception& e)
00337                 {
00338                     // the point is not in ANY element, store the nearest element and corresponding weights
00339                     elem_index = mrFineMesh.GetNearestElementIndexFromTestElements(rPoint,test_element_indices);
00340                     weight = mrFineMesh.GetElement(elem_index)->CalculateInterpolationWeights(rPoint);
00341 
00342                     mNotInMesh.push_back(index);
00343                     mNotInMeshNearestElementWeights.push_back(weight);
00344                     mStatisticsCounters[1]++;
00345                 }
00346             }
00347             else
00348             {
00349                 assert(test_element_indices.size()>0); // boxes probably too small if this fails
00350                 
00351                 // immediately assume it isn't in the rest of the mesh - this should be the 
00352                 // case assuming the box width was chosen suitably.
00353                 // Store the nearest element and corresponding weights
00354                 elem_index = mrFineMesh.GetNearestElementIndexFromTestElements(rPoint,test_element_indices);
00355                 weight = mrFineMesh.GetElement(elem_index)->CalculateInterpolationWeights(rPoint);
00356 
00357                 mNotInMesh.push_back(index);
00358                 mNotInMeshNearestElementWeights.push_back(weight);
00359                 mStatisticsCounters[1]++;
00360             }
00361         }
00362     }
00363 
00364     mFineMeshElementsAndWeights[index].ElementNum = elem_index;
00365     mFineMeshElementsAndWeights[index].Weights = weight;    
00366 }
00367 
00369 //
00370 //
00371 // ComputeCoarseElementsForFineNodes
00372 // and 
00373 // ComputeCoarseElementsForFineElementCentroids
00374 // and
00375 // common method
00376 //
00377 //
00379 
00380 
00381 
00382 template<unsigned DIM>
00383 void FineCoarseMeshPair<DIM>::ComputeCoarseElementsForFineNodes(bool safeMode)
00384 {
00385     if(mpCoarseMeshBoxCollection==NULL)
00386     {
00387         EXCEPTION("Call SetUpBoxesOnCoarseMesh() before ComputeCoarseElementsForFineNodes()");
00388     }
00389 
00390     #ifdef FINECOARSEMESHPAIR_VERBOSE
00391     std::cout << "\nComputing coarse elements for fine nodes\n";
00392     #endif
00393 
00394     mCoarseElementsForFineNodes.resize(mrFineMesh.GetNumNodes());
00395     
00396     ResetStatisticsVariables();
00397     for(unsigned i=0; i<mCoarseElementsForFineNodes.size(); i++)
00398     {
00399         #ifdef FINECOARSEMESHPAIR_VERBOSE
00400         std::cout << "\r " << i << " of " << mCoarseElementsForFineNodes.size() << std::flush;
00401         #endif
00402         
00403         ChastePoint<DIM> point = mrFineMesh.GetNode(i)->GetPoint();
00404         // get the box this point is in
00405         unsigned box_for_this_point = mpCoarseMeshBoxCollection->CalculateContainingBox( mrFineMesh.GetNode(i)->rGetModifiableLocation() );
00406 
00407         mCoarseElementsForFineNodes[i] = ComputeCoarseElementForGivenPoint(point, safeMode, box_for_this_point);
00408     }
00409 }
00410 
00411 
00412 template<unsigned DIM>
00413 void FineCoarseMeshPair<DIM>::ComputeCoarseElementsForFineElementCentroids(bool safeMode)
00414 {
00415     if(mpCoarseMeshBoxCollection==NULL)
00416     {
00417         EXCEPTION("Call SetUpBoxesOnCoarseMesh() before ComputeCoarseElementsForFineElementCentroids()");
00418     }    
00419 
00420     #ifdef FINECOARSEMESHPAIR_VERBOSE
00421     std::cout << "\nComputing coarse elements for fine element centroids\n";
00422     #endif
00423     
00424     mCoarseElementsForFineElementCentroids.resize(mrFineMesh.GetNumElements());
00425     
00426     ResetStatisticsVariables();
00427     for(unsigned i=0; i<mrFineMesh.GetNumElements(); i++)
00428     {
00429         #ifdef FINECOARSEMESHPAIR_VERBOSE
00430         std::cout << "\r " << i << " of " << mrFineMesh.GetNumElements() << std::flush;
00431         #endif
00432         
00433         c_vector<double,DIM> point_cvec = mrFineMesh.GetElement(i)->CalculateCentroid();
00434         ChastePoint<DIM> point(point_cvec);
00435 
00436         // get the box this point is in
00437         unsigned box_for_this_point = mpCoarseMeshBoxCollection->CalculateContainingBox( point_cvec );
00438         
00439         mCoarseElementsForFineElementCentroids[i] = ComputeCoarseElementForGivenPoint(point, safeMode, box_for_this_point);
00440     }
00441 }
00442 
00443 
00448 template<unsigned DIM>
00449 unsigned FineCoarseMeshPair<DIM>::ComputeCoarseElementForGivenPoint(ChastePoint<DIM>& rPoint, 
00450                                                                     bool safeMode,
00451                                                                     unsigned boxForThisPoint)
00452 {    
00453     std::set<unsigned> test_element_indices;
00454     CollectElementsInContainingBox(mpCoarseMeshBoxCollection, boxForThisPoint, test_element_indices);
00455 
00456     unsigned elem_index;
00457 
00458     try
00459     {
00460         elem_index = mrCoarseMesh.GetContainingElementIndex(rPoint,
00461                                                             false,
00462                                                             test_element_indices,
00463                                                             true /* quit if not in test_elements */);
00464 
00465         mStatisticsCounters[0]++;
00466     }
00467     catch(Exception& e)
00468     {
00469         // now try all the elements, but trying the elements contained in the boxes local to this
00470         // element first
00471         std::set<unsigned> test_element_indices;
00472         CollectElementsInLocalBoxes(mpCoarseMeshBoxCollection, boxForThisPoint, test_element_indices);
00473 
00474         try
00475         {
00476             elem_index = mrCoarseMesh.GetContainingElementIndex(rPoint,
00477                                                                 false,
00478                                                                 test_element_indices,
00479                                                                 true);
00480             
00481             mStatisticsCounters[0]++;
00482         }
00483         catch(Exception& e)
00484         {
00485             if(safeMode)
00486             {
00487                 // try the remaining elements
00488                 try
00489                 {
00490                     elem_index = mrCoarseMesh.GetContainingElementIndex(rPoint,
00491                                                                         false);
00492 
00493                     mStatisticsCounters[0]++;
00494                 }
00495                 catch (Exception& e)
00496                 {
00497                     // the point is not in ANY element, store the nearest element and corresponding weights
00498                     elem_index = mrCoarseMesh.GetNearestElementIndexFromTestElements(rPoint,test_element_indices);
00499                     mStatisticsCounters[1]++;
00500                 }
00501             }
00502             else
00503             {
00504                 assert(test_element_indices.size()>0); // boxes probably too small if this fails
00505                 
00506                 // immediately assume it isn't in the rest of the mesh - this should be the 
00507                 // case assuming the box width was chosen suitably.
00508                 // Store the nearest element and corresponding weights
00509                 elem_index = mrCoarseMesh.GetNearestElementIndexFromTestElements(rPoint,test_element_indices);
00510                 mStatisticsCounters[1]++;
00511             }
00512         }
00513     }
00514 
00515     return elem_index;
00516 }
00517 
00518 
00519 
00520 
00522 // 
00523 // helper methods for code
00524 //
00526 
00527 template<unsigned DIM>
00528 void FineCoarseMeshPair<DIM>::CollectElementsInContainingBox(BoxCollection<DIM>*& rpBoxCollection, 
00529                                                              unsigned boxIndex,
00530                                                              std::set<unsigned>& rElementIndices)
00531 {
00532     for(typename std::set<Element<DIM,DIM>*>::iterator elem_iter = rpBoxCollection->rGetBox(boxIndex).rGetElementsContained().begin();
00533         elem_iter != rpBoxCollection->rGetBox(boxIndex).rGetElementsContained().end();
00534         ++elem_iter)
00535     {
00536         rElementIndices.insert((*elem_iter)->GetIndex());
00537     }
00538 }
00539 
00540 
00541 template<unsigned DIM>
00542 void FineCoarseMeshPair<DIM>::CollectElementsInLocalBoxes(BoxCollection<DIM>*& rpBoxCollection, 
00543                                                           unsigned boxIndex,
00544                                                           std::set<unsigned>& rElementIndices)
00545 {
00546     std::set<unsigned> local_boxes = rpBoxCollection->GetLocalBoxes(boxIndex);
00547     for(std::set<unsigned>::iterator local_box_iter = local_boxes.begin();
00548         local_box_iter != local_boxes.end();
00549          ++local_box_iter)
00550     {
00551         for(typename std::set<Element<DIM,DIM>*>::iterator elem_iter = rpBoxCollection->rGetBox(*local_box_iter).rGetElementsContained().begin();
00552             elem_iter != rpBoxCollection->rGetBox(*local_box_iter).rGetElementsContained().end();
00553             ++elem_iter)
00554         {
00555             rElementIndices.insert((*elem_iter)->GetIndex());
00556         }
00557     }
00558 }
00559 
00560 
00562 // 
00563 // statistics related methods
00564 //
00566 
00567 
00568 template<unsigned DIM> 
00569 void FineCoarseMeshPair<DIM>::ResetStatisticsVariables()
00570 {
00571     mNotInMesh.clear();
00572     mNotInMeshNearestElementWeights.clear();
00573     mStatisticsCounters.resize(2, 0u);
00574 }
00575 
00576 template<unsigned DIM>
00577 void FineCoarseMeshPair<DIM>::PrintStatistics()
00578 {
00579     std::cout << "\nFineCoarseMeshPair statistics for the last-called method:\n";
00580 
00581 //    std::cout << "\tNum points for which containing element was found, using box containing that point = " << mStatisticsCounters[0] << "\n";
00582 //    std::cout << "\tNum points for which containing element was in local box = " << mStatisticsCounters[1] << "\n";
00583 //    std::cout << "\tNum points for which containing element was in an element in a non-local box = " << mStatisticsCounters[2] << "\n";
00584 //    std::cout << "\tNum points for which no containing element was found = " << mStatisticsCounters[3] << "\n";
00585 
00586     std::cout << "\tNum points for which containing element was found: " << mStatisticsCounters[0] << "\n";
00587     std::cout << "\tNum points for which no containing element was found = " << mStatisticsCounters[1] << "\n";
00588 
00589     if(mNotInMesh.size()>0)
00590     {
00591         std::cout << "\tIndices and weights for points (nodes/quad points) for which no containing element was found:\n";
00592         for(unsigned i=0; i<mNotInMesh.size(); i++)
00593         {
00594             std::cout << "\t\t" << mNotInMesh[i] << ", " << mNotInMeshNearestElementWeights[i] << "\n";
00595         }
00596     }
00597 }
00598 
00600 // Explicit instantiation
00602 
00603 
00604 template class FineCoarseMeshPair<1>;
00605 template class FineCoarseMeshPair<2>;
00606 template class FineCoarseMeshPair<3>;

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