1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 
set_p_level(unsigned int p)41 void Elem::set_p_level(unsigned int p)
42 {
43   // Maintain the parent's p level as the minimum of it's children
44   if (this->parent() != nullptr)
45     {
46       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       if (parent_p_level > p)
50         {
51           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           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       else if (parent_p_level == _p_level && _p_level < p)
60         {
61           _p_level = cast_int<unsigned char>(p);
62           parent_p_level = cast_int<unsigned char>(p);
63           for (auto & c : this->parent()->child_ref_range())
64             parent_p_level = std::min(parent_p_level,
65                                       c.p_level());
66 
67           // When its children all have a higher p level, the parent's
68           // should rise
69           if (parent_p_level > this->parent()->p_level())
70             {
71               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               this->parent()->set_p_refinement_flag(Elem::JUST_REFINED);
76             }
77 
78           return;
79         }
80     }
81 
82   _p_level = cast_int<unsigned char>(p);
83 }
84 
85 
86 
refine(MeshRefinement & mesh_refinement)87 void Elem::refine (MeshRefinement & mesh_refinement)
88 {
89   libmesh_assert_equal_to (this->refinement_flag(), Elem::REFINE);
90   libmesh_assert (this->active());
91 
92   const unsigned int nc = this->n_children();
93 
94   // Create my children if necessary
95   if (!_children)
96     {
97       _children = new Elem *[nc];
98 
99       unsigned int parent_p_level = this->p_level();
100       const unsigned int nei = this->n_extra_integers();
101       for (unsigned int c = 0; c != nc; c++)
102         {
103           auto current_child = Elem::build(this->type(), this);
104           _children[c] = current_child.get();
105 
106           current_child->set_refinement_flag(Elem::JUST_REFINED);
107           current_child->set_p_level(parent_p_level);
108           current_child->set_p_refinement_flag(this->p_refinement_flag());
109 
110           for (auto cnode : current_child->node_index_range())
111             {
112               Node * node =
113                 mesh_refinement.add_node(*this, c, cnode,
114                                          current_child->processor_id());
115               node->set_n_systems (this->n_systems());
116               current_child->set_node(cnode) = node;
117             }
118 
119           Elem * added_child = mesh_refinement.add_elem (std::move(current_child));
120           added_child->set_n_systems(this->n_systems());
121           libmesh_assert_equal_to (added_child->n_extra_integers(),
122                                    this->n_extra_integers());
123           for (unsigned int i=0; i != nei; ++i)
124             added_child->set_extra_integer(i, this->get_extra_integer(i));
125         }
126     }
127   else
128     {
129       unsigned int parent_p_level = this->p_level();
130       for (unsigned int c = 0; c != nc; c++)
131         {
132           Elem * current_child = this->child_ptr(c);
133           if (current_child != remote_elem)
134             {
135               libmesh_assert(current_child->subactive());
136               current_child->set_refinement_flag(Elem::JUST_REFINED);
137               current_child->set_p_level(parent_p_level);
138               current_child->set_p_refinement_flag(this->p_refinement_flag());
139             }
140         }
141     }
142 
143   // Un-set my refinement flag now
144   this->set_refinement_flag(Elem::INACTIVE);
145 
146   // Leave the p refinement flag set - we will need that later to get
147   // projection operations correct
148   // this->set_p_refinement_flag(Elem::INACTIVE);
149 
150 #ifndef NDEBUG
151   for (unsigned int c = 0; c != nc; c++)
152     if (this->child_ptr(c) != remote_elem)
153       {
154         libmesh_assert_equal_to (this->child_ptr(c)->parent(), this);
155         libmesh_assert(this->child_ptr(c)->active());
156       }
157 #endif
158   libmesh_assert (this->ancestor());
159 }
160 
161 
162 
coarsen()163 void Elem::coarsen()
164 {
165   libmesh_assert_equal_to (this->refinement_flag(), Elem::COARSEN_INACTIVE);
166   libmesh_assert (!this->active());
167 
168   // We no longer delete children until MeshRefinement::contract()
169   // delete [] _children;
170   // _children = nullptr;
171 
172   unsigned int parent_p_level = 0;
173 
174   const unsigned int n_n = this->n_nodes();
175 
176   // re-compute hanging node nodal locations
177   for (unsigned int c = 0, nc = this->n_children(); c != nc; ++c)
178     {
179       Elem * mychild = this->child_ptr(c);
180       if (mychild == remote_elem)
181         continue;
182       for (auto cnode : mychild->node_index_range())
183         {
184           Point new_pos;
185           bool calculated_new_pos = false;
186 
187           for (unsigned int n=0; n<n_n; n++)
188             {
189               // The value from the embedding matrix
190               const float em_val = this->embedding_matrix(c,cnode,n);
191 
192               // The node location is somewhere between existing vertices
193               if ((em_val != 0.) && (em_val != 1.))
194                 {
195                   new_pos.add_scaled (this->point(n), em_val);
196                   calculated_new_pos = true;
197                 }
198             }
199 
200           if (calculated_new_pos)
201             {
202               //Move the existing node back into it's original location
203               for (unsigned int i=0; i<LIBMESH_DIM; i++)
204                 {
205                   Point & child_node = mychild->point(cnode);
206                   child_node(i)=new_pos(i);
207                 }
208             }
209         }
210     }
211 
212   for (auto & mychild : this->child_ref_range())
213     {
214       if (&mychild == remote_elem)
215         continue;
216       libmesh_assert_equal_to (mychild.refinement_flag(), Elem::COARSEN);
217       mychild.set_refinement_flag(Elem::INACTIVE);
218       if (mychild.p_level() > parent_p_level)
219         parent_p_level = mychild.p_level();
220     }
221 
222   this->set_refinement_flag(Elem::JUST_COARSENED);
223   this->set_p_level(parent_p_level);
224 
225   libmesh_assert (this->active());
226 }
227 
228 
229 
contract()230 void Elem::contract()
231 {
232   // Subactive elements get deleted entirely, not contracted
233   libmesh_assert (this->active());
234 
235   // Active contracted elements no longer can have children
236   delete [] _children;
237   _children = nullptr;
238 
239   if (this->refinement_flag() == Elem::JUST_COARSENED)
240     this->set_refinement_flag(Elem::DO_NOTHING);
241 }
242 
243 #endif // #ifdef LIBMESH_ENABLE_AMR
244 
245 
246 } // namespace libMesh
247