refineable_line_spectral_element.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #ifndef OOMPH_REFINEABLE_LINE_SPECTRAL_ELEMENT_HEADER
27 #define OOMPH_REFINEABLE_LINE_SPECTRAL_ELEMENT_HEADER
28 
29 
30 // Config header generated by autoconfig
31 #ifdef HAVE_CONFIG_H
32 #include <oomph-lib-config.h>
33 #endif
34 
35 
36 // oomph-lib headers
37 #include "mesh.h"
38 #include "algebraic_elements.h"
40 #include "Qspectral_elements.h"
42 
43 namespace oomph
44 {
45  //==========================================================================
46  /// Refineable version of LineElements that add functionality for spectral
47  /// Elements.
48  //==========================================================================
49  template<>
50  class RefineableQSpectralElement<1> : public virtual RefineableQElement<1>
51  {
52  public:
53  /// Constructor
55  {
56 #ifdef LEAK_CHECK
57  LeakCheckNames::RefineableQSpectralElement<1> _build += 1;
58 #endif
59  }
60 
61  /// Broken copy constructor
63  delete;
64 
65  /// Broken assignment operator
66  // Commented out broken assignment operator because this can lead to a
67  // conflict warning when used in the virtual inheritence hierarchy.
68  // Essentially the compiler doesn't realise that two separate
69  // implementations of the broken function are the same and so, quite
70  // rightly, it shouts.
71  /*void operator=(const RefineableQSpectralElement<1>&) = delete;*/
72 
73  /// Destructor
75  {
76 #ifdef LEAK_CHECK
77  LeakCheckNames::RefineableQSpectralElement<1> _build -= 1;
78 #endif
79  }
80 
81  /// The only thing to add is rebuild from sons
82  void rebuild_from_sons(Mesh*& mesh_pt)
83  {
84  // The timestepper should be the same for all nodes and node 0 should
85  // never be deleted.
86  if (this->node_pt(0) == 0)
87  {
88  throw OomphLibError("The vertex node (0) does not exist",
89  OOMPH_CURRENT_FUNCTION,
90  OOMPH_EXCEPTION_LOCATION);
91  }
92 
93  TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
94 
95  // Determine number of history values stored
96  const unsigned ntstorage = time_stepper_pt->ntstorage();
97 
98  // Allocate storage for local coordinates
99  Vector<double> s_fraction(1), s(1);
100 
101  // Determine the number of nodes in the element
102  const unsigned n_node = this->nnode_1d();
103 
104  // Loop over the nodes in the element
105  for (unsigned n = 0; n < n_node; n++)
106  {
107  // Get the fractional position of the node in the direction of s[0]
108  s_fraction[0] = this->local_one_d_fraction_of_node(n, 0);
109 
110  // Determine the local coordinate in the father element
111  s[0] = -1.0 + 2.0 * s_fraction[0];
112 
113  // If the node has not been built
114  if (this->node_pt(n) == 0)
115  {
116  // Has the node been created by one of its neighbours?
117  bool is_periodic = false;
118  Node* created_node_pt =
119  this->node_created_by_neighbour(s_fraction, is_periodic);
120 
121  // If it has, set the pointer
122  if (created_node_pt != 0)
123  {
124  // If the node is periodic
125  if (is_periodic)
126  {
127  throw OomphLibError(
128  "Cannot handle periodic nodes in refineable spectral elements",
129  OOMPH_CURRENT_FUNCTION,
130  OOMPH_EXCEPTION_LOCATION);
131  }
132  // Non-periodic case, just set the pointer
133  else
134  {
135  this->node_pt(n) = created_node_pt;
136  }
137  }
138  // Otherwise, we need to build it
139  else
140  {
141  // First, find the son element in which the node should live
142 
143  // Find coordinate in the son
144  Vector<double> s_in_son(1);
145  using namespace BinaryTreeNames;
146  int son = -10;
147  // If s_fraction is between 0 and 0.5, we are in the left son
148  if (s_fraction[0] < 0.5)
149  {
150  son = L;
151  s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
152  }
153  // Otherwise we are in the right son
154  else
155  {
156  son = R;
157  s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
158  }
159 
160  // Get the pointer to the son element in which the new node
161  // would live
162  RefineableQSpectralElement<1>* son_el_pt =
163  dynamic_cast<RefineableQSpectralElement<1>*>(
164  this->tree_pt()->son_pt(son)->object_pt());
165 
166  // In 1D we should never be rebuilding an element's vertex nodes
167  // (since they will never be deleted), so throw an error if we
168  // appear to be doing so
169 #ifdef PARANOID
170  if (n == 0 || n == n_node - 1)
171  {
172  std::string error_message =
173  "I am trying to rebuild one of the (two) vertex nodes in\n";
174  error_message +=
175  "this 1D element. It should not have been possible to delete\n";
176  error_message += "either of these!\n";
177 
178  throw OomphLibError(error_message,
179  OOMPH_CURRENT_FUNCTION,
180  OOMPH_EXCEPTION_LOCATION);
181  }
182 #endif
183 
184  // With this in mind we will always be creating normal "bulk" nodes
185  this->node_pt(n) = this->construct_node(n, time_stepper_pt);
186 
187  // Now we set the position and values at the newly created node
188 
189  // In the first instance use macro element or FE representation
190  // to create past and present nodal positions.
191  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC ELEMENTS AS NOT
192  // ALL OF THEM NECESSARILY IMPLEMENT NONTRIVIAL NODE UPDATE
193  // FUNCTIONS. CALLING THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL
194  // LEAVE THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
195  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL NOT ASSIGN SENSIBLE
196  // INITIAL POSITIONS!)
197 
198  // Loop over history values
199  for (unsigned t = 0; t < ntstorage; t++)
200  {
201  // Allocate storage for the previous position of the node
202  Vector<double> x_prev(1);
203 
204  // Get position from son element -- this uses the macro element
205  // representation if appropriate
206  son_el_pt->get_x(t, s_in_son, x_prev);
207 
208  // Set the previous position of the new node
209  this->node_pt(n)->x(t, 0) = x_prev[0];
210 
211  // Allocate storage for the previous values at the node
212  // NOTE: the size of this vector is equal to the number of values
213  // (e.g. 3 velocity components and 1 pressure, say)
214  // associated with each node and NOT the number of history values
215  // which the node stores!
216  Vector<double> prev_values;
217 
218  // Get values from son element
219  // Note: get_interpolated_values() sets Vector size itself.
220  son_el_pt->get_interpolated_values(t, s_in_son, prev_values);
221 
222  // Determine the number of values at the new node
223  const unsigned n_value = this->node_pt(n)->nvalue();
224 
225  // Loop over all values and set the previous values
226  for (unsigned v = 0; v < n_value; v++)
227  {
228  this->node_pt(n)->set_value(t, v, prev_values[v]);
229  }
230  } // End of loop over history values
231 
232  // Add new node to mesh
233  mesh_pt->add_node_pt(this->node_pt(n));
234 
235  } // End of case where we build the node ourselves
236 
237  // Check if the element is an algebraic element
238  AlgebraicElementBase* alg_el_pt =
239  dynamic_cast<AlgebraicElementBase*>(this);
240 
241  // If so, throw error
242  if (alg_el_pt != 0)
243  {
244  std::string error_message =
245  "Have not implemented rebuilding from sons for";
246  error_message += "Algebraic Spectral elements yet\n";
247 
248  throw OomphLibError(
249  error_message,
250  "RefineableQSpectralElement<1>::rebuild_from_sons()",
251  OOMPH_EXCEPTION_LOCATION);
252  }
253 
254  } // End of if this node has not been built
255  } // End of loop over nodes in element
256  }
257 
258 
259  /// Overload the nodes built function
260  virtual bool nodes_built()
261  {
262  const unsigned n_node = this->nnode();
263  for (unsigned n = 0; n < n_node; n++)
264  {
265  if (node_pt(n) == 0)
266  {
267  return false;
268  }
269  }
270  // If we get to here, OK
271  return true;
272  }
273  };
274 
275 } // namespace oomph
276 
277 #endif
static char t char * s
Definition: cfortran.h:568
char t
Definition: cfortran.h:568
////////////////////////////////////////////////////////////////////
void get_x(const Vector< double > &s, Vector< double > &x) const
Global coordinates as function of local coordinates. Either via FE representation or via macro-elemen...
Definition: elements.h:1885
A general mesh class.
Definition: mesh.h:67
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:611
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
An OomphLibError object which should be thrown when an run-time error is encountered....
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
virtual void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get all continously interpolated function values in this element as a Vector. Note: Vector sets is ow...
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
Refineable version of LineElements that add functionality for spectral Elements.
virtual ~RefineableQSpectralElement()
Broken assignment operator.
void rebuild_from_sons(Mesh *&mesh_pt)
The only thing to add is rebuild from sons.
RefineableQSpectralElement(const RefineableQSpectralElement< 1 > &dummy)=delete
Broken copy constructor.
virtual bool nodes_built()
Overload the nodes built function.
A class that is used to template the refineable Q spectral elements by dimension. It's really nothing...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
long RefineableQElement< 2 > _build
//////////////////////////////////////////////////////////////////// ////////////////////////////////...