refineable_brick_spectral_element.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#ifndef OOMPH_REFINEABLE_BRICK_SPECTRAL_ELEMENT_HEADER
27#define OOMPH_REFINEABLE_BRICK_SPECTRAL_ELEMENT_HEADER
28
29
30// Config header generated by autoconfig
31#ifdef HAVE_CONFIG_H
32#include <oomph-lib-config.h>
33#endif
34
35
36// oomph-lib headers
38
39namespace oomph
40{
41 //=======================================================================
42 /// Refineable version of QuadElements that add functionality for spectral
43 /// Elements.
44 //=======================================================================
45 template<>
47 {
48 public:
49 /// Constructor
51 {
52#ifdef LEAK_CHECK
53 LeakCheckNames::RefineableQSpectralElement<3> _build += 1;
54#endif
55 }
56
57
58 /// Broken copy constructor
60 delete;
61
62 /// Broken assignment operator
63 // Commented out broken assignment operator because this can lead to a
64 // conflict warning when used in the virtual inheritence hierarchy.
65 // Essentially the compiler doesn't realise that two separate
66 // implementations of the broken function are the same and so, quite
67 // rightly, it shouts.
68 /*void operator=(const RefineableQSpectralElement<3>&) = delete;*/
69
70 /// Destructor
72 {
73#ifdef LEAK_CHECK
74 LeakCheckNames::RefineableQSpectralElement<3> _build -= 1;
75#endif
76 }
77
78 /// The only thing to add is rebuild from sons
79 void rebuild_from_sons(Mesh*& mesh_pt)
80 {
81 // The timestepper should be the same for all nodes and node 0 should
82 // never be deleted.
83 if (this->node_pt(0) == 0)
84 {
85 throw OomphLibError("The Corner node (0) does not exist",
86 OOMPH_CURRENT_FUNCTION,
87 OOMPH_EXCEPTION_LOCATION);
88 }
89
90 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
91 unsigned ntstorage = time_stepper_pt->ntstorage();
92
93 unsigned jnod = 0;
94 Vector<double> s_fraction(3), s(3);
95 // Loop over the nodes in the element
96 unsigned n_p = this->nnode_1d();
97 for (unsigned i0 = 0; i0 < n_p; i0++)
98 {
99 // Get the fractional position of the node
100 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
101 // Local coordinate
102 s[0] = -1.0 + 2.0 * s_fraction[0];
103
104 for (unsigned i1 = 0; i1 < n_p; i1++)
105 {
106 // Get the fractional position of the node in the direction of s[1]
107 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
108 // Local coordinate in father element
109 s[1] = -1.0 + 2.0 * s_fraction[1];
110
111 for (unsigned i2 = 0; i2 < n_p; i2++)
112 {
113 // Get the fractional position of the node in the direction of s[1]
114 s_fraction[2] = this->local_one_d_fraction_of_node(i2, 2);
115 // Local coordinate in father element
116 s[2] = -1.0 + 2.0 * s_fraction[2];
117
118 // Set the local node number
119 jnod = i0 + n_p * i1 + n_p * n_p * i2;
120
121 // If the node has not been built
122 if (this->node_pt(jnod) == 0)
123 {
124 // Has the node been created by one of its neighbours
125 bool is_periodic = false;
126 Node* created_node_pt =
127 this->node_created_by_neighbour(s_fraction, is_periodic);
128
129 // If it has set the pointer
130 if (created_node_pt != 0)
131 {
132 // If the node is periodic
133 if (is_periodic)
134 {
135 throw OomphLibError("Cannot handle periodic nodes in "
136 "refineable spectral elements yet",
137 OOMPH_CURRENT_FUNCTION,
138 OOMPH_EXCEPTION_LOCATION);
139 }
140 // Non-periodic case, just set the pointer
141 else
142 {
143 this->node_pt(jnod) = created_node_pt;
144 }
145 }
146 // Otherwise, we need to build it
147 else
148 {
149 // First we need to find the pointer to the son that
150 // would contain the node
151
152 // Find coordinates in the sons
153 Vector<double> s_in_son(3);
154 using namespace OcTreeNames;
155 int son = -10;
156 // If less than one-half on the left side
157 if (s_fraction[0] < 0.5)
158 {
159 // On the down side
160 if (s_fraction[1] < 0.5)
161 {
162 // On the back side
163 if (s_fraction[2] < 0.5)
164 {
165 // It's the left down back son
166 son = LDB;
167 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
168 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
169 s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
170 }
171 // On the front side
172 else
173 {
174 // It's the left down front son
175 son = LDF;
176 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
177 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
178 s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
179 }
180 } // End of on down side
181 // else its on the top
182 else
183 {
184 // On the back
185 if (s_fraction[2] < 0.5)
186 {
187 // It's the left up back son
188 son = LUB;
189 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
190 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
191 s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
192 }
193 // On the front side
194 else
195 {
196 // It's the left up front son
197 son = LUF;
198 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
199 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
200 s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
201 }
202 } // End of on top
203 } // End of on left
204 // Otherwise its on the right
205 else
206 {
207 // On the down side
208 if (s_fraction[1] < 0.5)
209 {
210 // On the back side
211 if (s_fraction[2] < 0.5)
212 {
213 // It's the right down back son
214 son = RDB;
215 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
216 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
217 s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
218 }
219 // On the front side
220 else
221 {
222 // It's the right down front son
223 son = RDF;
224 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
225 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
226 s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
227 }
228 } // End of on down side
229 // else its on the top
230 else
231 {
232 // On the back
233 if (s_fraction[2] < 0.5)
234 {
235 // It's the right up back son
236 son = RUB;
237 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
238 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
239 s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
240 }
241 // On the front side
242 else
243 {
244 // It's the right up front son
245 son = RUF;
246 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
247 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
248 s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
249 }
250 }
251 }
252
253 // Get the pointer to the son element
255 dynamic_cast<RefineableQSpectralElement<3>*>(
256 this->tree_pt()->son_pt(son)->object_pt());
257
258 // If we are rebuilding, then worry about the boundary
259 // conditions Find the boundary of the node Initially none
260 int boundary = Tree::OMEGA;
261 // If we are on the left boundary
262 if (i0 == 0)
263 {
264 boundary = L;
265 }
266 // If we are on the right boundary
267 else if (i0 == n_p - 1)
268 {
269 boundary = R;
270 }
271
272 // If we are on the lower boundary
273 if (i1 == 0)
274 {
275 // If we already have already set the boundary, we're on an
276 // edge
277 switch (boundary)
278 {
279 case L:
280 boundary = LD;
281 break;
282 case R:
283 boundary = RD;
284 break;
285 // Boundary not set
286 default:
287 boundary = D;
288 break;
289 }
290 }
291 // If we are the northern bounadry
292 else if (i1 == n_p - 1)
293 {
294 // If we already have a boundary
295 switch (boundary)
296 {
297 case L:
298 boundary = LU;
299 break;
300 case R:
301 boundary = RU;
302 break;
303 default:
304 boundary = U;
305 break;
306 }
307 }
308
309 // If we are on the back face
310 if (i2 == 0)
311 {
312 // If we already have boundaries
313 switch (boundary)
314 {
315 case L:
316 boundary = LB;
317 break;
318 case R:
319 boundary = RB;
320 break;
321 case U:
322 boundary = UB;
323 break;
324 case D:
325 boundary = DB;
326 break;
327 case LD:
328 boundary = LDB;
329 break;
330 case RD:
331 boundary = RDB;
332 break;
333 case LU:
334 boundary = LUB;
335 break;
336 case RU:
337 boundary = RUB;
338 break;
339 default:
340 boundary = B;
341 break;
342 }
343 }
344 else if (i2 == n_p - 1)
345 {
346 // If we already have boundaries
347 switch (boundary)
348 {
349 case L:
350 boundary = LF;
351 break;
352 case R:
353 boundary = RF;
354 break;
355 case U:
356 boundary = UF;
357 break;
358 case D:
359 boundary = DF;
360 break;
361 case LD:
362 boundary = LDF;
363 break;
364 case RD:
365 boundary = RDF;
366 break;
367 case LU:
368 boundary = LUF;
369 break;
370 case RU:
371 boundary = RUF;
372 break;
373 default:
374 boundary = F;
375 break;
376 }
377 }
378
379 // set of boundaries that this edge in the son lives on
380 std::set<unsigned> boundaries;
381 // If we are on a boundary of the Element, find the
382 // mesh boundaries on which we live
383 // The boundaries will be common to the son because there can be
384 // no rotations here
385 if (boundary != Tree::OMEGA)
386 {
387 son_el_pt->get_boundaries(boundary, boundaries);
388 }
389
390 // If the node lives on a boundary:
391 // construct a boundary node,
392 // get boundary conditions and
393 // update both lookup schemes
394 if (boundaries.size() > 0)
395 {
396 // Construct the new node
397 this->node_pt(jnod) =
398 this->construct_boundary_node(jnod, time_stepper_pt);
399
400 // Get the boundary conditions from the son
401 Vector<int> bound_cons(ncont_interpolated_values());
402 son_el_pt->get_bcs(boundary, bound_cons);
403
404 // Loop over the values and pin if necessary
405 unsigned nval = this->node_pt(jnod)->nvalue();
406 for (unsigned k = 0; k < nval; k++)
407 {
408 if (bound_cons[k])
409 {
410 this->node_pt(jnod)->pin(k);
411 }
412 }
413
414 // Solid node? If so, deal with the positional boundary
415 // conditions:
416 SolidNode* solid_node_pt =
417 dynamic_cast<SolidNode*>(this->node_pt(jnod));
418 if (solid_node_pt != 0)
419 {
420 // Get the positional boundary conditions from the father:
421 unsigned n_dim = this->node_pt(jnod)->ndim();
422 Vector<int> solid_bound_cons(n_dim);
423 RefineableSolidQElement<3>* son_solid_el_pt =
424 dynamic_cast<RefineableSolidQElement<3>*>(son_el_pt);
425#ifdef PARANOID
426 if (son_solid_el_pt == 0)
427 {
428 std::string error_message = "We have a SolidNode outside "
429 "a refineable SolidElement\n";
430 error_message +=
431 "during mesh refinement -- this doesn't make sense\n";
432
433 throw OomphLibError(error_message,
434 OOMPH_CURRENT_FUNCTION,
435 OOMPH_EXCEPTION_LOCATION);
436 }
437#endif
438 son_solid_el_pt->get_solid_bcs(boundary, solid_bound_cons);
439
440 // Loop over the positions and pin, if necessary
441 for (unsigned k = 0; k < n_dim; k++)
442 {
443 if (solid_bound_cons[k])
444 {
445 solid_node_pt->pin_position(k);
446 }
447 }
448 }
449
450 // Next we update the boundary look-up schemes
451 // Loop over the boundaries stored in the set
452 for (std::set<unsigned>::iterator it = boundaries.begin();
453 it != boundaries.end();
454 ++it)
455 {
456 // Add the node to the boundary
457 mesh_pt->add_boundary_node(*it, this->node_pt(jnod));
458
459 // If we have set an intrinsic coordinate on this
460 // mesh boundary then it must also be interpolated on
461 // the new node
462 // Now interpolate the intrinsic boundary coordinate
463 if (mesh_pt->boundary_coordinate_exists(*it) == true)
464 {
465 Vector<double> zeta(2);
466 son_el_pt->interpolated_zeta_on_face(
467 *it, boundary, s_in_son, zeta);
468
469 this->node_pt(jnod)->set_coordinates_on_boundary(*it,
470 zeta);
471 }
472 }
473 }
474 // Otherwise we construct a normal "bulk" node
475 else
476 {
477 // Construct the new node
478 this->node_pt(jnod) =
479 this->construct_node(jnod, time_stepper_pt);
480 }
481
482 // Now we set the position and values at the newly created node
483
484 // In the first instance use macro element or FE representation
485 // to create past and present nodal positions.
486 // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
487 // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
488 // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
489 // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
490 // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
491 // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
492 // NOT ASSIGN SENSIBLE INITIAL POSITONS!
493
494 // Loop over # of history values
495 for (unsigned t = 0; t < ntstorage; t++)
496 {
497 // Get the position from the son
498 Vector<double> x_prev(3);
499
500 // Now let's fill in the value
501 son_el_pt->get_x(t, s_in_son, x_prev);
502 for (unsigned i = 0; i < 3; i++)
503 {
504 this->node_pt(jnod)->x(t, i) = x_prev[i];
505 }
506 }
507
508 // Now set the values
509 // Loop over all history values
510 for (unsigned t = 0; t < ntstorage; t++)
511 {
512 // Get values from father element
513 // Note: get_interpolated_values() sets Vector size itself.
514 Vector<double> prev_values;
515 son_el_pt->get_interpolated_values(t, s_in_son, prev_values);
516
517 // Initialise the values at the new node
518 for (unsigned k = 0; k < this->node_pt(jnod)->nvalue(); k++)
519 {
520 this->node_pt(jnod)->set_value(t, k, prev_values[k]);
521 }
522 }
523
524 // Add the node to the mesh
525 mesh_pt->add_node_pt(this->node_pt(jnod));
526
527 } // Node has been constructed
528
529 // Algebraic stuff here
530 // Check whether the element is an algebraic element
531 AlgebraicElementBase* alg_el_pt =
532 dynamic_cast<AlgebraicElementBase*>(this);
533
534 // If we do have an algebraic element
535 if (alg_el_pt != 0)
536 {
537 std::string error_message =
538 "Have not implemented rebuilding from sons for";
539 error_message += "Algebraic Spectral elements yet\n";
540
541 throw OomphLibError(
542 error_message,
543 "RefineableQSpectralElement::rebuild_from_sons()",
544 OOMPH_EXCEPTION_LOCATION);
545 }
546
547 } // End of the case when the node was not built
548 }
549 }
550 }
551 }
552
553 /// Overload the nodes built function
554 virtual bool nodes_built()
555 {
556 unsigned n_node = this->nnode();
557 for (unsigned n = 0; n < n_node; n++)
558 {
559 if (node_pt(n) == 0)
560 {
561 return false;
562 }
563 }
564 // If we get to here, OK
565 return true;
566 }
567 };
568
569} // namespace oomph
570
571#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
////////////////////////////////////////////////////////////////////
void get_x(const Vector< double > &s, Vector< double > &x) const
Global coordinates as function of local coordinates. Either via FE representation or via macro-elemen...
Definition: elements.h:1885
A general mesh class.
Definition: mesh.h:67
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:243
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:611
bool boundary_coordinate_exists(const unsigned &i) const
Indicate whether the i-th boundary has an intrinsic coordinate.
Definition: mesh.h:565
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
An OomphLibError object which should be thrown when an run-time error is encountered....
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
virtual void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get all continously interpolated function values in this element as a Vector. Note: Vector sets is ow...
void interpolated_zeta_on_face(const unsigned &boundary, const int &face, const Vector< double > &s, Vector< double > &zeta)
Return the value of the intrinsic boundary coordinate interpolated along the face.
void get_bcs(int bound, Vector< int > &bound_cons) const
Determine Vector of boundary conditions along the element's boundary (or vertex) bound (S/W/N/E/SW/SE...
void get_boundaries(const int &edge, std::set< unsigned > &boundaries) const
Given an element edge/vertex, return a Vector which contains all the (mesh-)boundary numbers that thi...
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
Refineable version of QuadElements that add functionality for spectral Elements.
virtual ~RefineableQSpectralElement()
Broken assignment operator.
void rebuild_from_sons(Mesh *&mesh_pt)
The only thing to add is rebuild from sons.
RefineableQSpectralElement(const RefineableQSpectralElement< 3 > &dummy)=delete
Broken copy constructor.
virtual bool nodes_built()
Overload the nodes built function.
A class that is used to template the refineable Q spectral elements by dimension. It's really nothing...
Refineable version of Solid brick elements.
void get_solid_bcs(int bound, Vector< int > &solid_bound_cons) const
Determine vector of solid (positional) boundary conditions along edge (or on vertex) bound (S/W/N/E/S...
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1816
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:262
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
long RefineableQElement< 2 > _build
//////////////////////////////////////////////////////////////////// ////////////////////////////////...