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