hp_refineable_elements.cc
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// Non-inline member functions for hp-refineable elements
27
28// oomph-lib includes
29#include "algebraic_elements.h"
32//#include "shape.h"
33
34namespace oomph
35{
36 /// /////////////////////////////////////////////////////////////
37 // 1D PRefineableQElements
38 /// /////////////////////////////////////////////////////////////
39
40 /// Get local coordinates of node j in the element; vector sets its own size
41 template<unsigned INITIAL_NNODE_1D>
43 const unsigned& n, Vector<double>& s) const
44 {
45 s.resize(1);
46
47 switch (this->nnode_1d())
48 {
49 case 2:
52 break;
53 case 3:
56 break;
57 case 4:
60 break;
61 case 5:
64 break;
65 case 6:
68 break;
69 case 7:
72 break;
73 default:
74 oomph_info << "\n ERROR: Exceeded maximum polynomial order for";
75 oomph_info << "\n shape functions." << std::endl;
76 break;
77 }
78 }
79
80 /// Get the local fractino of node j in the element
81 template<unsigned INITIAL_NNODE_1D>
83 const unsigned& n, Vector<double>& s_fraction)
84 {
85 this->local_coordinate_of_node(n, s_fraction);
86 s_fraction[0] = 0.5 * (s_fraction[0] + 1.0);
87 }
88
89 template<unsigned INITIAL_NNODE_1D>
91 const unsigned& n1d, const unsigned& i)
92 {
93 switch (this->nnode_1d())
94 {
95 case 2:
97 return 0.5 *
99 case 3:
101 return 0.5 *
103 case 4:
105 return 0.5 *
107 case 5:
109 return 0.5 *
111 case 6:
113 return 0.5 *
115 case 7:
117 return 0.5 *
119 default:
120 std::ostringstream error_message;
121 error_message << "\nERROR: Exceeded maximum polynomial order for";
122 error_message << "\n shape functions.\n";
123 throw OomphLibError(error_message.str(),
124 OOMPH_CURRENT_FUNCTION,
125 OOMPH_EXCEPTION_LOCATION);
126 return 0.0;
127 }
128 }
129
130
131 //==================================================================
132 /// Return the node at the specified local coordinate
133 //==================================================================
134 template<unsigned INITIAL_NNODE_1D>
136 const Vector<double>& s) const
137 {
138 // Load the tolerance into a local variable
140 // There is one possible index.
141 Vector<int> index(1);
142
143 // Determine the index
144 // -------------------
145
146 // If we are at the lower limit, the index is zero
147 if (std::fabs(s[0] + 1.0) < tol)
148 {
149 index[0] = 0;
150 }
151 // If we are at the upper limit, the index is the number of nodes minus 1
152 else if (std::fabs(s[0] - 1.0) < tol)
153 {
154 index[0] = this->nnode_1d() - 1;
155 }
156 // Otherwise, we have to calculate the index in general
157 else
158 {
159 // Compute Gauss-Lobatto-Legendre node positions
161 Orthpoly::gll_nodes(this->nnode_1d(), z);
162 // Loop over possible internal nodal positions
163 for (unsigned n = 1; n < this->nnode_1d() - 1; n++)
164 {
165 if (std::fabs(z[n] - s[0]) < tol)
166 {
167 index[0] = n;
168 break;
169 }
170 }
171 // No matching nodes
172 return 0;
173 }
174 // If we've got here we have a node, so let's return a pointer to it
175 return this->node_pt(index[0]);
176 }
177
178 //===================================================================
179 /// If a neighbouring element's son has already created a node at
180 /// a position corresponding to the local fractional position within the
181 /// present element, s_fraction, return
182 /// a pointer to that node. If not, return NULL (0). If the node is
183 /// periodic the flag is_periodic will be true
184 //===================================================================
185 template<unsigned INITIAL_NNODE_1D>
188 bool& is_periodic)
189 {
190 // Not possible in 1D case, so return null pointer
191 return 0;
192 }
193
194 //==================================================================
195 /// Set the correct p-order of the element based on that of its
196 /// father. Then construct an integration scheme of the correct order.
197 /// If an adopted father is specified, information from this is
198 /// used instead of using the father found from the tree.
199 //==================================================================
200 template<unsigned INITIAL_NNODE_1D>
202 Tree* const& adopted_father_pt, const unsigned& initial_p_order)
203 {
204 // Storage for pointer to my father (in binarytree impersonation)
205 BinaryTree* father_pt = 0;
206
207 // Check if an adopted father has been specified
208 if (adopted_father_pt != 0)
209 {
210 // Get pointer to my father (in binarytree impersonation)
211 father_pt = dynamic_cast<BinaryTree*>(adopted_father_pt);
212 }
213 // Check if element is in a tree
214 else if (Tree_pt != 0)
215 {
216 // Get pointer to my father (in binarytree impersonation)
217 father_pt = dynamic_cast<BinaryTree*>(binary_tree_pt()->father_pt());
218 }
219 // else
220 // {
221 // throw OomphLibError(
222 // "Element not in a tree, and no adopted father has been
223 // specified!", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
224 // }
225
226 // Check if element has father
227 if (father_pt != 0 || initial_p_order != 0)
228 {
229 if (father_pt != 0)
230 {
233 father_pt->object_pt());
234 if (father_el_pt != 0)
235 {
236 unsigned father_p_order = father_el_pt->p_order();
237 // Set p-order to that of father
238 P_order = father_p_order;
239 }
240 }
241 else
242 {
243 P_order = initial_p_order;
244 }
245
246 // Now sort out the element...
247 // (has p nodes)
248 unsigned new_n_node = P_order;
249
250 // Allocate new space for Nodes (at the element level)
251 this->set_n_node(new_n_node);
252
253 // Set integration scheme
254 delete this->integral_pt();
255 switch (P_order)
256 {
257 case 2:
258 this->set_integration_scheme(new GaussLobattoLegendre<1, 2>);
259 break;
260 case 3:
261 this->set_integration_scheme(new GaussLobattoLegendre<1, 3>);
262 break;
263 case 4:
264 this->set_integration_scheme(new GaussLobattoLegendre<1, 4>);
265 break;
266 case 5:
267 this->set_integration_scheme(new GaussLobattoLegendre<1, 5>);
268 break;
269 case 6:
270 this->set_integration_scheme(new GaussLobattoLegendre<1, 6>);
271 break;
272 case 7:
273 this->set_integration_scheme(new GaussLobattoLegendre<1, 7>);
274 break;
275 default:
276 std::ostringstream error_message;
277 error_message << "\nERROR: Exceeded maximum polynomial order for";
278 error_message << "\n integration scheme.\n";
279 throw OomphLibError(error_message.str(),
280 OOMPH_CURRENT_FUNCTION,
281 OOMPH_EXCEPTION_LOCATION);
282 }
283 }
284 }
285
286 //==================================================================
287 /// Check the father element for any required nodes which
288 /// already exist
289 //==================================================================
290 template<unsigned INITIAL_NNODE_1D>
292 Mesh*& mesh_pt, Vector<Node*>& new_node_pt)
293 {
294 /*
295 //Pointer to my father (in binarytree impersonation)
296 BinaryTree* father_pt =
297 dynamic_cast<BinaryTree*>(binary_tree_pt()->father_pt());
298
299 // Check if element has father
300 if (father_pt!=0)
301 {
302 PRefineableQElement<1>* father_el_pt =
303 dynamic_cast<PRefineableQElement<1>*>
304 (this->tree_pt()->father_pt()->object_pt());
305 if (father_el_pt!=0)
306 {
307 // Pre-build actions
308 //??
309 }
310 else
311 {
312 std::ostringstream error_message;
313 error_message <<"\nERROR: Dynamic cast failed!\n";
314 throw OomphLibError(error_message.str(),
315 OOMPH_CURRENT_FUNCTION,
316 OOMPH_EXCEPTION_LOCATION);
317 }
318 }
319 */
320 }
321
322 //==================================================================
323 /// p-refine the element inc times. (If inc<0 then p-unrefinement
324 /// is performed.)
325 //==================================================================
326 template<unsigned INITIAL_NNODE_1D>
328 const int& inc, Mesh* const& mesh_pt, GeneralisedElement* const& clone_pt)
329 {
330 // Cast clone to correct type
332 dynamic_cast<PRefineableQElement<1, INITIAL_NNODE_1D>*>(clone_pt);
333
334 // Check if we can use it
335 if (clone_el_pt == 0)
336 {
337 throw OomphLibError(
338 "Cloned copy must be a PRefineableQElement<1,INITIAL_NNODE_1D>!",
339 OOMPH_CURRENT_FUNCTION,
340 OOMPH_EXCEPTION_LOCATION);
341 }
342#ifdef PARANOID
343 // Clone exists, so check that it is infact a copy of me
344 else
345 {
346 // Flag to keep track of check
347 bool clone_is_ok = true;
348
349 // Does the clone have the correct p-order?
350 clone_is_ok = clone_is_ok && (clone_el_pt->p_order() == this->p_order());
351
352 if (!clone_is_ok)
353 {
354 std::ostringstream err_stream;
355 err_stream << "Clone element has a different p-order from me,"
356 << std::endl
357 << "but it is supposed to be a copy!" << std::endl;
358
359 throw OomphLibError(
360 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
361 }
362
363 // Does the clone have the same number of nodes as me?
364 clone_is_ok = clone_is_ok && (clone_el_pt->nnode() == this->nnode());
365
366 if (!clone_is_ok)
367 {
368 std::ostringstream err_stream;
369 err_stream << "Clone element has a different number of nodes from me,"
370 << std::endl
371 << "but it is supposed to be a copy!" << std::endl;
372
373 throw OomphLibError(
374 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
375 }
376
377 // Does the clone have the same nodes as me?
378 for (unsigned n = 0; n < this->nnode(); n++)
379 {
380 clone_is_ok =
381 clone_is_ok && (clone_el_pt->node_pt(n) == this->node_pt(n));
382 }
383
384 if (!clone_is_ok)
385 {
386 std::ostringstream err_stream;
387 err_stream << "The nodes belonging to the clone element are different"
388 << std::endl
389 << "from mine, but it is supposed to be a copy!"
390 << std::endl;
391
392 throw OomphLibError(
393 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
394 }
395
396 // If we get to here then the clone has all the information we require
397 }
398#endif
399
400 // Currently we can't handle the case of generalised coordinates
401 // since we haven't established how they should be interpolated.
402 // Buffer this case:
403 if (clone_el_pt->node_pt(0)->nposition_type() != 1)
404 {
405 throw OomphLibError("Can't handle generalised nodal positions (yet).",
406 OOMPH_CURRENT_FUNCTION,
407 OOMPH_EXCEPTION_LOCATION);
408 }
409
410 // Timestepper should be the same for all nodes -- use it
411 // to create timesteppers for new nodes
412 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
413
414 // Get number of history values (incl. present)
415 unsigned ntstorage = time_stepper_pt->ntstorage();
416
417 // Increment p-order of the element
418 P_order += inc;
419
420 // Change integration scheme
421 delete this->integral_pt();
422 switch (P_order)
423 {
424 case 2:
425 this->set_integration_scheme(new GaussLobattoLegendre<1, 2>);
426 break;
427 case 3:
428 this->set_integration_scheme(new GaussLobattoLegendre<1, 3>);
429 break;
430 case 4:
431 this->set_integration_scheme(new GaussLobattoLegendre<1, 4>);
432 break;
433 case 5:
434 this->set_integration_scheme(new GaussLobattoLegendre<1, 5>);
435 break;
436 case 6:
437 this->set_integration_scheme(new GaussLobattoLegendre<1, 6>);
438 break;
439 case 7:
440 this->set_integration_scheme(new GaussLobattoLegendre<1, 7>);
441 break;
442 default:
443 std::ostringstream error_message;
444 error_message << "\nERROR: Exceeded maximum polynomial order for";
445 error_message << "\n integration scheme.\n";
446 throw OomphLibError(error_message.str(),
447 OOMPH_CURRENT_FUNCTION,
448 OOMPH_EXCEPTION_LOCATION);
449 }
450
451 // Allocate new space for Nodes (at the element level)
452 this->set_n_node(P_order);
453
454 // Copy vertex nodes and create new internal nodes
455 //------------------------------------------------
456
457 // Setup vertex coordinates in element:
458 //-------------------------------------
459 Vector<double> s_left(1);
460 Vector<double> s_right(1);
461 s_left[0] = -1.0;
462 s_right[0] = 1.0;
463
464 // Local coordinate in element
465 Vector<double> s(1);
466
467 Vector<double> s_fraction(1);
468 // Loop over all nodes in the element
469 for (unsigned n = 0; n < P_order; n++)
470 {
471 // Get the fractional position (in the current element) of the node
472 // in the direction of s[0]
473 s_fraction[0] = this->local_one_d_fraction_of_node(n, 0);
474
475 // Evaluate the local coordinate of the node in the father element
476 s[0] = s_left[0] + (s_right[0] - s_left[0]) * s_fraction[0];
477
478 // Initialise flag: So far, this node hasn't been built or copied yet
479 bool node_done = false;
480
481 // Get the pointer to the node in this element (or rather, its clone),
482 // returns NULL if there is not node
483 Node* created_node_pt = clone_el_pt->get_node_at_local_coordinate(s);
484
485 // Does this node already exist in this element?
486 //----------------------------------------------
487 if (created_node_pt != 0)
488 {
489 // Copy node across
490 this->node_pt(n) = created_node_pt;
491
492 // Make sure that we update the values at the node so that they
493 // are consistent with the present representation. This is only
494 // needed for mixed interpolation where the value at the father
495 // could now become active.
496
497 // Loop over all history values
498 for (unsigned t = 0; t < ntstorage; t++)
499 {
500 // Get values from father element
501 // Note: get_interpolated_values() sets Vector size itself
502 Vector<double> prev_values;
503 clone_el_pt->get_interpolated_values(t, s, prev_values);
504 // Find the minimum number of values
505 //(either those stored at the node, or those returned by
506 // the function)
507 unsigned n_val_at_node = created_node_pt->nvalue();
508 unsigned n_val_from_function = prev_values.size();
509 // Use the ternary conditional operator here
510 unsigned n_var = n_val_at_node < n_val_from_function ?
511 n_val_at_node :
512 n_val_from_function;
513 // Assign the values that we can
514 for (unsigned k = 0; k < n_var; k++)
515 {
516 created_node_pt->set_value(t, k, prev_values[k]);
517 }
518 }
519
520 // Indicate that node has been created by copy
521 node_done = true;
522 }
523
524 // Node does not exist in this element
525 //------------------------------------
526
527 // If the node has not been built anywhere ---> build it here
528 if (!node_done)
529 {
530 // In a 1D mesh any node which lies on the boundary must exist in
531 // the initial (coarse) mesh, so any newly-built nodes cannot be
532 // boundary nodes. Therefore we always create a normal "bulk" node.
533
534 // Create node and set the pointer to it from the element
535 created_node_pt = this->construct_node(n, time_stepper_pt);
536
537 // Now we set the position and values at the newly created node
538
539 // Now we set the position and values at the newly created node.
540 // In the first instance use macro element or FE representation
541 // to create past and present nodal positions.
542 // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC ELEMENTS AS NOT
543 // ALL OF THEM NECESSARILY IMPLEMENT NONTRIVIAL NODE UPDATE
544 // FUNCTIONS. CALLING THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL
545 // LEAVE THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
546 // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL NOT ASSIGN SENSIBLE
547 // INITIAL POSITIONS!)
548
549 // Loop over history values
550 for (unsigned t = 0; t < ntstorage; t++)
551 {
552 // Allocate storage for the previous position of the node
553 Vector<double> x_prev(1);
554
555 // Get position from father element -- this uses the macro
556 // element representation if appropriate.
557 clone_el_pt->get_x(t, s, x_prev);
558
559 // Set the previous position of the new node
560 created_node_pt->x(t, 0) = x_prev[0];
561
562 // Allocate storage for the previous values at the node
563 // NOTE: the size of this vector is equal to the number of values
564 // (e.g. 3 velocity components and 1 pressure, say)
565 // associated with each node and NOT the number of history values
566 // which the node stores!
567 Vector<double> prev_values;
568
569 // Get values from father element
570 // Note: get_interpolated_values() sets Vector size itself.
571 clone_el_pt->get_interpolated_values(t, s, prev_values);
572
573 // Determine the number of values at the new node
574 const unsigned n_value = created_node_pt->nvalue();
575
576 // Loop over all values and set the previous values
577 for (unsigned v = 0; v < n_value; v++)
578 {
579 created_node_pt->set_value(t, v, prev_values[v]);
580 }
581 } // End of loop over history values
582
583 // Add new node to mesh (if requested)
584 if (mesh_pt != 0)
585 {
586 mesh_pt->add_node_pt(created_node_pt);
587 }
588
589 AlgebraicElementBase* alg_el_pt =
590 dynamic_cast<AlgebraicElementBase*>(this);
591
592 // If we do have an algebraic element
593 if (alg_el_pt != 0)
594 {
595 std::string error_message = "Have not implemented p-refinement for";
596 error_message += "Algebraic p-refineable elements yet\n";
597
598 throw OomphLibError(
599 error_message,
600 "PRefineableQElement<1,INITIAL_NNODE_1D>::p_refine()",
601 OOMPH_EXCEPTION_LOCATION);
602 }
603
604 } // End of case where we build the node ourselves
605
606 // Check if the element is an algebraic element
607 AlgebraicElementBase* alg_el_pt =
608 dynamic_cast<AlgebraicElementBase*>(this);
609
610 // If the element is an algebraic element, setup
611 // node position (past and present) from algebraic node update
612 // function. This over-writes previous assingments that
613 // were made based on the macro-element/FE representation.
614 // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
615 // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
616 // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
617 // INFO FOR *ALL* ROOT ELEMENTS!
618 if (alg_el_pt != 0)
619 {
620 // Build algebraic node update info for new node
621 // This sets up the node update data for all node update
622 // functions that are shared by all nodes in the father
623 // element
625 this->node_pt(n), s, clone_el_pt);
626 }
627
628 } // End of loop over all nodes in element
629
630
631 // If the element is a MacroElementNodeUpdateElement, set
632 // the update parameters for the current element's nodes --
633 // all this needs is the vector of (pointers to the)
634 // geometric objects that affect the MacroElement-based
635 // node update -- this needs to be done to set the node
636 // update info for newly created nodes
637 MacroElementNodeUpdateElementBase* clone_m_el_pt =
638 dynamic_cast<MacroElementNodeUpdateElementBase*>(clone_el_pt);
639 if (clone_m_el_pt != 0)
640 {
641 // Get vector of geometric objects from father (construct vector
642 // via copy operation)
643 Vector<GeomObject*> geom_object_pt(clone_m_el_pt->geom_object_pt());
644
645 // Cast current element to MacroElementNodeUpdateElement:
647 dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
648
649#ifdef PARANOID
650 if (m_el_pt == 0)
651 {
652 std::string error_message =
653 "Failed to cast to MacroElementNodeUpdateElementBase*\n";
654 error_message +=
655 "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
656 error_message += "then I should be too....\n";
657
658 throw OomphLibError(
659 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
660 }
661#endif
662 // Build update info by passing vector of geometric objects:
663 // This sets the current element to be the update element
664 // for all of the element's nodes -- this is reversed
665 // if the element is ever un-refined in the father element's
666 // rebuild_from_sons() function which overwrites this
667 // assignment to avoid nasty segmentation faults that occur
668 // when a node tries to update itself via an element that no
669 // longer exists...
670 m_el_pt->set_node_update_info(geom_object_pt);
671 }
672
673
674 // Loop over all nodes in element again, to re-set the positions
675 // This must be done using the new element's macro-element
676 // representation, rather than the old version which may be
677 // of a different p-order!
678 for (unsigned n = 0; n < P_order; n++)
679 {
680 // Get the fractional position of the node in the direction of s[0]
681 s_fraction[0] = this->local_one_d_fraction_of_node(n, 0);
682 // Local coordinate
683 s[0] = s_left[0] + (s_right[0] - s_left[0]) * s_fraction[0];
684
685 // Loop over # of history values
686 for (unsigned t = 0; t < ntstorage; t++)
687 {
688 // Get position from father element -- this uses the macro
689 // element representation if appropriate. If the node
690 // turns out to be a hanging node later on, then
691 // its position gets adjusted in line with its
692 // hanging node interpolation.
693 Vector<double> x_prev(1);
694 this->get_x(t, s, x_prev);
695
696 // Set previous positions of the new node
697 this->node_pt(n)->x(t, 0) = x_prev[0];
698 }
699 }
700
701 // Not necessary to delete the old nodes since all original nodes are in the
702 // current mesh and so will be pruned as part of the mesh adaption process.
703
704 // Do any further-build required
705 this->further_build();
706 }
707
708 //=======================================================================
709 /// Shape functions for PRefineableQElement<DIM>
710 //=======================================================================
711 template<unsigned INITIAL_NNODE_1D>
713 Shape& psi) const
714 {
715 switch (p_order())
716 {
717 case 2:
718 {
719 // Calculate nodal positions
721 // Create one-dim shape functions
723 // Now let's loop over the nodal points in the element
724 // and copy the values back in
725 for (unsigned i = 0; i < p_order(); i++)
726 {
727 psi(i) = psi1[i];
728 }
729 break;
730 }
731 case 3:
732 {
733 // Calculate nodal positions
735 // Create one-dim shape functions
737 // Now let's loop over the nodal points in the element
738 // and copy the values back in
739 for (unsigned i = 0; i < p_order(); i++)
740 {
741 psi(i) = psi1[i];
742 }
743 break;
744 }
745 case 4:
746 {
747 // Calculate nodal positions
749 // Create one-dim shape functions
751 // Now let's loop over the nodal points in the element
752 // and copy the values back in
753 for (unsigned i = 0; i < p_order(); i++)
754 {
755 psi(i) = psi1[i];
756 }
757 break;
758 }
759 case 5:
760 {
761 // Calculate nodal positions
763 // Create one-dim shape functions
765 // Now let's loop over the nodal points in the element
766 // and copy the values back in
767 for (unsigned i = 0; i < p_order(); i++)
768 {
769 psi(i) = psi1[i];
770 }
771 break;
772 }
773 case 6:
774 {
775 // Calculate nodal positions
777 // Create one-dim shape functions
779 // Now let's loop over the nodal points in the element
780 // and copy the values back in
781 for (unsigned i = 0; i < p_order(); i++)
782 {
783 psi(i) = psi1[i];
784 }
785 break;
786 }
787 case 7:
788 {
789 // Calculate nodal positions
791 // Create one-dim shape functions
793 // Now let's loop over the nodal points in the element
794 // and copy the values back in
795 for (unsigned i = 0; i < p_order(); i++)
796 {
797 psi(i) = psi1[i];
798 }
799 break;
800 }
801 default:
802 oomph_info << "\n ERROR: PRefineableQElement::shape() exceeded maximum";
803 oomph_info << "\n polynomial order for shape functions."
804 << std::endl;
805 }
806 }
807
808 //=======================================================================
809 /// Derivatives of shape functions for PRefineableQElement<DIM>
810 //=======================================================================
811 template<unsigned INITIAL_NNODE_1D>
813 const Vector<double>& s, Shape& psi, DShape& dpsi) const
814 {
815 switch (p_order())
816 {
817 case 2:
818 {
819 // Calculate nodal positions
821 // Call the shape functions and derivatives
824 // Loop over shapes and copy across
825 for (unsigned i = 0; i < p_order(); i++)
826 {
827 psi(i) = psi1[i];
828 dpsi(i, 0) = dpsi1ds[i];
829 }
830
831 break;
832 }
833 case 3:
834 {
835 // Calculate nodal positions
837 // Call the shape functions and derivatives
840 // Loop over shapes and copy across
841 for (unsigned i = 0; i < p_order(); i++)
842 {
843 psi(i) = psi1[i];
844 dpsi(i, 0) = dpsi1ds[i];
845 }
846 break;
847 }
848 case 4:
849 {
850 // Calculate nodal positions
852 // Call the shape functions and derivatives
855 // Loop over shapes and copy across
856 for (unsigned i = 0; i < p_order(); i++)
857 {
858 psi(i) = psi1[i];
859 dpsi(i, 0) = dpsi1ds[i];
860 }
861 break;
862 }
863 case 5:
864 {
865 // Calculate nodal positions
867 // Call the shape functions and derivatives
870 // Loop over shapes and copy across
871 for (unsigned i = 0; i < p_order(); i++)
872 {
873 psi(i) = psi1[i];
874 dpsi(i, 0) = dpsi1ds[i];
875 }
876 break;
877 }
878 case 6:
879 {
880 // Calculate nodal positions
882 // Call the shape functions and derivatives
885 // Loop over shapes and copy across
886 for (unsigned i = 0; i < p_order(); i++)
887 {
888 psi(i) = psi1[i];
889 dpsi(i, 0) = dpsi1ds[i];
890 }
891 break;
892 }
893 case 7:
894 {
895 // Calculate nodal positions
897 // Call the shape functions and derivatives
900 // Loop over shapes and copy across
901 for (unsigned i = 0; i < p_order(); i++)
902 {
903 psi(i) = psi1[i];
904 dpsi(i, 0) = dpsi1ds[i];
905 }
906 break;
907 }
908 default:
910 << "\n ERROR: PRefineableQElement::dshape_local() exceeded maximum";
911 oomph_info << "\n polynomial order for shape functions."
912 << std::endl;
913 }
914 }
915
916 //=======================================================================
917 /// Second derivatives of shape functions for PRefineableQElement<DIM>
918 /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
919 //=======================================================================
920 template<unsigned INITIAL_NNODE_1D>
922 const Vector<double>& s, Shape& psi, DShape& dpsids, DShape& d2psids) const
923 {
924 std::ostringstream error_message;
925 error_message
926 << "\nd2shape_local currently not implemented for this element\n";
927 throw OomphLibError(
928 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
929 }
930
931 //=================================================================
932 /// Internal function to set up the hanging nodes on a particular
933 /// edge of the element. (Not required in 1D.)
934 //=================================================================
935 template<unsigned INITIAL_NNODE_1D>
937 const int& value_id, const int& my_edge, std::ofstream& output_hangfile)
938 {
939 }
940
941 //=======================================================================
942 /// Rebuild the element from nodes found in its sons
943 /// Adjusts its p-order to be the maximum of its sons' p-orders
944 //=======================================================================
945 template<unsigned INITIAL_NNODE_1D>
947 Mesh*& mesh_pt)
948 {
949 // Get p-orders of sons
950 unsigned n_sons = this->tree_pt()->nsons();
951 Vector<unsigned> son_p_order(n_sons);
952 unsigned max_son_p_order = 0;
953 for (unsigned ison = 0; ison < n_sons; ison++)
954 {
955 PRefineableElement* el_pt = dynamic_cast<PRefineableElement*>(
956 this->tree_pt()->son_pt(ison)->object_pt());
957 son_p_order[ison] = el_pt->p_order();
958 if (son_p_order[ison] > max_son_p_order)
959 max_son_p_order = son_p_order[ison];
960 }
961
962 unsigned old_Nnode = this->nnode();
963 unsigned old_P_order = this->p_order();
964 // Set p-order of the element
965 this->p_order() = max_son_p_order;
966
967 // Change integration scheme
968 delete this->integral_pt();
969 switch (this->p_order())
970 {
971 case 2:
972 this->set_integration_scheme(new GaussLobattoLegendre<1, 2>);
973 break;
974 case 3:
975 this->set_integration_scheme(new GaussLobattoLegendre<1, 3>);
976 break;
977 case 4:
978 this->set_integration_scheme(new GaussLobattoLegendre<1, 4>);
979 break;
980 case 5:
981 this->set_integration_scheme(new GaussLobattoLegendre<1, 5>);
982 break;
983 case 6:
984 this->set_integration_scheme(new GaussLobattoLegendre<1, 6>);
985 break;
986 case 7:
987 this->set_integration_scheme(new GaussLobattoLegendre<1, 7>);
988 break;
989 default:
990 std::ostringstream error_message;
991 error_message << "\nERROR: Exceeded maximum polynomial order for";
992 error_message << "\n integration scheme.\n";
993 throw OomphLibError(error_message.str(),
994 OOMPH_CURRENT_FUNCTION,
995 OOMPH_EXCEPTION_LOCATION);
996 }
997
998 // Back up pointers to old nodes before they are lost
999 Vector<Node*> old_node_pt(old_Nnode);
1000 for (unsigned n = 0; n < old_Nnode; n++)
1001 {
1002 old_node_pt[n] = this->node_pt(n);
1003 }
1004
1005 // Allocate new space for Nodes (at the element level)
1006 this->set_n_node(this->p_order());
1007
1008 // Copy vertex nodes and create new edge and internal nodes
1009 //---------------------------------------------------------
1010
1011 // Copy vertex nodes
1012 this->node_pt(0) = old_node_pt[0];
1013 this->node_pt(this->p_order() - 1) = old_node_pt[old_P_order - 1];
1014
1015
1016 //=============================================================
1017 // Below this line is copied from RefineableQSpectralElement<2>
1018
1019 // The timestepper should be the same for all nodes and node 0 should
1020 // never be deleted.
1021 if (this->node_pt(0) == 0)
1022 {
1023 throw OomphLibError("The vertex node (0) does not exist",
1024 OOMPH_CURRENT_FUNCTION,
1025 OOMPH_EXCEPTION_LOCATION);
1026 }
1027
1028 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
1029
1030 // Determine number of history values stored
1031 const unsigned ntstorage = time_stepper_pt->ntstorage();
1032
1033 // Allocate storage for local coordinates
1034 Vector<double> s_fraction(1), s(1);
1035
1036 // Determine the number of nodes in the element
1037 const unsigned n_node = this->nnode_1d();
1038
1039 // Loop over the nodes in the element
1040 for (unsigned n = 0; n < n_node; n++)
1041 {
1042 // Get the fractional position of the node in the direction of s[0]
1043 s_fraction[0] = this->local_one_d_fraction_of_node(n, 0);
1044
1045 // Determine the local coordinate in the father element
1046 s[0] = -1.0 + 2.0 * s_fraction[0];
1047
1048 // If the node has not been built
1049 if (this->node_pt(n) == 0)
1050 {
1051 // Has the node been created by one of its neighbours?
1052 bool is_periodic = false;
1053 Node* created_node_pt =
1054 this->node_created_by_neighbour(s_fraction, is_periodic);
1055
1056 // If it has, set the pointer
1057 if (created_node_pt != 0)
1058 {
1059 // If the node is periodic
1060 if (is_periodic)
1061 {
1062 throw OomphLibError("Cannot handle periodic nodes yet",
1063 OOMPH_CURRENT_FUNCTION,
1064 OOMPH_EXCEPTION_LOCATION);
1065 }
1066 // Non-periodic case, just set the pointer
1067 else
1068 {
1069 this->node_pt(n) = created_node_pt;
1070 }
1071 }
1072 // Otherwise, we need to build it
1073 else
1074 {
1075 // First, find the son element in which the node should live
1076
1077 // Find coordinate in the son
1078 Vector<double> s_in_son(1);
1079 using namespace BinaryTreeNames;
1080 int son = -10;
1081 // If s_fraction is between 0 and 0.5, we are in the left son
1082 if (s_fraction[0] < 0.5)
1083 {
1084 son = L;
1085 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
1086 }
1087 // Otherwise we are in the right son
1088 else
1089 {
1090 son = R;
1091 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
1092 }
1093
1094 // Get the pointer to the son element in which the new node
1095 // would live
1098 this->tree_pt()->son_pt(son)->object_pt());
1099
1100 // In 1D we should never be rebuilding an element's vertex nodes
1101 // (since they will never be deleted), so throw an error if we
1102 // appear to be doing so
1103#ifdef PARANOID
1104 if (n == 0 || n == n_node - 1)
1105 {
1106 std::string error_message =
1107 "I am trying to rebuild one of the (two) vertex nodes in\n";
1108 error_message +=
1109 "this 1D element. It should not have been possible to delete\n";
1110 error_message += "either of these!\n";
1111
1112 throw OomphLibError(
1113 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1114 }
1115#endif
1116
1117 // With this in mind we will always be creating normal "bulk" nodes
1118 this->node_pt(n) = this->construct_node(n, time_stepper_pt);
1119
1120 // Now we set the position and values at the newly created node
1121
1122 // In the first instance use macro element or FE representation
1123 // to create past and present nodal positions.
1124 // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC ELEMENTS AS NOT
1125 // ALL OF THEM NECESSARILY IMPLEMENT NONTRIVIAL NODE UPDATE
1126 // FUNCTIONS. CALLING THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL
1127 // LEAVE THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
1128 // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL NOT ASSIGN SENSIBLE
1129 // INITIAL POSITIONS!)
1130
1131 // Loop over history values
1132 for (unsigned t = 0; t < ntstorage; t++)
1133 {
1134 // Allocate storage for the previous position of the node
1135 Vector<double> x_prev(1);
1136
1137 // Get position from son element -- this uses the macro element
1138 // representation if appropriate
1139 son_el_pt->get_x(t, s_in_son, x_prev);
1140
1141 // Set the previous position of the new node
1142 this->node_pt(n)->x(t, 0) = x_prev[0];
1143
1144 // Allocate storage for the previous values at the node
1145 // NOTE: the size of this vector is equal to the number of values
1146 // (e.g. 3 velocity components and 1 pressure, say)
1147 // associated with each node and NOT the number of history values
1148 // which the node stores!
1149 Vector<double> prev_values;
1150
1151 // Get values from son element
1152 // Note: get_interpolated_values() sets Vector size itself.
1153 son_el_pt->get_interpolated_values(t, s_in_son, prev_values);
1154
1155 // Determine the number of values at the new node
1156 const unsigned n_value = this->node_pt(n)->nvalue();
1157
1158 // Loop over all values and set the previous values
1159 for (unsigned v = 0; v < n_value; v++)
1160 {
1161 this->node_pt(n)->set_value(t, v, prev_values[v]);
1162 }
1163 } // End of loop over history values
1164
1165 // Add new node to mesh
1166 mesh_pt->add_node_pt(this->node_pt(n));
1167
1168 } // End of case where we build the node ourselves
1169
1170 } // End of if this node has not been built
1171 } // End of loop over nodes in element
1172
1173 // Check if the element is an algebraic element
1174 // This is done on all nodes in the element after reconstruction
1175 // rather than as the nodes are built
1176 AlgebraicElementBase* alg_el_pt = dynamic_cast<AlgebraicElementBase*>(this);
1177
1178 // If so, throw error
1179 if (alg_el_pt != 0)
1180 {
1181 std::string error_message =
1182 "Have not implemented rebuilding from sons for";
1183 error_message += "Algebraic p-refineable elements yet\n";
1184
1185 throw OomphLibError(
1186 error_message,
1187 "PRefineableQElement<1,INITIAL_NNODE_1D>::rebuild_from_sons()",
1188 OOMPH_EXCEPTION_LOCATION);
1189 }
1190 }
1191
1192 //=================================================================
1193 /// Check inter-element continuity of
1194 /// - nodal positions
1195 /// - (nodally) interpolated function values
1196 //====================================================================
1197 template<unsigned INITIAL_NNODE_1D>
1199 double& max_error)
1200 {
1202 }
1203
1204 /// /////////////////////////////////////////////////////////////
1205 // 2D PRefineableQElements
1206 /// /////////////////////////////////////////////////////////////
1207
1208 /// Get local coordinates of node j in the element; vector sets its own size
1209 template<unsigned INITIAL_NNODE_1D>
1211 const unsigned& n, Vector<double>& s) const
1212 {
1213 s.resize(2);
1214 unsigned Nnode_1d = this->nnode_1d();
1215 unsigned n0 = n % Nnode_1d;
1216 unsigned n1 = n / Nnode_1d;
1217
1218 switch (Nnode_1d)
1219 {
1220 case 2:
1224 break;
1225 case 3:
1229 break;
1230 case 4:
1234 break;
1235 case 5:
1239 break;
1240 case 6:
1244 break;
1245 case 7:
1249 break;
1250 default:
1251 std::ostringstream error_message;
1252 error_message << "\nERROR: Exceeded maximum polynomial order for";
1253 error_message << "\n shape functions.\n";
1254 throw OomphLibError(error_message.str(),
1255 OOMPH_CURRENT_FUNCTION,
1256 OOMPH_EXCEPTION_LOCATION);
1257 break;
1258 }
1259 }
1260
1261 /// Get the local fractino of node j in the element
1262 template<unsigned INITIAL_NNODE_1D>
1264 const unsigned& n, Vector<double>& s_fraction)
1265 {
1266 this->local_coordinate_of_node(n, s_fraction);
1267 s_fraction[0] = 0.5 * (s_fraction[0] + 1.0);
1268 s_fraction[1] = 0.5 * (s_fraction[1] + 1.0);
1269 }
1270
1271 template<unsigned INITIAL_NNODE_1D>
1273 const unsigned& n1d, const unsigned& i)
1274 {
1275 switch (this->nnode_1d())
1276 {
1277 case 2:
1279 return 0.5 *
1281 case 3:
1283 return 0.5 *
1285 case 4:
1287 return 0.5 *
1289 case 5:
1291 return 0.5 *
1293 case 6:
1295 return 0.5 *
1297 case 7:
1299 return 0.5 *
1301 default:
1302 std::ostringstream error_message;
1303 error_message << "\nERROR: Exceeded maximum polynomial order for";
1304 error_message << "\n shape functions.\n";
1305 throw OomphLibError(error_message.str(),
1306 OOMPH_CURRENT_FUNCTION,
1307 OOMPH_EXCEPTION_LOCATION);
1308 return 0.0;
1309 }
1310 }
1311
1312
1313 //==================================================================
1314 /// Return the node at the specified local coordinate
1315 //==================================================================
1316 template<unsigned INITIAL_NNODE_1D>
1318 const Vector<double>& s) const
1319 {
1320 unsigned Nnode_1d = this->nnode_1d();
1321 // Load the tolerance into a local variable
1323 // There are two possible indices.
1324 Vector<int> index(2);
1325
1326 // Loop over indices
1327 for (unsigned i = 0; i < 2; i++)
1328 {
1329 // Determine the index
1330 // -------------------
1331
1332 bool is_found = false;
1333
1334 // If we are at the lower limit, the index is zero
1335 if (std::fabs(s[i] + 1.0) < tol)
1336 {
1337 index[i] = 0;
1338 is_found = true;
1339 }
1340 // If we are at the upper limit, the index is the number of nodes minus 1
1341 else if (std::fabs(s[i] - 1.0) < tol)
1342 {
1343 index[i] = Nnode_1d - 1;
1344 is_found = true;
1345 }
1346 // Otherwise, we have to calculate the index in general
1347 else
1348 {
1349 // Compute Gauss-Lobatto-Legendre node positions
1351 Orthpoly::gll_nodes(Nnode_1d, z);
1352 // Loop over possible internal nodal positions
1353 for (unsigned n = 1; n < Nnode_1d - 1; n++)
1354 {
1355 if (std::fabs(z[n] - s[i]) < tol)
1356 {
1357 index[i] = n;
1358 is_found = true;
1359 break;
1360 }
1361 }
1362 }
1363
1364 if (!is_found)
1365 {
1366 // No matching nodes
1367 return 0;
1368 }
1369 }
1370 // If we've got here we have a node, so let's return a pointer to it
1371 return this->node_pt(index[0] + Nnode_1d * index[1]);
1372 }
1373
1374 //===================================================================
1375 /// If a neighbouring element has already created a node at
1376 /// a position corresponding to the local fractional position within the
1377 /// present element, s_fraction, return
1378 /// a pointer to that node. If not, return NULL (0). If the node is
1379 /// periodic the flag is_periodic will be true
1380 //===================================================================
1381 template<unsigned INITIAL_NNODE_1D>
1383 const Vector<double>& s_fraction, bool& is_periodic)
1384 {
1385 using namespace QuadTreeNames;
1386
1387 // Calculate the edges on which the node lies
1388 Vector<int> edges;
1389 if (s_fraction[0] == 0.0)
1390 {
1391 edges.push_back(W);
1392 }
1393 if (s_fraction[0] == 1.0)
1394 {
1395 edges.push_back(E);
1396 }
1397 if (s_fraction[1] == 0.0)
1398 {
1399 edges.push_back(S);
1400 }
1401 if (s_fraction[1] == 1.0)
1402 {
1403 edges.push_back(N);
1404 }
1405
1406 // Find the number of edges
1407 unsigned n_size = edges.size();
1408 // If there are no edges, then there is no neighbour, return 0
1409 if (n_size == 0)
1410 {
1411 return 0;
1412 }
1413
1414 Vector<unsigned> translate_s(2);
1415 Vector<double> s_lo_neigh(2);
1416 Vector<double> s_hi_neigh(2);
1417 Vector<double> s(2);
1418
1419 int neigh_edge, diff_level;
1420 bool in_neighbouring_tree;
1421
1422 // Loop over the edges
1423 for (unsigned j = 0; j < n_size; j++)
1424 {
1425 // Find pointer to neighbouring element along edge
1426 QuadTree* neigh_pt;
1427 neigh_pt = quadtree_pt()->gteq_edge_neighbour(edges[j],
1428 translate_s,
1429 s_lo_neigh,
1430 s_hi_neigh,
1431 neigh_edge,
1432 diff_level,
1433 in_neighbouring_tree);
1434
1435 // Neighbour exists
1436 if (neigh_pt != 0)
1437 {
1438 // Have any of its vertex nodes been created yet?
1439 // (Must look in incomplete neighbours because after the
1440 // pre-build they may contain pointers to the required nodes. e.g.
1441 // h-refinement of neighbouring linear and quadratic elements)
1442 bool a_vertex_node_is_built = false;
1443 QElement<2, INITIAL_NNODE_1D>* neigh_obj_pt =
1444 dynamic_cast<QElement<2, INITIAL_NNODE_1D>*>(neigh_pt->object_pt());
1445 if (neigh_obj_pt == 0)
1446 {
1447 throw OomphLibError("Not a quad element!",
1448 "PRefineableQElement<2,INITIAL_NNODE_1D>::node_"
1449 "created_by_neighbour()",
1450 OOMPH_EXCEPTION_LOCATION);
1451 }
1452 for (unsigned vnode = 0; vnode < neigh_obj_pt->nvertex_node(); vnode++)
1453 {
1454 if (neigh_obj_pt->vertex_node_pt(vnode) != 0)
1455 a_vertex_node_is_built = true;
1456 break;
1457 }
1458 if (a_vertex_node_is_built)
1459 {
1460 // We now need to translate the nodal location
1461 // as defined in terms of the fractional coordinates of the present
1462 // element into those of its neighbour
1463
1464 // Calculate the local coordinate in the neighbour
1465 // Note that we need to use the translation scheme in case
1466 // the local coordinates are swapped in the neighbour.
1467 for (unsigned i = 0; i < 2; i++)
1468 {
1469 s[i] = s_lo_neigh[i] +
1470 s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
1471 }
1472
1473 // Find the node in the neighbour
1474 Node* neighbour_node_pt =
1476
1477 // If there is a node, return it
1478 if (neighbour_node_pt != 0)
1479 {
1480 // Now work out whether it's a periodic boundary
1481 // only possible if we have moved into a neighbouring tree
1482 if (in_neighbouring_tree)
1483 {
1484 // Return whether the neighbour is periodic
1485 is_periodic =
1486 quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
1487 }
1488 // Return the pointer to the neighbouring node
1489 return neighbour_node_pt;
1490 }
1491 }
1492 }
1493 }
1494 // Node not found, return null
1495 return 0;
1496 }
1497
1498 //===================================================================
1499 /// If a neighbouring element's son has already created a node at
1500 /// a position corresponding to the local fractional position within the
1501 /// present element, s_fraction, return
1502 /// a pointer to that node. If not, return NULL (0). If the node is
1503 /// periodic the flag is_periodic will be true
1504 //===================================================================
1505 template<unsigned INITIAL_NNODE_1D>
1508 bool& is_periodic)
1509 {
1510 using namespace QuadTreeNames;
1511
1512 // Calculate the edges on which the node lies
1513 Vector<int> edges;
1514 if (s_fraction[0] == 0.0)
1515 {
1516 edges.push_back(W);
1517 }
1518 if (s_fraction[0] == 1.0)
1519 {
1520 edges.push_back(E);
1521 }
1522 if (s_fraction[1] == 0.0)
1523 {
1524 edges.push_back(S);
1525 }
1526 if (s_fraction[1] == 1.0)
1527 {
1528 edges.push_back(N);
1529 }
1530
1531 // Find the number of edges
1532 unsigned n_size = edges.size();
1533 // If there are no edges, then there is no neighbour, return 0
1534 if (n_size == 0)
1535 {
1536 return 0;
1537 }
1538
1539 Vector<unsigned> translate_s(2);
1540 Vector<double> s_lo_neigh(2);
1541 Vector<double> s_hi_neigh(2);
1542 Vector<double> s(2);
1543
1544 int neigh_edge, diff_level;
1545 bool in_neighbouring_tree;
1546
1547 // Loop over the edges
1548 for (unsigned j = 0; j < n_size; j++)
1549 {
1550 // Find pointer to neighbouring element along edge
1551 QuadTree* neigh_pt;
1552 neigh_pt = quadtree_pt()->gteq_edge_neighbour(edges[j],
1553 translate_s,
1554 s_lo_neigh,
1555 s_hi_neigh,
1556 neigh_edge,
1557 diff_level,
1558 in_neighbouring_tree);
1559
1560 // Neighbour exists
1561 if (neigh_pt != 0)
1562 {
1563 // Have its nodes been created yet?
1564 // (Must look in sons of unfinished neighbours too!!!)
1565 if (1)
1566 {
1567 // We now need to translate the nodal location
1568 // as defined in terms of the fractional coordinates of the present
1569 // element into those of its neighbour
1570
1571 // Calculate the local coordinate in the neighbour
1572 // Note that we need to use the translation scheme in case
1573 // the local coordinates are swapped in the neighbour.
1574 for (unsigned i = 0; i < 2; i++)
1575 {
1576 s[i] = s_lo_neigh[i] +
1577 s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
1578 }
1579
1580 // Check if the element has sons
1581 if (neigh_pt->nsons() != 0)
1582 {
1583 // First, find the son element in which the node should live
1584
1585 // Find coordinates in the sons
1586 Vector<double> s_in_son(2);
1587 int son = -10;
1588 // If negative on the west side
1589 if (s[0] < 0.0)
1590 {
1591 // On the south side
1592 if (s[1] < 0.0)
1593 {
1594 // It's the southwest son
1595 son = SW;
1596 s_in_son[0] = 1.0 + 2.0 * s[0];
1597 s_in_son[1] = 1.0 + 2.0 * s[1];
1598 }
1599 // On the north side
1600 else
1601 {
1602 // It's the northwest son
1603 son = NW;
1604 s_in_son[0] = 1.0 + 2.0 * s[0];
1605 s_in_son[1] = -1.0 + 2.0 * s[1];
1606 }
1607 }
1608 else
1609 {
1610 // On the south side
1611 if (s[1] < 0.0)
1612 {
1613 // It's the southeast son
1614 son = SE;
1615 s_in_son[0] = -1.0 + 2.0 * s[0];
1616 s_in_son[1] = 1.0 + 2.0 * s[1];
1617 }
1618 // On the north side
1619 else
1620 {
1621 // It's the northeast son
1622 son = NE;
1623 s_in_son[0] = -1.0 + 2.0 * s[0];
1624 s_in_son[1] = -1.0 + 2.0 * s[1];
1625 }
1626 }
1627
1628 // Find the node in the neighbour's son
1629 Node* neighbour_son_node_pt =
1631 s_in_son);
1632
1633 // If there is a node, return it
1634 if (neighbour_son_node_pt != 0)
1635 {
1636 // Now work out whether it's a periodic boundary
1637 // only possible if we have moved into a neighbouring tree
1638 if (in_neighbouring_tree)
1639 {
1640 // Return whether the neighbour is periodic
1641 is_periodic =
1642 quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
1643 }
1644 // Return the pointer to the neighbouring node
1645 return neighbour_son_node_pt;
1646 }
1647 }
1648 else
1649 {
1650 // No sons to search in, so no node can be found
1651 return 0;
1652 }
1653 }
1654 }
1655 }
1656 // Node not found, return null
1657 return 0;
1658 }
1659
1660 //==================================================================
1661 /// Set the correct p-order of the element based on that of its
1662 /// father. Then construct an integration scheme of the correct order.
1663 /// If an adopted father is specified, information from this is
1664 /// used instead of using the father found from the tree.
1665 //==================================================================
1666 template<unsigned INITIAL_NNODE_1D>
1668 Tree* const& adopted_father_pt, const unsigned& initial_p_order)
1669 {
1670 // Storage for pointer to my father (in quadtree impersonation)
1671 QuadTree* father_pt = 0;
1672
1673 // Check if an adopted father has been specified
1674 if (adopted_father_pt != 0)
1675 {
1676 // Get pointer to my father (in quadtree impersonation)
1677 father_pt = dynamic_cast<QuadTree*>(adopted_father_pt);
1678 }
1679 // Check if element is in a tree
1680 else if (Tree_pt != 0)
1681 {
1682 // Get pointer to my father (in quadtree impersonation)
1683 father_pt = dynamic_cast<QuadTree*>(quadtree_pt()->father_pt());
1684 }
1685 // else
1686 // {
1687 // throw OomphLibError(
1688 // "Element not in a tree, and no adopted father has been
1689 // specified!", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1690 // }
1691
1692 // Check if we found a father
1693 if (father_pt != 0 || initial_p_order != 0)
1694 {
1695 if (father_pt != 0)
1696 {
1699 father_pt->object_pt());
1700 if (father_el_pt != 0)
1701 {
1702 unsigned father_p_order = father_el_pt->p_order();
1703 // Set p-order to that of father
1704 P_order = father_p_order;
1705 }
1706 }
1707 else
1708 {
1709 P_order = initial_p_order;
1710 }
1711
1712 // Now sort out the element...
1713 // (has p^2 nodes)
1714 unsigned new_n_node = P_order * P_order;
1715
1716 // Allocate new space for Nodes (at the element level)
1717 this->set_n_node(new_n_node);
1718
1719 // Set integration scheme
1720 delete this->integral_pt();
1721 switch (P_order)
1722 {
1723 case 2:
1724 this->set_integration_scheme(new GaussLobattoLegendre<2, 2>);
1725 break;
1726 case 3:
1727 this->set_integration_scheme(new GaussLobattoLegendre<2, 3>);
1728 break;
1729 case 4:
1730 this->set_integration_scheme(new GaussLobattoLegendre<2, 4>);
1731 break;
1732 case 5:
1733 this->set_integration_scheme(new GaussLobattoLegendre<2, 5>);
1734 break;
1735 case 6:
1736 this->set_integration_scheme(new GaussLobattoLegendre<2, 6>);
1737 break;
1738 case 7:
1739 this->set_integration_scheme(new GaussLobattoLegendre<2, 7>);
1740 break;
1741 default:
1742 std::ostringstream error_message;
1743 error_message << "\nERROR: Exceeded maximum polynomial order for";
1744 error_message << "\n integration scheme.\n";
1745 throw OomphLibError(error_message.str(),
1746 OOMPH_CURRENT_FUNCTION,
1747 OOMPH_EXCEPTION_LOCATION);
1748 }
1749 }
1750 }
1751
1752 //==================================================================
1753 /// Check the father element for any required nodes which
1754 /// already exist
1755 //==================================================================
1756 template<unsigned INITIAL_NNODE_1D>
1758 Mesh*& mesh_pt, Vector<Node*>& new_node_pt)
1759 {
1760 using namespace QuadTreeNames;
1761
1762 // Get the number of 1d nodes
1763 unsigned n_p = nnode_1d();
1764
1765 // Check whether static father_bound needs to be created
1766 if (Father_bound[n_p].nrow() == 0)
1767 {
1768 setup_father_bounds();
1769 }
1770
1771 // Pointer to my father (in quadtree impersonation)
1772 QuadTree* father_pt = dynamic_cast<QuadTree*>(quadtree_pt()->father_pt());
1773
1774 // What type of son am I? Ask my quadtree representation...
1775 int son_type = Tree_pt->son_type();
1776
1777 // Has somebody build me already? (If any nodes have not been built)
1778 if (!nodes_built())
1779 {
1780#ifdef PARANOID
1781 if (father_pt == 0)
1782 {
1783 std::string error_message =
1784 "Something fishy here: I have no father and yet \n";
1785 error_message += "I have no nodes. Who has created me then?!\n";
1786
1787 throw OomphLibError(
1788 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1789 }
1790#endif
1791
1792 // Return pointer to father element
1793 RefineableQElement<2>* father_el_pt =
1794 dynamic_cast<RefineableQElement<2>*>(father_pt->object_pt());
1795
1796 // Timestepper should be the same for all nodes in father
1797 // element -- use it create timesteppers for new nodes
1798 TimeStepper* time_stepper_pt =
1799 father_el_pt->node_pt(0)->time_stepper_pt();
1800
1801 // Number of history values (incl. present)
1802 unsigned ntstorage = time_stepper_pt->ntstorage();
1803
1804 Vector<double> s_lo(2);
1805 Vector<double> s_hi(2);
1806 Vector<double> s(2);
1807 Vector<double> x(2);
1808
1809 // Setup vertex coordinates in father element:
1810 //--------------------------------------------
1811 switch (son_type)
1812 {
1813 case SW:
1814 s_lo[0] = -1.0;
1815 s_hi[0] = 0.0;
1816 s_lo[1] = -1.0;
1817 s_hi[1] = 0.0;
1818 break;
1819
1820 case SE:
1821 s_lo[0] = 0.0;
1822 s_hi[0] = 1.0;
1823 s_lo[1] = -1.0;
1824 s_hi[1] = 0.0;
1825 break;
1826
1827 case NE:
1828 s_lo[0] = 0.0;
1829 s_hi[0] = 1.0;
1830 s_lo[1] = 0.0;
1831 s_hi[1] = 1.0;
1832 break;
1833
1834 case NW:
1835 s_lo[0] = -1.0;
1836 s_hi[0] = 0.0;
1837 s_lo[1] = 0.0;
1838 s_hi[1] = 1.0;
1839 break;
1840 }
1841
1842 // Pass macro element pointer on to sons and
1843 // set coordinates in macro element
1844 // hierher why can I see this?
1845 if (father_el_pt->macro_elem_pt() != 0)
1846 {
1847 set_macro_elem_pt(father_el_pt->macro_elem_pt());
1848 for (unsigned i = 0; i < 2; i++)
1849 {
1850 s_macro_ll(i) =
1851 father_el_pt->s_macro_ll(i) +
1852 0.5 * (s_lo[i] + 1.0) *
1853 (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
1854 s_macro_ur(i) =
1855 father_el_pt->s_macro_ll(i) +
1856 0.5 * (s_hi[i] + 1.0) *
1857 (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
1858 }
1859 }
1860
1861
1862 // If the father element hasn't been generated yet, we're stuck...
1863 if (father_el_pt->node_pt(0) == 0)
1864 {
1865 throw OomphLibError(
1866 "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
1867 OOMPH_CURRENT_FUNCTION,
1868 OOMPH_EXCEPTION_LOCATION);
1869 }
1870 else
1871 {
1872 unsigned jnod = 0;
1873 Vector<double> x_small(2);
1874 Vector<double> x_large(2);
1875
1876 Vector<double> s_fraction(2);
1877 // Loop over nodes in element
1878 for (unsigned i0 = 0; i0 < n_p; i0++)
1879 {
1880 // Get the fractional position of the node in the direction of s[0]
1881 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
1882 // Local coordinate in father element
1883 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
1884
1885 for (unsigned i1 = 0; i1 < n_p; i1++)
1886 {
1887 // Get the fractional position of the node in the direction of s[1]
1888 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
1889 // Local coordinate in father element
1890 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
1891
1892 // Local node number
1893 jnod = i0 + n_p * i1;
1894
1895 // Check whether the father's node is periodic if so, complain
1896 /* {
1897 Node* father_node_pt = father_el_pt->node_pt(jnod);
1898 if((father_node_pt->is_a_copy()) ||
1899 (father_node_pt->position_is_a_copy()))
1900 {
1901 throw OomphLibError(
1902 "Can't handle periodic nodes (yet).",
1903 OOMPH_CURRENT_FUNCTION,
1904 OOMPH_EXCEPTION_LOCATION);
1905 }
1906 }*/
1907
1908 // Initialise flag: So far, this node hasn't been built
1909 // or copied yet
1910 // bool node_done=false;
1911
1912 // Get the pointer to the node in the father, returns NULL
1913 // if there is not node
1914 Node* created_node_pt =
1915 father_el_pt->get_node_at_local_coordinate(s);
1916
1917 // Does this node already exist in father element?
1918 //------------------------------------------------
1919 if (created_node_pt != 0)
1920 {
1921 // Copy node across
1922 this->node_pt(jnod) = created_node_pt;
1923
1924 // Make sure that we update the values at the node so that
1925 // they are consistent with the present representation.
1926 // This is only need for mixed interpolation where the value
1927 // at the father could now become active.
1928
1929 // Loop over all history values
1930 for (unsigned t = 0; t < ntstorage; t++)
1931 {
1932 // Get values from father element
1933 // Note: get_interpolated_values() sets Vector size itself.
1934 Vector<double> prev_values;
1935 father_el_pt->get_interpolated_values(t, s, prev_values);
1936 // Find the minimum number of values
1937 //(either those stored at the node, or those returned by
1938 // the function)
1939 unsigned n_val_at_node = created_node_pt->nvalue();
1940 unsigned n_val_from_function = prev_values.size();
1941 // Use the ternary conditional operator here
1942 unsigned n_var = n_val_at_node < n_val_from_function ?
1943 n_val_at_node :
1944 n_val_from_function;
1945 // Assign the values that we can
1946 for (unsigned k = 0; k < n_var; k++)
1947 {
1948 created_node_pt->set_value(t, k, prev_values[k]);
1949 }
1950 }
1951
1952 // Node has been created by copy
1953 // node_done=true;
1954 }
1955 } // End of vertical loop over nodes in element
1956 } // End of horizontal loop over nodes in element
1957 } // Sanity check: Father element has been generated
1958
1959 } // End of nodes not built
1960 else
1961 {
1962 // If this is element updates by macro element node update then we need
1963 // to set the correct info in the nodes here since this won't get done
1964 // later in build() because we already have all our nodes created.
1966 dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
1967 if (m_el_pt != 0)
1968 {
1969 // Get vector of geometric objects
1970 Vector<GeomObject*> geom_object_pt = m_el_pt->geom_object_pt();
1971
1972 /// / Build update info by passing vector of geometric objects:
1973 /// / This sets the current element to be the update element
1974 /// / for all of the element's nodes -- this is reversed
1975 /// / if the element is ever un-refined in the father element's
1976 /// / rebuild_from_sons() function which overwrites this
1977 /// / assignment to avoid nasty segmentation faults that occur
1978 /// / when a node tries to update itself via an element that no
1979 /// / longer exists...
1980 m_el_pt->set_node_update_info(geom_object_pt);
1981 }
1982 }
1983 }
1984
1985 //==================================================================
1986 /// p-refine the element inc times. (If inc<0 then p-unrefinement
1987 /// is performed.)
1988 //==================================================================
1989 template<unsigned INITIAL_NNODE_1D>
1991 const int& inc, Mesh* const& mesh_pt, GeneralisedElement* const& clone_pt)
1992 {
1993 // Cast clone to correct type
1995 dynamic_cast<PRefineableQElement<2, INITIAL_NNODE_1D>*>(clone_pt);
1996
1997 // If it is a MacroElement then we need to copy the geometric data too.
1999 dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
2000 if (m_el_pt != 0)
2001 {
2002 MacroElementNodeUpdateElementBase* clone_m_el_pt =
2003 dynamic_cast<MacroElementNodeUpdateElementBase*>(clone_pt);
2004
2005 // Store local copy of geom object vector, so it can be passed on
2006 // to son elements (and their nodes) during refinement
2007 unsigned ngeom_object = m_el_pt->ngeom_object();
2008 clone_m_el_pt->geom_object_pt().resize(ngeom_object);
2009 for (unsigned i = 0; i < ngeom_object; i++)
2010 {
2011 clone_m_el_pt->geom_object_pt()[i] = m_el_pt->geom_object_pt(i);
2012 }
2013
2014 // clone_m_el_pt->set_node_update_info(m_el_pt->geom_object_pt());
2015 }
2016
2017 // Check if we can use it
2018 if (clone_el_pt == 0)
2019 {
2020 throw OomphLibError(
2021 "Cloned copy must be a PRefineableQElement<2,INITIAL_NNODE_1D>!",
2022 OOMPH_CURRENT_FUNCTION,
2023 OOMPH_EXCEPTION_LOCATION);
2024 }
2025#ifdef PARANOID
2026 // Clone exists, so check that it is infact a copy of me
2027 else
2028 {
2029 // Flag to keep track of check
2030 bool clone_is_ok = true;
2031
2032 // Does the clone have the correct p-order?
2033 clone_is_ok = clone_is_ok && (clone_el_pt->p_order() == this->p_order());
2034
2035 if (!clone_is_ok)
2036 {
2037 std::ostringstream err_stream;
2038 err_stream << "Clone element has a different p-order from me,"
2039 << std::endl
2040 << "but it is supposed to be a copy!" << std::endl;
2041
2042 throw OomphLibError(
2043 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2044 }
2045
2046 // Does the clone have the same number of nodes as me?
2047 clone_is_ok = clone_is_ok && (clone_el_pt->nnode() == this->nnode());
2048
2049 if (!clone_is_ok)
2050 {
2051 std::ostringstream err_stream;
2052 err_stream << "Clone element has a different number of nodes from me,"
2053 << std::endl
2054 << "but it is supposed to be a copy!" << std::endl;
2055
2056 throw OomphLibError(
2057 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2058 }
2059
2060 // Does the clone have the same nodes as me?
2061 for (unsigned n = 0; n < this->nnode(); n++)
2062 {
2063 clone_is_ok =
2064 clone_is_ok && (clone_el_pt->node_pt(n) == this->node_pt(n));
2065 }
2066
2067 if (!clone_is_ok)
2068 {
2069 std::ostringstream err_stream;
2070 err_stream << "The nodes belonging to the clone element are different"
2071 << std::endl
2072 << "from mine, but it is supposed to be a copy!"
2073 << std::endl;
2074
2075 throw OomphLibError(
2076 err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2077 }
2078
2079 // Is this a macro element?
2081 dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
2082 if (m_el_pt != 0)
2083 {
2084 // Get macro element version of clone
2085 MacroElementNodeUpdateElementBase* clone_m_el_pt =
2086 dynamic_cast<MacroElementNodeUpdateElementBase*>(clone_el_pt);
2087
2088 // Does the clone have the same geometric objects?
2089 Vector<GeomObject*> clone_geom_object_pt(
2090 clone_m_el_pt->geom_object_pt());
2091 Vector<GeomObject*> geom_object_pt(m_el_pt->geom_object_pt());
2092
2093 clone_is_ok =
2094 clone_is_ok && (geom_object_pt.size() == clone_geom_object_pt.size());
2095 for (unsigned n = 0;
2096 n < std::min(geom_object_pt.size(), clone_geom_object_pt.size());
2097 n++)
2098 {
2099 clone_is_ok =
2100 clone_is_ok && (clone_geom_object_pt[n] == geom_object_pt[n]);
2101 }
2102
2103 if (!clone_is_ok)
2104 {
2105 std::ostringstream err_stream;
2106 err_stream << "The clone element has different geometric objects"
2107 << std::endl
2108 << "from mine, but it is supposed to be a copy!"
2109 << std::endl;
2110
2111 throw OomphLibError(
2112 err_stream.str(),
2113 "PRefineableQElement<2,INITIAL_NNODE_1D>::p_refine()",
2114 OOMPH_EXCEPTION_LOCATION);
2115 }
2116 }
2117
2118 // If we get to here then the clone has all the information we require
2119 }
2120#endif
2121
2122 // Currently we can't handle the case of generalised coordinates
2123 // since we haven't established how they should be interpolated.
2124 // Buffer this case:
2125 if (clone_el_pt->node_pt(0)->nposition_type() != 1)
2126 {
2127 throw OomphLibError("Can't handle generalised nodal positions (yet).",
2128 OOMPH_CURRENT_FUNCTION,
2129 OOMPH_EXCEPTION_LOCATION);
2130 }
2131
2132 // Timestepper should be the same for all nodes -- use it
2133 // to create timesteppers for new nodes
2134 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
2135
2136 // Get number of history values (incl. present)
2137 unsigned ntstorage = time_stepper_pt->ntstorage();
2138
2139 // Increment p-order of the element
2140 P_order += inc;
2141
2142 // Change integration scheme
2143 delete this->integral_pt();
2144 switch (P_order)
2145 {
2146 case 2:
2147 this->set_integration_scheme(new GaussLobattoLegendre<2, 2>);
2148 break;
2149 case 3:
2150 this->set_integration_scheme(new GaussLobattoLegendre<2, 3>);
2151 break;
2152 case 4:
2153 this->set_integration_scheme(new GaussLobattoLegendre<2, 4>);
2154 break;
2155 case 5:
2156 this->set_integration_scheme(new GaussLobattoLegendre<2, 5>);
2157 break;
2158 case 6:
2159 this->set_integration_scheme(new GaussLobattoLegendre<2, 6>);
2160 break;
2161 case 7:
2162 this->set_integration_scheme(new GaussLobattoLegendre<2, 7>);
2163 break;
2164 default:
2165 std::ostringstream error_message;
2166 error_message << "\nERROR: Exceeded maximum polynomial order for";
2167 error_message << "\n integration scheme.\n";
2168 throw OomphLibError(error_message.str(),
2169 OOMPH_CURRENT_FUNCTION,
2170 OOMPH_EXCEPTION_LOCATION);
2171 }
2172
2173 // Allocate new space for Nodes (at the element level)
2174 this->set_n_node(P_order * P_order);
2175
2176 // Copy vertex nodes and create new edge and internal nodes
2177 //---------------------------------------------------------
2178
2179 // Setup vertex coordinates in element:
2180 //-------------------------------------
2181 Vector<double> s_lo(2);
2182 Vector<double> s_hi(2);
2183 s_lo[0] = -1.0;
2184 s_hi[0] = 1.0;
2185 s_lo[1] = -1.0;
2186 s_hi[1] = 1.0;
2187
2188 // Local coordinate in element
2189 Vector<double> s(2);
2190
2191 unsigned jnod = 0;
2192
2193 Vector<double> s_fraction(2);
2194 // Loop over nodes in element
2195 for (unsigned i0 = 0; i0 < P_order; i0++)
2196 {
2197 // Get the fractional position of the node in the direction of s[0]
2198 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
2199 // Local coordinate
2200 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
2201
2202 for (unsigned i1 = 0; i1 < P_order; i1++)
2203 {
2204 // Get the fractional position of the node in the direction of s[1]
2205 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
2206 // Local coordinate
2207 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
2208
2209 // Local node number
2210 jnod = i0 + P_order * i1;
2211
2212 // Initialise flag: So far, this node hasn't been built
2213 // or copied yet
2214 bool node_done = false;
2215
2216 // Get the pointer to the node in this element (or rather, its clone),
2217 // returns NULL if there is not node
2218 Node* created_node_pt = clone_el_pt->get_node_at_local_coordinate(s);
2219
2220 // Does this node already exist in this element?
2221 //----------------------------------------------
2222 if (created_node_pt != 0)
2223 {
2224 // Copy node across
2225 this->node_pt(jnod) = created_node_pt;
2226
2227 // Make sure that we update the values at the node so that
2228 // they are consistent with the present representation.
2229 // This is only need for mixed interpolation where the value
2230 // at the father could now become active.
2231
2232 // Loop over all history values
2233 for (unsigned t = 0; t < ntstorage; t++)
2234 {
2235 // Get values from father element
2236 // Note: get_interpolated_values() sets Vector size itself.
2237 Vector<double> prev_values;
2238 clone_el_pt->get_interpolated_values(t, s, prev_values);
2239 // Find the minimum number of values
2240 //(either those stored at the node, or those returned by
2241 // the function)
2242 unsigned n_val_at_node = created_node_pt->nvalue();
2243 unsigned n_val_from_function = prev_values.size();
2244 // Use the ternary conditional operator here
2245 unsigned n_var = n_val_at_node < n_val_from_function ?
2246 n_val_at_node :
2247 n_val_from_function;
2248 // Assign the values that we can
2249 for (unsigned k = 0; k < n_var; k++)
2250 {
2251 created_node_pt->set_value(t, k, prev_values[k]);
2252 }
2253 }
2254
2255 // Node has been created by copy
2256 node_done = true;
2257 }
2258 // Node does not exist in this element but might already
2259 //------------------------------------------------------
2260 // have been created by neighbouring elements
2261 //-------------------------------------------
2262 else
2263 {
2264 // Was the node created by one of its neighbours
2265 // Whether or not the node lies on an edge can be calculated
2266 // by from the fractional position
2267 bool is_periodic = false;
2268 created_node_pt = node_created_by_neighbour(s_fraction, is_periodic);
2269
2270 // If the node was so created, assign the pointers
2271 if (created_node_pt != 0)
2272 {
2273 // If the node is periodic
2274 if (is_periodic)
2275 {
2276 // Now the node must be on a boundary, but we don't know which
2277 // one
2278 // The returned created_node_pt is actually the neighbouring
2279 // periodic node
2280 Node* neighbour_node_pt = created_node_pt;
2281
2282 // Determine the edge on which the new node will live
2283 //(cannot be a vertex node)
2284 int my_bound = Tree::OMEGA;
2285 if (s_fraction[0] == 0.0) my_bound = QuadTreeNames::W;
2286 else if (s_fraction[0] == 1.0)
2287 my_bound = QuadTreeNames::E;
2288 else if (s_fraction[1] == 0.0)
2289 my_bound = QuadTreeNames::S;
2290 else if (s_fraction[1] == 1.0)
2291 my_bound = QuadTreeNames::N;
2292
2293 // Storage for the set of Mesh boundaries on which the
2294 // appropriate edge lives.
2295 // [New nodes should always be mid-edge nodes and therefore
2296 // only live on one boundary but just to play it safe...]
2297 std::set<unsigned> boundaries;
2298 // Only get the boundaries if we are at the edge of
2299 // an element. Nodes in the centre of an element cannot be
2300 // on Mesh boundaries
2301 if (my_bound != Tree::OMEGA)
2302 {
2303 clone_el_pt->get_boundaries(my_bound, boundaries);
2304 }
2305
2306#ifdef PARANOID
2307 // Case where a new node lives on more than one boundary
2308 // seems fishy enough to flag
2309 if (boundaries.size() > 1)
2310 {
2311 throw OomphLibError(
2312 "boundaries.size()!=1 seems a bit strange..\n",
2313 OOMPH_CURRENT_FUNCTION,
2314 OOMPH_EXCEPTION_LOCATION);
2315 }
2316
2317 // Case when there are no boundaries, we are in big trouble
2318 if (boundaries.size() == 0)
2319 {
2320 std::ostringstream error_stream;
2321 error_stream << "Periodic node is not on a boundary...\n"
2322 << "Coordinates: " << created_node_pt->x(0) << " "
2323 << created_node_pt->x(1) << "\n";
2324 throw OomphLibError(error_stream.str(),
2325 OOMPH_CURRENT_FUNCTION,
2326 OOMPH_EXCEPTION_LOCATION);
2327 }
2328#endif
2329
2330 // Create node and set the pointer to it from the element
2331 created_node_pt =
2332 this->construct_boundary_node(jnod, time_stepper_pt);
2333 // Make the node periodic from the neighbour
2334 created_node_pt->make_periodic(neighbour_node_pt);
2335
2336 // Loop over # of history values
2337 for (unsigned t = 0; t < ntstorage; t++)
2338 {
2339 // Get position from father element -- this uses the macro
2340 // element representation if appropriate. If the node
2341 // turns out to be a hanging node later on, then
2342 // its position gets adjusted in line with its
2343 // hanging node interpolation.
2344 Vector<double> x_prev(2);
2345 clone_el_pt->get_x(t, s, x_prev);
2346 // Set previous positions of the new node
2347 for (unsigned i = 0; i < 2; i++)
2348 {
2349 created_node_pt->x(t, i) = x_prev[i];
2350 }
2351 }
2352
2353 // Check if we need to add nodes to the mesh
2354 if (mesh_pt != 0)
2355 {
2356 // Next, we Update the boundary lookup schemes
2357 // Loop over the boundaries stored in the set
2358 for (std::set<unsigned>::iterator it = boundaries.begin();
2359 it != boundaries.end();
2360 ++it)
2361 {
2362 // Add the node to the boundary
2363 mesh_pt->add_boundary_node(*it, created_node_pt);
2364
2365 // If we have set an intrinsic coordinate on this
2366 // mesh boundary then it must also be interpolated on
2367 // the new node
2368 // Now interpolate the intrinsic boundary coordinate
2369 if (mesh_pt->boundary_coordinate_exists(*it) == true)
2370 {
2371 Vector<double> zeta(1);
2372 clone_el_pt->interpolated_zeta_on_edge(
2373 *it, my_bound, s, zeta);
2374
2375 created_node_pt->set_coordinates_on_boundary(*it, zeta);
2376 }
2377 }
2378
2379 // Make sure that we add the node to the mesh
2380 mesh_pt->add_node_pt(created_node_pt);
2381 }
2382 } // End of periodic case
2383 // Otherwise the node is not periodic, so just set the
2384 // pointer to the neighbours node
2385 else
2386 {
2387 this->node_pt(jnod) = created_node_pt;
2388 }
2389 // Node has been created
2390 node_done = true;
2391 }
2392 // Node does not exist in neighbour element but might already
2393 //-----------------------------------------------------------
2394 // have been created by a son of a neighbouring element
2395 //-----------------------------------------------------
2396 else
2397 {
2398 // Was the node created by one of its neighbours' sons
2399 // Whether or not the node lies on an edge can be calculated
2400 // by from the fractional position
2401 bool is_periodic = false;
2402 created_node_pt =
2403 node_created_by_son_of_neighbour(s_fraction, is_periodic);
2404
2405 // If the node was so created, assign the pointers
2406 if (created_node_pt != 0)
2407 {
2408 // If the node is periodic
2409 if (is_periodic)
2410 {
2411 // Now the node must be on a boundary, but we don't know which
2412 // one
2413 // The returned created_node_pt is actually the neighbouring
2414 // periodic node
2415 Node* neighbour_node_pt = created_node_pt;
2416
2417 // Determine the edge on which the new node will live
2418 //(cannot be a vertex node)
2419 int my_bound = Tree::OMEGA;
2420 if (s_fraction[0] == 0.0) my_bound = QuadTreeNames::W;
2421 else if (s_fraction[0] == 1.0)
2422 my_bound = QuadTreeNames::E;
2423 else if (s_fraction[1] == 0.0)
2424 my_bound = QuadTreeNames::S;
2425 else if (s_fraction[1] == 1.0)
2426 my_bound = QuadTreeNames::N;
2427
2428 // Storage for the set of Mesh boundaries on which the
2429 // appropriate edge lives.
2430 // [New nodes should always be mid-edge nodes and therefore
2431 // only live on one boundary but just to play it safe...]
2432 std::set<unsigned> boundaries;
2433 // Only get the boundaries if we are at the edge of
2434 // an element. Nodes in the centre of an element cannot be
2435 // on Mesh boundaries
2436 if (my_bound != Tree::OMEGA)
2437 {
2438 clone_el_pt->get_boundaries(my_bound, boundaries);
2439 }
2440
2441#ifdef PARANOID
2442 // Case where a new node lives on more than one boundary
2443 // seems fishy enough to flag
2444 if (boundaries.size() > 1)
2445 {
2446 throw OomphLibError(
2447 "boundaries.size()!=1 seems a bit strange..\n",
2448 OOMPH_CURRENT_FUNCTION,
2449 OOMPH_EXCEPTION_LOCATION);
2450 }
2451
2452 // Case when there are no boundaries, we are in big trouble
2453 if (boundaries.size() == 0)
2454 {
2455 std::ostringstream error_stream;
2456 error_stream << "Periodic node is not on a boundary...\n"
2457 << "Coordinates: " << created_node_pt->x(0)
2458 << " " << created_node_pt->x(1) << "\n";
2459 throw OomphLibError(error_stream.str(),
2460 OOMPH_CURRENT_FUNCTION,
2461 OOMPH_EXCEPTION_LOCATION);
2462 }
2463#endif
2464
2465 // Create node and set the pointer to it from the element
2466 created_node_pt =
2467 this->construct_boundary_node(jnod, time_stepper_pt);
2468 // Make the node periodic from the neighbour
2469 created_node_pt->make_periodic(neighbour_node_pt);
2470
2471 // Loop over # of history values
2472 for (unsigned t = 0; t < ntstorage; t++)
2473 {
2474 // Get position from father element -- this uses the macro
2475 // element representation if appropriate. If the node
2476 // turns out to be a hanging node later on, then
2477 // its position gets adjusted in line with its
2478 // hanging node interpolation.
2479 Vector<double> x_prev(2);
2480 clone_el_pt->get_x(t, s, x_prev);
2481 // Set previous positions of the new node
2482 for (unsigned i = 0; i < 2; i++)
2483 {
2484 created_node_pt->x(t, i) = x_prev[i];
2485 }
2486 }
2487
2488 // Check if we need to add nodes to the mesh
2489 if (mesh_pt != 0)
2490 {
2491 // Next, we Update the boundary lookup schemes
2492 // Loop over the boundaries stored in the set
2493 for (std::set<unsigned>::iterator it = boundaries.begin();
2494 it != boundaries.end();
2495 ++it)
2496 {
2497 // Add the node to the boundary
2498 mesh_pt->add_boundary_node(*it, created_node_pt);
2499
2500 // If we have set an intrinsic coordinate on this
2501 // mesh boundary then it must also be interpolated on
2502 // the new node
2503 // Now interpolate the intrinsic boundary coordinate
2504 if (mesh_pt->boundary_coordinate_exists(*it) == true)
2505 {
2506 Vector<double> zeta(1);
2507 clone_el_pt->interpolated_zeta_on_edge(
2508 *it, my_bound, s, zeta);
2509
2510 created_node_pt->set_coordinates_on_boundary(*it, zeta);
2511 }
2512 }
2513
2514 // Make sure that we add the node to the mesh
2515 mesh_pt->add_node_pt(created_node_pt);
2516 }
2517 } // End of periodic case
2518 // Otherwise the node is not periodic, so just set the
2519 // pointer to the neighbours node
2520 else
2521 {
2522 this->node_pt(jnod) = created_node_pt;
2523 }
2524 // Node has been created
2525 node_done = true;
2526 } // Node does not exist in son of neighbouring element
2527 } // Node does not exist in neighbouring element
2528 } // Node does not exist in this element
2529
2530 // Node has not been built anywhere ---> build it here
2531 if (!node_done)
2532 {
2533 // Firstly, we need to determine whether or not a node lies
2534 // on the boundary before building it, because
2535 // we actually assign a different type of node on boundaries.
2536
2537 // Determine the edge on which the new node will live
2538 //(cannot be a vertex node)
2539 int my_bound = Tree::OMEGA;
2540 if (s_fraction[0] == 0.0) my_bound = QuadTreeNames::W;
2541 else if (s_fraction[0] == 1.0)
2542 my_bound = QuadTreeNames::E;
2543 else if (s_fraction[1] == 0.0)
2544 my_bound = QuadTreeNames::S;
2545 else if (s_fraction[1] == 1.0)
2546 my_bound = QuadTreeNames::N;
2547
2548 // Storage for the set of Mesh boundaries on which the
2549 // appropriate edge lives.
2550 // [New nodes should always be mid-edge nodes and therefore
2551 // only live on one boundary but just to play it safe...]
2552 std::set<unsigned> boundaries;
2553 // Only get the boundaries if we are at the edge of
2554 // an element. Nodes in the centre of an element cannot be
2555 // on Mesh boundaries
2556 if (my_bound != Tree::OMEGA)
2557 {
2558 clone_el_pt->get_boundaries(my_bound, boundaries);
2559 }
2560
2561#ifdef PARANOID
2562 // Case where a new node lives on more than one boundary
2563 // seems fishy enough to flag
2564 if (boundaries.size() > 1)
2565 {
2566 throw OomphLibError("boundaries.size()!=1 seems a bit strange..\n",
2567 OOMPH_CURRENT_FUNCTION,
2568 OOMPH_EXCEPTION_LOCATION);
2569 }
2570#endif
2571
2572 // If the node lives on a mesh boundary,
2573 // then we need to create a boundary node
2574 if (boundaries.size() > 0)
2575 {
2576 // Create node and set the pointer to it from the element
2577 created_node_pt =
2578 this->construct_boundary_node(jnod, time_stepper_pt);
2579
2580 // Now we need to work out whether to pin the values at
2581 // the new node based on the boundary conditions applied at
2582 // its Mesh boundary
2583
2584 // Get the boundary conditions from the father
2585 Vector<int> bound_cons(this->ncont_interpolated_values());
2586 clone_el_pt->get_bcs(my_bound, bound_cons);
2587
2588 // Loop over the values and pin, if necessary
2589 unsigned n_value = created_node_pt->nvalue();
2590 for (unsigned k = 0; k < n_value; k++)
2591 {
2592 if (bound_cons[k])
2593 {
2594 created_node_pt->pin(k);
2595 }
2596 }
2597
2598 // Solid node? If so, deal with the positional boundary
2599 // conditions:
2600 SolidNode* solid_node_pt =
2601 dynamic_cast<SolidNode*>(created_node_pt);
2602 if (solid_node_pt != 0)
2603 {
2604 // Get the positional boundary conditions from the father:
2605 unsigned n_dim = created_node_pt->ndim();
2606 Vector<int> solid_bound_cons(n_dim);
2607 RefineableSolidQElement<2>* clone_solid_el_pt =
2608 dynamic_cast<RefineableSolidQElement<2>*>(clone_el_pt);
2609#ifdef PARANOID
2610 if (clone_solid_el_pt == 0)
2611 {
2612 std::string error_message =
2613 "We have a SolidNode outside a refineable SolidElement\n";
2614 error_message +=
2615 "during mesh refinement -- this doesn't make sense";
2616
2617 throw OomphLibError(error_message,
2618 OOMPH_CURRENT_FUNCTION,
2619 OOMPH_EXCEPTION_LOCATION);
2620 }
2621#endif
2622 clone_solid_el_pt->get_solid_bcs(my_bound, solid_bound_cons);
2623
2624 // Loop over the positions and pin, if necessary
2625 for (unsigned k = 0; k < n_dim; k++)
2626 {
2627 if (solid_bound_cons[k])
2628 {
2629 solid_node_pt->pin_position(k);
2630 }
2631 }
2632 } // End of if solid_node_pt
2633
2634
2635 // Check if we need to add nodes to the mesh
2636 if (mesh_pt != 0)
2637 {
2638 // Next, we Update the boundary lookup schemes
2639 // Loop over the boundaries stored in the set
2640 for (std::set<unsigned>::iterator it = boundaries.begin();
2641 it != boundaries.end();
2642 ++it)
2643 {
2644 // Add the node to the boundary
2645 mesh_pt->add_boundary_node(*it, created_node_pt);
2646
2647 // If we have set an intrinsic coordinate on this
2648 // mesh boundary then it must also be interpolated on
2649 // the new node
2650 // Now interpolate the intrinsic boundary coordinate
2651 if (mesh_pt->boundary_coordinate_exists(*it) == true)
2652 {
2653 Vector<double> zeta(1);
2654 clone_el_pt->interpolated_zeta_on_edge(
2655 *it, my_bound, s, zeta);
2656
2657 created_node_pt->set_coordinates_on_boundary(*it, zeta);
2658 }
2659 }
2660 }
2661 }
2662 // Otherwise the node is not on a Mesh boundary and
2663 // we create a normal "bulk" node
2664 else
2665 {
2666 // Create node and set the pointer to it from the element
2667 created_node_pt = this->construct_node(jnod, time_stepper_pt);
2668 }
2669
2670 // Now we set the position and values at the newly created node
2671
2672 // In the first instance use macro element or FE representation
2673 // to create past and present nodal positions.
2674 // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
2675 // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
2676 // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
2677 // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
2678 // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
2679 // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
2680 // NOT ASSIGN SENSIBLE INITIAL POSITONS!
2681
2682 // Loop over # of history values
2683 for (unsigned t = 0; t < ntstorage; t++)
2684 {
2685 // Get position from father element -- this uses the macro
2686 // element representation if appropriate. If the node
2687 // turns out to be a hanging node later on, then
2688 // its position gets adjusted in line with its
2689 // hanging node interpolation.
2690 Vector<double> x_prev(2);
2691 clone_el_pt->get_x(t, s, x_prev);
2692
2693 // Set previous positions of the new node
2694 for (unsigned i = 0; i < 2; i++)
2695 {
2696 created_node_pt->x(t, i) = x_prev[i];
2697 }
2698 }
2699
2700 // Loop over all history values
2701 for (unsigned t = 0; t < ntstorage; t++)
2702 {
2703 // Get values from father element
2704 // Note: get_interpolated_values() sets Vector size itself.
2705 Vector<double> prev_values;
2706 clone_el_pt->get_interpolated_values(t, s, prev_values);
2707 // Initialise the values at the new node
2708 unsigned n_value = created_node_pt->nvalue();
2709 for (unsigned k = 0; k < n_value; k++)
2710 {
2711 created_node_pt->set_value(t, k, prev_values[k]);
2712 }
2713 }
2714
2715 // Add new node to mesh (if requested)
2716 if (mesh_pt != 0)
2717 {
2718 mesh_pt->add_node_pt(created_node_pt);
2719 }
2720
2721 AlgebraicElementBase* alg_el_pt =
2722 dynamic_cast<AlgebraicElementBase*>(this);
2723
2724 // If we do have an algebraic element
2725 if (alg_el_pt != 0)
2726 {
2727 std::string error_message = "Have not implemented p-refinement for";
2728 error_message += "Algebraic p-refineable elements yet\n";
2729
2730 throw OomphLibError(
2731 error_message,
2732 "PRefineableQElement<2,INITIAL_NNODE_1D>::p_refine()",
2733 OOMPH_EXCEPTION_LOCATION);
2734 }
2735
2736 } // End of case when we build the node ourselves
2737
2738 // Check if the element is an algebraic element
2739 AlgebraicElementBase* alg_el_pt =
2740 dynamic_cast<AlgebraicElementBase*>(this);
2741
2742 // If the element is an algebraic element, setup
2743 // node position (past and present) from algebraic node update
2744 // function. This over-writes previous assingments that
2745 // were made based on the macro-element/FE representation.
2746 // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
2747 // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
2748 // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
2749 // INFO FOR *ALL* ROOT ELEMENTS!
2750 if (alg_el_pt != 0)
2751 {
2752 // Build algebraic node update info for new node
2753 // This sets up the node update data for all node update
2754 // functions that are shared by all nodes in the father
2755 // element
2756 alg_el_pt->setup_algebraic_node_update(
2757 this->node_pt(jnod), s, clone_el_pt);
2758 }
2759
2760 } // End of vertical loop over nodes in element
2761
2762 } // End of horizontal loop over nodes in element
2763
2764
2765 // If the element is a MacroElementNodeUpdateElement, set
2766 // the update parameters for the current element's nodes --
2767 // all this needs is the vector of (pointers to the)
2768 // geometric objects that affect the MacroElement-based
2769 // node update -- this needs to be done to set the node
2770 // update info for newly created nodes
2771 MacroElementNodeUpdateElementBase* clone_m_el_pt =
2772 dynamic_cast<MacroElementNodeUpdateElementBase*>(clone_el_pt);
2773 if (clone_m_el_pt != 0)
2774 {
2775 // Get vector of geometric objects from father (construct vector
2776 // via copy operation)
2777 Vector<GeomObject*> geom_object_pt(clone_m_el_pt->geom_object_pt());
2778
2779 // Cast current element to MacroElementNodeUpdateElement:
2781 dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
2782
2783#ifdef PARANOID
2784 if (m_el_pt == 0)
2785 {
2786 std::string error_message =
2787 "Failed to cast to MacroElementNodeUpdateElementBase*\n";
2788 error_message +=
2789 "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
2790 error_message += "then I should be too....\n";
2791
2792 throw OomphLibError(
2793 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2794 }
2795#endif
2796 // Build update info by passing vector of geometric objects:
2797 // This sets the current element to be the update element
2798 // for all of the element's nodes -- this is reversed
2799 // if the element is ever un-refined in the father element's
2800 // rebuild_from_sons() function which overwrites this
2801 // assignment to avoid nasty segmentation faults that occur
2802 // when a node tries to update itself via an element that no
2803 // longer exists...
2804 m_el_pt->set_node_update_info(geom_object_pt);
2805
2806 /// / Now loop over nodes in element
2807 // unsigned n_node = this->nnode();
2808 // for (unsigned j=0;j<n_node;j++)
2809 // {
2810 // // Get local coordinate in element (Vector sets its own size)
2811 // Vector<double> s_in_node_update_element;
2812 // this->local_coordinate_of_node(j,s_in_node_update_element);
2813 //
2814 // // Pass the lot to the node
2815 // static_cast<MacroElementNodeUpdateNode*>(this->node_pt(j))->
2816 // set_node_update_info(this,s_in_node_update_element,m_el_pt->geom_object_pt());
2817 // }
2818
2819 /// /BENFLAG:
2820 // std::cout << "Checking that all the nodes have this as their update
2821 // element..." << std::endl;
2822 /// /std::cout << "this = " << this << std::endl;
2823 // for(unsigned j=0; j<this->nnode(); j++)
2824 // {
2825 // //std::cout << this->node_pt(j) << ": [" << this->node_pt(j)->x(0)
2826 // << "," << this->node_pt(j)->x(1) << "] update element: " <<
2827 // dynamic_cast<MacroElementNodeUpdateNode*>(this->node_pt(j))->node_update_element_pt()
2828 // << std::endl; MacroElementNodeUpdateNode* mac_nod_pt =
2829 // dynamic_cast<MacroElementNodeUpdateNode*>(this->node_pt(j));
2830 // if(mac_nod_pt->node_update_element_pt()!=this)
2831 // {
2832 // std::cout << "Something's not right! Update element is wrong..." <<
2833 // std::endl;
2834 // }
2835 // FiniteElement* up_el_pt =
2836 // dynamic_cast<FiniteElement*>(mac_nod_pt->node_update_element_pt());
2837 // bool not_good = true;
2838 // for(unsigned l=0; l<up_el_pt->nnode(); l++)
2839 // {
2840 // if(up_el_pt->node_pt(l)==mac_nod_pt)
2841 // {
2842 // not_good = false;
2843 // break;
2844 // }
2845 // }
2846 // if(not_good==true)
2847 // {
2848 // std::cout << "Macro node doesn't belong to its update element!" <<
2849 // std::endl;
2850 // }
2851 // }
2852
2853
2854 // Loop over all nodes in element again, to re-set the positions
2855 // This must be done using the new element's macro-element
2856 // representation, rather than the old version which may be
2857 // of a different p-order!
2858 for (unsigned i0 = 0; i0 < P_order; i0++)
2859 {
2860 // Get the fractional position of the node in the direction of s[0]
2861 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
2862 // Local coordinate
2863 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
2864
2865 for (unsigned i1 = 0; i1 < P_order; i1++)
2866 {
2867 // Get the fractional position of the node in the direction of s[1]
2868 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
2869 // Local coordinate
2870 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
2871
2872 // Local node number
2873 jnod = i0 + P_order * i1;
2874
2875 // Loop over # of history values
2876 for (unsigned t = 0; t < ntstorage; t++)
2877 {
2878 // Get position from father element -- this uses the macro
2879 // element representation if appropriate. If the node
2880 // turns out to be a hanging node later on, then
2881 // its position gets adjusted in line with its
2882 // hanging node interpolation.
2883 Vector<double> x_prev(2);
2884 this->get_x(t, s, x_prev);
2885
2886 // Set previous positions of the new node
2887 for (unsigned i = 0; i < 2; i++)
2888 {
2889 this->node_pt(jnod)->x(t, i) = x_prev[i];
2890 }
2891 }
2892 }
2893 }
2894 }
2895
2896 // Not necessary to delete the old nodes since all original nodes are in the
2897 // current mesh and so will be pruned as part of the mesh adaption process.
2898
2899 // Do any further-build required
2900 this->further_build();
2901 }
2902
2903 //=======================================================================
2904 /// Shape functions for PRefineableQElement<DIM>
2905 //=======================================================================
2906 template<unsigned INITIAL_NNODE_1D>
2908 Shape& psi) const
2909 {
2910 switch (p_order())
2911 {
2912 case 2:
2913 {
2914 // Call the OneDimensional Shape functions
2916 OneDimensionalLegendreShape<2> psi1(s[0]), psi2(s[1]);
2917
2918 // Now let's loop over the nodal points in the element
2919 // and copy the values back in
2920 for (unsigned i = 0; i < 2; i++)
2921 {
2922 for (unsigned j = 0; j < 2; j++)
2923 {
2924 psi(2 * i + j) = psi2[i] * psi1[j];
2925 }
2926 }
2927 break;
2928 }
2929 case 3:
2930 {
2931 // Call the OneDimensional Shape functions
2933 OneDimensionalLegendreShape<3> psi1(s[0]), psi2(s[1]);
2934
2935 // Now let's loop over the nodal points in the element
2936 // and copy the values back in
2937 for (unsigned i = 0; i < 3; i++)
2938 {
2939 for (unsigned j = 0; j < 3; j++)
2940 {
2941 psi(3 * i + j) = psi2[i] * psi1[j];
2942 }
2943 }
2944 break;
2945 }
2946 case 4:
2947 {
2948 // Call the OneDimensional Shape functions
2950 OneDimensionalLegendreShape<4> psi1(s[0]), psi2(s[1]);
2951
2952 // Now let's loop over the nodal points in the element
2953 // and copy the values back in
2954 for (unsigned i = 0; i < 4; i++)
2955 {
2956 for (unsigned j = 0; j < 4; j++)
2957 {
2958 psi(4 * i + j) = psi2[i] * psi1[j];
2959 }
2960 }
2961 break;
2962 }
2963 case 5:
2964 {
2965 // Call the OneDimensional Shape functions
2967 OneDimensionalLegendreShape<5> psi1(s[0]), psi2(s[1]);
2968
2969 // Now let's loop over the nodal points in the element
2970 // and copy the values back in
2971 for (unsigned i = 0; i < 5; i++)
2972 {
2973 for (unsigned j = 0; j < 5; j++)
2974 {
2975 psi(5 * i + j) = psi2[i] * psi1[j];
2976 }
2977 }
2978 break;
2979 }
2980 case 6:
2981 {
2982 // Call the OneDimensional Shape functions
2984 OneDimensionalLegendreShape<6> psi1(s[0]), psi2(s[1]);
2985
2986 // Now let's loop over the nodal points in the element
2987 // and copy the values back in
2988 for (unsigned i = 0; i < 6; i++)
2989 {
2990 for (unsigned j = 0; j < 6; j++)
2991 {
2992 psi(6 * i + j) = psi2[i] * psi1[j];
2993 }
2994 }
2995 break;
2996 }
2997 case 7:
2998 {
2999 // Call the OneDimensional Shape functions
3001 OneDimensionalLegendreShape<7> psi1(s[0]), psi2(s[1]);
3002
3003 // Now let's loop over the nodal points in the element
3004 // and copy the values back in
3005 for (unsigned i = 0; i < 7; i++)
3006 {
3007 for (unsigned j = 0; j < 7; j++)
3008 {
3009 psi(7 * i + j) = psi2[i] * psi1[j];
3010 }
3011 }
3012 break;
3013 }
3014 default:
3015 std::ostringstream error_message;
3016 error_message << "\nERROR: Exceeded maximum polynomial order for";
3017 error_message << "\n polynomial order for shape functions.\n";
3018 throw OomphLibError(error_message.str(),
3019 OOMPH_CURRENT_FUNCTION,
3020 OOMPH_EXCEPTION_LOCATION);
3021 }
3022 }
3023
3024 //=======================================================================
3025 /// Derivatives of shape functions for PRefineableQElement<DIM>
3026 //=======================================================================
3027 template<unsigned INITIAL_NNODE_1D>
3029 const Vector<double>& s, Shape& psi, DShape& dpsids) const
3030 {
3031 switch (p_order())
3032 {
3033 case 2:
3034 {
3035 // Call the shape functions and derivatives
3037 OneDimensionalLegendreShape<2> psi1(s[0]), psi2(s[1]);
3038 OneDimensionalLegendreDShape<2> dpsi1ds(s[0]), dpsi2ds(s[1]);
3039
3040 // Index for the shape functions
3041 unsigned index = 0;
3042 // Loop over shape functions in element
3043 for (unsigned i = 0; i < 2; i++)
3044 {
3045 for (unsigned j = 0; j < 2; j++)
3046 {
3047 // Assign the values
3048 dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3049 dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3050 psi[index] = psi2[i] * psi1[j];
3051 // Increment the index
3052 ++index;
3053 }
3054 }
3055 break;
3056 }
3057 case 3:
3058 {
3059 // Call the shape functions and derivatives
3061 OneDimensionalLegendreShape<3> psi1(s[0]), psi2(s[1]);
3062 OneDimensionalLegendreDShape<3> dpsi1ds(s[0]), dpsi2ds(s[1]);
3063
3064 // Index for the shape functions
3065 unsigned index = 0;
3066 // Loop over shape functions in element
3067 for (unsigned i = 0; i < 3; i++)
3068 {
3069 for (unsigned j = 0; j < 3; j++)
3070 {
3071 // Assign the values
3072 dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3073 dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3074 psi[index] = psi2[i] * psi1[j];
3075 // Increment the index
3076 ++index;
3077 }
3078 }
3079 break;
3080 }
3081 case 4:
3082 {
3083 // Call the shape functions and derivatives
3085 OneDimensionalLegendreShape<4> psi1(s[0]), psi2(s[1]);
3086 OneDimensionalLegendreDShape<4> dpsi1ds(s[0]), dpsi2ds(s[1]);
3087
3088 // Index for the shape functions
3089 unsigned index = 0;
3090 // Loop over shape functions in element
3091 for (unsigned i = 0; i < 4; i++)
3092 {
3093 for (unsigned j = 0; j < 4; j++)
3094 {
3095 // Assign the values
3096 dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3097 dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3098 psi[index] = psi2[i] * psi1[j];
3099 // Increment the index
3100 ++index;
3101 }
3102 }
3103 break;
3104 }
3105 case 5:
3106 {
3107 // Call the shape functions and derivatives
3109 OneDimensionalLegendreShape<5> psi1(s[0]), psi2(s[1]);
3110 OneDimensionalLegendreDShape<5> dpsi1ds(s[0]), dpsi2ds(s[1]);
3111
3112 // Index for the shape functions
3113 unsigned index = 0;
3114 // Loop over shape functions in element
3115 for (unsigned i = 0; i < 5; i++)
3116 {
3117 for (unsigned j = 0; j < 5; j++)
3118 {
3119 // Assign the values
3120 dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3121 dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3122 psi[index] = psi2[i] * psi1[j];
3123 // Increment the index
3124 ++index;
3125 }
3126 }
3127 break;
3128 }
3129 case 6:
3130 {
3131 // Call the shape functions and derivatives
3133 OneDimensionalLegendreShape<6> psi1(s[0]), psi2(s[1]);
3134 OneDimensionalLegendreDShape<6> dpsi1ds(s[0]), dpsi2ds(s[1]);
3135
3136 // Index for the shape functions
3137 unsigned index = 0;
3138 // Loop over shape functions in element
3139 for (unsigned i = 0; i < 6; i++)
3140 {
3141 for (unsigned j = 0; j < 6; j++)
3142 {
3143 // Assign the values
3144 dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3145 dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3146 psi[index] = psi2[i] * psi1[j];
3147 // Increment the index
3148 ++index;
3149 }
3150 }
3151 break;
3152 }
3153 case 7:
3154 {
3155 // Call the shape functions and derivatives
3157 OneDimensionalLegendreShape<7> psi1(s[0]), psi2(s[1]);
3158 OneDimensionalLegendreDShape<7> dpsi1ds(s[0]), dpsi2ds(s[1]);
3159
3160 // Index for the shape functions
3161 unsigned index = 0;
3162 // Loop over shape functions in element
3163 for (unsigned i = 0; i < 7; i++)
3164 {
3165 for (unsigned j = 0; j < 7; j++)
3166 {
3167 // Assign the values
3168 dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3169 dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3170 psi[index] = psi2[i] * psi1[j];
3171 // Increment the index
3172 ++index;
3173 }
3174 }
3175 break;
3176 }
3177 default:
3178 std::ostringstream error_message;
3179 error_message << "\nERROR: Exceeded maximum polynomial order for";
3180 error_message << "\n polynomial order for shape functions.\n";
3181 throw OomphLibError(error_message.str(),
3182 OOMPH_CURRENT_FUNCTION,
3183 OOMPH_EXCEPTION_LOCATION);
3184 }
3185 }
3186
3187 //=======================================================================
3188 /// Second derivatives of shape functions for PRefineableQElement<DIM>
3189 /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
3190 //=======================================================================
3191 template<unsigned INITIAL_NNODE_1D>
3193 const Vector<double>& s, Shape& psi, DShape& dpsids, DShape& d2psids) const
3194 {
3195 std::ostringstream error_message;
3196 error_message
3197 << "\nd2shape_local currently not implemented for this element\n";
3198 throw OomphLibError(
3199 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3200 }
3201
3202 //=======================================================================
3203 /// Rebuild the element from nodes found in its sons
3204 /// Adjusts its p-order to be the maximum of its sons' p-orders
3205 //=======================================================================
3206 template<unsigned INITIAL_NNODE_1D>
3208 Mesh*& mesh_pt)
3209 {
3210 using namespace QuadTreeNames;
3211
3212 // Get p-orders of sons
3213 unsigned n_sons = this->tree_pt()->nsons();
3214 Vector<unsigned> son_p_order(n_sons);
3215 unsigned max_son_p_order = 0;
3216 for (unsigned ison = 0; ison < n_sons; ison++)
3217 {
3218 PRefineableElement* el_pt = dynamic_cast<PRefineableElement*>(
3219 this->tree_pt()->son_pt(ison)->object_pt());
3220 son_p_order[ison] = el_pt->p_order();
3221 if (son_p_order[ison] > max_son_p_order)
3222 max_son_p_order = son_p_order[ison];
3223 }
3224
3225 unsigned old_Nnode = this->nnode();
3226 unsigned old_P_order = this->p_order();
3227 // Set p-order of the element
3228 this->p_order() = max_son_p_order;
3229
3230 // Change integration scheme
3231 delete this->integral_pt();
3232 switch (this->p_order())
3233 {
3234 case 2:
3235 this->set_integration_scheme(new GaussLobattoLegendre<2, 2>);
3236 break;
3237 case 3:
3238 this->set_integration_scheme(new GaussLobattoLegendre<2, 3>);
3239 break;
3240 case 4:
3241 this->set_integration_scheme(new GaussLobattoLegendre<2, 4>);
3242 break;
3243 case 5:
3244 this->set_integration_scheme(new GaussLobattoLegendre<2, 5>);
3245 break;
3246 case 6:
3247 this->set_integration_scheme(new GaussLobattoLegendre<2, 6>);
3248 break;
3249 case 7:
3250 this->set_integration_scheme(new GaussLobattoLegendre<2, 7>);
3251 break;
3252 default:
3253 std::ostringstream error_message;
3254 error_message << "\nERROR: Exceeded maximum polynomial order for";
3255 error_message << "\n integration scheme.\n";
3256 throw OomphLibError(error_message.str(),
3257 OOMPH_CURRENT_FUNCTION,
3258 OOMPH_EXCEPTION_LOCATION);
3259 }
3260
3261 // Back up pointers to old nodes before they are lost
3262 Vector<Node*> old_node_pt(old_Nnode);
3263 for (unsigned n = 0; n < old_Nnode; n++)
3264 {
3265 old_node_pt[n] = this->node_pt(n);
3266 }
3267
3268 // Allocate new space for Nodes (at the element level)
3269 this->set_n_node(this->p_order() * this->p_order());
3270
3271 // Copy vertex nodes which were populated in the pre-build
3272 this->node_pt(0) = old_node_pt[0];
3273 this->node_pt(this->p_order() - 1) = old_node_pt[old_P_order - 1];
3274 this->node_pt(this->p_order() * (this->p_order() - 1)) =
3275 old_node_pt[(old_P_order) * (old_P_order - 1)];
3276 this->node_pt(this->p_order() * this->p_order() - 1) =
3277 old_node_pt[(old_P_order) * (old_P_order)-1];
3278
3279 // Copy midpoint nodes from sons if new p-order is odd
3280 if (this->p_order() % 2 == 1)
3281 {
3282 // Work out which is midpoint node
3283 unsigned midpoint = (this->p_order() - 1) / 2;
3284
3285 // Bottom edge
3286 this->node_pt(midpoint) = dynamic_cast<RefineableQElement<2>*>(
3287 quadtree_pt()->son_pt(SW)->object_pt())
3288 ->vertex_node_pt(1);
3289 // Left edge
3290 this->node_pt(midpoint * this->p_order()) =
3291 dynamic_cast<RefineableQElement<2>*>(
3292 quadtree_pt()->son_pt(SW)->object_pt())
3293 ->vertex_node_pt(2);
3294 // Top edge
3295 this->node_pt((this->p_order() - 1) * this->p_order() + midpoint) =
3296 dynamic_cast<RefineableQElement<2>*>(
3297 quadtree_pt()->son_pt(NE)->object_pt())
3298 ->vertex_node_pt(2);
3299 // Right edge
3300 this->node_pt((midpoint + 1) * this->p_order() - 1) =
3301 dynamic_cast<RefineableQElement<2>*>(
3302 quadtree_pt()->son_pt(NE)->object_pt())
3303 ->vertex_node_pt(1);
3304 }
3305
3306
3307 // The timestepper should be the same for all nodes and node 0 should
3308 // never be deleted.
3309 if (this->node_pt(0) == 0)
3310 {
3311 throw OomphLibError("The Corner node (0) does not exist",
3312 OOMPH_CURRENT_FUNCTION,
3313 OOMPH_EXCEPTION_LOCATION);
3314 }
3315
3316 TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
3317 unsigned ntstorage = time_stepper_pt->ntstorage();
3318
3319 unsigned jnod = 0;
3320 Vector<double> s_fraction(2), s(2);
3321 // Loop over the nodes in the element
3322 unsigned n_p = this->nnode_1d();
3323 for (unsigned i0 = 0; i0 < n_p; i0++)
3324 {
3325 // Get the fractional position of the node
3326 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3327 // Local coordinate
3328 s[0] = -1.0 + 2.0 * s_fraction[0];
3329
3330 for (unsigned i1 = 0; i1 < n_p; i1++)
3331 {
3332 // Get the fractional position of the node in the direction of s[1]
3333 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
3334 // Local coordinate in father element
3335 s[1] = -1.0 + 2.0 * s_fraction[1];
3336
3337 // Set the local node number
3338 jnod = i0 + n_p * i1;
3339
3340 // Initialise flag: So far, this node hasn't been built
3341 // or copied yet
3342 bool node_done = false;
3343
3344 // Get the pointer to the node in this element, returns NULL
3345 // if there is not node
3346 Node* created_node_pt = this->get_node_at_local_coordinate(s);
3347
3348 // Does this node already exist in this element?
3349 //----------------------------------------------
3350 if (created_node_pt != 0)
3351 {
3352 // Copy node across
3353 this->node_pt(jnod) = created_node_pt;
3354
3355 // Node has been created by copy
3356 node_done = true;
3357 }
3358 // Node does not exist in this element but might already
3359 //------------------------------------------------------
3360 // have been created by neighbouring elements
3361 //-------------------------------------------
3362 else
3363 {
3364 // Was the node created by one of its neighbours
3365 // Whether or not the node lies on an edge can be calculated
3366 // by from the fractional position
3367 bool is_periodic = false;
3368 created_node_pt = node_created_by_neighbour(s_fraction, is_periodic);
3369
3370 // If the node was so created, assign the pointers
3371 if (created_node_pt != 0)
3372 {
3373 // If the node is periodic
3374 if (is_periodic)
3375 {
3376 throw OomphLibError("Cannot handle periodic nodes yet",
3377 OOMPH_CURRENT_FUNCTION,
3378 OOMPH_EXCEPTION_LOCATION);
3379 }
3380 // Non-periodic case, just set the pointer
3381 else
3382 {
3383 this->node_pt(jnod) = created_node_pt;
3384 }
3385 // Node has been created
3386 node_done = true;
3387 }
3388 } // Node does not exist in this element
3389
3390 // Node has not been built anywhere ---> build it here
3391 if (!node_done)
3392 {
3393 // First, find the son element in which the node should live
3394
3395 // Find coordinates in the sons
3396 Vector<double> s_in_son(2);
3397 using namespace QuadTreeNames;
3398 int son = -10;
3399 // If negative on the west side
3400 if (s_fraction[0] < 0.5)
3401 {
3402 // On the south side
3403 if (s_fraction[1] < 0.5)
3404 {
3405 // It's the southwest son
3406 son = SW;
3407 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
3408 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
3409 }
3410 // On the north side
3411 else
3412 {
3413 // It's the northwest son
3414 son = NW;
3415 s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
3416 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
3417 }
3418 }
3419 else
3420 {
3421 // On the south side
3422 if (s_fraction[1] < 0.5)
3423 {
3424 // It's the southeast son
3425 son = SE;
3426 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
3427 s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
3428 }
3429 // On the north side
3430 else
3431 {
3432 // It's the northeast son
3433 son = NE;
3434 s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
3435 s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
3436 }
3437 }
3438
3439 // Get the pointer to the son element in which the new node
3440 // would live
3443 this->tree_pt()->son_pt(son)->object_pt());
3444
3445 // If we are rebuilding, then worry about the boundary conditions
3446 // Find the boundary of the node
3447 // Initially none
3448 int boundary = Tree::OMEGA;
3449 // If we are on the western boundary
3450 if (i0 == 0)
3451 {
3452 boundary = W;
3453 }
3454 // If we are on the eastern boundary
3455 else if (i0 == n_p - 1)
3456 {
3457 boundary = E;
3458 }
3459
3460 // If we are on the southern boundary
3461 if (i1 == 0)
3462 {
3463 // If we already have already set the boundary, we're on a corner
3464 switch (boundary)
3465 {
3466 case W:
3467 boundary = SW;
3468 break;
3469 case E:
3470 boundary = SE;
3471 break;
3472 // Boundary not set
3473 default:
3474 boundary = S;
3475 break;
3476 }
3477 }
3478 // If we are the northern bounadry
3479 else if (i1 == n_p - 1)
3480 {
3481 // If we already have a boundary
3482 switch (boundary)
3483 {
3484 case W:
3485 boundary = NW;
3486 break;
3487 case E:
3488 boundary = NE;
3489 break;
3490 default:
3491 boundary = N;
3492 break;
3493 }
3494 }
3495
3496 // set of boundaries that this edge in the son lives on
3497 std::set<unsigned> boundaries;
3498
3499 // Now get the boundary conditions from the son
3500 // The boundaries will be common to the son because there can be
3501 // no rotations here
3502 if (boundary != Tree::OMEGA)
3503 {
3504 son_el_pt->get_boundaries(boundary, boundaries);
3505 }
3506
3507 // If the node lives on a boundary:
3508 // Construct a boundary node,
3509 // Get boundary conditions and
3510 // update all lookup schemes
3511 if (boundaries.size() > 0)
3512 {
3513 // Construct the new node
3514 created_node_pt =
3515 this->construct_boundary_node(jnod, time_stepper_pt);
3516
3517 // Get the boundary conditions from the son
3518 Vector<int> bound_cons(this->ncont_interpolated_values());
3519 son_el_pt->get_bcs(boundary, bound_cons);
3520
3521 // Loop over the values and pin if necessary
3522 unsigned nval = created_node_pt->nvalue();
3523 for (unsigned k = 0; k < nval; k++)
3524 {
3525 if (bound_cons[k])
3526 {
3527 created_node_pt->pin(k);
3528 }
3529 }
3530
3531 // Solid node? If so, deal with the positional boundary
3532 // conditions:
3533 SolidNode* solid_node_pt =
3534 dynamic_cast<SolidNode*>(created_node_pt);
3535 if (solid_node_pt != 0)
3536 {
3537 // Get the positional boundary conditions from the father:
3538 unsigned n_dim = created_node_pt->ndim();
3539 Vector<int> solid_bound_cons(n_dim);
3540 RefineableSolidQElement<2>* son_solid_el_pt =
3541 dynamic_cast<RefineableSolidQElement<2>*>(son_el_pt);
3542#ifdef PARANOID
3543 if (son_solid_el_pt == 0)
3544 {
3545 std::string error_message =
3546 "We have a SolidNode outside a refineable SolidElement\n";
3547 error_message +=
3548 "during mesh refinement -- this doesn't make sense\n";
3549
3550 throw OomphLibError(error_message,
3551 OOMPH_CURRENT_FUNCTION,
3552 OOMPH_EXCEPTION_LOCATION);
3553 }
3554#endif
3555 son_solid_el_pt->get_solid_bcs(boundary, solid_bound_cons);
3556
3557 // Loop over the positions and pin, if necessary
3558 for (unsigned k = 0; k < n_dim; k++)
3559 {
3560 if (solid_bound_cons[k])
3561 {
3562 solid_node_pt->pin_position(k);
3563 }
3564 }
3565 } // End of if solid_node_pt
3566
3567
3568 // Next we update the boundary look-up schemes
3569 // Loop over the boundaries stored in the set
3570 for (std::set<unsigned>::iterator it = boundaries.begin();
3571 it != boundaries.end();
3572 ++it)
3573 {
3574 // Add the node to the boundary
3575 mesh_pt->add_boundary_node(*it, created_node_pt);
3576
3577 // If we have set an intrinsic coordinate on this
3578 // mesh boundary then it must also be interpolated on
3579 // the new node
3580 // Now interpolate the intrinsic boundary coordinate
3581 if (mesh_pt->boundary_coordinate_exists(*it) == true)
3582 {
3583 Vector<double> zeta(1);
3584 son_el_pt->interpolated_zeta_on_edge(
3585 *it, boundary, s_in_son, zeta);
3586
3587 created_node_pt->set_coordinates_on_boundary(*it, zeta);
3588 }
3589 }
3590 }
3591 // Otherwise the node is not on a Mesh boundary
3592 // and we create a normal "bulk" node
3593 else
3594 {
3595 // Construct the new node
3596 created_node_pt = this->construct_node(jnod, time_stepper_pt);
3597 }
3598
3599 // Now we set the position and values at the newly created node
3600
3601 // In the first instance use macro element or FE representation
3602 // to create past and present nodal positions.
3603 // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
3604 // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
3605 // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
3606 // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
3607 // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
3608 // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
3609 // NOT ASSIGN SENSIBLE INITIAL POSITONS!
3610
3611 // Loop over # of history values
3612 // Loop over # of history values
3613 for (unsigned t = 0; t < ntstorage; t++)
3614 {
3615 using namespace QuadTreeNames;
3616 // Get the position from the son
3617 Vector<double> x_prev(2);
3618
3619 // Now let's fill in the value
3620 son_el_pt->get_x(t, s_in_son, x_prev);
3621 for (unsigned i = 0; i < 2; i++)
3622 {
3623 created_node_pt->x(t, i) = x_prev[i];
3624 }
3625 }
3626
3627 // Now set up the values
3628 // Loop over all history values
3629 for (unsigned t = 0; t < ntstorage; t++)
3630 {
3631 // Get values from father element
3632 // Note: get_interpolated_values() sets Vector size itself.
3633 Vector<double> prev_values;
3634 son_el_pt->get_interpolated_values(t, s_in_son, prev_values);
3635
3636 // Initialise the values at the new node
3637 for (unsigned k = 0; k < created_node_pt->nvalue(); k++)
3638 {
3639 created_node_pt->set_value(t, k, prev_values[k]);
3640 }
3641 }
3642
3643 // Add the node to the mesh
3644 mesh_pt->add_node_pt(created_node_pt);
3645
3646 // Check if the element is an algebraic element
3647 AlgebraicElementBase* alg_el_pt =
3648 dynamic_cast<AlgebraicElementBase*>(this);
3649
3650 // If we do have an algebraic element
3651 if (alg_el_pt != 0)
3652 {
3653 std::string error_message =
3654 "Have not implemented rebuilding from sons for";
3655 error_message += "Algebraic p-refineable elements yet\n";
3656
3657 throw OomphLibError(
3658 error_message,
3659 "PRefineableQElement<2,INITIAL_NNODE_1D>::rebuild_from_sons()",
3660 OOMPH_EXCEPTION_LOCATION);
3661 }
3662
3663 } // End of the case when we build the node ourselves
3664 }
3665 } // End of loop over all nodes in element
3666
3667
3668 // If the element is a MacroElementNodeUpdateElement, set the update
3669 // parameters for the current element's nodes. These need to be reset
3670 // (as in MacroElementNodeUpdateElement<ELEMENT>::rebuild_from_sons())
3671 // because the nodes in this element have changed
3673 dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
3674 if (m_el_pt != 0)
3675 {
3676 // Loop over the nodes
3677 for (unsigned j = 0; j < this->nnode(); j++)
3678 {
3679 // Get local coordinate in element (Vector sets its own size)
3680 Vector<double> s_in_node_update_element;
3681 this->local_coordinate_of_node(j, s_in_node_update_element);
3682
3683 // Get vector of geometric objects
3684 Vector<GeomObject*> geom_object_pt(m_el_pt->geom_object_pt());
3685
3686 // Pass the lot to the node
3687 static_cast<MacroElementNodeUpdateNode*>(this->node_pt(j))
3688 ->set_node_update_info(
3689 this, s_in_node_update_element, geom_object_pt);
3690 }
3691 }
3692
3693 // MacroElementNodeUpdateElementBase* m_el_pt=dynamic_cast<
3694 // MacroElementNodeUpdateElementBase*>(this);
3695 // if(m_el_pt!=0)
3696 // {
3697 // Loop over all nodes in element again, to re-set the positions
3698 // This must be done using the new element's macro-element
3699 // representation, rather than the old version which may be
3700 // of a different p-order!
3701 for (unsigned i0 = 0; i0 < n_p; i0++)
3702 {
3703 // Get the fractional position of the node
3704 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3705 // Local coordinate
3706 s[0] = -1.0 + 2.0 * s_fraction[0];
3707
3708 for (unsigned i1 = 0; i1 < n_p; i1++)
3709 {
3710 // Get the fractional position of the node in the direction of s[1]
3711 s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
3712 // Local coordinate in father element
3713 s[1] = -1.0 + 2.0 * s_fraction[1];
3714
3715 // Set the local node number
3716 jnod = i0 + n_p * i1;
3717
3718 // Loop over # of history values
3719 for (unsigned t = 0; t < ntstorage; t++)
3720 {
3721 // Get position from father element -- this uses the macro
3722 // element representation if appropriate. If the node
3723 // turns out to be a hanging node later on, then
3724 // its position gets adjusted in line with its
3725 // hanging node interpolation.
3726 Vector<double> x_prev(2);
3727 this->get_x(t, s, x_prev);
3728
3729 // Set previous positions of the new node
3730 for (unsigned i = 0; i < 2; i++)
3731 {
3732 this->node_pt(jnod)->x(t, i) = x_prev[i];
3733 }
3734 }
3735 }
3736 }
3737 // }
3738 }
3739
3740 //=================================================================
3741 /// Check inter-element continuity of
3742 /// - nodal positions
3743 /// - (nodally) interpolated function values
3744 /// Overloaded to not check differences in the value. Mortaring
3745 /// doesn't enforce strong continuity between elements.
3746 //====================================================================
3747 template<unsigned INITIAL_NNODE_1D>
3749 double& max_error)
3750 {
3751 // Overloaded to *not* check for continuity in value of interpolated
3752 // variables. This is necessary because mortaring does not ensure continuity
3753 // across element boundaries. It therefore makes no sense to test for this.
3754
3755 // Dummy set max_error to 0
3756 max_error = 0.0;
3757
3758 // With macro-elements, (strong) continuity in position is nolonger
3759 // guaranteed either, so we don't check for this either. In fact, we do
3760 // nothing at all.
3761 if (this->macro_elem_pt() != 0)
3762 {
3763 // We have a macro element, so do nothing!
3764 return;
3765 }
3766
3767 using namespace QuadTreeNames;
3768
3769 // Number of nodes along edge
3770 unsigned n_p = nnode_1d();
3771
3772 // Number of timesteps (incl. present) for which continuity is
3773 // to be checked.
3774 unsigned n_time = 1;
3775
3776 // Initialise errors
3777 max_error = 0.0;
3778 Vector<double> max_error_x(2, 0.0);
3779 double max_error_val = 0.0;
3780
3781 Vector<int> edges(4);
3782 edges[0] = S;
3783 edges[1] = N;
3784 edges[2] = W;
3785 edges[3] = E;
3786
3787 // Loop over the edges
3788 for (unsigned edge_counter = 0; edge_counter < 4; edge_counter++)
3789 {
3790 Vector<unsigned> translate_s(2);
3791 Vector<double> s(2), s_lo_neigh(2), s_hi_neigh(2), s_fraction(2);
3792 int neigh_edge, diff_level;
3793 bool in_neighbouring_tree;
3794
3795 // Find pointer to neighbour in this direction
3796 QuadTree* neigh_pt;
3797 neigh_pt = quadtree_pt()->gteq_edge_neighbour(edges[edge_counter],
3798 translate_s,
3799 s_lo_neigh,
3800 s_hi_neigh,
3801 neigh_edge,
3802 diff_level,
3803 in_neighbouring_tree);
3804
3805 // Neighbour exists and has existing nodes
3806 if ((neigh_pt != 0) && (neigh_pt->object_pt()->nodes_built()))
3807 {
3808 // Need to exclude periodic nodes from this check
3809 // There are only periodic nodes if we are in a neighbouring tree
3810 bool is_periodic = false;
3811 if (in_neighbouring_tree)
3812 {
3813 // Is it periodic
3814 is_periodic = this->tree_pt()->root_pt()->is_neighbour_periodic(
3815 edges[edge_counter]);
3816 }
3817
3818 // We also need to exclude edges which may have hanging nodes
3819 // because mortaring does not guarantee (strong) continuity
3820 // in position or in value at nonconforming element boundaries
3821 bool exclude_this_edge = false;
3822 if (diff_level != 0)
3823 {
3824 // h-type nonconformity (dependent)
3825 exclude_this_edge = true;
3826 }
3827 else if (neigh_pt->nsons() != 0)
3828 {
3829 // h-type nonconformity (master)
3830 exclude_this_edge = true;
3831 }
3832 else
3833 {
3834 unsigned my_p_order = this->p_order();
3835 unsigned neigh_p_order =
3836 dynamic_cast<PRefineableQElement*>(neigh_pt->object_pt())
3837 ->p_order();
3838 if (my_p_order != neigh_p_order)
3839 {
3840 // p-type nonconformity
3841 exclude_this_edge = true;
3842 }
3843 }
3844
3845 // With macro-elements, (strong) continuity in position is nolonger
3846 // guaranteed either, so we don't check for this either. In fact, we do
3847 // nothing at all.
3848 if (dynamic_cast<FiniteElement*>(neigh_pt->object_pt())
3849 ->macro_elem_pt() != 0)
3850 {
3851 // We have a macro element, so do nothing!
3852 break;
3853 }
3854
3855 // Check conforming edges
3856 if (!exclude_this_edge)
3857 {
3858 // Loop over nodes along the edge
3859 for (unsigned i0 = 0; i0 < n_p; i0++)
3860 {
3861 // Storage for pointer to the local node
3862 Node* local_node_pt = 0;
3863
3864 switch (edge_counter)
3865 {
3866 case 0:
3867 // Local fraction of node
3868 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3869 s_fraction[1] = 0.0;
3870 // Get pointer to local node
3871 local_node_pt = this->node_pt(i0);
3872 break;
3873
3874 case 1:
3875 // Local fraction of node
3876 s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3877 s_fraction[1] = 1.0;
3878 // Get pointer to local node
3879 local_node_pt = this->node_pt(i0 + n_p * (n_p - 1));
3880 break;
3881
3882 case 2:
3883 // Local fraction of node
3884 s_fraction[0] = 0.0;
3885 s_fraction[1] = this->local_one_d_fraction_of_node(i0, 1);
3886 // Get pointer to local node
3887 local_node_pt = this->node_pt(n_p * i0);
3888 break;
3889
3890 case 3:
3891 // Local fraction of node
3892 s_fraction[0] = 1.0;
3893 s_fraction[1] = this->local_one_d_fraction_of_node(i0, 1);
3894 // Get pointer to local node
3895 local_node_pt = this->node_pt(n_p - 1 + n_p * i0);
3896 break;
3897 }
3898
3899 // Calculate the local coordinate and the local coordinate as viewed
3900 // from the neighbour
3901 Vector<double> s_in_neighb(2);
3902 for (unsigned i = 0; i < 2; i++)
3903 {
3904 // Local coordinate in this element
3905 s[i] = -1.0 + 2.0 * s_fraction[i];
3906 // Local coordinate in the neighbour
3907 s_in_neighb[i] =
3908 s_lo_neigh[i] +
3909 s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
3910 }
3911
3912 // Loop over timesteps
3913 for (unsigned t = 0; t < n_time; t++)
3914 {
3915 // Get the nodal position from neighbour element
3916 Vector<double> x_in_neighb(2);
3917 neigh_pt->object_pt()->interpolated_x(
3918 t, s_in_neighb, x_in_neighb);
3919
3920 // Check error only if the node is NOT periodic
3921 if (is_periodic == false)
3922 {
3923 for (int i = 0; i < 2; i++)
3924 {
3925 // Find the spatial error
3926 double err =
3927 std::fabs(local_node_pt->x(t, i) - x_in_neighb[i]);
3928
3929 // If it's bigger than our tolerance, say so
3930 if (err > 1e-9)
3931 {
3932 oomph_info << "errx " << err << " " << t << " "
3933 << local_node_pt->x(t, i) << " "
3934 << x_in_neighb[i] << std::endl;
3935
3936 oomph_info << "at " << local_node_pt->x(0) << " "
3937 << local_node_pt->x(1) << std::endl;
3938 }
3939
3940 // If it's bigger than the previous max error, it is the
3941 // new max error!
3942 if (err > max_error_x[i])
3943 {
3944 max_error_x[i] = err;
3945 }
3946 }
3947 }
3948
3949 // Get the values from neighbour element. Note: # of values
3950 // gets set by routine (because in general we don't know
3951 // how many interpolated values a certain element has
3952 Vector<double> values_in_neighb;
3954 t, s_in_neighb, values_in_neighb);
3955
3956 // Get the values in current element.
3957 Vector<double> values;
3958 this->get_interpolated_values(t, s, values);
3959
3960 // Now figure out how many continuously interpolated values there
3961 // are
3962 unsigned num_val =
3963 neigh_pt->object_pt()->ncont_interpolated_values();
3964
3965 // Check error
3966 for (unsigned ival = 0; ival < num_val; ival++)
3967 {
3968 double err = std::fabs(values[ival] - values_in_neighb[ival]);
3969
3970 if (err > 1.0e-10)
3971 {
3972 oomph_info << local_node_pt->x(0) << " "
3973 << local_node_pt->x(1) << " \n# "
3974 << "erru (S)" << err << " " << ival << " "
3975 << this->get_node_number(local_node_pt) << " "
3976 << values[ival] << " " << values_in_neighb[ival]
3977 << std::endl;
3978 }
3979
3980 if (err > max_error_val)
3981 {
3982 max_error_val = err;
3983 }
3984 }
3985 }
3986 }
3987 }
3988 }
3989 }
3990
3991 max_error = max_error_x[0];
3992 if (max_error_x[1] > max_error) max_error = max_error_x[1];
3993 if (max_error_val > max_error) max_error = max_error_val;
3994
3995 if (max_error > 1e-9)
3996 {
3997 oomph_info << "\n#------------------------------------ \n#Max error ";
3998 oomph_info << max_error_x[0] << " " << max_error_x[1] << " "
3999 << max_error_val << std::endl;
4000 oomph_info << "#------------------------------------ \n " << std::endl;
4001 }
4002 }
4003
4004 //=================================================================
4005 /// Internal function to set up the hanging nodes on a particular
4006 /// edge of the element.
4007 /// Implements the mortarting method to enforce continuity weakly
4008 /// across non-conforming element boundaries \f$\Gamma\f$ using an
4009 /// integral matching condition
4010 /// \f[ \int_\Gamma (u_{\mbox{S}} - u_{\mbox{M}}) \psi \mbox{d} s = 0 \f]
4011 /// for all polynomials \f$\psi\f$ on \f$\Gamma\f$ of degree at most
4012 /// p-2 (where p is the spectral-order of the dependent element) and a
4013 /// vertex matching condition
4014 /// \f[ (u_{\mbox{S}} - u_{\mbox{M}})\big\vert_{\partial\Gamma} = 0.\f]
4015 ///
4016 /// The algorithm works as follows:
4017 /// - First the element determines if its edge my_edge is on the
4018 /// master or dependent side of the non-conformity.
4019 /// At h-type non-conformities
4020 /// we choose long edges to be masters, and at p-type nonconformities the
4021 /// edge with lower p-order is the master.
4022 /// - Mortaring is performed by the dependent side.
4023 /// - If a vertex node of the mortar is shared between master and dependent
4024 /// element then the vertex matching condition is enforced automatically.
4025 /// Otherwise it must be imposed by constraining its value to that at on
4026 /// the master side.
4027 /// - The integral matching condition is discretised and the mortar test
4028 /// functions \f$ \psi \f$ are chosen to be derivatives of Legendre
4029 /// polynomials of degree p-1.
4030 /// - The mortar mass matrix M is constructed. Its entries are the
4031 /// mortar test functions evaluated at the dependent nodal positions,
4032 /// so it is diagonal.
4033 /// - The local projection matrix is constructed for the master element by
4034 /// applying the discretised integral matching condition along the mortar
4035 /// using the appropriate quadrature order.
4036 /// - The global projection matrix is then assembled by subtracting
4037 /// contributions from the mortar vertex nodes.
4038 /// - The mortar system \f$ M\xi^s = P\hat{\xi^m} \f$ is constructed,
4039 /// where \f$ \xi^m \f$ and \f$ \xi^s \f$ are the nodal values at
4040 /// the master and dependent nodes respectively.
4041 /// - The conformity matrix \f$ C = M^{-1}P \f$ is computed. This is
4042 /// straightforward since the mass matrix is diagonal.
4043 /// - Finally, the master nodes and weights for each dependent node
4044 /// are read from
4045 /// the conformity matrix and stored in the dependents' hanging schemes.
4046 ///
4047 /// The positions of the dependent nodes are set to be consistent with their
4048 /// hanging schemes.
4049 //=================================================================
4050 template<unsigned INITIAL_NNODE_1D>
4052 const int& value_id, const int& my_edge, std::ofstream& output_hangfile)
4053 {
4054 using namespace QuadTreeNames;
4055
4056 Vector<unsigned> translate_s(2);
4057 Vector<double> s_lo_neigh(2);
4058 Vector<double> s_hi_neigh(2);
4059 int neigh_edge, diff_level;
4060 bool in_neighbouring_tree;
4061
4062 // Find pointer to neighbour in this direction
4063 QuadTree* neigh_pt;
4064 neigh_pt = this->quadtree_pt()->gteq_edge_neighbour(my_edge,
4065 translate_s,
4066 s_lo_neigh,
4067 s_hi_neigh,
4068 neigh_edge,
4069 diff_level,
4070 in_neighbouring_tree);
4071
4072 // Work out master/dependent edges
4073 //----------------------------
4074
4075 // Set up booleans
4076 // bool h_type_master = false;
4077 bool h_type_dependent = false;
4078 // bool p_type_master = false;
4079 bool p_type_dependent = false;
4080
4081 // Neighbour exists and all nodes have been created
4082 if (neigh_pt != 0)
4083 {
4084 // Check if neighbour is bigger than me
4085 if (diff_level != 0)
4086 {
4087 // Dependent at h-type non-conformity
4088 h_type_dependent = true;
4089 }
4090 // Check if neighbour is the same size as me
4091 else if (neigh_pt->nsons() == 0)
4092 {
4093 // Neighbour is same size as me
4094 // Find p-orders of each element
4095 unsigned my_p_order =
4096 dynamic_cast<PRefineableQElement<2, INITIAL_NNODE_1D>*>(this)
4097 ->p_order();
4098 unsigned neigh_p_order =
4100 neigh_pt->object_pt())
4101 ->p_order();
4102
4103 // Check for p-type non-conformity
4104 if (neigh_p_order == my_p_order)
4105 {
4106 // At a conforming interface
4107 }
4108 else if (neigh_p_order < my_p_order)
4109 {
4110 // Dependent at p-type non-conformity
4111 p_type_dependent = true;
4112 }
4113 else
4114 {
4115 // Master at p-type non-conformity
4116 // p_type_master = true;
4117 }
4118 }
4119 // Neighbour must be smaller than me
4120 else
4121 {
4122 // Master at h-type non-conformity
4123 // h_type_master = true;
4124 }
4125 }
4126 else
4127 {
4128 // Edge is on a boundary
4129 }
4130
4131 // Now do the mortaring
4132 //---------------------
4133 if (h_type_dependent || p_type_dependent)
4134 {
4135 // Compute the active coordinate index along the this side of mortar
4136 unsigned active_coord_index;
4137 if (my_edge == N || my_edge == S) active_coord_index = 0;
4138 else if (my_edge == E || my_edge == W)
4139 active_coord_index = 1;
4140 else
4141 {
4142 throw OomphLibError("Cannot transform coordinates",
4143 OOMPH_CURRENT_FUNCTION,
4144 OOMPH_EXCEPTION_LOCATION);
4145 }
4146
4147 // Get pointer to neighbouring master element (in p-refineable form)
4149 neigh_obj_pt = dynamic_cast<PRefineableQElement<2, INITIAL_NNODE_1D>*>(
4150 neigh_pt->object_pt());
4151
4152 // Create vector of master and dependent nodes
4153 //----------------------------------------
4154 Vector<Node*> master_node_pt, dependent_node_pt;
4155 Vector<unsigned> master_node_number, dependent_node_number;
4156 Vector<Vector<double>> dependent_node_s_fraction;
4157
4158 // Number of nodes in one dimension
4159 const unsigned my_n_p = this->ninterpolating_node_1d(value_id);
4160 const unsigned neigh_n_p = neigh_obj_pt->ninterpolating_node_1d(value_id);
4161
4162 // Test for the periodic node case
4163 // Are we crossing a periodic boundary
4164 bool is_periodic = false;
4165 if (in_neighbouring_tree)
4166 {
4167 is_periodic =
4168 this->tree_pt()->root_pt()->is_neighbour_periodic(my_edge);
4169 }
4170
4171 // If it is periodic we actually need to get the node in
4172 // the neighbour of the neighbour (which will be a parent of
4173 // the present element) so that the "fixed" coordinate is
4174 // correctly calculated.
4175 // The idea is to replace the neigh_pt and associated data
4176 // with those of the neighbour of the neighbour
4177 if (is_periodic)
4178 {
4179 throw OomphLibError(
4180 "Cannot do mortaring with periodic hanging nodes yet! (Fix me!)",
4181 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4182 OOMPH_EXCEPTION_LOCATION);
4183 } // End of special treatment for periodic hanging nodes
4184
4185 // Storage for pointers to the nodes and their numbers along the master
4186 // edge
4187 unsigned neighbour_node_number = 0;
4188 Node* neighbour_node_pt = 0;
4189
4190 // Loop over nodes along the edge
4191 for (unsigned i0 = 0; i0 < neigh_n_p; i0++)
4192 {
4193 // Find the neighbour's node
4194 switch (neigh_edge)
4195 {
4196 case N:
4197 neighbour_node_number = i0 + neigh_n_p * (neigh_n_p - 1);
4198 neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
4199 neighbour_node_number, value_id);
4200 break;
4201
4202 case S:
4203 neighbour_node_number = i0;
4204 neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
4205 neighbour_node_number, value_id);
4206 break;
4207
4208 case E:
4209 neighbour_node_number = neigh_n_p - 1 + neigh_n_p * i0;
4210 neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
4211 neighbour_node_number, value_id);
4212 break;
4213
4214 case W:
4215 neighbour_node_number = neigh_n_p * i0;
4216 neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
4217 neighbour_node_number, value_id);
4218 break;
4219
4220 default:
4221 throw OomphLibError("my_edge not N, S, W, E\n",
4222 OOMPH_CURRENT_FUNCTION,
4223 OOMPH_EXCEPTION_LOCATION);
4224 }
4225
4226 // Set node as master node
4227 master_node_number.push_back(neighbour_node_number);
4228 master_node_pt.push_back(neighbour_node_pt);
4229 }
4230
4231 // Storage for pointers to the local nodes and their numbers along my edge
4232 unsigned local_node_number = 0;
4233 Node* local_node_pt = 0;
4234
4235 // Loop over the nodes along my edge
4236 for (unsigned i0 = 0; i0 < my_n_p; i0++)
4237 {
4238 // Storage for the fractional position of the node in the element
4239 Vector<double> s_fraction(2);
4240
4241 // Find the local node and the fractional position of the node
4242 // which depends on the edge, of course
4243 switch (my_edge)
4244 {
4245 case N:
4246 s_fraction[0] =
4247 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
4248 s_fraction[1] = 1.0;
4249 local_node_number = i0 + my_n_p * (my_n_p - 1);
4250 local_node_pt =
4251 this->interpolating_node_pt(local_node_number, value_id);
4252 break;
4253
4254 case S:
4255 s_fraction[0] =
4256 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
4257 s_fraction[1] = 0.0;
4258 local_node_number = i0;
4259 local_node_pt =
4260 this->interpolating_node_pt(local_node_number, value_id);
4261 break;
4262
4263 case E:
4264 s_fraction[0] = 1.0;
4265 s_fraction[1] =
4266 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
4267 local_node_number = my_n_p - 1 + my_n_p * i0;
4268 local_node_pt =
4269 this->interpolating_node_pt(local_node_number, value_id);
4270 break;
4271
4272 case W:
4273 s_fraction[0] = 0.0;
4274 s_fraction[1] =
4275 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
4276 local_node_number = my_n_p * i0;
4277 local_node_pt =
4278 this->interpolating_node_pt(local_node_number, value_id);
4279 break;
4280
4281 default:
4282 throw OomphLibError("my_edge not N, S, W, E\n",
4283 OOMPH_CURRENT_FUNCTION,
4284 OOMPH_EXCEPTION_LOCATION);
4285 }
4286
4287 // Add node to vector of dependent element nodes
4288 dependent_node_number.push_back(local_node_number);
4289 dependent_node_pt.push_back(local_node_pt);
4290
4291 // Store node's local fraction
4292 dependent_node_s_fraction.push_back(s_fraction);
4293 }
4294
4295 // Store the number of dependent and master nodes for use later
4296 const unsigned n_dependent_nodes = dependent_node_pt.size();
4297 const unsigned n_master_nodes = master_node_pt.size();
4298 const unsigned dependent_element_nnode_1d = my_n_p;
4299 const unsigned master_element_nnode_1d = neigh_n_p;
4300
4301 // Storage for master shapes
4302 Shape master_shapes(neigh_obj_pt->ninterpolating_node(value_id));
4303
4304 // Get master and dependent nodal positions
4305 //-------------------------------------
4306 Vector<double> dependent_nodal_position;
4307 Vector<double> dependent_weight(dependent_element_nnode_1d);
4309 dependent_element_nnode_1d, dependent_nodal_position, dependent_weight);
4310 Vector<double> master_nodal_position;
4311 Vector<double> master_weight(master_element_nnode_1d);
4313 master_element_nnode_1d, master_nodal_position, master_weight);
4314
4315 // Apply the vertex matching condition
4316 //------------------------------------
4317 // Vertiex matching is ensured automatically in cases where there is a
4318 // node at each end of the mortar that is shared between the master and
4319 // dependent elements. Where this is not the case, the vertex matching
4320 // condition must be enforced by constraining the value of the unknown at
4321 // the node on the dependent side to be the same as the value at that
4322 // location in the master.
4323
4324 // Store positions of the mortar vertex/non-vertex nodes in the dependent
4325 // element
4326 const unsigned n_mortar_vertices = 2;
4327 Vector<unsigned> vertex_pos(n_mortar_vertices);
4328 vertex_pos[0] = 0;
4329 vertex_pos[1] = this->ninterpolating_node_1d(value_id) - 1;
4330 Vector<unsigned> non_vertex_pos(my_n_p - n_mortar_vertices);
4331 for (unsigned i = 0; i < my_n_p - n_mortar_vertices; i++)
4332 {
4333 non_vertex_pos[i] = i + 1;
4334 }
4335
4336 // Check if the mortar vertices are shared
4337 for (unsigned v = 0; v < n_mortar_vertices; v++)
4338 {
4339 // Search master node storage for the node
4340 bool node_is_shared = false;
4341 for (unsigned i = 0; i < master_node_pt.size(); i++)
4342 {
4343 if (dependent_node_pt[vertex_pos[v]] == master_node_pt[i])
4344 {
4345 node_is_shared = true;
4346 break;
4347 }
4348 }
4349
4350 // If the node is not shared then we must constrain its value by setting
4351 // up a hanging scheme
4352 if (!node_is_shared)
4353 {
4354 // Calculate weights. These are just the master shapes evaluated at
4355 // this dependent node's position
4356
4357 // Work out this node's location in the master
4358 Vector<double> s_in_neigh(2);
4359 for (unsigned i = 0; i < 2; i++)
4360 {
4361 s_in_neigh[i] =
4362 s_lo_neigh[i] +
4363 dependent_node_s_fraction[vertex_pos[v]][translate_s[i]] *
4364 (s_hi_neigh[i] - s_lo_neigh[i]);
4365 }
4366
4367 // Get master shapes at dependent nodal positions
4368 neigh_obj_pt->interpolating_basis(
4369 s_in_neigh, master_shapes, value_id);
4370
4371 // Remove any existing hanging node info
4372 // (This may be a bit wasteful, but guarantees correctness)
4373 dependent_node_pt[vertex_pos[v]]->set_nonhanging();
4374
4375 // Set up hanging scheme for this node
4376 HangInfo* explicit_hang_pt = new HangInfo(n_master_nodes);
4377 for (unsigned m = 0; m < n_master_nodes; m++)
4378 {
4379 explicit_hang_pt->set_master_node_pt(
4380 m, master_node_pt[m], master_shapes[master_node_number[m]]);
4381 }
4382
4383 // Make node hang
4384 dependent_node_pt[vertex_pos[v]]->set_hanging_pt(explicit_hang_pt,
4385 -1);
4386 }
4387 }
4388
4389 // Check that there are actually dependent nodes for which we still need
4390 // to construct a hanging scheme. If not then there is nothing more to do.
4391 if (n_dependent_nodes - n_mortar_vertices > 0)
4392 {
4393 // Assemble mass matrix for mortar
4394 //--------------------------------
4395 Vector<double> psi(n_dependent_nodes - n_mortar_vertices);
4396 Vector<double> diag_M(n_dependent_nodes - n_mortar_vertices);
4397 Vector<Vector<double>> shared_node_M(n_mortar_vertices);
4398 for (unsigned i = 0; i < shared_node_M.size(); i++)
4399 {
4400 shared_node_M[i].resize(n_dependent_nodes - n_mortar_vertices);
4401 }
4402
4403 // Fill in part corresponding to dependent nodal positions (unknown)
4404 for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4405 {
4406 // Use L'Hosptal's rule:
4407 psi[i] =
4408 pow(-1.0, int((dependent_element_nnode_1d - 1) - i - 1)) *
4409 -Orthpoly::ddlegendre(dependent_element_nnode_1d - 1,
4410 dependent_nodal_position[non_vertex_pos[i]]);
4411 // Put in contribution
4412 diag_M[i] = psi[i] * dependent_weight[non_vertex_pos[i]];
4413 }
4414
4415 // Fill in part corresponding to dependent element vertices (known)
4416 for (unsigned v = 0; v < shared_node_M.size(); v++)
4417 {
4418 for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4419 {
4420 // Check if denominator is zero
4421 if (std::fabs(dependent_nodal_position[non_vertex_pos[i]] -
4422 dependent_nodal_position[vertex_pos[v]]) >= 1.0e-8)
4423 {
4424 // We're ok
4425 psi[i] =
4426 pow(-1.0, int((dependent_element_nnode_1d - 1) - i - 1)) *
4427 Orthpoly::dlegendre(dependent_element_nnode_1d - 1,
4428 dependent_nodal_position[vertex_pos[v]]) /
4429 (dependent_nodal_position[non_vertex_pos[i]] -
4430 dependent_nodal_position[vertex_pos[v]]);
4431 }
4432 // Check if numerator is zero
4433 else if (std::fabs(Orthpoly::dlegendre(
4434 dependent_element_nnode_1d - 1,
4435 dependent_nodal_position[vertex_pos[v]])) < 1.0e-8)
4436 {
4437 // We can use l'hopital's rule
4438 psi[i] =
4439 pow(-1.0, int((dependent_element_nnode_1d - 1) - i - 1)) *
4441 dependent_element_nnode_1d - 1,
4442 dependent_nodal_position[non_vertex_pos[i]]);
4443 }
4444 else
4445 {
4446 // We can't use l'hopital's rule
4447 throw OomphLibError(
4448 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
4449 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4450 OOMPH_EXCEPTION_LOCATION);
4451 }
4452 // Put in contribution
4453 shared_node_M[v][i] = psi[i] * dependent_weight[vertex_pos[v]];
4454 }
4455 }
4456
4457 // Assemble local projection matrix for mortar
4458 //--------------------------------------------
4459
4460 // Have only one local projection matrix because there is just one
4461 // master
4462 Vector<Vector<double>> P(n_dependent_nodes - n_mortar_vertices);
4463 for (unsigned i = 0; i < P.size(); i++)
4464 {
4465 P[i].resize(n_master_nodes, 0.0);
4466 }
4467
4468 // Storage for local coordinate
4469 Vector<double> s(2);
4470
4471 // Sum contributions from master element shapes (quadrature).
4472 // The order of quadrature must be high enough to evaluate a polynomial
4473 // of order N_s + N_m - 3 exactly, where N_s = n_dependent_nodes, N_m =
4474 // n_master_nodes.
4475 // (Use pointers for the quadrature knots and weights so that
4476 // data is not unnecessarily copied)
4477 // unsigned quadrature_order =
4478 // std::max(dependent_element_nnode_1d,master_element_nnode_1d);
4479 Vector<double>*quadrature_knot, *quadrature_weight;
4480 if (dependent_element_nnode_1d >= master_element_nnode_1d)
4481 {
4482 // Use the same quadrature order as the dependent element (me)
4483 quadrature_knot = &dependent_nodal_position;
4484 quadrature_weight = &dependent_weight;
4485 }
4486 else
4487 {
4488 // Use the same quadrature order as the master element (neighbour)
4489 quadrature_knot = &master_nodal_position;
4490 quadrature_weight = &master_weight;
4491 }
4492
4493 // Quadrature loop
4494 for (unsigned q = 0; q < (*quadrature_weight).size(); q++)
4495 {
4496 // Evaluate mortar test functions at the quadrature knots in the
4497 // dependent
4498 // s[active_coord_index] = (*quadrature_knot)[q];
4499 Vector<double> s_on_mortar(1);
4500 s_on_mortar[0] = (*quadrature_knot)[q];
4501
4502 // Get psi
4503 for (unsigned k = 0; k < n_dependent_nodes - n_mortar_vertices; k++)
4504 {
4505 // Check if denominator is zero
4506 if (std::fabs(dependent_nodal_position[non_vertex_pos[k]] -
4507 s_on_mortar[0]) >= 1.0e-08)
4508 {
4509 // We're ok
4510 psi[k] =
4511 pow(-1.0, int((dependent_element_nnode_1d - 1) - k - 1)) *
4512 Orthpoly::dlegendre(dependent_element_nnode_1d - 1,
4513 s_on_mortar[0]) /
4514 (dependent_nodal_position[non_vertex_pos[k]] - s_on_mortar[0]);
4515 }
4516 // Check if numerator is zero
4517 else if (std::fabs(Orthpoly::dlegendre(
4518 dependent_element_nnode_1d - 1, s_on_mortar[0])) <
4519 1.0e-8)
4520 {
4521 // We can use l'Hopital's rule
4522 psi[k] =
4523 pow(-1.0, int((dependent_element_nnode_1d - 1) - k - 1)) *
4524 -Orthpoly::ddlegendre(dependent_element_nnode_1d - 1,
4525 s_on_mortar[0]);
4526 }
4527 else
4528 {
4529 // We can't use l'hopital's rule
4530 throw OomphLibError(
4531 "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
4532 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4533 OOMPH_EXCEPTION_LOCATION);
4534 }
4535 }
4536
4537 // Convert coordinate on mortar to local fraction in dependent element
4538 Vector<double> s_fraction(2);
4539 for (unsigned i = 0; i < 2; i++)
4540 {
4541 s_fraction[i] = (i == active_coord_index) ?
4542 0.5 * (s_on_mortar[0] + 1.0) :
4543 dependent_node_s_fraction[vertex_pos[0]][i];
4544 }
4545
4546 // Project active coordinate into master element
4547 Vector<double> s_in_neigh(2);
4548 for (unsigned i = 0; i < 2; i++)
4549 {
4550 s_in_neigh[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
4551 (s_hi_neigh[i] - s_lo_neigh[i]);
4552 }
4553
4554 // Evaluate master shapes at projections of local quadrature knots
4555 neigh_obj_pt->interpolating_basis(
4556 s_in_neigh, master_shapes, value_id);
4557
4558 // Populate local projection matrix
4559 for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4560 {
4561 for (unsigned j = 0; j < n_master_nodes; j++)
4562 {
4563 P[i][j] += master_shapes[master_node_number[j]] * psi[i] *
4564 (*quadrature_weight)[q];
4565 }
4566 }
4567 }
4568
4569 // Assemble global projection matrices for mortar
4570 //-----------------------------------------------
4571 // Need to subtract contributions from the "known unknowns"
4572 // corresponding to the nodes at the vertices of the mortar
4573
4574 // Assemble contributions from mortar vertex nodes
4575 for (unsigned v = 0; v < n_mortar_vertices; v++)
4576 {
4577 // Convert coordinate on mortar to local fraction in dependent element
4578 Vector<double> s_fraction(2);
4579 for (unsigned i = 0; i < 2; i++)
4580 {
4581 s_fraction[i] =
4582 (i == active_coord_index) ?
4583 0.5 * (dependent_nodal_position[vertex_pos[v]] + 1.0) :
4584 dependent_node_s_fraction[vertex_pos[0]][i];
4585 }
4586
4587 // Project active coordinate into master element
4588 Vector<double> s_in_neigh(2);
4589 for (unsigned i = 0; i < 2; i++)
4590 {
4591 s_in_neigh[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
4592 (s_hi_neigh[i] - s_lo_neigh[i]);
4593 }
4594
4595 // Get master shapes at dependent nodal positions
4596 neigh_obj_pt->interpolating_basis(
4597 s_in_neigh, master_shapes, value_id);
4598
4599 for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4600 {
4601 for (unsigned k = 0; k < n_master_nodes; k++)
4602 {
4603 P[i][k] -=
4604 master_shapes[master_node_number[k]] * shared_node_M[v][i];
4605 }
4606 }
4607 }
4608
4609 // Solve mortar system
4610 //--------------------
4611 for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4612 {
4613 for (unsigned j = 0; j < n_master_nodes; j++)
4614 {
4615 P[i][j] /= diag_M[i];
4616 }
4617 }
4618
4619 // Create structures to hold the hanging info
4620 //-------------------------------------------
4621 Vector<HangInfo*> hang_info_pt(n_dependent_nodes);
4622 for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4623 {
4624 hang_info_pt[i] = new HangInfo(n_master_nodes);
4625 }
4626
4627 // Copy information to hanging nodes
4628 //----------------------------------
4629 for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4630 {
4631 for (unsigned j = 0; j < n_master_nodes; j++)
4632 {
4633 hang_info_pt[i]->set_master_node_pt(j, master_node_pt[j], P[i][j]);
4634 }
4635 }
4636
4637 // Set pointers to hanging info
4638 //-----------------------------
4639 for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4640 {
4641 // Check that the node shoule actually hang.
4642 // That is, if the polynomial orders of the elements at a p-type
4643 // non-conormity are both odd then the middle node on the edge is
4644 // shared but a hanging scheme has been constructed for it.
4645 bool node_is_really_shared = false;
4646 for (unsigned m = 0; m < hang_info_pt[i]->nmaster(); m++)
4647 {
4648 // We can simply check if the hanging scheme lists itself as a
4649 // master
4650 if (hang_info_pt[i]->master_node_pt(m) ==
4651 dependent_node_pt[non_vertex_pos[i]])
4652 {
4653 node_is_really_shared = true;
4654
4655#ifdef PARANOID
4656 // Also check the corresponding weight: it should be 1
4657 if (std::fabs(hang_info_pt[i]->master_weight(m) - 1.0) > 1.0e-06)
4658 {
4659 throw OomphLibError(
4660 "Something fishy here -- with shared node at a mortar vertex",
4661 "PRefineableQElemen<2,INITIAL_NNODE_1D>t::quad_hang_helper()",
4662 OOMPH_EXCEPTION_LOCATION);
4663 }
4664#endif
4665 }
4666 }
4667
4668 // Now we can make the node hang, if it isn't a shared node
4669 if (!node_is_really_shared)
4670 {
4671 dependent_node_pt[non_vertex_pos[i]]->set_hanging_pt(
4672 hang_info_pt[i], -1);
4673 }
4674 }
4675
4676 } // End of case where there are still dependent nodes
4677
4678#ifdef PARANOID
4679 // Check all dependent nodes, hanging or otherwise
4680 for (unsigned i = 0; i < n_dependent_nodes; i++)
4681 {
4682 // Check that weights sum to 1 for those that hang
4683 if (dependent_node_pt[i]->is_hanging())
4684 {
4685 // Check that weights sum to 1
4686 double weight_sum = 0.0;
4687 for (unsigned m = 0;
4688 m < dependent_node_pt[i]->hanging_pt()->nmaster();
4689 m++)
4690 {
4691 weight_sum += dependent_node_pt[i]->hanging_pt()->master_weight(m);
4692 }
4693
4694 // Warn if not
4695 if (std::fabs(weight_sum - 1.0) > 1.0e-08)
4696 {
4697 oomph_info << "Sum of master weights: " << weight_sum << std::endl;
4699 "Weights in hanging scheme do not sum to 1",
4700 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4701 OOMPH_EXCEPTION_LOCATION);
4702 }
4703 }
4704 else
4705 {
4706 // Check that this node is shared with the master element if it
4707 // isn't hanging
4708 bool is_master = false;
4709 for (unsigned n = 0; n < n_master_nodes; n++)
4710 {
4711 if (dependent_node_pt[i] == master_node_pt[n])
4712 {
4713 // Node is a master
4714 is_master = true;
4715 break;
4716 }
4717 }
4718
4719 if (!is_master)
4720 {
4721 // Throw error
4722 std::ostringstream error_string;
4723 error_string << "This node in the dependent element is neither"
4724 << std::endl
4725 << "hanging or shared with a master element."
4726 << std::endl;
4727
4728 throw OomphLibError(
4729 error_string.str(),
4730 "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4731 OOMPH_EXCEPTION_LOCATION);
4732 }
4733 }
4734 }
4735#endif
4736
4737 // Finally, Loop over all dependent nodes and fine-tune their positions
4738 //-----------------------------------------------------------------
4739 // Here we simply set the node's positions to be consistent
4740 // with the hanging scheme. This is not strictly necessary
4741 // because it is done in the mesh adaptation before the node
4742 // becomes non-hanging later on. We make no attempt to ensure
4743 // (strong) continuity in the position across the mortar.
4744 for (unsigned i = 0; i < n_dependent_nodes; i++)
4745 {
4746 // Only fine-tune hanging nodes
4747 if (dependent_node_pt[i]->is_hanging())
4748 {
4749 // If we are doing the position, then
4750 if (value_id == -1)
4751 {
4752 // Get the local coordinate of this dependent node
4753 Vector<double> s_local(2);
4754 this->local_coordinate_of_node(dependent_node_number[i], s_local);
4755
4756 // Get the position from interpolation in this element via
4757 // the hanging scheme
4758 Vector<double> x_in_neighb(2);
4759 this->interpolated_x(s_local, x_in_neighb);
4760
4761 // Fine adjust the coordinates (macro map will pick up boundary
4762 // accurately but will lead to different element edges)
4763 dependent_node_pt[i]->x(0) = x_in_neighb[0];
4764 dependent_node_pt[i]->x(1) = x_in_neighb[1];
4765 }
4766 }
4767 }
4768 } // End of case where this interface is to be mortared
4769 }
4770
4771 /// /////////////////////////////////////////////////////////////
4772 // 3D PRefineableQElements
4773 /// /////////////////////////////////////////////////////////////
4774
4775 /// Get local coordinates of node j in the element; vector sets its own size
4776 template<unsigned INITIAL_NNODE_1D>
4778 const unsigned& n, Vector<double>& s) const
4779 {
4780 s.resize(3);
4781 unsigned Nnode_1d = this->nnode_1d();
4782 unsigned n0 = n % Nnode_1d;
4783 unsigned n1 = unsigned(double(n) / double(Nnode_1d)) % Nnode_1d;
4784 unsigned n2 = unsigned(double(n) / double(Nnode_1d * Nnode_1d));
4785
4786 switch (Nnode_1d)
4787 {
4788 case 2:
4793 break;
4794 case 3:
4799 break;
4800 case 4:
4805 break;
4806 case 5:
4811 break;
4812 case 6:
4817 break;
4818 case 7:
4823 break;
4824 default:
4825 std::ostringstream error_message;
4826 error_message << "\nERROR: Exceeded maximum polynomial order for";
4827 error_message << "\n shape functions.\n";
4828 throw OomphLibError(error_message.str(),
4829 OOMPH_CURRENT_FUNCTION,
4830 OOMPH_EXCEPTION_LOCATION);
4831 break;
4832 }
4833 }
4834
4835 /// Get the local fractino of node j in the element
4836 template<unsigned INITIAL_NNODE_1D>
4838 const unsigned& n, Vector<double>& s_fraction)
4839 {
4840 this->local_coordinate_of_node(n, s_fraction);
4841 s_fraction[0] = 0.5 * (s_fraction[0] + 1.0);
4842 s_fraction[1] = 0.5 * (s_fraction[1] + 1.0);
4843 s_fraction[2] = 0.5 * (s_fraction[2] + 1.0);
4844 }
4845
4846 template<unsigned INITIAL_NNODE_1D>
4848 const unsigned& n1d, const unsigned& i)
4849 {
4850 switch (this->nnode_1d())
4851 {
4852 case 2:
4854 return 0.5 *
4856 case 3:
4858 return 0.5 *
4860 case 4:
4862 return 0.5 *
4864 case 5:
4866 return 0.5 *
4868 case 6:
4870 return 0.5 *
4872 case 7:
4874 return 0.5 *
4876 default:
4877 std::ostringstream error_message;
4878 error_message << "\nERROR: Exceeded maximum polynomial order for";
4879 error_message << "\n shape functions.\n";
4880 throw OomphLibError(error_message.str(),
4881 OOMPH_CURRENT_FUNCTION,
4882 OOMPH_EXCEPTION_LOCATION);
4883 return 0.0;
4884 }
4885 }
4886
4887 //==================================================================
4888 /// Return the node at the specified local coordinate
4889 //==================================================================
4890 template<unsigned INITIAL_NNODE_1D>
4892 const Vector<double>& s) const
4893 {
4894 unsigned Nnode_1d = this->nnode_1d();
4895 // Load the tolerance into a local variable
4897 // There are two possible indices.
4898 Vector<int> index(3);
4899
4900 // Loop over indices
4901 for (unsigned i = 0; i < 3; i++)
4902 {
4903 // Determine the index
4904 // -------------------
4905
4906 bool is_found = false;
4907
4908 // If we are at the lower limit, the index is zero
4909 if (std::fabs(s[i] + 1.0) < tol)
4910 {
4911 index[i] = 0;
4912 is_found = true;
4913 }
4914 // If we are at the upper limit, the index is the number of nodes minus 1
4915 else if (std::fabs(s[i] - 1.0) < tol)
4916 {
4917 index[i] = Nnode_1d - 1;
4918 is_found = true;
4919 }
4920 // Otherwise, we have to calculate the index in general
4921 else
4922 {
4923 // Compute Gauss-Lobatto-Legendre node positions
4925 Orthpoly::gll_nodes(Nnode_1d, z);
4926 // Loop over possible internal nodal positions
4927 for (unsigned n = 1; n < Nnode_1d - 1; n++)
4928 {
4929 if (std::fabs(z[n] - s[i]) < tol)
4930 {
4931 index[i] = n;
4932 is_found = true;
4933 break;
4934 }
4935 }
4936 }
4937
4938 if (!is_found)
4939 {
4940 // No matching nodes
4941 return 0;
4942 }
4943 }
4944 // If we've got here we have a node, so let's return a pointer to it
4945 return this->node_pt(index[0] + Nnode_1d * index[1] +
4946 Nnode_1d * Nnode_1d * index[2]);
4947 }
4948
4949 //===================================================================
4950 /// If a neighbouring element has already created a node at
4951 /// a position corresponding to the local fractional position within the
4952 /// present element, s_fraction, return
4953 /// a pointer to that node. If not, return NULL (0).
4954 //===================================================================
4955 template<unsigned INITIAL_NNODE_1D>
4957 const Vector<double>& s_fraction, bool& is_periodic)
4958 {
4959 using namespace OcTreeNames;
4960
4961 // Calculate the faces/edges on which the node lies
4962 Vector<int> faces;
4963 Vector<int> edges;
4964
4965 if (s_fraction[0] == 0.0)
4966 {
4967 faces.push_back(L);
4968 if (s_fraction[1] == 0.0)
4969 {
4970 edges.push_back(LD);
4971 }
4972 if (s_fraction[2] == 0.0)
4973 {
4974 edges.push_back(LB);
4975 }
4976 if (s_fraction[1] == 1.0)
4977 {
4978 edges.push_back(LU);
4979 }
4980 if (s_fraction[2] == 1.0)
4981 {
4982 edges.push_back(LF);
4983 }
4984 }
4985
4986 if (s_fraction[0] == 1.0)
4987 {
4988 faces.push_back(R);
4989 if (s_fraction[1] == 0.0)
4990 {
4991 edges.push_back(RD);
4992 }
4993 if (s_fraction[2] == 0.0)
4994 {
4995 edges.push_back(RB);
4996 }
4997 if (s_fraction[1] == 1.0)
4998 {
4999 edges.push_back(RU);
5000 }
5001 if (s_fraction[2] == 1.0)
5002 {
5003 edges.push_back(RF);
5004 }
5005 }
5006
5007 if (s_fraction[1] == 0.0)
5008 {
5009 faces.push_back(D);
5010 if (s_fraction[2] == 0.0)
5011 {
5012 edges.push_back(DB);
5013 }
5014 if (s_fraction[2] == 1.0)
5015 {
5016 edges.push_back(DF);
5017 }
5018 }
5019
5020 if (s_fraction[1] == 1.0)
5021 {
5022 faces.push_back(U);
5023 if (s_fraction[2] == 0.0)
5024 {
5025 edges.push_back(UB);
5026 }
5027 if (s_fraction[2] == 1.0)
5028 {
5029 edges.push_back(UF);
5030 }
5031 }
5032
5033 if (s_fraction[2] == 0.0)
5034 {
5035 faces.push_back(B);
5036 }
5037
5038 if (s_fraction[2] == 1.0)
5039 {
5040 faces.push_back(F);
5041 }
5042
5043 // Find the number of faces
5044 unsigned n_face = faces.size();
5045
5046 // Find the number of edges
5047 unsigned n_edge = edges.size();
5048
5049 Vector<unsigned> translate_s(3);
5050 Vector<double> s_lo_neigh(3);
5051 Vector<double> s_hi_neigh(3);
5052 Vector<double> s(3);
5053
5054 int neigh_face, diff_level;
5055 bool in_neighbouring_tree;
5056
5057 // Loop over the faces on which the node lies
5058 //------------------------------------------
5059 for (unsigned j = 0; j < n_face; j++)
5060 {
5061 // Find pointer to neighbouring element along face
5062 OcTree* neigh_pt;
5063 neigh_pt = octree_pt()->gteq_face_neighbour(faces[j],
5064 translate_s,
5065 s_lo_neigh,
5066 s_hi_neigh,
5067 neigh_face,
5068 diff_level,
5069 in_neighbouring_tree);
5070
5071 // Neighbour exists
5072 if (neigh_pt != 0)
5073 {
5074 // Have any of its vertex nodes been created yet?
5075 // (Must look in incomplete neighbours because after the
5076 // pre-build they may contain pointers to the required nodes. e.g.
5077 // h-refinement of neighbouring linear and quadratic elements)
5078 bool a_vertex_node_is_built = false;
5079 QElement<3, INITIAL_NNODE_1D>* neigh_obj_pt =
5080 dynamic_cast<QElement<3, INITIAL_NNODE_1D>*>(neigh_pt->object_pt());
5081 if (neigh_obj_pt == 0)
5082 {
5083 throw OomphLibError("Not a quad element!",
5084 OOMPH_CURRENT_FUNCTION,
5085 OOMPH_EXCEPTION_LOCATION);
5086 }
5087 for (unsigned vnode = 0; vnode < neigh_obj_pt->nvertex_node(); vnode++)
5088 {
5089 if (neigh_obj_pt->vertex_node_pt(vnode) != 0)
5090 a_vertex_node_is_built = true;
5091 break;
5092 }
5093 if (a_vertex_node_is_built)
5094 {
5095 // We now need to translate the nodal location, defined in terms
5096 // of the fractional coordinates of the present element into
5097 // those of its neighbour. For this we use the information returned
5098 // to use from the octree function.
5099
5100 // Calculate the local coordinate in the neighbour
5101 // Note that we need to use the translation scheme in case
5102 // the local coordinates are swapped in the neighbour.
5103 for (unsigned i = 0; i < 3; i++)
5104 {
5105 s[i] = s_lo_neigh[i] +
5106 s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
5107 }
5108
5109 // Find the node in the neighbour
5110 Node* neighbour_node_pt =
5112
5113 // If there is a node, return it
5114 if (neighbour_node_pt != 0)
5115 {
5116 // Now work out whether it's a periodic boundary
5117 // only possible if we have moved into a neighbouring tree
5118 if (in_neighbouring_tree)
5119 {
5120 // Return whether the neighbour is periodic
5121 is_periodic =
5122 octree_pt()->root_pt()->is_neighbour_periodic(faces[j]);
5123 }
5124
5125 return neighbour_node_pt;
5126 }
5127 }
5128 }
5129 } // End of loop over faces
5130
5131
5132 // Loop over the edges on which the node lies
5133 //------------------------------------------
5134 for (unsigned j = 0; j < n_edge; j++)
5135 {
5136 // Even if we restrict ourselves to true edge neighbours (i.e.
5137 // elements that are not also face neighbours) there may be multiple
5138 // edge neighbours across the edges between multiple root octrees.
5139 // When making the first call to OcTree::gteq_true_edge_neighbour(...)
5140 // we simply return the first one of these multiple edge neighbours
5141 // (if there are any at all, of course) and also return the total number
5142 // of true edge neighbours. If the node in question already exists
5143 // on the first edge neighbour we're done. If it doesn't it may exist
5144 // on other edge neighbours so we repeat the process over all
5145 // other edge neighbours (bailing out if a node is found, of course).
5146
5147 // Initially return the zero-th true edge neighbour
5148 unsigned i_root_edge_neighbour = 0;
5149
5150 // Initialise the total number of true edge neighbours
5151 unsigned nroot_edge_neighbour = 0;
5152
5153 // Keep searching until we've found the node or until we've checked
5154 // all available edge neighbours
5155 bool keep_searching = true;
5156 while (keep_searching)
5157 {
5158 // Find pointer to neighbouring element along edge
5159 OcTree* neigh_pt;
5160 neigh_pt = octree_pt()->gteq_true_edge_neighbour(edges[j],
5161 i_root_edge_neighbour,
5162 nroot_edge_neighbour,
5163 translate_s,
5164 s_lo_neigh,
5165 s_hi_neigh,
5166 neigh_face,
5167 diff_level);
5168
5169 // Neighbour exists
5170 if (neigh_pt != 0)
5171 {
5172 // Have any of its vertex nodes been created yet?
5173 // (Must look in incomplete neighbours because after the
5174 // pre-build they may contain pointers to the required nodes. e.g.
5175 // h-refinement of neighbouring linear and quadratic elements)
5176 bool a_vertex_node_is_built = false;
5177 QElement<3, INITIAL_NNODE_1D>* neigh_obj_pt =
5178 dynamic_cast<QElement<3, INITIAL_NNODE_1D>*>(neigh_pt->object_pt());
5179 if (neigh_obj_pt == 0)
5180 {
5181 throw OomphLibError("Not a quad element!",
5182 OOMPH_CURRENT_FUNCTION,
5183 OOMPH_EXCEPTION_LOCATION);
5184 }
5185 for (unsigned vnode = 0; vnode < neigh_obj_pt->nvertex_node();
5186 vnode++)
5187 {
5188 if (neigh_obj_pt->vertex_node_pt(vnode) != 0)
5189 a_vertex_node_is_built = true;
5190 break;
5191 }
5192 if (a_vertex_node_is_built)
5193 {
5194 // We now need to translate the nodal location, defined in terms
5195 // of the fractional coordinates of the present element into
5196 // those of its neighbour. For this we use the information returned
5197 // to use from the octree function.
5198
5199 // Calculate the local coordinate in the neighbour
5200 // Note that we need to use the translation scheme in case
5201 // the local coordinates are swapped in the neighbour.
5202 for (unsigned i = 0; i < 3; i++)
5203 {
5204 s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
5205 (s_hi_neigh[i] - s_lo_neigh[i]);
5206 }
5207
5208 // Find the node in the neighbour
5209 Node* neighbour_node_pt =
5211
5212 // If there is a node, return it
5213 if (neighbour_node_pt != 0)
5214 {
5215 // Get the faces on which the edge lies
5216 Vector<int> faces_attached_to_edge =
5218
5219 // Get the number of entries in the vector
5220 unsigned n_faces_attached_to_edge = faces_attached_to_edge.size();
5221
5222 // Loop over the faces
5223 for (unsigned i_face = 0; i_face < n_faces_attached_to_edge;
5224 i_face++)
5225 {
5226 // Is the node periodic in the face direction?
5227 is_periodic = octree_pt()->root_pt()->is_neighbour_periodic(
5228 faces_attached_to_edge[i_face]);
5229
5230 // Check if the edge is periodic in the i_face-th face direction
5231 if (is_periodic)
5232 {
5233 // We're done!
5234 break;
5235 }
5236 } // for (unsigned
5237 // i_face=0;i_face<n_faces_attached_to_edge;i_face++)
5238
5239 return neighbour_node_pt;
5240 }
5241 }
5242 }
5243
5244 // Keep searching, but only if there are further edge neighbours
5245 // Try next root edge neighbour
5246 i_root_edge_neighbour++;
5247
5248 // Have we exhausted the search
5249 if (i_root_edge_neighbour >= nroot_edge_neighbour)
5250 {
5251 keep_searching = false;
5252 }
5253
5254 } // End of while keep searching over all true edge neighbours
5255
5256 } // End of loop over edges
5257
5258 // Node not found, return null
5259 return 0;
5260 }
5261
5262 //===================================================================
5263 /// If a neighbouring element's son has already created a node at
5264 /// a position corresponding to the local fractional position within the
5265 /// present element, s_fraction, return
5266 /// a pointer to that node. If not, return NULL (0). If the node is
5267 /// periodic the flag is_periodic will be true
5268 //===================================================================
5269 template<unsigned INITIAL_NNODE_1D>
5272 bool& is_periodic)
5273 {
5274 using namespace OcTreeNames;
5275
5276 // Calculate the faces/edges on which the node lies
5277 Vector<int> faces;
5278 Vector<int> edges;
5279
5280 if (s_fraction[0] == 0.0)
5281 {
5282 faces.push_back(L);
5283 if (s_fraction[1] == 0.0)
5284 {
5285 edges.push_back(LD);
5286 }
5287 if (s_fraction[2] == 0.0)
5288 {
5289 edges.push_back(LB);
5290 }
5291 if (s_fraction[1] == 1.0)
5292 {
5293 edges.push_back(LU);
5294 }
5295 if (s_fraction[2] == 1.0)
5296 {
5297 edges.push_back(LF);
5298 }
5299 }
5300
5301 if (s_fraction[0] == 1.0)
5302 {
5303 faces.push_back(R);
5304 if (s_fraction[1] == 0.0)
5305 {
5306 edges.push_back(RD);
5307 }
5308 if (s_fraction[2] == 0.0)
5309 {
5310 edges.push_back(RB);
5311 }
5312 if (s_fraction[1] == 1.0)
5313 {
5314 edges.push_back(RU);
5315 }
5316 if (s_fraction[2] == 1.0)
5317 {
5318 edges.push_back(RF);
5319 }
5320 }
5321
5322 if (s_fraction[1] == 0.0)
5323 {
5324 faces.push_back(D);
5325 if (s_fraction[2] == 0.0)
5326 {
5327 edges.push_back(DB);
5328 }
5329 if (s_fraction[2] == 1.0)
5330 {
5331 edges.push_back(DF);
5332 }
5333 }
5334
5335 if (s_fraction[1] == 1.0)
5336 {
5337 faces.push_back(U);
5338 if (s_fraction[2] == 0.0)
5339 {
5340 edges.push_back(UB);
5341 }
5342 if (s_fraction[2] == 1.0)
5343 {
5344 edges.push_back(UF);
5345 }
5346 }
5347
5348 if (s_fraction[2] == 0.0)
5349 {
5350 faces.push_back(B);
5351 }
5352
5353 if (s_fraction[2] == 1.0)
5354 {
5355 faces.push_back(F);
5356 }
5357
5358 // Find the number of faces
5359 unsigned n_face = faces.size();
5360
5361 // Find the number of edges
5362 unsigned n_edge = edges.size();
5363
5364 Vector<unsigned> translate_s(3);
5365 Vector<double> s_lo_neigh(3);
5366 Vector<double> s_hi_neigh(3);
5367 Vector<double> s(3);
5368
5369 int neigh_face, diff_level;
5370 bool in_neighbouring_tree;
5371
5372 // Loop over the faces on which the node lies
5373 //------------------------------------------
5374 for (unsigned j = 0; j < n_face; j++)
5375 {
5376 // Find pointer to neighbouring element along face
5377 OcTree* neigh_pt;
5378 neigh_pt = octree_pt()->gteq_face_neighbour(faces[j],
5379 translate_s,
5380 s_lo_neigh,
5381 s_hi_neigh,
5382 neigh_face,
5383 diff_level,
5384 in_neighbouring_tree);
5385
5386 // Neighbour exists
5387 if (neigh_pt != 0)
5388 {
5389 // Have its nodes been created yet?
5390 // (Must look in sons of unfinished neighbours too!!!)
5391 if (true)
5392 {
5393 // We now need to translate the nodal location, defined in terms
5394 // of the fractional coordinates of the present element into
5395 // those of its neighbour. For this we use the information returned
5396 // to use from the octree function.
5397
5398 // Calculate the local coordinate in the neighbour
5399 // Note that we need to use the translation scheme in case
5400 // the local coordinates are swapped in the neighbour.
5401 for (unsigned i = 0; i < 3; i++)
5402 {
5403 s[i] = s_lo_neigh[i] +
5404 s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
5405 }
5406
5407 // Check if the element has sons
5408 if (neigh_pt->nsons() != 0)
5409 {
5410 // First, find the son element in which the node should live
5411
5412 // Find coordinates in the sons
5413 Vector<double> s_in_son(3);
5414 int son = -10;
5415 // On the left
5416 if (s[0] < 0.0)
5417 {
5418 // On the bottom
5419 if (s[1] < 0.0)
5420 {
5421 // On the back
5422 if (s[2] < 0.0)
5423 {
5424 // It's the LDB son
5425 son = OcTreeNames::LDB;
5426 s_in_son[0] = 1.0 + 2.0 * s[0];
5427 s_in_son[1] = 1.0 + 2.0 * s[1];
5428 s_in_son[2] = 1.0 + 2.0 * s[2];
5429 }
5430 // On the front
5431 else
5432 {
5433 // It's the LDF son
5434 son = OcTreeNames::LDF;
5435 s_in_son[0] = 1.0 + 2.0 * s[0];
5436 s_in_son[1] = 1.0 + 2.0 * s[1];
5437 s_in_son[2] = -1.0 + 2.0 * s[2];
5438 }
5439 }
5440 // On the top
5441 else
5442 {
5443 // On the back
5444 if (s[2] < 0.0)
5445 {
5446 // It's the LUB son
5447 son = OcTreeNames::LUB;
5448 s_in_son[0] = 1.0 + 2.0 * s[0];
5449 s_in_son[1] = -1.0 + 2.0 * s[1];
5450 s_in_son[2] = 1.0 + 2.0 * s[2];
5451 }
5452 // On the front
5453 else
5454 {
5455 // It's the LUF son
5456 son = OcTreeNames::LUF;
5457 s_in_son[0] = 1.0 + 2.0 * s[0];
5458 s_in_son[1] = -1.0 + 2.0 * s[1];
5459 s_in_son[2] = -1.0 + 2.0 * s[2];
5460 }
5461 }
5462 }
5463 // On the right
5464 else
5465 {
5466 // On the bottom
5467 if (s[1] < 0.0)
5468 {
5469 // On the back
5470 if (s[2] < 0.0)
5471 {
5472 // It's the RDB son
5473 son = OcTreeNames::RDB;
5474 s_in_son[0] = -1.0 + 2.0 * s[0];
5475 s_in_son[1] = 1.0 + 2.0 * s[1];
5476 s_in_son[2] = 1.0 + 2.0 * s[2];
5477 }
5478 // On the front
5479 else
5480 {
5481 // It's the RDF son
5482 son = OcTreeNames::RDF;
5483 s_in_son[0] = -1.0 + 2.0 * s[0];
5484 s_in_son[1] = 1.0 + 2.0 * s[1];
5485 s_in_son[2] = -1.0 + 2.0 * s[2];
5486 }
5487 }
5488 // On the top
5489 else
5490 {
5491 // On the back
5492 if (s[2] < 0.0)
5493 {
5494 // It's the RUB son
5495 son = OcTreeNames::RUB;
5496 s_in_son[0] = -1.0 + 2.0 * s[0];
5497 s_in_son[1] = -1.0 + 2.0 * s[1];
5498 s_in_son[2] = 1.0 + 2.0 * s[2];
5499 }
5500 // On the front
5501 else
5502 {
5503 // It's the RUF son
5504 son = OcTreeNames::RUF;
5505 s_in_son[0] = -1.0 + 2.0 * s[0];
5506 s_in_son[1] = -1.0 + 2.0 * s[1];
5507 s_in_son[2] = -1.0 + 2.0 * s[2];
5508 }
5509 }
5510 }
5511
5512 // Find the node in the neighbour's son
5513 Node* neighbour_son_node_pt =
5515 s_in_son);
5516
5517 // If there is a node, return it
5518 if (neighbour_son_node_pt != 0)
5519 {
5520 // Now work out whether it's a periodic boundary
5521 // only possible if we have moved into a neighbouring tree
5522 if (in_neighbouring_tree)
5523 {
5524 // Return whether the neighbour is periodic
5525 is_periodic =
5526 octree_pt()->root_pt()->is_neighbour_periodic(faces[j]);
5527 }
5528
5529 // Return the pointer to the neighbouring node
5530 return neighbour_son_node_pt;
5531 }
5532 }
5533 }
5534 }
5535 } // End of loop over faces
5536
5537
5538 // Loop over the edges on which the node lies
5539 //------------------------------------------
5540 for (unsigned j = 0; j < n_edge; j++)
5541 {
5542 // Even if we restrict ourselves to true edge neighbours (i.e.
5543 // elements that are not also face neighbours) there may be multiple
5544 // edge neighbours across the edges between multiple root octrees.
5545 // When making the first call to OcTree::gteq_true_edge_neighbour(...)
5546 // we simply return the first one of these multiple edge neighbours
5547 // (if there are any at all, of course) and also return the total number
5548 // of true edge neighbours. If the node in question already exists
5549 // on the first edge neighbour we're done. If it doesn't it may exist
5550 // on other edge neighbours so we repeat the process over all
5551 // other edge neighbours (bailing out if a node is found, of course).
5552
5553 // Initially return the zero-th true edge neighbour
5554 unsigned i_root_edge_neighbour = 0;
5555
5556 // Initialise the total number of true edge neighbours
5557 unsigned nroot_edge_neighbour = 0;
5558
5559 // Keep searching until we've found the node or until we've checked
5560 // all available edge neighbours
5561 bool keep_searching = true;
5562 while (keep_searching)
5563 {
5564 // Find pointer to neighbouring element along edge
5565 OcTree* neigh_pt;
5566 neigh_pt = octree_pt()->gteq_true_edge_neighbour(edges[j],
5567 i_root_edge_neighbour,
5568 nroot_edge_neighbour,
5569 translate_s,
5570 s_lo_neigh,
5571 s_hi_neigh,
5572 neigh_face,
5573 diff_level);
5574
5575 // Neighbour exists
5576 if (neigh_pt != 0)
5577 {
5578 // Have its nodes been created yet?
5579 // (Must look in sons of unfinished neighbours too!!!)
5580 if (true)
5581 {
5582 // We now need to translate the nodal location, defined in terms
5583 // of the fractional coordinates of the present element into
5584 // those of its neighbour. For this we use the information returned
5585 // to use from the octree function.
5586
5587 // Calculate the local coordinate in the neighbour
5588 // Note that we need to use the translation scheme in case
5589 // the local coordinates are swapped in the neighbour.
5590 for (unsigned i = 0; i < 3; i++)
5591 {
5592 s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
5593 (s_hi_neigh[i] - s_lo_neigh[i]);
5594 }
5595
5596 // Check if the element has sons
5597 if (neigh_pt->nsons() != 0)
5598 {
5599 // First, find the son element in which the node should live
5600
5601 // Find coordinates in the sons
5602 Vector<double> s_in_son(3);
5603 int son = -10;
5604 // On the left
5605 if (s[0] < 0.0)
5606 {
5607 // On the bottom
5608 if (s[1] < 0.0)
5609 {
5610 // On the back
5611 if (s[2] < 0.0)
5612 {
5613 // It's the LDB son
5614 son = OcTreeNames::LDB;
5615 s_in_son[0] = 1.0 + 2.0 * s[0];
5616 s_in_son[1] = 1.0 + 2.0 * s[1];
5617 s_in_son[2] = 1.0 + 2.0 * s[2];
5618 }
5619 // On the front
5620 else
5621 {
5622 // It's the LDF son
5623 son = OcTreeNames::LDF;
5624 s_in_son[0] = 1.0 + 2.0 * s[0];
5625 s_in_son[1] = 1.0 + 2.0 * s[1];
5626 s_in_son[2] = -1.0 + 2.0 * s[2];
5627 }
5628 }
5629 // On the top
5630 else
5631 {
5632 // On the back
5633 if (s[2] < 0.0)
5634 {
5635 // It's the LUB son
5636 son = OcTreeNames::LUB;
5637 s_in_son[0] = 1.0 + 2.0 * s[0];
5638 s_in_son[1] = -1.0 + 2.0 * s[1];
5639 s_in_son[2] = 1.0 + 2.0 * s[2];
5640 }
5641 // On the front
5642 else
5643 {
5644 // It's the LUF son
5645 son = OcTreeNames::LUF;
5646 s_in_son[0] = 1.0 + 2.0 * s[0];
5647 s_in_son[1] = -1.0 + 2.0 * s[1];
5648 s_in_son[2] = -1.0 + 2.0 * s[2];
5649 }
5650 }
5651 }
5652 // On the right
5653 else
5654 {
5655 // On the bottom
5656 if (s[1] < 0.0)
5657 {
5658 // On the back
5659 if (s[2] < 0.0)
5660 {
5661 // It's the RDB son
5662 son = OcTreeNames::RDB;
5663 s_in_son[0] = -1.0 + 2.0 * s[0];
5664 s_in_son[1] = 1.0 + 2.0 * s[1];
5665 s_in_son[2] = 1.0 + 2.0 * s[2];
5666 }
5667 // On the front
5668 else
5669 {
5670 // It's the RDF son
5671 son = OcTreeNames::RDF;
5672 s_in_son[0] = -1.0 + 2.0 * s[0];
5673 s_in_son[1] = 1.0 + 2.0 * s[1];
5674 s_in_son[2] = -1.0 + 2.0 * s[2];
5675 }
5676 }
5677 // On the top
5678 else
5679 {
5680 // On the back
5681 if (s[2] < 0.0)
5682 {
5683 // It's the RUB son
5684 son = OcTreeNames::RUB;
5685 s_in_son[0] = -1.0 + 2.0 * s[0];
5686 s_in_son[1] = -1.0 + 2.0 * s[1];
5687 s_in_son[2] = 1.0 + 2.0 * s[2];
5688 }
5689 // On the front
5690 else
5691 {
5692 // It's the RUF son
5693 son = OcTreeNames::RUF;
5694 s_in_son[0] = -1.0 + 2.0 * s[0];
5695 s_in_son[1] = -1.0 + 2.0 * s[1];
5696 s_in_son[2] = -1.0 + 2.0 * s[2];
5697 }
5698 }
5699 }
5700
5701 // Find the node in the neighbour's son
5702 Node* neighbour_son_node_pt =
5703 neigh_pt->son_pt(son)
5704 ->object_pt()
5705 ->get_node_at_local_coordinate(s_in_son);
5706
5707 // If there is a node, return it
5708 if (neighbour_son_node_pt != 0)
5709 {
5710 // Get the faces on which the edge lies
5711 Vector<int> faces_attached_to_edge =
5713
5714 // Get the number of entries in the vector
5715 unsigned n_faces_attached_to_edge =
5716 faces_attached_to_edge.size();
5717
5718 // Loop over the faces
5719 for (unsigned i_face = 0; i_face < n_faces_attached_to_edge;
5720 i_face++)
5721 {
5722 // Is the node periodic in the face direction?
5723 is_periodic = octree_pt()->root_pt()->is_neighbour_periodic(
5724 faces_attached_to_edge[i_face]);
5725
5726 // Check if the edge is periodic in the i_face-th face
5727 // direction
5728 if (is_periodic)
5729 {
5730 // We're done!
5731 break;
5732 }
5733 } // for (unsigned
5734 // i_face=0;i_face<n_faces_attached_to_edge;i_face++)
5735
5736 // Return the pointer to the neighbouring node
5737 return neighbour_son_node_pt;
5738 }
5739 }
5740 }
5741 }
5742
5743 // Keep searching, but only if there are further edge neighbours
5744 // Try next root edge neighbour
5745 i_root_edge_neighbour++;
5746
5747 // Have we exhausted the search
5748 if (i_root_edge_neighbour >= nroot_edge_neighbour)
5749 {
5750 keep_searching = false;
5751 }
5752
5753 } // End of while keep searching over all true edge neighbours
5754
5755 } // End of loop over edges
5756
5757 // Node not found, return null
5758 return 0;
5759 }
5760
5761 //==================================================================
5762 /// Set the correct p-order of the element based on that of its
5763 /// father. Then construct an integration scheme of the correct order.
5764 /// If an adopted father is specified, information from this is
5765 /// used instead of using the father found from the tree.
5766 //==================================================================
5767 template<unsigned INITIAL_NNODE_1D>
5769 Tree* const& adopted_father_pt, const unsigned& initial_p_order)
5770 {
5771 // Storage for pointer to my father (in octree impersonation)
5772 OcTree* father_pt = 0;
5773
5774 // Check if an adopted father has been specified
5775 if (adopted_father_pt != 0)
5776 {
5777 // Get pointer to my father (in octree impersonation)
5778 father_pt = dynamic_cast<OcTree*>(adopted_father_pt);
5779 }
5780 // Check if element is in a tree
5781 else if (Tree_pt != 0)
5782 {
5783 // Get pointer to my father (in octree impersonation)
5784 father_pt = dynamic_cast<OcTree*>(octree_pt()->father_pt());
5785 }
5786 else
5787 {
5788 // throw OomphLibError(
5789 // "Element not in a tree, and no adopted father has been
5790 // specified!",
5791 // "PRefineableQElement<2,INITIAL_NNODE_1D>::initial_setup()",
5792 // OOMPH_EXCEPTION_LOCATION);
5793 }
5794
5795 // Check if element has father
5796 if (father_pt != 0 || initial_p_order != 0)
5797 {
5798 if (father_pt != 0)
5799 {
5802 father_pt->object_pt());
5803 if (father_el_pt != 0)
5804 {
5805 unsigned father_p_order = father_el_pt->p_order();
5806 // Set p-order to that of father
5807 P_order = father_p_order;
5808 }
5809 }
5810 else
5811 {
5812 P_order = initial_p_order;
5813 }
5814
5815 // Now sort out the element...
5816 // (has p^3 nodes)
5817 unsigned new_n_node = P_order * P_order * P_order;
5818
5819 // Allocate new space for Nodes (at the element level)
5820 this->set_n_node(new_n_node);
5821
5822 // Set integration scheme
5823 delete this->integral_pt();
5824 switch (P_order)
5825 {
5826 case 2:
5827 this->set_integration_scheme(new GaussLobattoLegendre<3, 2>);
5828 break;
5829 case 3:
5830 this->set_integration_scheme(new GaussLobattoLegendre<3, 3>);
5831 break;
5832 case 4:
5833 this->set_integration_scheme(new GaussLobattoLegendre<3, 4>);
5834 break;
5835 case 5:
5836 this->set_integration_scheme(new GaussLobattoLegendre<3, 5>);
5837 break;
5838 case 6:
5839 this->set_integration_scheme(new GaussLobattoLegendre<3, 6>);
5840 break;
5841 case 7:
5842 this->set_integration_scheme(new GaussLobattoLegendre<3, 7>);
5843 break;
5844 default:
5845 std::ostringstream error_message;
5846 error_message << "\nERROR: Exceeded maximum polynomial order for";
5847 error_message << "\n integration scheme.\n";
5848 throw OomphLibError(error_message.str(),
5849 OOMPH_CURRENT_FUNCTION,
5850 OOMPH_EXCEPTION_LOCATION);
5851 }
5852 }
5853 }
5854
5855 //==================================================================
5856 /// Check the father element for any required nodes which
5857 /// already exist
5858 //==================================================================
5859 template<unsigned INITIAL_NNODE_1D>
5861 Mesh*& mesh_pt, Vector<Node*>& new_node_pt)
5862 {
5863 using namespace OcTreeNames;
5864
5865 // Get the number of 1d nodes
5866 unsigned n_p = nnode_1d();
5867
5868 // Check whether static father_bound needs to be created
5869 if (Father_bound[n_p].nrow() == 0)
5870 {
5871 setup_father_bounds();
5872 }
5873
5874 // Pointer to my father (in quadtree impersonation)
5875 OcTree* father_pt = dynamic_cast<OcTree*>(octree_pt()->father_pt());
5876
5877 // What type of son am I? Ask my quadtree representation...
5878 int son_type = Tree_pt->son_type();
5879
5880 // Has somebody build me already? (If any nodes have not been built)
5881 if (!nodes_built())
5882 {
5883#ifdef PARANOID
5884 if (father_pt == 0)
5885 {
5886 std::string error_message =
5887 "Something fishy here: I have no father and yet \n";
5888 error_message += "I have no nodes. Who has created me then?!\n";
5889
5890 throw OomphLibError(
5891 error_message,
5892 "PRefineableQElement<3,INITIAL_NNODE_1D>::pre_build()",
5893 OOMPH_EXCEPTION_LOCATION);
5894 }
5895#endif
5896
5897 // Return pointer to father element
5898 RefineableQElement<3>* father_el_pt =
5899 dynamic_cast<RefineableQElement<3>*>(father_pt->object_pt());
5900
5901 // Timestepper should be the same for all nodes in father
5902 // element -- use it create timesteppers for new nodes
5903 TimeStepper* time_stepper_pt =
5904 father_el_pt->node_pt(0)->time_stepper_pt();
5905
5906 // Number of history values (incl. present)
5907 unsigned ntstorage = time_stepper_pt->ntstorage();
5908
5909 /// / Pass pointer to time object:
5910 // this->time_pt()=father_el_pt->time_pt();
5911
5912 Vector<double> s_lo(3);
5913 Vector<double> s_hi(3);
5914 Vector<double> s(3);
5915 Vector<double> x(3);
5916
5917 // Setup vertex coordinates in father element:
5918 //--------------------------------------------
5919 switch (son_type)
5920 {
5921 case LDB:
5922 s_lo[0] = -1.0;
5923 s_hi[0] = 0.0;
5924 s_lo[1] = -1.0;
5925 s_hi[1] = 0.0;
5926 s_lo[2] = -1.0;
5927 s_hi[2] = 0.0;
5928 break;
5929
5930 case LDF:
5931 s_lo[0] = -1.0;
5932 s_hi[0] = 0.0;
5933 s_lo[1] = -1.0;
5934 s_hi[1] = 0.0;
5935 s_lo[2] = 0.0;
5936 s_hi[2] = 1.0;
5937 break;
5938
5939 case LUB:
5940 s_lo[0] = -1.0;
5941 s_hi[0] = 0.0;
5942 s_lo[1] = 0.0;
5943 s_hi[1] = 1.0;
5944 s_lo[2] = -1.0;
5945 s_hi[2] = 0.0;
5946 break;
5947
5948 case LUF:
5949 s_lo[0] = -1.0;
5950 s_hi[0] = 0.0;
5951 s_lo[1] = 0.0;
5952 s_hi[1] = 1.0;
5953 s_lo[2] = 0.0;
5954 s_hi[2] = 1.0;
5955 break;
5956
5957 case RDB:
5958 s_lo[0] = 0.0;
5959 s_hi[0] = 1.0;
5960 s_lo[1] = -1.0;
5961 s_hi[1] = 0.0;
5962 s_lo[2] = -1.0;
5963 s_hi[2] = 0.0;
5964 break;
5965
5966 case RDF:
5967 s_lo[0] = 0.0;
5968 s_hi[0] = 1.0;
5969 s_lo[1] = -1.0;
5970 s_hi[1] = 0.0;
5971 s_lo[2] = 0.0;
5972 s_hi[2] = 1.0;
5973 break;
5974
5975 case RUB:
5976 s_lo[0] = 0.0;
5977 s_hi[0] = 1.0;
5978 s_lo[1] = 0.0;
5979 s_hi[1] = 1.0;
5980 s_lo[2] = -1.0;
5981 s_hi[2] = 0.0;
5982 break;
5983
5984 case RUF:
5985 s_lo[0] = 0.0;
5986 s_hi[0] = 1.0;
5987 s_lo[1] = 0.0;
5988 s_hi[1] = 1.0;
5989 s_lo[2] = 0.0;
5990 s_hi[2] = 1.0;
5991 break;
5992 }
5993
5994 /// / Pass macro element pointer on to sons and
5995 /// / set coordinates in macro element
5996 /// / hierher why can I see this?
5997 // if(father_el_pt->macro_elem_pt()!=0)
5998 // {
5999 // set_macro_elem_pt(father_el_pt->macro_elem_pt());
6000 // for(unsigned i=0;i<2;i++)
6001 // {
6002 // s_macro_ll(i)= father_el_pt->s_macro_ll(i)+
6003 // 0.5*(s_lo[i]+1.0)*(father_el_pt->s_macro_ur(i)-
6004 // father_el_pt->s_macro_ll(i));
6005 // s_macro_ur(i)= father_el_pt->s_macro_ll(i)+
6006 // 0.5*(s_hi[i]+1.0)*(father_el_pt->s_macro_ur(i)-
6007 // father_el_pt->s_macro_ll(i));
6008 // }
6009 // }
6010
6011
6012 // If the father element hasn't been generated yet, we're stuck...
6013 if (father_el_pt->node_pt(0) == 0)
6014 {