refineable_quad_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// Version 1.0; svn revision $LastChangedRevision: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #ifndef OOMPH_REFINEABLE_QUAD_SPECTRAL_ELEMENT_HEADER
31 #define OOMPH_REFINEABLE_QUAD_SPECTRAL_ELEMENT_HEADER
32 
33 
34 // Config header generated by autoconfig
35 #ifdef HAVE_CONFIG_H
36  #include <oomph-lib-config.h>
37 #endif
38 
39 
40 //oomph-lib headers
42 
43 namespace oomph
44 {
45 
46 //=======================================================================
47 /// Refineable version of QuadElements that add functionality for spectral
48 /// Elements.
49 //=======================================================================
50 template<>
51 class RefineableQSpectralElement<2> : public virtual RefineableQElement<2>
52 {
53 
54 public:
55 
56  /// Constructor
58  {
59 #ifdef LEAK_CHECK
60  LeakCheckNames::RefineableQSpectralElement<2>_build+=1;
61 #endif
62  }
63 
64 
65  /// Broken copy constructor
67  {
68  BrokenCopy::broken_copy("RefineableQSpectralElement<2>");
69  }
70 
71  /// Broken assignment operator
72 //Commented out broken assignment operator because this can lead to a conflict warning
73 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
74 //realise that two separate implementations of the broken function are the same and so,
75 //quite rightly, it shouts.
76  /*void operator=(const RefineableQSpectralElement<2>&)
77  {
78  BrokenCopy::broken_assign("RefineableQSpecralElement<2>");
79  }*/
80 
81  /// Destructor
83  {
84 #ifdef LEAK_CHECK
85  LeakCheckNames::RefineableQSpectralElement<2>_build-=1;
86 #endif
87  }
88 
89  /// The only thing to add is rebuild from sons
90  void rebuild_from_sons(Mesh* &mesh_pt)
91  {
92  //The timestepper should be the same for all nodes and node 0 should
93  //never be deleted.
94  if(this->node_pt(0)==0)
95  {
96  throw OomphLibError("The Corner node (0) does not exist",
97  OOMPH_CURRENT_FUNCTION,
98  OOMPH_EXCEPTION_LOCATION);
99  }
100 
101  TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
102  unsigned ntstorage = time_stepper_pt->ntstorage();
103 
104  unsigned jnod=0;
105  Vector<double> s_fraction(2), s(2);
106  //Loop over the nodes in the element
107  unsigned n_p = this->nnode_1d();
108  for(unsigned i0=0;i0<n_p;i0++)
109  {
110  //Get the fractional position of the node
111  s_fraction[0] = this->local_one_d_fraction_of_node(i0,0);
112  //Local coordinate
113  s[0] = -1.0 + 2.0*s_fraction[0];
114 
115  for(unsigned i1=0;i1<n_p;i1++)
116  {
117  //Get the fractional position of the node in the direction of s[1]
118  s_fraction[1] = this->local_one_d_fraction_of_node(i1,1);
119  // Local coordinate in father element
120  s[1] = -1.0 + 2.0*s_fraction[1];
121 
122  //Set the local node number
123  jnod = i0 + n_p*i1;
124 
125  //If the node has not been built
126  if(this->node_pt(jnod)==0)
127  {
128  //Has the node been created by one of its neighbours
129  bool is_periodic = false;
130  Node* created_node_pt =
131  this->node_created_by_neighbour(s_fraction,is_periodic);
132 
133  //If it has set the pointer
134  if(created_node_pt!=0)
135  {
136  //If the node is periodic
137  if(is_periodic)
138  {
139  throw OomphLibError(
140  "Cannot handle periodic nodes in refineable spectral elements yet",
141  OOMPH_CURRENT_FUNCTION,
142  OOMPH_EXCEPTION_LOCATION);
143  }
144  //Non-periodic case, just set the pointer
145  else
146  {
147  this->node_pt(jnod) = created_node_pt;
148  }
149  }
150  //Otherwise, we need to build it
151  else
152  {
153  //First, find the son element in which the node should live
154 
155  //Find coordinates in the sons
156  Vector<double> s_in_son(2);
157  using namespace QuadTreeNames;
158  int son=-10;
159  //If negative on the west side
160  if(s_fraction[0] < 0.5)
161  {
162  //On the south side
163  if(s_fraction[1] < 0.5)
164  {
165  //It's the southwest son
166  son = SW;
167  s_in_son[0] = -1.0 + 4.0*s_fraction[0];
168  s_in_son[1] = -1.0 + 4.0*s_fraction[1];
169  }
170  //On the north side
171  else
172  {
173  //It's the northwest son
174  son = NW;
175  s_in_son[0] = -1.0 + 4.0*s_fraction[0];
176  s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
177  }
178  }
179  else
180  {
181  //On the south side
182  if(s_fraction[1] < 0.5)
183  {
184  //It's the southeast son
185  son = SE;
186  s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
187  s_in_son[1] = -1.0 + 4.0*s_fraction[1];
188  }
189  //On the north side
190  else
191  {
192  //It's the northeast son
193  son = NE;
194  s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
195  s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
196  }
197  }
198 
199  //Get the pointer to the son element in which the new node
200  //would live
201  RefineableQSpectralElement<2>* son_el_pt =
202  dynamic_cast<RefineableQSpectralElement<2>*>(
203  this->tree_pt()->son_pt(son)->object_pt());
204 
205  //If we are rebuilding, then worry about the boundary conditions
206  //Find the boundary of the node
207  //Initially none
208  int boundary=Tree::OMEGA;
209  //If we are on the western boundary
210  if(i0==0) {boundary = W;}
211  //If we are on the eastern boundary
212  else if(i0==n_p-1) {boundary = E;}
213 
214  //If we are on the southern boundary
215  if(i1==0)
216  {
217  //If we already have already set the boundary, we're on a corner
218  switch(boundary)
219  {
220  case W:
221  boundary = SW;
222  break;
223  case E:
224  boundary = SE;
225  break;
226  //Boundary not set
227  default:
228  boundary = S;
229  break;
230  }
231  }
232  //If we are the northern bounadry
233  else if(i1==n_p-1)
234  {
235  //If we already have a boundary
236  switch(boundary)
237  {
238  case W:
239  boundary = NW;
240  break;
241  case E:
242  boundary = NE;
243  break;
244  default:
245  boundary = N;
246  break;
247  }
248  }
249 
250  // set of boundaries that this edge in the son lives on
251  std::set<unsigned> boundaries;
252 
253  //Now get the boundary conditions from the son
254  //The boundaries will be common to the son because there can be
255  //no rotations here
256  if(boundary!=Tree::OMEGA)
257  {
258  son_el_pt->get_boundaries(boundary,boundaries);
259  }
260 
261  // If the node lives on a boundary:
262  // Construct a boundary node,
263  // Get boundary conditions and
264  // update all lookup schemes
265  if(boundaries.size()>0)
266  {
267  //Construct the new node
268  this->node_pt(jnod) =
269  this->construct_boundary_node(jnod,time_stepper_pt);
270 
271  //Get the boundary conditions from the son
272  Vector<int> bound_cons(ncont_interpolated_values());
273  son_el_pt->get_bcs(boundary,bound_cons);
274 
275  //Loop over the values and pin if necessary
276  unsigned nval = this->node_pt(jnod)->nvalue();
277  for(unsigned k=0;k<nval;k++)
278  {
279  if(bound_cons[k]) {this->node_pt(jnod)->pin(k);}
280  }
281 
282  // Solid node? If so, deal with the positional boundary
283  // conditions:
284  SolidNode* solid_node_pt =
285  dynamic_cast<SolidNode*>(this->node_pt(jnod));
286  if (solid_node_pt!=0)
287  {
288  //Get the positional boundary conditions from the father:
289  unsigned n_dim = this->node_pt(jnod)->ndim();
290  Vector<int> solid_bound_cons(n_dim);
291  RefineableSolidQElement<2>* son_solid_el_pt=
292  dynamic_cast<RefineableSolidQElement<2>*>(son_el_pt);
293 #ifdef PARANOID
294  if (son_solid_el_pt==0)
295  {
296  std::string error_message =
297  "We have a SolidNode outside a refineable SolidElement\n";
298  error_message +=
299  "during mesh refinement -- this doesn't make sense\n";
300 
301  throw OomphLibError(
302  error_message,
303  OOMPH_CURRENT_FUNCTION,
304  OOMPH_EXCEPTION_LOCATION);
305  }
306 #endif
307  son_solid_el_pt->get_solid_bcs(boundary,solid_bound_cons);
308 
309  //Loop over the positions and pin, if necessary
310  for(unsigned k=0;k<n_dim;k++)
311  {
312  if (solid_bound_cons[k]) {solid_node_pt->pin_position(k);}
313  }
314  }
315 
316  //Next we update the boundary look-up schemes
317  //Loop over the boundaries stored in the set
318  for(std::set<unsigned>::iterator it = boundaries.begin();
319  it != boundaries.end(); ++it)
320  {
321  //Add the node to the boundary
322  mesh_pt->add_boundary_node(*it,this->node_pt(jnod));
323 
324  //If we have set an intrinsic coordinate on this
325  //mesh boundary then it must also be interpolated on
326  //the new node
327  //Now interpolate the intrinsic boundary coordinate
328  if(mesh_pt->boundary_coordinate_exists(*it)==true)
329  {
330  Vector<double> zeta(1);
331  son_el_pt->interpolated_zeta_on_edge(*it,boundary,
332  s_in_son,zeta);
333 
334  this->node_pt(jnod)
335  ->set_coordinates_on_boundary(*it,zeta);
336  }
337  }
338  }
339  //Otherwise the node is not on a Mesh boundary
340  //and we create a normal "bulk" node
341  else
342  {
343  //Construct the new node
344  this->node_pt(jnod) = this->construct_node(jnod,time_stepper_pt);
345  }
346 
347  //Now we set the position and values at the newly created node
348 
349  // In the first instance use macro element or FE representation
350  // to create past and present nodal positions.
351  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
352  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
353  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
354  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
355  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
356  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
357  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
358 
359  // Loop over # of history values
360  //Loop over # of history values
361  for(unsigned t=0;t<ntstorage;t++)
362  {
363  using namespace QuadTreeNames;
364  //Get the position from the son
365  Vector<double> x_prev(2);
366 
367  //Now let's fill in the value
368  son_el_pt->get_x(t,s_in_son,x_prev);
369  for(unsigned i=0;i<2;i++)
370  {
371  this->node_pt(jnod)->x(t,i) = x_prev[i];
372  }
373  }
374 
375  // Now set up the values
376  // Loop over all history values
377  for(unsigned t=0;t<ntstorage;t++)
378  {
379  // Get values from father element
380  // Note: get_interpolated_values() sets Vector size itself.
381  Vector<double> prev_values;
382  son_el_pt->get_interpolated_values(t,s_in_son,prev_values);
383 
384  //Initialise the values at the new node
385  for(unsigned k=0;k<this->node_pt(jnod)->nvalue();k++)
386  {
387  this->node_pt(jnod)->set_value(t,k,prev_values[k]);
388  }
389  }
390 
391  //Add the node to the mesh
392  mesh_pt->add_node_pt(this->node_pt(jnod));
393 
394  } //End of the case when we build the node ourselvesx
395 
396  //Algebraic stuff here
397  //Check whether the element is an algebraic element
398  AlgebraicElementBase* alg_el_pt =
399  dynamic_cast<AlgebraicElementBase*>(this);
400 
401  //If we do have an algebraic element
402  if(alg_el_pt!=0)
403  {
404  std::string error_message =
405  "Have not implemented rebuilding from sons for";
406  error_message +=
407  "Algebraic Spectral elements yet\n";
408 
409  throw
410  OomphLibError(error_message,
411  "RefineableQSpectralElement::rebuild_from_sons()",
412  OOMPH_EXCEPTION_LOCATION);
413  }
414 
415  }
416  }
417  }
418  }
419 
420 
421  /// Overload the nodes built function
422  virtual bool nodes_built()
423  {
424  unsigned n_node = this->nnode();
425  for(unsigned n=0;n<n_node;n++)
426  {
427  if(node_pt(n)==0) {return false;}
428  }
429  //If we get to here, OK
430  return true;
431  }
432 
433 };
434 
435 }
436 
437 #endif
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
virtual ~RefineableQSpectralElement()
Broken assignment operator.
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:246
cstr elem_len * i
Definition: cfortran.h:607
void get_solid_bcs(int bound, Vector< int > &solid_bound_cons) const
Determine vector of solid (positional) boundary conditions along edge (or on vertex) bound (S/W/N/E/S...
char t
Definition: cfortran.h:572
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:602
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1693
static char t char * s
Definition: cfortran.h:572
bool boundary_coordinate_exists(const unsigned &i) const
Indicate whether the i-th boundary has an intrinsic coordinate.
Definition: mesh.h:565
void rebuild_from_sons(Mesh *&mesh_pt)
The only thing to add is rebuild from sons.
virtual bool nodes_built()
Overload the nodes built function.
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:241
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
Definition: nodes.h:1568
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
Refineable version of Solid quad elements.
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
long RefineableQElement< 2 > _build
RefineableQSpectralElement(const RefineableQSpectralElement< 2 > &dummy)
Broken copy constructor.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
A general mesh class.
Definition: mesh.h:74