refineable_brick_spectral_element.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// 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_BRICK_SPECTRAL_ELEMENT_HEADER
31 #define OOMPH_REFINEABLE_BRICK_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<3> : public virtual RefineableQElement<3>
52 {
53 
54 public:
55 
56  /// Constructor
58  {
59 #ifdef LEAK_CHECK
60  LeakCheckNames::RefineableQSpectralElement<3>_build+=1;
61 #endif
62  }
63 
64 
65  /// Broken copy constructor
67  {
68  BrokenCopy::broken_copy("RefineableQSpectralElement<3>");
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<3>&)
77  {
78  BrokenCopy::broken_assign("RefineableQSpecralElement<3>");
79  }*/
80 
81  /// Destructor
83  {
84 #ifdef LEAK_CHECK
85  LeakCheckNames::RefineableQSpectralElement<3>_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(3), s(3);
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  for(unsigned i2=0;i2<n_p;i2++)
123  {
124  //Get the fractional position of the node in the direction of s[1]
125  s_fraction[2] = this->local_one_d_fraction_of_node(i2,2);
126  // Local coordinate in father element
127  s[2] = -1.0 + 2.0*s_fraction[2];
128 
129  //Set the local node number
130  jnod = i0 + n_p*i1 + n_p*n_p*i2;
131 
132  //If the node has not been built
133  if(this->node_pt(jnod)==0)
134  {
135  //Has the node been created by one of its neighbours
136  Node* created_node_pt = this->node_created_by_neighbour(s_fraction);
137 
138  //If it has set the pointer
139  if(created_node_pt!=0)
140  {
141  this->node_pt(jnod) = created_node_pt;
142  }
143  //Otherwise, we need to build it
144  else
145  {
146  //First we need to find the pointer to the son that
147  //would contain the node
148 
149  //Find coordinates in the sons
150  Vector<double> s_in_son(3);
151  using namespace OcTreeNames;
152  int son=-10;
153  //If less than one-half on the left side
154  if(s_fraction[0] < 0.5)
155  {
156  //On the down side
157  if(s_fraction[1] < 0.5)
158  {
159  //On the back side
160  if(s_fraction[2] < 0.5)
161  {
162  //It's the left down back son
163  son = LDB;
164  s_in_son[0] = -1.0 + 4.0*s_fraction[0];
165  s_in_son[1] = -1.0 + 4.0*s_fraction[1];
166  s_in_son[2] = -1.0 + 4.0*s_fraction[2];
167  }
168  //On the front side
169  else
170  {
171  //It's the left down front son
172  son = LDF;
173  s_in_son[0] = -1.0 + 4.0*s_fraction[0];
174  s_in_son[1] = -1.0 + 4.0*s_fraction[1];
175  s_in_son[2] = -1.0 + 4.0*(s_fraction[2]-0.5);
176  }
177  } //End of on down side
178  //else its on the top
179  else
180  {
181  //On the back
182  if(s_fraction[2] < 0.5)
183  {
184  //It's the left up back son
185  son = LUB;
186  s_in_son[0] = -1.0 + 4.0*s_fraction[0];
187  s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
188  s_in_son[2] = -1.0 + 4.0*s_fraction[2];
189  }
190  //On the front side
191  else
192  {
193  //It's the left up front son
194  son = LUF;
195  s_in_son[0] = -1.0 + 4.0*s_fraction[0];
196  s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
197  s_in_son[2] = -1.0 + 4.0*(s_fraction[2]-0.5);
198  }
199  } //End of on top
200  } //End of on left
201  //Otherwise its on the right
202  else
203  {
204  //On the down side
205  if(s_fraction[1] < 0.5)
206  {
207  //On the back side
208  if(s_fraction[2] < 0.5)
209  {
210  //It's the right down back son
211  son = RDB;
212  s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
213  s_in_son[1] = -1.0 + 4.0*s_fraction[1];
214  s_in_son[2] = -1.0 + 4.0*s_fraction[2];
215  }
216  //On the front side
217  else
218  {
219  //It's the right down front son
220  son = RDF;
221  s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
222  s_in_son[1] = -1.0 + 4.0*s_fraction[1];
223  s_in_son[2] = -1.0 + 4.0*(s_fraction[2]-0.5);
224  }
225  } //End of on down side
226  //else its on the top
227  else
228  {
229  //On the back
230  if(s_fraction[2] < 0.5)
231  {
232  //It's the right up back son
233  son = RUB;
234  s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
235  s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
236  s_in_son[2] = -1.0 + 4.0*s_fraction[2];
237  }
238  //On the front side
239  else
240  {
241  //It's the right up front son
242  son = RUF;
243  s_in_son[0] = -1.0 + 4.0*(s_fraction[0]-0.5);
244  s_in_son[1] = -1.0 + 4.0*(s_fraction[1]-0.5);
245  s_in_son[2] = -1.0 + 4.0*(s_fraction[2]-0.5);
246  }
247  }
248  }
249 
250  //Get the pointer to the son element
251  RefineableQSpectralElement<3>* son_el_pt =
252  dynamic_cast<RefineableQSpectralElement<3>*>(
253  this->tree_pt()->son_pt(son)->object_pt());
254 
255  //If we are rebuilding, then worry about the boundary conditions
256  //Find the boundary of the node
257  //Initially none
258  int boundary=Tree::OMEGA;
259  //If we are on the left boundary
260  if(i0==0) {boundary = L;}
261  //If we are on the right boundary
262  else if(i0==n_p-1) {boundary = R;}
263 
264  //If we are on the lower boundary
265  if(i1==0)
266  {
267  //If we already have already set the boundary, we're on an edge
268  switch(boundary)
269  {
270  case L:
271  boundary = LD;
272  break;
273  case R:
274  boundary = RD;
275  break;
276  //Boundary not set
277  default:
278  boundary = D;
279  break;
280  }
281  }
282  //If we are the northern bounadry
283  else if(i1==n_p-1)
284  {
285  //If we already have a boundary
286  switch(boundary)
287  {
288  case L:
289  boundary = LU;
290  break;
291  case R:
292  boundary = RU;
293  break;
294  default:
295  boundary = U;
296  break;
297  }
298  }
299 
300  //If we are on the back face
301  if(i2==0)
302  {
303  //If we already have boundaries
304  switch(boundary)
305  {
306  case L:
307  boundary = LB;
308  break;
309  case R:
310  boundary = RB;
311  break;
312  case U:
313  boundary = UB;
314  break;
315  case D:
316  boundary = DB;
317  break;
318  case LD:
319  boundary = LDB;
320  break;
321  case RD:
322  boundary = RDB;
323  break;
324  case LU:
325  boundary = LUB;
326  break;
327  case RU:
328  boundary = RUB;
329  break;
330  default:
331  boundary = B;
332  break;
333  }
334  }
335  else if(i2==n_p-1)
336  {
337  //If we already have boundaries
338  switch(boundary)
339  {
340  case L:
341  boundary = LF;
342  break;
343  case R:
344  boundary = RF;
345  break;
346  case U:
347  boundary = UF;
348  break;
349  case D:
350  boundary = DF;
351  break;
352  case LD:
353  boundary = LDF;
354  break;
355  case RD:
356  boundary = RDF;
357  break;
358  case LU:
359  boundary = LUF;
360  break;
361  case RU:
362  boundary = RUF;
363  break;
364  default:
365  boundary = F;
366  break;
367  }
368  }
369 
370  // set of boundaries that this edge in the son lives on
371  std::set<unsigned> boundaries;
372  //If we are on a boundary of the Element, find the
373  //mesh boundaries on which we live
374  //The boundaries will be common to the son because there can be
375  //no rotations here
376  if(boundary!=Tree::OMEGA)
377  {
378  son_el_pt->get_boundaries(boundary,boundaries);
379  }
380 
381  // If the node lives on a boundary:
382  // construct a boundary node,
383  // get boundary conditions and
384  // update both lookup schemes
385  if (boundaries.size()>0)
386  {
387  //Construct the new node
388  this->node_pt(jnod) =
389  this->construct_boundary_node(jnod,time_stepper_pt);
390 
391  //Get the boundary conditions from the son
392  Vector<int> bound_cons(ncont_interpolated_values());
393  son_el_pt->get_bcs(boundary,bound_cons);
394 
395  //Loop over the values and pin if necessary
396  unsigned nval = this->node_pt(jnod)->nvalue();
397  for(unsigned k=0;k<nval;k++)
398  {
399  if(bound_cons[k]) {this->node_pt(jnod)->pin(k);}
400  }
401 
402  // Solid node? If so, deal with the positional boundary
403  // conditions:
404  SolidNode* solid_node_pt =
405  dynamic_cast<SolidNode*>(this->node_pt(jnod));
406  if (solid_node_pt!=0)
407  {
408  //Get the positional boundary conditions from the father:
409  unsigned n_dim = this->node_pt(jnod)->ndim();
410  Vector<int> solid_bound_cons(n_dim);
411  RefineableSolidQElement<3>* son_solid_el_pt=
412  dynamic_cast<RefineableSolidQElement<3>*>(son_el_pt);
413 #ifdef PARANOID
414  if (son_solid_el_pt==0)
415  {
416  std::string error_message =
417  "We have a SolidNode outside a refineable SolidElement\n";
418  error_message +=
419  "during mesh refinement -- this doesn't make sense\n";
420 
421  throw OomphLibError(
422  error_message,
423  OOMPH_CURRENT_FUNCTION,
424  OOMPH_EXCEPTION_LOCATION);
425  }
426 #endif
427  son_solid_el_pt->get_solid_bcs(boundary,solid_bound_cons);
428 
429  //Loop over the positions and pin, if necessary
430  for(unsigned k=0;k<n_dim;k++)
431  {
432  if (solid_bound_cons[k]) {solid_node_pt->pin_position(k);}
433  }
434  }
435 
436  //Next we update the boundary look-up schemes
437  //Loop over the boundaries stored in the set
438  for(std::set<unsigned>::iterator it = boundaries.begin();
439  it != boundaries.end(); ++it)
440  {
441  //Add the node to the boundary
442  mesh_pt->add_boundary_node(*it,this->node_pt(jnod));
443 
444  //If we have set an intrinsic coordinate on this
445  //mesh boundary then it must also be interpolated on
446  //the new node
447  //Now interpolate the intrinsic boundary coordinate
448  if(mesh_pt->boundary_coordinate_exists(*it)==true)
449  {
450  Vector<double> zeta(2);
451  son_el_pt->interpolated_zeta_on_face(*it,boundary,
452  s_in_son,zeta);
453 
454  this->node_pt(jnod)->set_coordinates_on_boundary(*it,zeta);
455  }
456  }
457  }
458  //Otherwise we construct a normal "bulk" node
459  else
460  {
461  //Construct the new node
462  this->node_pt(jnod) =
463  this->construct_node(jnod,time_stepper_pt);
464  }
465 
466  //Now we set the position and values at the newly created node
467 
468  // In the first instance use macro element or FE representation
469  // to create past and present nodal positions.
470  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
471  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
472  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
473  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
474  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
475  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
476  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
477 
478  //Loop over # of history values
479  for(unsigned t=0;t<ntstorage;t++)
480  {
481  //Get the position from the son
482  Vector<double> x_prev(3);
483 
484  //Now let's fill in the value
485  son_el_pt->get_x(t,s_in_son,x_prev);
486  for(unsigned i=0;i<3;i++)
487  {
488  this->node_pt(jnod)->x(t,i) = x_prev[i];
489  }
490  }
491 
492  // Now set the values
493  // Loop over all history values
494  for(unsigned t=0;t<ntstorage;t++)
495  {
496  // Get values from father element
497  // Note: get_interpolated_values() sets Vector size itself.
498  Vector<double> prev_values;
499  son_el_pt->get_interpolated_values(t,s_in_son,prev_values);
500 
501  //Initialise the values at the new node
502  for(unsigned k=0;k<this->node_pt(jnod)->nvalue();k++)
503  {
504  this->node_pt(jnod)->set_value(t,k,prev_values[k]);
505  }
506  }
507 
508  //Add the node to the mesh
509  mesh_pt->add_node_pt(this->node_pt(jnod));
510 
511  } //Node has been constructed
512 
513  //Algebraic stuff here
514  //Check whether the element is an algebraic element
515  AlgebraicElementBase* alg_el_pt =
516  dynamic_cast<AlgebraicElementBase*>(this);
517 
518  //If we do have an algebraic element
519  if(alg_el_pt!=0)
520  {
521  std::string error_message =
522  "Have not implemented rebuilding from sons for";
523  error_message +=
524  "Algebraic Spectral elements yet\n";
525 
526  throw
527  OomphLibError(error_message,
528  "RefineableQSpectralElement::rebuild_from_sons()",
529  OOMPH_EXCEPTION_LOCATION);
530  }
531 
532  } //End of the case when the node was not built
533  }
534  }
535  }
536  }
537 
538  /// Overload the nodes built function
539  virtual bool nodes_built()
540  {
541  unsigned n_node = this->nnode();
542  for(unsigned n=0;n<n_node;n++)
543  {
544  if(node_pt(n)==0) {return false;}
545  }
546  //If we get to here, OK
547  return true;
548  }
549 
550 
551 };
552 
553 }
554 
555 #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
virtual bool nodes_built()
Overload the nodes built function.
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
void rebuild_from_sons(Mesh *&mesh_pt)
The only thing to add is rebuild from sons.
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
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...
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:602
RefineableQSpectralElement(const RefineableQSpectralElement< 3 > &dummy)
Broken copy constructor.
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1693
Refineable version of Solid brick elements.
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
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.
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
long RefineableQElement< 2 > _build
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