Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : // C++ includes
      20             : #include <algorithm> // for std::fill
      21             : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
      22             : #include <cmath>     // for std::sqrt std::pow std::abs
      23             : 
      24             : 
      25             : // Local Includes
      26             : #include "libmesh/libmesh_common.h"
      27             : #include "libmesh/patch.h"
      28             : #include "libmesh/elem.h"
      29             : 
      30             : namespace libMesh
      31             : {
      32             : 
      33             : 
      34             : 
      35             : //-----------------------------------------------------------------
      36             : // Patch implementations
      37      301353 : void Patch::find_face_neighbors(std::set<const Elem *> & new_neighbors)
      38             : {
      39             :   // Loop over all the elements in the patch
      40     2203521 :   for (const auto & elem : *this)
      41     9800982 :       for (auto neighbor : elem->neighbor_ptr_range())
      42     7608672 :         if (neighbor != nullptr)        // we have a neighbor on this side
      43             :           {
      44             : #ifdef LIBMESH_ENABLE_AMR
      45      875293 :             if (!neighbor->active())          // the neighbor is *not* active,
      46             :               {                               // so add *all* neighboring
      47             :                                               // active children to the patch
      48      135400 :                 std::vector<const Elem *> active_neighbor_children;
      49             : 
      50             :                 neighbor->active_family_tree_by_neighbor
      51      427724 :                   (active_neighbor_children, elem);
      52             : 
      53     1283172 :                 for (const auto & child : active_neighbor_children)
      54      720092 :                   new_neighbors.insert(child);
      55             :               }
      56             :             else
      57             : #endif // #ifdef LIBMESH_ENABLE_AMR
      58     4552609 :               new_neighbors.insert (neighbor); // add active neighbors
      59             :           }
      60      301353 : }
      61             : 
      62             : 
      63             : 
      64           0 : void Patch::add_face_neighbors()
      65             : {
      66           0 :   std::set<const Elem *> new_neighbors;
      67             : 
      68           0 :   this->find_face_neighbors(new_neighbors);
      69             : 
      70           0 :   this->insert(new_neighbors.begin(), new_neighbors.end());
      71           0 : }
      72             : 
      73             : 
      74             : 
      75      301353 : void Patch::add_local_face_neighbors()
      76             : {
      77       72272 :   std::set<const Elem *> new_neighbors;
      78             : 
      79      301353 :   this->find_face_neighbors(new_neighbors);
      80             : 
      81             :   // if neighbor belongs to this processor then add it to the patch
      82     3492971 :   for (const auto & neighbor : new_neighbors)
      83     3191618 :     if (neighbor->processor_id() == _my_procid)
      84     2167656 :       this->insert (neighbor);
      85      301353 : }
      86             : 
      87             : 
      88             : 
      89           0 : void Patch::add_semilocal_face_neighbors()
      90             : {
      91           0 :   std::set<const Elem *> new_neighbors;
      92             : 
      93           0 :   this->find_face_neighbors(new_neighbors);
      94             : 
      95           0 :   for (const auto & neighbor : new_neighbors)
      96           0 :     if (neighbor->is_semilocal(_my_procid))
      97           0 :       this->insert (neighbor);
      98           0 : }
      99             : 
     100             : 
     101             : 
     102           0 : void Patch::find_point_neighbors(std::set<const Elem *> & new_neighbors)
     103             : {
     104             :   // Loop over all the elements in the patch
     105           0 :   for (const auto & elem : *this)
     106             :     {
     107           0 :       std::set<const Elem *> elem_point_neighbors;
     108             : 
     109           0 :       elem->find_point_neighbors(elem_point_neighbors);
     110             : 
     111           0 :       new_neighbors.insert(elem_point_neighbors.begin(),
     112             :                            elem_point_neighbors.end());
     113             :     }
     114           0 : }
     115             : 
     116             : 
     117             : 
     118           0 : void Patch::add_point_neighbors()
     119             : {
     120           0 :   std::set<const Elem *> new_neighbors;
     121             : 
     122           0 :   this->find_point_neighbors(new_neighbors);
     123             : 
     124           0 :   this->insert(new_neighbors.begin(), new_neighbors.end());
     125           0 : }
     126             : 
     127             : 
     128             : 
     129           0 : void Patch::add_local_point_neighbors()
     130             : {
     131           0 :   std::set<const Elem *> new_neighbors;
     132             : 
     133           0 :   this->find_point_neighbors(new_neighbors);
     134             : 
     135             :   // if neighbor belongs to this processor then add it to the patch
     136           0 :   for (const auto & neighbor : new_neighbors)
     137           0 :     if (neighbor->processor_id() == _my_procid)
     138           0 :       this->insert (neighbor);
     139           0 : }
     140             : 
     141             : 
     142             : 
     143           0 : void Patch::add_semilocal_point_neighbors()
     144             : {
     145           0 :   std::set<const Elem *> new_neighbors;
     146             : 
     147           0 :   this->find_point_neighbors(new_neighbors);
     148             : 
     149           0 :   for (const auto & neighbor : new_neighbors)
     150           0 :     if (neighbor->is_semilocal(_my_procid))
     151           0 :       this->insert (neighbor);
     152           0 : }
     153             : 
     154             : 
     155             : 
     156      194067 : void Patch::build_around_element (const Elem * e0,
     157             :                                   const unsigned int target_patch_size,
     158             :                                   PMF patchtype)
     159             : {
     160             : 
     161             :   // Make sure we are building a patch for an active element.
     162       17403 :   libmesh_assert(e0);
     163       17403 :   libmesh_assert (e0->active());
     164             :   // Make sure we are either starting with a local element or
     165             :   // requesting a nonlocal patch
     166       17403 :   libmesh_assert ((patchtype != &Patch::add_local_face_neighbors &&
     167             :                    patchtype != &Patch::add_local_point_neighbors) ||
     168             :                   e0->processor_id() == _my_procid);
     169             : 
     170             :   // First clear the current set, then add the element of interest.
     171       17403 :   this->clear();
     172      176668 :   this->insert (e0);
     173             : 
     174             :   // Repeatedly add the neighbors of the elements in the patch until
     175             :   // the target patch size is met
     176      445855 :   while (this->size() < target_patch_size)
     177             :     {
     178             :       // It is possible that the target patch size is larger than the number
     179             :       // of elements that can be added to the patch.  Since we don't
     180             :       // have access to the Mesh object here, the only way we can
     181             :       // detect this case is by detecting a "stagnant patch," i.e. a
     182             :       // patch whose size does not increase after adding face neighbors
     183       36136 :       const std::size_t old_patch_size = this->size();
     184             : 
     185             :       // We profile the patch-extending functions separately
     186      301353 :       (this->*patchtype)();
     187             : 
     188             :       // Check for a "stagnant" patch
     189      301353 :       if (this->size() == old_patch_size)
     190             :         {
     191       49827 :           libmesh_do_once(libMesh::err <<
     192             :                           "WARNING: stagnant patch of " << this->size() << " elements."
     193             :                           << std::endl <<
     194             :                           "Does the target patch size exceed the number of local elements?"
     195             :                           << std::endl;
     196             :                           libmesh_here(););
     197        2418 :           break;
     198             :         }
     199             :     } // end while loop
     200             : 
     201             : 
     202             :   // make sure all the elements in the patch are active and local
     203             :   // if we are in debug mode
     204             : #ifdef DEBUG
     205      145139 :     for (const auto & elem : *this)
     206             :       {
     207      127736 :         libmesh_assert (elem->active());
     208      127736 :         if ((patchtype == &Patch::add_local_face_neighbors ||
     209           0 :              patchtype == &Patch::add_local_point_neighbors))
     210      127736 :           libmesh_assert_equal_to (elem->processor_id(), _my_procid);
     211             :       }
     212             : #endif
     213             : 
     214      194067 : }
     215             : 
     216             : } // namespace libMesh