Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2026 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             : 
      20             : // C++ includes
      21             : 
      22             : // Local includes
      23             : #include "libmesh/elem.h"
      24             : #include "libmesh/mesh_refinement.h"
      25             : #include "libmesh/remote_elem.h"
      26             : 
      27             : namespace libMesh
      28             : {
      29             : 
      30             : 
      31             : //--------------------------------------------------------------------
      32             : // Elem methods
      33             : 
      34             : /**
      35             :  * The following functions only apply when
      36             :  * AMR is enabled and thus are not present
      37             :  * otherwise.
      38             :  */
      39             : #ifdef LIBMESH_ENABLE_AMR
      40             : 
      41   321932567 : void Elem::set_p_level(unsigned int p)
      42             : {
      43             :   // Maintain the parent's p level as the minimum of it's children
      44   338375545 :   if (this->parent() != nullptr)
      45             :     {
      46    65186973 :       unsigned int parent_p_level = this->parent()->p_level();
      47             : 
      48             :       // If our new p level is less than our parents, our parents drops
      49    65186973 :       if (parent_p_level > p)
      50             :         {
      51         242 :           this->parent()->set_p_level(p);
      52             : 
      53             :           // And we should keep track of the drop, in case we need to
      54             :           // do a projection later.
      55          32 :           this->parent()->set_p_refinement_flag(Elem::JUST_COARSENED);
      56             :         }
      57             :       // If we are the lowest p level and it increases, so might
      58             :       // our parent's, but we have to check every other child to see
      59    65186731 :       else if (parent_p_level == _p_level && _p_level < p)
      60             :         {
      61       24145 :           _p_level = cast_int<unsigned char>(p);
      62       24145 :           parent_p_level = cast_int<unsigned char>(p);
      63      176941 :           for (auto & c : this->parent()->child_ref_range())
      64      152796 :             parent_p_level = std::min(parent_p_level,
      65      184838 :                                       c.p_level());
      66             : 
      67             :           // When its children all have a higher p level, the parent's
      68             :           // should rise
      69       27045 :           if (parent_p_level > this->parent()->p_level())
      70             :             {
      71        1809 :               this->parent()->set_p_level(parent_p_level);
      72             : 
      73             :               // And we should keep track of the rise, in case we need to
      74             :               // do a projection later.
      75         382 :               this->parent()->set_p_refinement_flag(Elem::JUST_REFINED);
      76             :             }
      77             : 
      78       24145 :           return;
      79             :         }
      80             :     }
      81             : 
      82    48233544 :   this->hack_p_level(p);
      83             : }
      84             : 
      85             : 
      86             : 
      87     1248720 : void Elem::refine (MeshRefinement & mesh_refinement)
      88             : {
      89       74944 :   libmesh_assert_equal_to (this->refinement_flag(), Elem::REFINE);
      90       74944 :   libmesh_assert (this->active());
      91             : 
      92     1248720 :   const unsigned int nc = this->n_children();
      93             : 
      94             :   // Create my children if necessary
      95     1248720 :   if (!_children)
      96             :     {
      97     2315626 :       _children = std::make_unique<Elem *[]>(nc);
      98             : 
      99      148730 :       unsigned int parent_p_level = this->p_level();
     100     1232178 :       const unsigned int nei = this->n_extra_integers();
     101     7766255 :       for (unsigned int c = 0; c != nc; c++)
     102             :         {
     103     6923863 :           auto current_child = Elem::build(this->type(), this);
     104     6923863 :           _children[c] = current_child.get();
     105             : 
     106      389786 :           current_child->set_refinement_flag(Elem::JUST_REFINED);
     107     6534077 :           current_child->set_p_level(parent_p_level);
     108      779572 :           current_child->set_p_refinement_flag(this->p_refinement_flag());
     109             : 
     110    53981206 :           for (auto cnode : current_child->node_index_range())
     111             :             {
     112             :               Node * node =
     113    50118409 :                 mesh_refinement.add_node(*this, c, cnode,
     114     2671280 :                                          current_child->processor_id());
     115    59735375 :               node->set_n_systems (this->n_systems());
     116    47447129 :               current_child->set_node(cnode, node);
     117             :             }
     118             : 
     119     6534077 :           Elem * added_child = mesh_refinement.add_elem (std::move(current_child));
     120     8409203 :           added_child->set_n_systems(this->n_systems());
     121      389786 :           libmesh_assert_equal_to (added_child->n_extra_integers(),
     122             :                                    this->n_extra_integers());
     123     6547919 :           for (unsigned int i=0; i != nei; ++i)
     124       12734 :             added_child->set_extra_integer(i, this->get_extra_integer(i));
     125     5754505 :         }
     126             :     }
     127             :   else
     128             :     {
     129        1156 :       unsigned int parent_p_level = this->p_level();
     130       82710 :       for (unsigned int c = 0; c != nc; c++)
     131             :         {
     132        2312 :           Elem * current_child = this->child_ptr(c);
     133       66168 :           if (current_child != remote_elem)
     134             :             {
     135        2240 :               libmesh_assert(current_child->subactive());
     136        2240 :               current_child->set_refinement_flag(Elem::JUST_REFINED);
     137       62121 :               current_child->set_p_level(parent_p_level);
     138        4480 :               current_child->set_p_refinement_flag(this->p_refinement_flag());
     139             :             }
     140             :         }
     141             :     }
     142             : 
     143             :   // Get any new remote neighbor links right; even find_neighbors
     144             :   // relies on those
     145     6040941 :   for (unsigned int s : this->side_index_range())
     146             :     {
     147     5088791 :       if (this->neighbor_ptr(s) != remote_elem)
     148     4334695 :         continue;
     149             : 
     150     1175312 :       for (unsigned int c = 0; c != nc; c++)
     151             :         {
     152        1036 :           Elem * current_child = this->child_ptr(c);
     153     2024241 :           if (current_child != remote_elem &&
     154     1010067 :               this->is_child_on_side(c, s))
     155             :             current_child->set_neighbor
     156      504077 :               (s, const_cast<RemoteElem *>(remote_elem));
     157             :         }
     158             :     }
     159             : 
     160             :   // Un-set my refinement flag now
     161       74944 :   this->set_refinement_flag(Elem::INACTIVE);
     162             : 
     163             :   // Leave the p refinement flag set - we will need that later to get
     164             :   // projection operations correct
     165             :   // this->set_p_refinement_flag(Elem::INACTIVE);
     166             : 
     167             : #ifndef NDEBUG
     168      467042 :   for (unsigned int c = 0; c != nc; c++)
     169      392098 :     if (this->child_ptr(c) != remote_elem)
     170             :       {
     171      392026 :         libmesh_assert_equal_to (this->child_ptr(c)->parent(), this);
     172      392026 :         libmesh_assert(this->child_ptr(c)->active());
     173             :       }
     174             : #endif
     175       74944 :   libmesh_assert (this->ancestor());
     176     1248720 : }
     177             : 
     178             : 
     179             : 
     180      396682 : void Elem::coarsen()
     181             : {
     182       18476 :   libmesh_assert_equal_to (this->refinement_flag(), Elem::COARSEN_INACTIVE);
     183       18476 :   libmesh_assert (!this->active());
     184             : 
     185             :   // We no longer delete children until MeshRefinement::contract()
     186             : 
     187       18476 :   unsigned int parent_p_level = 0;
     188             : 
     189      396682 :   const unsigned int n_n = this->n_nodes();
     190             : 
     191             :   // re-compute hanging node nodal locations
     192     1986142 :   for (unsigned int c = 0, nc = this->n_children(); c != nc; ++c)
     193             :     {
     194       74096 :       Elem * mychild = this->child_ptr(c);
     195     1589460 :       if (mychild == remote_elem)
     196       29746 :         continue;
     197     8943873 :       for (auto cnode : mychild->node_index_range())
     198             :         {
     199      420936 :           Point new_pos;
     200      420936 :           bool calculated_new_pos = false;
     201             : 
     202    51395254 :           for (unsigned int n=0; n<n_n; n++)
     203             :             {
     204             :               // The value from the embedding matrix
     205    44011023 :               const Real em_val = this->embedding_matrix(c,cnode,n);
     206             : 
     207             :               // The node location is somewhere between existing vertices
     208    44011023 :               if ((em_val != 0.) && (em_val != 1.))
     209             :                 {
     210    14597961 :                   new_pos.add_scaled (this->point(n), em_val);
     211      894600 :                   calculated_new_pos = true;
     212             :                 }
     213             :             }
     214             : 
     215     7384231 :           if (calculated_new_pos)
     216             :             {
     217             :               //Move the existing node back into it's original location
     218    19038800 :               for (unsigned int i=0; i<LIBMESH_DIM; i++)
     219             :                 {
     220     1571904 :                   Point & child_node = mychild->point(cnode);
     221    14279100 :                   child_node(i)=new_pos(i);
     222             :                 }
     223             :             }
     224             :         }
     225             :     }
     226             : 
     227     1986142 :   for (auto & mychild : this->child_ref_range())
     228             :     {
     229     1589460 :       if (&mychild == remote_elem)
     230       29746 :         continue;
     231       74024 :       libmesh_assert_equal_to (mychild.refinement_flag(), Elem::COARSEN);
     232       74024 :       mychild.set_refinement_flag(Elem::INACTIVE);
     233     1633666 :       if (mychild.p_level() > parent_p_level)
     234           0 :         parent_p_level = mychild.p_level();
     235             :     }
     236             : 
     237       18476 :   this->set_refinement_flag(Elem::JUST_COARSENED);
     238      396682 :   this->set_p_level(parent_p_level);
     239             : 
     240       18476 :   libmesh_assert (this->active());
     241      396682 : }
     242             : 
     243             : 
     244             : 
     245     6776970 : void Elem::contract()
     246             : {
     247             :   // Subactive elements get deleted entirely, not contracted
     248      359394 :   libmesh_assert (this->active());
     249             : 
     250             :   // Active contracted elements no longer can have children
     251      359394 :   _children.reset(nullptr);
     252             : 
     253     6776970 :   if (this->refinement_flag() == Elem::JUST_COARSENED)
     254       17896 :     this->set_refinement_flag(Elem::DO_NOTHING);
     255     6776970 : }
     256             : 
     257             : #endif // #ifdef LIBMESH_ENABLE_AMR
     258             : 
     259             : 
     260             : } // namespace libMesh