refineable_elements.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for classes that define refineable element objects
27
28// Include guard to prevent multiple inclusions of the header
29#ifndef OOMPH_REFINEABLE_ELEMENTS_HEADER
30#define OOMPH_REFINEABLE_ELEMENTS_HEADER
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37#include "elements.h"
38#include "tree.h"
39
40namespace oomph
41{
42 class Mesh;
43
44 //=======================================================================
45 /// RefineableElements are FiniteElements that may be subdivided into
46 /// children to provide a better local approximation to the solution.
47 /// After non-uniform refinement adjacent elements need not necessarily have
48 /// nodes in common. A node that does not have a counterpart in its
49 /// neighbouring element is known as a hanging node and its position and
50 /// any data that it stores must be constrained to ensure inter-element
51 /// continuity.
52 ///
53 /// Generic data and function interfaces associated with refinement
54 /// are defined in this class.
55 ///
56 /// Additional data includes:
57 /// - a pointer to a general Tree object that is used to track the
58 /// refinement history,
59 /// - a refinement level (not necessarily the same as the level in the tree!),
60 /// - a flag indicating whether the element should be refined,
61 /// - a flag indicating whether the element should be de-refined,
62 /// - a global element number for plotting/validation purposes,
63 /// - storage for local equation numbers associated with hanging nodes.
64 ///
65 /// Additional functions perform the following generic tasks:
66 /// - provide access to additional data,
67 /// - setup local equation numbering for data associated with hanging nodes,
68 /// - generic finite-difference calculation of contributions to the
69 /// elemental jacobian from nodal data to include hanging nodes,
70 /// - split of the element into its sons,
71 /// - select and deselect the element for refinement,
72 /// - select and deselect the sons of the element for de-refinement (merging),
73 ///
74 /// In addition, there are a number of interfaces that specify
75 /// element-specific tasks. These should be overloaded in RefineableElements
76 /// of particular geometric types and perform the following tasks:
77 /// - return a pointer to the root and father elements in the tree structure,
78 /// - define the number of sons into which the element is divided,
79 /// - build the element: construct nodes, assign their positions,
80 /// values and any boundary conditions,
81 /// - recreate the element from its sons if they are merged,
82 /// - deactivate the element, perform any operations that are
83 /// required when the element is still in the tree, but no longer active
84 /// - set the number and provide access to the values interpolated
85 /// by the nodes,
86 /// - setup the hanging nodes
87 ///
88 /// In mixed element different sets of nodes are used to interpolate different
89 /// unknowns. Interfaces are provided for functions that can be used to find
90 /// the position of the nodes that interpolate the different unknowns. These
91 /// functions are used to setup hanging node information automatically in
92 /// particular elements, e.g. Taylor Hood Navier--Stokes.
93 /// The default implementation assumes that the elements are isoparameteric.
94 ///
95 //======================================================================
96 class RefineableElement : public virtual FiniteElement
97 {
98 protected:
99 /// A pointer to a general tree object
101
102 /// Refinement level
103 unsigned Refine_level;
104
105 /// Flag for refinement
107
108 /// Flag to indicate suppression of any refinement
110
111 /// Flag for unrefinement
113
114 /// Global element number -- for plotting/validation purposes
115 long Number;
116
117 /// Max. allowed discrepancy in element integrity check
119
120 /// Static helper function that is used to check that the value_id
121 /// is in range
122 static void check_value_id(const int& n_continuously_interpolated_values,
123 const int& value_id);
124
125
126 /// Assemble the jacobian matrix for the mapping from local
127 /// to Eulerian coordinates, given the derivatives of the shape function
128 /// w.r.t the local coordinates.
129 /// Overload the standard version to use the hanging information for
130 /// the Eulerian coordinates.
132 const DShape& dpsids, DenseMatrix<double>& jacobian) const;
133
134 /// Assemble the the "jacobian" matrix of second derivatives of the
135 /// mapping from local to Eulerian coordinates, given
136 /// the second derivatives of the shape functions w.r.t. local coordinates.
137 /// Overload the standard version to use the hanging information for
138 /// the Eulerian coordinates.
140 const DShape& d2psids, DenseMatrix<double>& jacobian2) const;
141
142 /// Assemble the covariant Eulerian base vectors, assuming that
143 /// the derivatives of the shape functions with respect to the local
144 /// coordinates have already been constructed.
145 /// Overload the standard version to account for hanging nodes.
147 const DShape& dpsids, DenseMatrix<double>& interpolated_G) const;
148
149 /// Calculate the mapping from local to Eulerian coordinates given
150 /// the derivatives of the shape functions w.r.t the local coordinates.
151 /// assuming that the coordinates are aligned in the direction of the local
152 /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
153 /// This funciton returns the determinant of the jacobian, the jacobian
154 /// and the inverse jacobian. Overload the standard version to take
155 /// hanging info into account.
157 const DShape& dpsids,
158 DenseMatrix<double>& jacobian,
159 DenseMatrix<double>& inverse_jacobian) const;
160
161 private:
162 /// Storage for local equation numbers of hanging node variables
163 /// (values stored at master nodes). It is
164 /// essential that these are indexed by a Node pointer because the Node
165 /// may be internal or external to the element.
166 /// local equation number = Local_hang_eqn(master_node_pt,ival)
167 std::map<Node*, int>* Local_hang_eqn;
168
169 /// Lookup scheme for unique number associated with any of the nodes
170 /// that actively control the shape of the element (i.e. they are either
171 /// non-hanging nodes of this element or master nodes of hanging nodes.
172 std::map<Node*, unsigned> Shape_controlling_node_lookup;
173
174 protected:
175 /// Assign the local equation numbers for hanging node variables
176 void assign_hanging_local_eqn_numbers(const bool& store_local_dof_pt);
177
178 /// Calculate the contributions to the jacobian from the nodal
179 /// degrees of freedom using finite differences.
180 /// This version is overloaded to take hanging node information into
181 /// account
183 Vector<double>& residuals, DenseMatrix<double>& jacobian);
184
185 public:
186 /// Constructor, calls the FiniteElement constructor and initialises
187 /// the member data
189 : FiniteElement(),
190 Tree_pt(0),
191 Refine_level(0),
192 To_be_refined(false),
195 Number(-1),
197 {
198 }
199
200 /// Destructor, delete the allocated storage for the hanging equations
201 // (The body is now in the cc file to keep the xlC compiler happy under AIX)
202 virtual ~RefineableElement();
203
204 /// Broken copy constructor
206
207 /// Broken assignment operator
208 void operator=(const RefineableElement&) = delete;
209
210 /// Access function: Pointer to quadtree representation of this element
212 {
213 return Tree_pt;
214 }
215
216 /// Set pointer to quadtree representation of this element
217 void set_tree_pt(Tree* my_tree_pt)
218 {
219 Tree_pt = my_tree_pt;
220 }
221
222 /// Set the number of sons that can be constructed by the element
223 /// The default is none
224 virtual unsigned required_nsons() const
225 {
226 return 0;
227 }
228
229
230 /// Flag to indicate suppression of any refinement
232 {
234 }
235
236 /// Suppress of any refinement for this element
238 {
239 Refinement_is_enabled = false;
240 }
241
242 /// Emnable refinement for this element
244 {
246 }
247
248 /// Split the element into the number of sons to be
249 /// constructed and return a
250 /// vector of pointers to the sons. Elements are allocated, but they are
251 /// not given any properties. The refinement level of the sons is one
252 /// higher than that of the father elemern.
253 template<class ELEMENT>
254 void split(Vector<ELEMENT*>& son_pt) const
255 {
256 // Increase refinement level
257 int son_refine_level = Refine_level + 1;
258
259 // How many sons are to be constructed
260 unsigned n_sons = required_nsons();
261 // Resize the son pointer
262 son_pt.resize(n_sons);
263
264 // Loop over the sons and construct
265 for (unsigned i = 0; i < n_sons; i++)
266 {
267 son_pt[i] = new ELEMENT;
268 // Set the refinement level of the newly constructed son.
269 son_pt[i]->set_refinement_level(son_refine_level);
270 }
271 }
272
273
274 /// Access function that returns the local equation number for the
275 /// hanging node variables (values stored at master nodes). The local
276 /// equation number corresponds to the i-th unknown stored at the node
277 /// addressed by node_pt
278 inline int local_hang_eqn(Node* const& node_pt, const unsigned& i)
279 {
280#ifdef RANGE_CHECKING
282 {
283 std::ostringstream error_message;
284 error_message << "Range Error: Value " << i
285 << " is not in the range (0,"
286 << ncont_interpolated_values() - 1 << ")";
287 throw OomphLibError(error_message.str(),
288 OOMPH_CURRENT_FUNCTION,
289 OOMPH_EXCEPTION_LOCATION);
290 }
291#endif
292
293 return Local_hang_eqn[i][node_pt];
294 }
295
296 /// Interface to function that builds the element: i.e. construct
297 /// the nodes, assign their positions, apply boundary conditions, etc. The
298 /// required procedures depend on the geometrical type of the element and
299 /// must be implemented in specific refineable elements. Any new nodes
300 /// created during the build process are returned in the vector
301 /// new_node_pt.
302 virtual void build(Mesh*& mesh_pt,
303 Vector<Node*>& new_node_pt,
304 bool& was_already_built,
305 std::ofstream& new_nodes_file) = 0;
306
307 /// Set the refinement level
308 void set_refinement_level(const int& refine_level)
309 {
310 Refine_level = refine_level;
311 }
312
313 /// Return the Refinement level
314 unsigned refinement_level() const
315 {
316 return Refine_level;
317 }
318
319 /// Select the element for refinement
321 {
322 To_be_refined = true;
323 }
324
325 /// Deselect the element for refinement
327 {
328 To_be_refined = false;
329 }
330
331 /// Unrefinement will be performed by merging the four sons of this element
333 {
335 }
336
337 /// No unrefinement will be performed by merging the four sons of this
338 /// element
340 {
341 Sons_to_be_unrefined = false;
342 }
343
344 /// Has the element been selected for refinement?
346 {
347 return To_be_refined;
348 }
349
350 /// Has the element been selected for unrefinement?
352 {
354 }
355
356 /// Rebuild the element, e.g. set internal values in line with
357 /// those of the sons that have now merged
358 virtual void rebuild_from_sons(Mesh*& mesh_pt) = 0;
359
360 /// Unbuild the element, i.e. mark the nodes that were created
361 /// during its creation for possible deletion
362 virtual void unbuild()
363 {
364 // Get pointer to father element
366 // If there is no father, nothing to do
367 if (father_pt == 0)
368 {
369 return;
370 }
371
372 // Loop over all the nodes
373 unsigned n_node = this->nnode();
374 for (unsigned n = 0; n < n_node; n++)
375 {
376 // If any node in this element is in the father, it can't be deleted
377 if (father_pt->get_node_number(this->node_pt(n)) >= 0)
378 {
380 }
381 }
382 }
383
384 /// Final operations that must be performed when the element is no
385 /// longer active in the mesh, but still resident in the QuadTree.
386 virtual void deactivate_element();
387
388 /// Return true if all the nodes have been built, false if not
389 virtual bool nodes_built()
390 {
391 return node_pt(0) != 0;
392 }
393
394 /// Element number (for debugging/plotting)
395 long number() const
396 {
397 return Number;
398 }
399
400 /// Set element number (for debugging/plotting)
401 void set_number(const long& mynumber)
402 {
403 Number = mynumber;
404 }
405
406 /// Number of continuously interpolated values. Note: We assume
407 /// that they are located at the beginning of the value_pt Vector!
408 /// (Used for interpolation to son elements, for integrity check
409 /// and post-processing -- we can only expect
410 /// the continously interpolated values to be continous across
411 /// element boundaries).
412 virtual unsigned ncont_interpolated_values() const = 0;
413
414 /// Get all continously interpolated function values in this
415 /// element as a Vector. Note: Vector sets is own size to ensure that
416 /// that this function can be used in black-box fashion.
418 Vector<double>& values)
419 {
420 get_interpolated_values(0, s, values);
421 }
422
423 /// Get all continously interpolated function values at previous
424 /// timestep in this element as a Vector. (t=0: present; t>0:
425 /// prev. timestep) Note: Vector sets is own size to ensure that that this
426 /// function can be used in black-box fashion.
427 virtual void get_interpolated_values(const unsigned& t,
428 const Vector<double>& s,
429 Vector<double>& values) = 0;
430
431 /// In mixed elements, different sets of nodes are used to
432 /// interpolate different unknowns. This function returns the n-th node that
433 /// interpolates the value_id-th unknown. Default implementation is that all
434 /// variables use the positional nodes, i.e. isoparametric elements. Note
435 /// that any overloaded versions of this function MUST provide a set
436 /// of nodes for the position, which always has the value_id -1.
437 virtual Node* interpolating_node_pt(const unsigned& n, const int& value_id)
438
439 {
440 return node_pt(n);
441 }
442
443 /// Return the local one dimensional fraction of the n1d-th node
444 /// in the direction of the local coordinate s[i] that is used to
445 /// interpolate the value_id-th continuously interpolated unknown. Default
446 /// assumes isoparametric interpolation for all unknowns
448 const unsigned& n1d, const unsigned& i, const int& value_id)
449 {
450 return local_one_d_fraction_of_node(n1d, i);
451 }
452
453 /// Return a pointer to the node that interpolates the value-id-th
454 /// unknown at local coordinate s. If there is not a node at that position,
455 /// then return 0.
457 const Vector<double>& s, const int& value_id)
458
459 {
461 }
462
463
464 /// Return the number of nodes that are used to interpolate the
465 /// value_id-th unknown. Default is to assume isoparametric elements.
466 virtual unsigned ninterpolating_node(const int& value_id)
467 {
468 return nnode();
469 }
470
471 /// Return the number of nodes in a one_d direction that are
472 /// used to interpolate the value_id-th unknown. Default is to assume
473 /// an isoparametric mapping.
474 virtual unsigned ninterpolating_node_1d(const int& value_id)
475 {
476 return nnode_1d();
477 }
478
479 /// Return the basis functions that are used to interpolate
480 /// the value_id-th unknown. By default assume isoparameteric interpolation
482 Shape& psi,
483 const int& value_id) const
484 {
485 shape(s, psi);
486 }
487
488 /// Check the integrity of the element: Continuity of positions
489 /// values, etc. Essentially, check that the approximation of the functions
490 /// is consistent when viewed from both sides of the element boundaries
491 /// Must be overloaded for each different geometric element
492 virtual void check_integrity(double& max_error) = 0;
493
494 /// Max. allowed discrepancy in element integrity check
495 static double& max_integrity_tolerance()
496 {
498 }
499
500 /// The purpose of this function is to identify all possible
501 /// Data that can affect the fields interpolated by the FiniteElement.
502 /// This must be overloaded to include data from any hanging nodes
503 /// correctly
505 std::set<std::pair<Data*, unsigned>>& paired_field_data);
506
507
508 /// Overload the function that assigns local equation numbers
509 /// for the Data stored at the nodes so that hanging data is taken
510 /// into account
511 inline void assign_nodal_local_eqn_numbers(const bool& store_local_dof_pt)
512 {
514 assign_hanging_local_eqn_numbers(store_local_dof_pt);
515 }
516
517 /// Pointer to the root element in refinement hierarchy (must be
518 /// implemented in specific elements that do refinement via
519 /// tree-like refinement structure. Here we provide a default
520 /// implementation that is appropriate for cases where tree-like
521 /// refinement doesn't exist or if the element doesn't have
522 /// root in that tree (i.e. if it's a root itself): We return
523 /// "this".
525 {
526 // If there is no tree -- the element is its own root
527 if (Tree_pt == 0)
528 {
529 return this;
530 }
531 // Otherwise it's the tree's root object
532 else
533 {
534 return Tree_pt->root_pt()->object_pt();
535 }
536 }
537
538 /// Return a pointer to the father element.
540 {
541 // If we have no tree, we have no father
542 if (Tree_pt == 0)
543 {
544 return 0;
545 }
546 else
547 {
548 // Otherwise get the father of the tree
549 Tree* father_pt = Tree_pt->father_pt();
550 // If the tree has no father then return null, no father
551 if (father_pt == 0)
552 {
553 return 0;
554 }
555 else
556 {
557 return father_pt->object_pt();
558 }
559 }
560 }
561
562 /// Return a pointer to the "father" element at the specified refinement
563 /// level
565 unsigned& refinement_level, RefineableElement*& father_at_reflevel_pt)
566 {
567 // Get the father in the tree (it shouldn't try to get a null Tree...)
568 Tree* father_pt = Tree_pt->father_pt();
569 // Get the refineable element associated with this father
570 RefineableElement* father_el_pt =
571 dynamic_cast<RefineableElement*>(father_pt->object_pt());
572 // Get the refinement level
573 unsigned level = father_el_pt->refinement_level();
574 // If the level matches the required one then return, if not call again
575 if (level == refinement_level)
576 {
577 father_at_reflevel_pt = father_el_pt;
578 }
579 else
580 {
581 // Recursive call
583 father_at_reflevel_pt);
584 }
585 }
586
587 /// Initial setup of the element: e.g. set the appropriate internal
588 /// p-order. If an adopted father is specified, information from this is
589 /// used instead of using the father found from the tree.
590 virtual void initial_setup(Tree* const& adopted_father_pt = 0,
591 const unsigned& initial_p_order = 0)
592 {
593 }
594
595 /// Pre-build the element
596 virtual void pre_build(Mesh*& mesh_pt, Vector<Node*>& new_node_pt) {}
597
598 /// Further build: e.g. deal with interpolation of internal values
599 virtual void further_build() {}
600
601 /// Mark up any hanging nodes that arise as a result of non-uniform
602 /// refinement. Any hanging nodes will be documented in files addressed by
603 /// the streams in the vector output_stream, if the streams are open.
604 virtual void setup_hanging_nodes(Vector<std::ofstream*>& output_stream) {}
605
606 /// Perform additional hanging node procedures for variables
607 /// that are not interpolated by all nodes (e.g. lower order interpolations
608 /// for the pressure in Taylor Hood).
610
611 /// Compute derivatives of elemental residual vector with respect
612 /// to nodal coordinates. Default implementation by FD can be overwritten
613 /// for specific elements.
614 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
615 /// This version is overloaded from the version in FiniteElement
616 /// and takes hanging nodes into account -- j in the above loop
617 /// loops over all the nodes that actively control the
618 /// shape of the element (i.e. they are non-hanging or master nodes of
619 /// hanging nodes in this element).
621 RankThreeTensor<double>& dresidual_dnodal_coordinates);
622
623
624 /// Number of shape-controlling nodes = the number
625 /// of non-hanging nodes plus the number of master nodes associated
626 /// with hanging nodes.
628 {
629 return Shape_controlling_node_lookup.size();
630 }
631
632 /// Return lookup scheme for unique number associated
633 /// with any of the nodes that actively control the shape of the
634 /// element (i.e. they are either non-hanging nodes of this element
635 /// or master nodes of hanging nodes.
636 std::map<Node*, unsigned> shape_controlling_node_lookup()
637 {
639 }
640 };
641
642
643 /// ///////////////////////////////////////////////////////////////////
644 /// ///////////////////////////////////////////////////////////////////
645 /// ///////////////////////////////////////////////////////////////////
646
647
648 //======================================================================
649 /// p-refineable version of RefineableElement
650 //======================================================================
652 {
653 protected:
654 /// The polynomial expansion order of the elemental basis functions
655 unsigned P_order;
656
657 /// Flag for p-refinement
659
660 /// Flag to indicate suppression of any refinement
662
663 /// Flag for unrefinement
665
666 public:
667 /// Constructor, calls the RefineableElement constructor
670 P_order(2),
671 To_be_p_refined(false),
673 To_be_p_unrefined(false)
674 {
675 }
676
677 /// Destructor, empty
679
680 /// Broken copy constructor
682
683 /// Broken assignment operator
684 void operator=(const PRefineableElement&) = delete;
685
686 /// Flag to indicate suppression of any refinement
688 {
690 }
691
692 /// Suppress of any refinement for this element
694 {
696 }
697
698 /// Emnable refinement for this element
700 {
702 }
703
704 /// Access function to P_order
705 unsigned& p_order()
706 {
707 return P_order;
708 }
709
710 /// Access function to P_order (const version)
711 unsigned p_order() const
712 {
713 return P_order;
714 }
715
716 /// Get the initial P_order
717 /// This is required so that elements which are constructed with a
718 /// higher p-order initially are not un-refined past this level (e.g.
719 /// in fluid problems where elements initially use quadratic velocity
720 /// and linear pressure)
721 /// Virtual because this needs to be set in templated derived classes
722 virtual unsigned initial_p_order() const = 0;
723
724 /// Select the element for p-refinement
726 {
727 To_be_p_refined = true;
728 }
729
730 /// Deselect the element for p-refinement
732 {
733 To_be_p_refined = false;
734 }
735
736 /// Select the element for p-unrefinement
738 {
739 To_be_p_unrefined = true;
740 }
741
742 /// Deselect the element for p-unrefinement
744 {
745 To_be_p_unrefined = false;
746 }
747
748 /// Has the element been selected for refinement?
750 {
751 return To_be_p_refined;
752 }
753
754 /// Has the element been selected for p-unrefinement?
756 {
757 return To_be_p_unrefined;
758 }
759
760 /// p-refine the element
761 virtual void p_refine(const int& inc,
762 Mesh* const& mesh_pt,
763 GeneralisedElement* const& clone_pt) = 0;
764
765 // Overload the nodes_built function to check every node
767 {
768 // Must check that EVERY node exists
769 unsigned n_node = this->nnode();
770 for (unsigned n = 0; n < n_node; n++)
771 {
772 if (this->node_pt(n) == 0)
773 {
774 return false;
775 }
776 }
777 // If we get here then all the nodes are built
778 return true;
779 }
780 };
781
782
783 /// ///////////////////////////////////////////////////////////////////
784 /// ///////////////////////////////////////////////////////////////////
785 /// ///////////////////////////////////////////////////////////////////
786
787
788 //=======================================================================
789 /// A base class for elements that can have hanging nodes
790 /// but are not refineable as such. This class is usually used as a
791 /// base class for FaceElements that are attached to refineable
792 /// bulk elements (and stripped out before adapting the bulk
793 /// mesh, so they don't participate in the refimenent process
794 /// itself). We therefore simply break the pure virtual functions
795 /// that don't make any sense for such elements
796 //======================================================================
798 {
799 public:
800 /// Broken build function -- shouldn't really be needed
801 void build(Mesh*& mesh_pt,
802 Vector<Node*>& new_node_pt,
803 bool& was_already_built,
804 std::ofstream& new_nodes_file)
805 {
806 std::ostringstream error_message;
807 error_message << "This function is broken as it's only needed/used \n"
808 << "during \"proper\" refinement\n";
809 throw OomphLibError(
810 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
811 }
812
813
814 /// Broken function -- this shouldn't really be needed.
816 Vector<double>& values)
817 {
818 std::ostringstream error_message;
819 error_message << "This function is broken as it's only needed/used \n"
820 << "during \"proper\" refinement\n";
821 throw OomphLibError(
822 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
823 }
824
825 /// Broken function -- this shouldn't really be needed.
826 virtual void get_interpolated_values(const unsigned& t,
827 const Vector<double>& s,
828 Vector<double>& values)
829 {
830 std::ostringstream error_message;
831 error_message << "This function is broken as it's only needed/used \n"
832 << "during \"proper\" refinement\n";
833 throw OomphLibError(
834 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
835 }
836
837
838 /// Broken function -- this shouldn't really be needed.
839 void check_integrity(double& max_error)
840 {
841 std::ostringstream error_message;
842 error_message << "This function is broken as it's only needed/used \n"
843 << "during \"proper\" refinement\n";
844 throw OomphLibError(
845 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
846 }
847
848 /// Broken function -- this shouldn't really be needed.
849 void rebuild_from_sons(Mesh*& mesh_pt)
850 {
851 std::ostringstream error_message;
852 error_message << "This function is broken as it's only needed/used \n"
853 << "during \"proper\" refinement\n";
854 throw OomphLibError(
855 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
856 }
857 };
858
859
860 /// ///////////////////////////////////////////////////////////////////
861 /// ///////////////////////////////////////////////////////////////////
862 /// ///////////////////////////////////////////////////////////////////
863
864
865 //=======================================================================
866 /// RefineableSolidElements are SolidFiniteElements that may
867 /// be subdivided into children to provide a better local approximation
868 /// to the solution. The distinction is required to keep a clean
869 /// separation between problems that alter nodal positions and others.
870 /// A number of procedures are generic and are included in this class.
871 //======================================================================
873 public virtual SolidFiniteElement
874 {
875 private:
876 /// Storage for local equation numbers of
877 /// hanging node variables associated with nodal positions.
878 /// local position equation number =
879 /// Local_position_hang_eqn(master_node_pt,ival)
880 std::map<Node*, DenseMatrix<int>> Local_position_hang_eqn;
881
882
883 /// Assign local equation numbers to the hanging values associated
884 /// with positions or additional solid values.
885 void assign_solid_hanging_local_eqn_numbers(const bool& store_local_dof_pt);
886
887
888 protected:
889 /// Flag deciding if the Lagrangian coordinates of newly-created
890 /// interior SolidNodes are to be determined by the father element's
891 /// undeformed MacroElement representation (if it has one). Default: False
892 /// as it means that, following a refinement an element is no longer in
893 /// equilbrium (as the Lagrangian coordinate is determined differently from
894 /// the Eulerian one). However, basing the Lagrangian coordinates on the
895 /// undeformed MacroElement can be more stable numerically and for steady
896 /// problems it's not a big deal either way as the difference between the
897 /// two formulations only matters at finite resolution so we have no right
898 /// to say that one is "more correct" than the other...
900
901 /// Assemble the jacobian matrix for the mapping from local
902 /// to lagrangian coordinates, given the derivatives of the shape function
903 /// Overload the standard version to use the hanging information for
904 /// the lagrangian coordinates.
906 const DShape& dpsids, DenseMatrix<double>& jacobian) const;
907
908 /// Assemble the the "jacobian" matrix of second derivatives, given
909 /// the second derivatives of the shape functions w.r.t. local coordinates
910 /// Overload the standard version to use the hanging information for
911 /// the lagrangian coordinates.
913 const DShape& d2psids, DenseMatrix<double>& jacobian2) const;
914
915 /// Calculate the mapping from local to Lagrangian coordinates given
916 /// the derivatives of the shape functions w.r.t the local coorindates.
917 /// assuming that the coordinates are aligned in the direction of the local
918 /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
919 /// This function returns the determinant of the jacobian, the jacobian
920 /// and the inverse jacobian.
922 const DShape& dpsids,
923 DenseMatrix<double>& jacobian,
924 DenseMatrix<double>& inverse_jacobian) const;
925
926 public:
927 /// Constructor
932 {
933 }
934
935 /// Virtual Destructor, delete any allocated storage
937
938 /// Overload the local equation numbers for Data stored as part
939 /// of solid nodes to include hanging node data
940 void assign_solid_local_eqn_numbers(const bool& store_local_dof_pt)
941 {
943 assign_solid_hanging_local_eqn_numbers(store_local_dof_pt);
944 }
945
946 /// The number of geometric data affecting a
947 /// RefineableSolidFiniteElement is the positional Data of all
948 /// non-hanging nodes plus the geometric Data of all distinct
949 /// master nodes. Recomputed on the fly.
950 unsigned ngeom_data() const;
951
952 /// Return pointer to the j-th Data item that the object's
953 /// shape depends on: Positional data of non-hanging nodes and
954 /// positional data of master nodes. Recomputed on the fly.
955 Data* geom_data_pt(const unsigned& j);
956
957 /// Specify Data that affects the geometry of the element
958 /// by adding the position Data to the set that's passed in.
959 /// (This functionality is required in FSI problems; set is used to
960 /// avoid double counting). Refineable version includes hanging nodes
961 void identify_geometric_data(std::set<Data*>& geometric_data_pt);
962
963 /// Compute element residual Vector and element Jacobian matrix
964 /// corresponding to the solid positions. Overloaded version to take
965 /// the hanging nodes into account
967 Vector<double>& residuals, DenseMatrix<double>& jacobian);
968
969 /// Return the flag deciding if the Lagrangian coordinates of
970 /// newly-created interior SolidNodes are to be determined by the father
971 /// element's undeformed MacroElement representation (if it has one).
972 /// Default: False as it means that, following a refinement an element
973 /// is no longer in equilbrium (as the Lagrangian coordinate is
974 /// determined differently from the Eulerian one). However, basing
975 /// the Lagrangian coordinates on the undeformed MacroElement can be
976 /// more stable numerically and for steady problems it's not a big deal
977 /// either way as the difference between the two formulations only matters
978 /// at finite resolution so we have no right to say that one is "more
979 /// correct" than the other...
981 {
983 }
984
985 /// Set the flag deciding if the Lagrangian coordinates of
986 /// newly-created interior SolidNodes are to be determined by the father
987 /// element's undeformed MacroElement representation (if it has one).
989 {
991 }
992
993 /// Unset the flag deciding if the Lagrangian coordinates of
994 /// newly-created interior SolidNodes are to be determined by the father
995 /// element's undeformed MacroElement representation (if it has one).
997 {
999 }
1000
1001 /// Access the local equation number of of hanging node variables
1002 /// associated with nodal positions. The function returns a dense
1003 /// matrix that contains all the local equation numbers corresponding to
1004 /// the positional degrees of freedom.
1006 {
1008 }
1009
1010 /// Further build: Pass the father's
1011 /// Use_undeformed_macro_element_for_new_lagrangian_coords
1012 /// flag down, then call the underlying RefineableElement's
1013 /// version.
1014 virtual void further_build()
1015 {
1019
1021 }
1022 };
1023
1024
1025 /// ///////////////////////////////////////////////////////////////////
1026 /// ///////////////////////////////////////////////////////////////////
1027 /// ///////////////////////////////////////////////////////////////////
1028
1029
1030 //=======================================================================
1031 /// A base class for SolidElements that can have hanging nodes
1032 /// but are not refineable as such. This class is usually used as a
1033 /// base class for FaceElements that are attached to refineable
1034 /// bulk elements (and stripped out before adapting the bulk
1035 /// mesh, so they don't participate in the refimenent process
1036 /// itself). We therefore simply break the pure virtual functions
1037 /// that don't make any sense for such elements
1038 //======================================================================
1040 : public virtual NonRefineableElementWithHangingNodes,
1041 public virtual RefineableSolidElement
1042 {
1043 };
1044
1045
1046 /// ///////////////////////////////////////////////////////////////////
1047 /// ///////////////////////////////////////////////////////////////////
1048 /// ///////////////////////////////////////////////////////////////////
1049
1050
1051} // namespace oomph
1052
1053#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
A general Finite Element class.
Definition: elements.h:1313
virtual Node * get_node_at_local_coordinate(const Vector< double > &s) const
If there is a node at this local coordinate, return the pointer to the node.
Definition: elements.cc:3882
virtual void assign_nodal_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers for Data stored at the nodes Virtual so that it can be overloaded b...
Definition: elements.cc:3544
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2218
virtual double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
Get the local fraction of any node in the n-th position in a one dimensional expansion along the i-th...
Definition: elements.h:1858
int get_node_number(Node *const &node_pt) const
Return the number of the node *node_pt if this node is in the element, else return -1;.
Definition: elements.cc:3814
A Generalised Element class.
Definition: elements.h:73
A general mesh class.
Definition: mesh.h:67
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
void set_non_obsolete()
Mark node as non-obsolete.
Definition: nodes.h:1442
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Broken function – this shouldn't really be needed.
void check_integrity(double &max_error)
Broken function – this shouldn't really be needed.
void rebuild_from_sons(Mesh *&mesh_pt)
Broken function – this shouldn't really be needed.
virtual void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Broken function – this shouldn't really be needed.
void build(Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)
Broken build function – shouldn't really be needed.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
bool to_be_p_unrefined()
Has the element been selected for p-unrefinement?
bool p_refinement_is_enabled()
Flag to indicate suppression of any refinement.
void enable_p_refinement()
Emnable refinement for this element.
bool nodes_built()
Return true if all the nodes have been built, false if not.
void operator=(const PRefineableElement &)=delete
Broken assignment operator.
PRefineableElement(const PRefineableElement &)=delete
Broken copy constructor.
virtual ~PRefineableElement()
Destructor, empty.
virtual unsigned initial_p_order() const =0
Get the initial P_order This is required so that elements which are constructed with a higher p-order...
virtual void p_refine(const int &inc, Mesh *const &mesh_pt, GeneralisedElement *const &clone_pt)=0
p-refine the element
void disable_p_refinement()
Suppress of any refinement for this element.
bool To_be_p_refined
Flag for p-refinement.
unsigned P_order
The polynomial expansion order of the elemental basis functions.
unsigned & p_order()
Access function to P_order.
void select_for_p_refinement()
Select the element for p-refinement.
bool To_be_p_unrefined
Flag for unrefinement.
void deselect_for_p_unrefinement()
Deselect the element for p-unrefinement.
bool to_be_p_refined()
Has the element been selected for refinement?
void select_for_p_unrefinement()
Select the element for p-unrefinement.
bool P_refinement_is_enabled
Flag to indicate suppression of any refinement.
unsigned p_order() const
Access function to P_order (const version)
void deselect_for_p_refinement()
Deselect the element for p-refinement.
PRefineableElement()
Constructor, calls the RefineableElement constructor.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
void set_number(const long &mynumber)
Set element number (for debugging/plotting)
void select_sons_for_unrefinement()
Unrefinement will be performed by merging the four sons of this element.
void assemble_local_to_eulerian_jacobian2(const DShape &d2psids, DenseMatrix< double > &jacobian2) const
Assemble the the "jacobian" matrix of second derivatives of the mapping from local to Eulerian coordi...
virtual void build(Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)=0
Interface to function that builds the element: i.e. construct the nodes, assign their positions,...
static double Max_integrity_tolerance
Max. allowed discrepancy in element integrity check.
virtual void further_build()
Further build: e.g. deal with interpolation of internal values.
void select_for_refinement()
Select the element for refinement.
virtual void rebuild_from_sons(Mesh *&mesh_pt)=0
Rebuild the element, e.g. set internal values in line with those of the sons that have now merged.
void assign_hanging_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers for hanging node variables.
virtual void pre_build(Mesh *&mesh_pt, Vector< Node * > &new_node_pt)
Pre-build the element.
void set_refinement_level(const int &refine_level)
Set the refinement level.
std::map< Node *, unsigned > shape_controlling_node_lookup()
Return lookup scheme for unique number associated with any of the nodes that actively control the sha...
virtual void unbuild()
Unbuild the element, i.e. mark the nodes that were created during its creation for possible deletion.
Tree * Tree_pt
A pointer to a general tree object.
virtual void deactivate_element()
Final operations that must be performed when the element is no longer active in the mesh,...
long number() const
Element number (for debugging/plotting)
Tree * tree_pt()
Access function: Pointer to quadtree representation of this element.
double local_to_eulerian_mapping_diagonal(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates given the derivatives of the shape functions...
virtual bool nodes_built()
Return true if all the nodes have been built, false if not.
virtual unsigned ncont_interpolated_values() const =0
Number of continuously interpolated values. Note: We assume that they are located at the beginning of...
virtual void setup_hanging_nodes(Vector< std::ofstream * > &output_stream)
Mark up any hanging nodes that arise as a result of non-uniform refinement. Any hanging nodes will be...
static double & max_integrity_tolerance()
Max. allowed discrepancy in element integrity check.
virtual void initial_setup(Tree *const &adopted_father_pt=0, const unsigned &initial_p_order=0)
Initial setup of the element: e.g. set the appropriate internal p-order. If an adopted father is spec...
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
void get_father_at_refinement_level(unsigned &refinement_level, RefineableElement *&father_at_reflevel_pt)
Return a pointer to the "father" element at the specified refinement level.
unsigned refinement_level() const
Return the Refinement level.
static void check_value_id(const int &n_continuously_interpolated_values, const int &value_id)
Static helper function that is used to check that the value_id is in range.
unsigned Refine_level
Refinement level.
virtual Node * interpolating_node_pt(const unsigned &n, const int &value_id)
In mixed elements, different sets of nodes are used to interpolate different unknowns....
void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
virtual void fill_in_jacobian_from_nodal_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the jacobian from the nodal degrees of freedom using finite difference...
virtual unsigned ninterpolating_node_1d(const int &value_id)
Return the number of nodes in a one_d direction that are used to interpolate the value_id-th unknown....
void disable_refinement()
Suppress of any refinement for this element.
virtual void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes (e....
long Number
Global element number – for plotting/validation purposes.
virtual unsigned ninterpolating_node(const int &value_id)
Return the number of nodes that are used to interpolate the value_id-th unknown. Default is to assume...
virtual RefineableElement * root_element_pt()
Pointer to the root element in refinement hierarchy (must be implemented in specific elements that do...
virtual void interpolating_basis(const Vector< double > &s, Shape &psi, const int &value_id) const
Return the basis functions that are used to interpolate the value_id-th unknown. By default assume is...
bool refinement_is_enabled()
Flag to indicate suppression of any refinement.
virtual unsigned required_nsons() const
Set the number of sons that can be constructed by the element The default is none.
virtual Node * get_interpolating_node_at_local_coordinate(const Vector< double > &s, const int &value_id)
Return a pointer to the node that interpolates the value-id-th unknown at local coordinate s....
void deselect_for_refinement()
Deselect the element for refinement.
void assemble_local_to_eulerian_jacobian(const DShape &dpsids, DenseMatrix< double > &jacobian) const
Assemble the jacobian matrix for the mapping from local to Eulerian coordinates, given the derivative...
unsigned nshape_controlling_nodes()
Number of shape-controlling nodes = the number of non-hanging nodes plus the number of master nodes a...
bool To_be_refined
Flag for refinement.
std::map< Node *, int > * Local_hang_eqn
Storage for local equation numbers of hanging node variables (values stored at master nodes)....
RefineableElement()
Constructor, calls the FiniteElement constructor and initialises the member data.
void assign_nodal_local_eqn_numbers(const bool &store_local_dof_pt)
Overload the function that assigns local equation numbers for the Data stored at the nodes so that ha...
void split(Vector< ELEMENT * > &son_pt) const
Split the element into the number of sons to be constructed and return a vector of pointers to the so...
bool to_be_refined()
Has the element been selected for refinement?
void enable_refinement()
Emnable refinement for this element.
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
void assemble_eulerian_base_vectors(const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
Assemble the covariant Eulerian base vectors, assuming that the derivatives of the shape functions wi...
virtual double local_one_d_fraction_of_interpolating_node(const unsigned &n1d, const unsigned &i, const int &value_id)
Return the local one dimensional fraction of the n1d-th node in the direction of the local coordinate...
void identify_field_data_for_interactions(std::set< std::pair< Data *, unsigned > > &paired_field_data)
The purpose of this function is to identify all possible Data that can affect the fields interpolated...
RefineableElement(const RefineableElement &)=delete
Broken copy constructor.
void set_tree_pt(Tree *my_tree_pt)
Set pointer to quadtree representation of this element.
bool Sons_to_be_unrefined
Flag for unrefinement.
virtual void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get all continously interpolated function values in this element as a Vector. Note: Vector sets is ow...
std::map< Node *, unsigned > Shape_controlling_node_lookup
Lookup scheme for unique number associated with any of the nodes that actively control the shape of t...
virtual void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)=0
Get all continously interpolated function values at previous timestep in this element as a Vector....
bool sons_to_be_unrefined()
Has the element been selected for unrefinement?
bool Refinement_is_enabled
Flag to indicate suppression of any refinement.
virtual ~RefineableElement()
Destructor, delete the allocated storage for the hanging equations.
virtual void check_integrity(double &max_error)=0
Check the integrity of the element: Continuity of positions values, etc. Essentially,...
void operator=(const RefineableElement &)=delete
Broken assignment operator.
void deselect_sons_for_unrefinement()
No unrefinement will be performed by merging the four sons of this element.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
virtual void further_build()
Further build: Pass the father's Use_undeformed_macro_element_for_new_lagrangian_coords flag down,...
void disable_use_of_undeformed_macro_element_for_new_lagrangian_coords()
Unset the flag deciding if the Lagrangian coordinates of newly-created interior SolidNodes are to be ...
bool is_undeformed_macro_element_used_for_new_lagrangian_coords() const
Return the flag deciding if the Lagrangian coordinates of newly-created interior SolidNodes are to be...
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object's shape depends on: Positional data of non-hangi...
double local_to_lagrangian_mapping_diagonal(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Lagrangian coordinates given the derivatives of the shape functio...
DenseMatrix< int > & local_position_hang_eqn(Node *const &node_pt)
Access the local equation number of of hanging node variables associated with nodal positions....
void assemble_local_to_lagrangian_jacobian2(const DShape &d2psids, DenseMatrix< double > &jacobian2) const
Assemble the the "jacobian" matrix of second derivatives, given the second derivatives of the shape f...
void enable_use_of_undeformed_macro_element_for_new_lagrangian_coords()
Set the flag deciding if the Lagrangian coordinates of newly-created interior SolidNodes are to be de...
unsigned ngeom_data() const
The number of geometric data affecting a RefineableSolidFiniteElement is the positional Data of all n...
void assign_solid_local_eqn_numbers(const bool &store_local_dof_pt)
Overload the local equation numbers for Data stored as part of solid nodes to include hanging node da...
bool Use_undeformed_macro_element_for_new_lagrangian_coords
Flag deciding if the Lagrangian coordinates of newly-created interior SolidNodes are to be determined...
virtual ~RefineableSolidElement()
Virtual Destructor, delete any allocated storage.
void assign_solid_hanging_local_eqn_numbers(const bool &store_local_dof_pt)
Assign local equation numbers to the hanging values associated with positions or additional solid val...
void fill_in_jacobian_from_solid_position_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute element residual Vector and element Jacobian matrix corresponding to the solid positions....
std::map< Node *, DenseMatrix< int > > Local_position_hang_eqn
Storage for local equation numbers of hanging node variables associated with nodal positions....
void identify_geometric_data(std::set< Data * > &geometric_data_pt)
Specify Data that affects the geometry of the element by adding the position Data to the set that's p...
void assemble_local_to_lagrangian_jacobian(const DShape &dpsids, DenseMatrix< double > &jacobian) const
Assemble the jacobian matrix for the mapping from local to lagrangian coordinates,...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: elements.h:3561
virtual void assign_solid_local_eqn_numbers(const bool &store_local_dof)
Assigns local equation numbers for the generic solid local equation numbering schemes....
Definition: elements.cc:6898
A generalised tree base class that abstracts the common functionality between the quad- and octrees u...
Definition: tree.h:74
TreeRoot *& root_pt()
Return pointer to root of the tree.
Definition: tree.h:141
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Definition: tree.h:88
Tree * father_pt() const
Return pointer to father: NULL if it's a root node.
Definition: tree.h:235
//////////////////////////////////////////////////////////////////// ////////////////////////////////...