nodes.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1297 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-08-21 13:56:24 +0100 (Mon, 21 Aug 2017) $
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 //Functions for the Node/Data/SolidNode classes
31 
32 #include<algorithm>
33 
34 //oomph-lib headers
35 #include "nodes.h"
36 #include "timesteppers.h"
37 
38 
39 namespace oomph
40 {
41 
42 ///////////////////////////////////////////////////////////////////
43 ///////////////////////////////////////////////////////////////////
44 //Functions for the Data class
45 ///////////////////////////////////////////////////////////////////
46 ///////////////////////////////////////////////////////////////////
47 
48  //=================================================================
49  /// \short Private function to check that the arguments are within
50  /// the range of the stored data values and timesteps.
51  //=================================================================
52  void Data::range_check(const unsigned &t, const unsigned &i) const
53  {
54  //If either the value or the time history value are out of range
55  if((i>= Nvalue) || (t >= ntstorage()))
56  {
57  std::ostringstream error_message;
58  //Value range check
59  if(i>= Nvalue)
60  {
61  error_message << "Range Error: Value " << i
62  << " is not in the range (0,"
63  << Nvalue-1 << ")";
64  }
65  //Time range check
66  if(t >= ntstorage())
67  {
68  error_message << "Range Error: Time Value " << t
69  << " is not in the range (0,"
70  << ntstorage() - 1 << ")";
71  }
72  //Throw the error
73  throw OomphLibError(error_message.str(),
74  OOMPH_CURRENT_FUNCTION,
75  OOMPH_EXCEPTION_LOCATION);
76  }
77  }
78 
79 
80 //====================================================================
81 /// Add the pointer data_pt to the internal storage used to keep track
82 /// of copies of the Data object.
83 //====================================================================
84  void Data::add_copy(Data* const &data_pt)
85  {
86  //Find the current number of copies
87  const unsigned n_copies = Ncopies;
88  //Allocate storage for the pointers to the new number of copies
89  Data** new_copy_of_data_pt = new Data*[n_copies+1];
90  //Copy over the exisiting pointers
91  for(unsigned i=0;i<n_copies;i++)
92  {new_copy_of_data_pt[i] = Copy_of_data_pt[i];}
93  //Add the new pointer to the end
94  new_copy_of_data_pt[n_copies] = data_pt;
95 
96  //Delete the old storage
97  delete[] Copy_of_data_pt;
98  //Allocate the new storage
99  Copy_of_data_pt = new_copy_of_data_pt;
100  //Increase the number of copies
101  ++Ncopies;
102  }
103 
104 //=====================================================================
105 /// Remove the pointer data_pt from the internal storage used to keep
106 /// track of copies
107 //=====================================================================
108  void Data::remove_copy(Data* const &data_pt)
109  {
110  //Find the current number of copies
111  const unsigned n_copies = Ncopies;
112  //Index of the copy
113  unsigned data_index = n_copies;
114  //Check that the existing data is actually a copy
115  for(unsigned i=0;i<n_copies;i++)
116  {
117  if(Copy_of_data_pt[i]==data_pt) {data_index = i; break;}
118  }
119 
120  //If we have not found an index throw an error
121  if(data_index==n_copies)
122  {
123  std::ostringstream error_stream;
124  error_stream << "Data pointer " << data_pt
125  << " is not stored as a copy of the data object " << this
126  << std::endl;
127  throw OomphLibError(error_stream.str(),
128  OOMPH_CURRENT_FUNCTION,
129  OOMPH_EXCEPTION_LOCATION);
130  }
131 
132  //If we still here remove the data
133  //Allocate storage for the pointers to the new number of copies
134  Data** new_copy_of_data_pt = new Data*[n_copies-1];
135 
136  unsigned index=0;
137  //Copy over the exisiting pointers
138  for(unsigned i=0;i<n_copies;i++)
139  {
140  //If we are not at the copied index
141  if(i!=data_index)
142  {
143  //Copy the data across
144  new_copy_of_data_pt[index] = Copy_of_data_pt[i];
145  //Increase the index
146  ++index;
147  }
148  }
149 
150  //Delete the old storage
151  delete[] Copy_of_data_pt;
152  //Allocate the new storage
153  Copy_of_data_pt = new_copy_of_data_pt;
154  //Set the new number of copies
155  --Ncopies;
156  }
157 
158 //================================================================
159 /// \short Helper function that should be overloaded in classes
160 /// that contain copies of Data. The function must
161 /// reset the internal pointers to the copied data. This is used
162 /// when resizing data to ensure that all the pointers remain valid.
163 /// The base Data class cannot be a copy, so throw an error
164  //==================================================================
166  {
167  throw OomphLibError("Data can never be a copy",
168  OOMPH_CURRENT_FUNCTION,
169  OOMPH_EXCEPTION_LOCATION);
170  }
171 
172 //=======================================================================
173 /// \short Helper function that should be overloaded classes
174 /// that contain copies of data. The function must
175 /// unset (NULL out) the internal pointers to the copied data.
176 /// This is used when destructing data to ensure that all pointers remain
177 /// valid.
178 //======================================================================
180  {
181  throw OomphLibError("Data can never be a copy",
182  OOMPH_CURRENT_FUNCTION,
183  OOMPH_EXCEPTION_LOCATION);
184  }
185 
186 //================================================================
187 /// Delete all the storage allocated by the Data object and
188 /// set its pointers to NULL
189 //================================================================
191  {
192  //If we have nulled out the storage already return immediately
193  if((Value==0) && (Eqn_number==0)) {return;}
194 
195  //Delete the double storage arrays at once (they were allocated at once)
196  delete[] Value[0];
197  //Delete the pointers to the arrays.
198  delete[] Value; delete[] Eqn_number;
199  //Null out the pointers
200  Value = 0; Eqn_number = 0;
201  }
202 
203 //================================================================
204 /// Default (steady) timestepper for steady Data
205 //================================================================
207 
208 //================================================================
209 /// Static "Magic number" to indicate pinned values
210 //================================================================
211  long Data::Is_pinned=-1;
212 
213 //================================================================
214 /// \short Static "Magic number" to indicate values that haven't
215 /// been classified as pinned or free
216 //================================================================
217 long Data::Is_unclassified=-10;
218 
219 //================================================================
220 /// Static "Magic number" to indicate that the value is constrained,
221 /// usually because is it associated with non-conforming data,
222 /// otherwise known as hanging nodes
223 //================================================================
224 long Data::Is_constrained=-2;
225 
226 /// \short Static "Magic number" used in place of the equation number to
227 /// indicate that the value is pinned, but only for the duration of a
228 /// segregated solve.
230 
231 
232 //================================================================
233 /// Default constructor.
234 //================================================================
235  Data::Data() : Value(0), Eqn_number(0),
236  Time_stepper_pt(Data::Default_static_time_stepper_pt),
237  Copy_of_data_pt(0),
238  Ncopies(0), Nvalue(0)
239 #ifdef OOMPH_HAS_MPI
240  , Non_halo_proc_ID(-1)
241 #endif
242 
243  {}
244 
245 //================================================================
246 /// Default constructor for steady problems. Memory is assigned for a given
247 /// number of values, which are assumed to be free (not pinned)
248 //================================================================
249  Data::Data(const unsigned &initial_n_value) :
250  Value(0), Eqn_number(0),
251  Time_stepper_pt(Data::Default_static_time_stepper_pt),
252  Copy_of_data_pt(0),
253  Ncopies(0),
254  Nvalue(initial_n_value)
255 #ifdef OOMPH_HAS_MPI
256  , Non_halo_proc_ID(-1)
257 #endif
258  {
259  //Only bother to do something if there are values
260  if(initial_n_value > 0)
261  {
262  //Allocate initial_n_value values in the value and equation number
263  //storage schemes.
264  Value = new double*[initial_n_value];
265  Eqn_number = new long[initial_n_value];
266 
267  //Allocate contiguous arrays of doubles and longs to
268  //hold the data values.
269  double *values = new double[initial_n_value];
270 
271  //Set the pointers to the data values and equation numbers
272  //and initialise the actual values.
273  for(unsigned i=0;i<initial_n_value;i++)
274  {
275  //Set the pointer from the address in the contiguous array
276  Value[i] = &values[i];
277  //Initialise the value to zero
278  Value[i][0] = 0.0;
279  //Initialise the equation number to Is_unclassified
281  }
282  }
283  }
284 
285 //================================================================
286 /// Constructor for unsteady problems. Memory is assigned for a given
287 /// number of values; and the additional storage required by the Timestepper.
288 /// The values are assumed to be free (not pinned).
289 //================================================================
290 Data::Data(TimeStepper* const &time_stepper_pt_,
291  const unsigned &initial_n_value,
292  const bool &allocate_storage)
293  :
294  Value(0), Eqn_number(0), Time_stepper_pt(time_stepper_pt_),
295  Copy_of_data_pt(0),
296  Ncopies(0),
297  Nvalue(initial_n_value)
298 #ifdef OOMPH_HAS_MPI
299  , Non_halo_proc_ID(-1)
300 #endif
301 {
302  //If we are in charge of allocating the storage,
303  //and there are data to allocate, do so
304  if((allocate_storage) && (initial_n_value > 0))
305  {
306  //Allocate storage for initial_n_value equation numbers
307  Eqn_number = new long[initial_n_value];
308 
309  //Locally cache the number of time history values
310  const unsigned n_tstorage = ntstorage();
311 
312  //There will be initial_n_value pointers each addressing and array
313  //of n_tstorage doubles.
314  Value = new double*[initial_n_value];
315 
316  //Allocate all the data values in one big array for data locality.
317  double *values = new double[initial_n_value*n_tstorage];
318 
319  //Set the pointers to the data values and equation numbers
320  for(unsigned i=0;i<initial_n_value;i++)
321  {
322  //Set the pointers to the start of the time history values
323  //allocated for each value.
324  Value[i] = &values[i*n_tstorage];
325  //Initialise all values to zero
326  for(unsigned t=0;t<n_tstorage;t++) {Value[i][t] = 0.0;}
327  //Initialise the equation number to be unclassified.
329  }
330  }
331 }
332 
333  /// Data output operator: output equation numbers and values at all times,
334  /// along with any extra information stored for the timestepper.
335  std::ostream& operator<< (std::ostream &out, const Data& d)
336  {
337  const unsigned nvalue = d.nvalue();
338  const unsigned nt = d.ntstorage();
339 
340  out << "Data: [" << std::endl;
341 
342  for(unsigned j=0; j<nvalue; j++)
343  {
344  out << "global eq " << d.eqn_number(j) << ": [";
345  for(unsigned t=0; t<nt-1; t++)
346  {
347  out << d.value(t, j) << ", ";
348  }
349  out << d.value(nt-1, j) << "]" << std::endl;
350  }
351  out << "]" << std::endl;
352 
353 
354  return out;
355  }
356 
357  /// Node output operator: output equation numbers and values at all times,
358  /// along with any extra information stored for the timestepper.
359  std::ostream& operator<< (std::ostream &out, const Node& nd)
360  {
361  const unsigned nt = nd.ntstorage();
362  const unsigned dim = nd.ndim();
363 
364  // Output position, only doing current value for now but add position
365  // history if you need it - David.
366  out << "Position: [";
367  for(unsigned j=0; j<dim; j++)
368  {
369  out << "dimension " << dim << ": [";
370  for(unsigned t=0; t<nt-1; t++)
371  {
372  out << nd.x(t, j) << ", ";
373  }
374  out << nd.x(nt-1, j) << "]" << std::endl;
375  }
376  out << "]" << std::endl;
377 
378  // Use the function for data to output the data (we can't use overloading
379  // here because operator<< is not a class function, so instead we
380  // typecast the node to a data).
381  out << dynamic_cast<const Data&>(nd);
382  return out;
383  }
384 
385 
386 //================================================================
387 /// Set a new TimeStepper be resizing the appropriate storage.
388 /// Equation numbering (if already performed) will be unaffected.
389 /// The current (zero) values will be unaffected, but all other entries
390 /// will be set to zero.
391 //================================================================
392  void Data::set_time_stepper(TimeStepper* const &time_stepper_pt,
393  const bool &preserve_existing_data)
394 {
395  //If the timestepper is unchanged do nothing
396  if(Time_stepper_pt==time_stepper_pt) {return;}
397 
398  //Find the amount of data to be preserved
399  //Default is just the current values
400  unsigned n_preserved_tstorage = 1;
401  if(preserve_existing_data) {n_preserved_tstorage = this->ntstorage();}
402 
403  //Set the new time stepper
405 
406  //If the data is a copy don't mess with it
407  if(this->is_a_copy()) {return;}
408 
409  //Find the current number of values
410  const unsigned n_value = nvalue();
411 
412  //IF there are data to allocate, do so
413  if(n_value > 0)
414  {
415  //Locally cache the new number of time storage values
416  const unsigned n_tstorage = time_stepper_pt->ntstorage();
417 
418  //Allocate all the data values in one big array for data locality.
419  double *values = new double[n_value*n_tstorage];
420 
421  //Copy the old "preserved" values into the new storage scheme
422  //Make sure that we limit the values to the level of storage
423  if(n_tstorage < n_preserved_tstorage) {n_preserved_tstorage = n_tstorage;}
424  for(unsigned i=0;i<n_value;i++)
425  {
426  for(unsigned t=0;t<n_preserved_tstorage;t++)
427  {
428  values[i*n_tstorage + t] = Value[i][t];
429  }
430  }
431 
432  //Now delete the old value storage
433  delete[] Value[0];
434 
435  //Reset the pointers to the new data values
436  for(unsigned i=0;i<n_value;i++)
437  {
438  Value[i] = &values[i*n_tstorage];
439  //Initialise all new time storage values to zero
440  for(unsigned t=n_preserved_tstorage;t<n_tstorage;t++) {Value[i][t] = 0.0;}
441  }
442 
443  //Update any pointers in any copies of this data
444  for(unsigned i=0;i<Ncopies;i++)
445  {
447  }
448  }
449 }
450 
451 //================================================================
452 ///Virtual destructor, deallocates memory assigned for data
453 //================================================================
455 {
456  //If we have any copies clear their pointers
457  for(unsigned i=0;i<Ncopies;i++)
458  {
460  }
461 
462  //Now delete the storage allocated for pointers to the copies
463  delete[] Copy_of_data_pt; Copy_of_data_pt=0;
464 
465  //Clean up the allocated storage
467 }
468 
469 
470 //==================================================================
471 /// Compute Vector of values (dofs or pinned) at this Data object
472 //==================================================================
473 void Data::value(Vector<double>& values) const
474 {
475  //Loop over all the values and set the appropriate value
476  const unsigned n_value = nvalue();
477  for(unsigned i=0;i<n_value;i++) {values[i] = value(i);}
478 }
479 
480 //==================================================================
481 /// Compute Vector of values (dofs or pinned) at this node
482 /// at time level t (t=0: present; t>0: previous)
483 //==================================================================
484 void Data::value(const unsigned& t, Vector<double>& values) const
485 {
486  //Loop over all the values and set the value at time level t
487  const unsigned n_value = nvalue();
488  for(unsigned i=0;i<n_value;i++) {values[i] = value(t,i);}
489 }
490 
491 
492 
493 //================================================================
494 /// Assign (global) equation number. Overloaded version for nodes.
495 /// Checks if a hanging value has a non-negative equation number
496 /// and if so shouts and then sets it to Is_constrained. Then drops
497 /// down to Data version which does the actual work
498 //================================================================
499 void Node::assign_eqn_numbers(unsigned long &global_number,
500  Vector<double *> &dof_pt)
501 {
502 
503  //Loop over the number of values
504  const unsigned eqn_number_range = nvalue();
505  for(unsigned i=0;i<eqn_number_range;i++)
506  {
507  // Is it hanging and not constrained or pinned? If so shout and constrain it
508  if((is_hanging(i)) && (!is_constrained(i)) && (!is_pinned(i)))
509  {
510 #ifdef PARANOID
511  std::ostringstream warn_message;
512  warn_message
513  << "Node::assign_eqn_numbers(...) noticed that " << i << " -th value\n"
514  << "is hanging but not constrained. This shouldn't happen and is\n"
515  << "probably because a hanging value was unpinned manually\n"
516  << "Rectifying this now...\n";
517  OomphLibWarning(warn_message.str(),
518  OOMPH_CURRENT_FUNCTION,
519  OOMPH_EXCEPTION_LOCATION);
520 #endif
521  constrain(i);
522  }
523  }
524  // Now descend
525  Data::assign_eqn_numbers(global_number, dof_pt);
526 }
527 
528 
529 //=====================================================================
530 /// If pointer parameter_pt addresses internal data values then return
531 /// return true, otherwise return false
532 //======================================================================
533 bool Data::does_pointer_correspond_to_value(double*const &parameter_pt)
534  {
535  //If there is no value data then return false
536  if(Value==0) {return false;}
537 
538  //Find the amount of data stored
539  const unsigned n_value = nvalue();
540  const unsigned n_time = ntstorage();
541  const unsigned n_storage = n_value*n_time;
542 
543  //Pointer to the local data
544  double *local_value_pt = Value[0];
545 
546  //Loop over data and if we find the pointer then return true
547  for(unsigned i=0;i<n_storage;++i)
548  {
549  if(parameter_pt==(local_value_pt+i)) {return true;}
550  }
551 
552  //If we get to here we haven't found the data
553  return false;
554  }
555 
556 
557 
558 //================================================================
559 /// Copy Data values from specified Data object
560 //================================================================
561 void Data::copy(Data* orig_data_pt)
562 {
563 
564  //Find the amount of data stored
565  const unsigned n_value = nvalue();
566  const unsigned n_time = ntstorage();
567 
568  // Check # of values:
569  const unsigned long n_value_orig=orig_data_pt->nvalue();
570  if (n_value!=n_value_orig)
571  {
572  std::ostringstream error_stream;
573  error_stream << "The number of values, " << n_value
574  << " is not the same of those in the original data "
575  << n_value_orig << std::endl;
576 
577  throw OomphLibError(error_stream.str(),
578  OOMPH_CURRENT_FUNCTION,
579  OOMPH_EXCEPTION_LOCATION);
580  }
581  const unsigned long n_time_orig=orig_data_pt->ntstorage();
582  if (n_time!=n_time_orig)
583  {
584  std::ostringstream error_stream;
585  error_stream << "The number of time history values, " << n_time
586  << " is not the same of those in the original data "
587  << n_time_orig << std::endl;
588 
589  throw OomphLibError(error_stream.str(),
590  OOMPH_CURRENT_FUNCTION,
591  OOMPH_EXCEPTION_LOCATION);
592  }
593 
594  // Read data
595  for(unsigned t=0;t<n_time;t++)
596  {
597  for(unsigned j=0;j<n_value;j++)
598  {
599  set_value(t,j,orig_data_pt->value(t,j));
600  }
601  }
602 }
603 
604 
605 //================================================================
606 ///Dump data object to a file
607 //================================================================
608 void Data::dump(std::ostream& dump_file) const
609 {
610  //Find the amount of storage used
611  const unsigned value_pt_range = nvalue();
612  const unsigned time_steps_range = ntstorage();
613 
614  //Only write data if there is some stored
615  if (value_pt_range*time_steps_range > 0)
616  {
617  dump_file << value_pt_range << " # number of data values" << std::endl;
618  dump_file << time_steps_range << " # number of doubles for time history"
619  << std::endl;
620 
621  // Write data
622  for(unsigned t=0;t<time_steps_range;t++)
623  {
624  for(unsigned j=0;j<value_pt_range;j++)
625  {
626  dump_file << value(t,j) << std::endl ;
627  }
628  }
629  }
630 }
631 
632 //================================================================
633 ///Read data object from file
634 //================================================================
635 void Data::read(std::ifstream& restart_file)
636 {
637  std::string input_string;
638  std::ostringstream error_stream;
639 
640  //Find the amount of data stored
641  const unsigned value_pt_range = nvalue();
642  const unsigned time_steps_range = ntstorage();
643 
644  //Only read in data if there is some storage available
645  if (value_pt_range*time_steps_range > 0)
646  {
647  // Read line up to termination sign
648  getline(restart_file,input_string,'#');
649  // Ignore rest of line
650  restart_file.ignore(80,'\n');
651  // Check # of values:
652  const unsigned long check_nvalues=atoi(input_string.c_str());
653  if (check_nvalues!=value_pt_range)
654  {
655  error_stream
656  << "Number of values stored in dump file is not equal to the amount "
657  << "of storage allocated in Data object "
658  << check_nvalues << " " << value_pt_range;
659  if (check_nvalues>value_pt_range)
660  {
661  error_stream << " [ignoring extra entries]";
662  }
663  error_stream << std::endl;
664  Node* nod_pt=dynamic_cast<Node*>(this);
665  if (nod_pt!=0)
666  {
667  unsigned n_dim=nod_pt->ndim();
668  error_stream << "Node coordinates: ";
669  for (unsigned i=0;i<n_dim;i++)
670  {
671  error_stream << nod_pt->x(i) << " ";
672  }
673  error_stream << nod_pt << " ";
674 #ifdef OOMPH_HAS_MPI
675  if (nod_pt->is_halo())
676  {
677  error_stream << " (halo)\n";
678  }
679  else
680  {
681  error_stream << " (not halo)\n";
682  }
683 #endif
684  }
685  if (check_nvalues<value_pt_range)
686  {
687  throw OomphLibError(error_stream.str(),
688  OOMPH_CURRENT_FUNCTION,
689  OOMPH_EXCEPTION_LOCATION);
690  }
691  }
692 
693  // Read line up to termination sign
694  getline(restart_file,input_string,'#');
695 
696  // Ignore rest of line
697  restart_file.ignore(80,'\n');
698 
699  // Check # of values:
700  const unsigned check_ntvalues=atoi(input_string.c_str());
701 
702  // Dynamic run restarted from steady run
703  if (check_ntvalues<time_steps_range)
704  {
705  std::ostringstream warning_stream;
706  warning_stream
707  << "Number of time history values in dump file is less "
708  << "than the storage allocated in Data object: "
709  << check_ntvalues << " " << time_steps_range << std::endl;
710  warning_stream
711  << "We're using steady data as initial data for unsteady \n"
712  << "run. I'll fill in the remaining history values with zeroes. \n"
713  << "If you don't like this \n"
714  << "you'll have to overwrite this yourself with whatever is \n "
715  << "appropriate for your timestepping scheme. \n";
716  //Issue the warning
717  OomphLibWarning(warning_stream.str(),
718  "Data::read()",
719  OOMPH_EXCEPTION_LOCATION);
720 
721  // Read data
722  for(unsigned t=0;t<time_steps_range;t++)
723  {
724  for(unsigned j=0;j<check_nvalues;j++)
725  {
726  if (t==0)
727  {
728  // Read line
729  getline(restart_file,input_string);
730 
731  // Transform to double
732  if (j<value_pt_range)
733  {
734  set_value(t,j,atof(input_string.c_str()));
735  }
736  else
737  {
738  error_stream
739  << "Not setting j=" << j
740  << " -th history restart value [t = " << t << " ] to "
741  << atof(input_string.c_str()) << " because Data "
742  << " hasn't been sufficiently resized\n";
743  }
744  }
745  else
746  {
747  if (j<value_pt_range)
748  {
749  set_value(t,j,0.0);
750  }
751  else
752  {
753  error_stream
754  << "Not setting j=" << j
755  << " -th restart history value [t = " << t << " ] to "
756  << 0.0 << " because Data "
757  << " hasn't been sufficiently resized\n";
758  }
759  }
760  }
761  }
762  }
763  // Static run restarted from unsteady run
764  else if (check_ntvalues>time_steps_range)
765  {
766  std::ostringstream warning_stream;
767  warning_stream
768  << "Warning: number of time history values in dump file is greater "
769  << "than the storage allocated in Data object: "
770  << check_ntvalues << " " << time_steps_range << std::endl;
771  warning_stream << "We're using the current values from an unsteady \n"
772  << "restart file to initialise a static run. \n";
773  //Issue the warning
774  OomphLibWarning(warning_stream.str(),
775  "Data::read()",
776  OOMPH_EXCEPTION_LOCATION);
777 
778  // Read data
779  for(unsigned t=0;t<check_ntvalues;t++)
780  {
781  for(unsigned j=0;j<check_nvalues;j++)
782  {
783  // Read line
784  getline(restart_file,input_string);
785  if (t==0)
786  {
787  // Transform to double
788  if (j<value_pt_range)
789  {
790  set_value(t,j,atof(input_string.c_str()));
791  }
792  else
793  {
794  error_stream
795  << "Not setting j=" << j
796  << " -th restart history value [t = " << t << " ] to "
797  << atof(input_string.c_str()) << " because Data "
798  << " hasn't been sufficiently resized\n";
799  }
800  }
801  }
802  }
803  }
804  // Proper dynamic restart
805  else
806  {
807  // Read data
808  for(unsigned t=0;t<time_steps_range;t++)
809  {
810  for(unsigned j=0;j<check_nvalues;j++)
811  {
812  // Read line
813  getline(restart_file,input_string);
814 
815  // Transform to double
816  if (j<value_pt_range)
817  {
818  set_value(t,j,atof(input_string.c_str()));
819  }
820  else
821  {
822  error_stream << "Not setting j=" << j
823  << " -th restart history value [t = " << t << " ] to "
824  << atof(input_string.c_str()) << " because Data "
825  << " hasn't been sufficiently resized\n";
826  }
827  }
828  }
829  }
830  if (check_nvalues>value_pt_range)
831  {
832  OomphLibWarning(error_stream.str(),
833  "Data::read()",
834  OOMPH_EXCEPTION_LOCATION);
835  }
836  }
837 
838 
839 }
840 
841 
842 //===================================================================
843 /// Return the total number of doubles stored per value to record
844 /// the time history of ecah value. The information is read from the
845 /// time stepper
846 //===================================================================
847 unsigned Data::ntstorage() const {return Time_stepper_pt->ntstorage();}
848 
849 //================================================================
850 ///Assign (global) equation number.
851 /// This function does NOT initialise the value because
852 /// if we're using things like node position as variables in the problem
853 /// they will have been set before the call to assign equation numbers
854 /// and setting it to zero will wipe it out :(.
855 ///
856 /// Pass:
857 /// - current number of global dofs global_number (which gets incremented)
858 /// - the Vector of pointers to global dofs (to which new dofs
859 /// get added)
860 //================================================================
861 void Data::assign_eqn_numbers(unsigned long &global_number,
862  Vector<double *> &dof_pt)
863 {
864 
865  //Loop over the number of variables
866  //Set temporary to hold range
867  const unsigned eqn_number_range = Nvalue;
868  for(unsigned i=0;i<eqn_number_range;i++)
869  {
870 #ifdef OOMPH_HAS_MPI
871  // Is the node a halo? If so, treat it as pinned for now
872  // This will be overwritten with the actual equation number
873  // during the synchronisation phase.
874  if (is_halo())
875  {
876  eqn_number(i) = Is_pinned;
877  }
878  else
879 #endif
880  {
881  //Boundary conditions test: if it's not a pinned or constrained variable,
882  //The assign a new global equation number
883  if((!is_pinned(i)) && (!is_constrained(i))
885  {
886  //Assign the equation number and increment global equation number
887  Eqn_number[i] = global_number++;
888  //Add pointer to global dof vector
889  dof_pt.push_back(value_pt(i));
890  }
891  }
892  }
893 }
894 
895 //================================================================
896 /// \short Function to describe the dofs of the Data. The ostream
897 /// specifies the output stream to which the description
898 /// is written; the string stores the currently
899 /// assembled output that is ultimately written to the
900 /// output stream by Data::describe_dofs(...); it is typically
901 /// built up incrementally as we descend through the
902 /// call hierarchy of this function when called from
903 /// Problem::describe_dofs(...)
904 //================================================================
905 void Data::describe_dofs(std::ostream& out,
906  const std::string& current_string) const
907 {
908  //Loop over the number of variables
909  const unsigned eqn_number_range = Nvalue;
910  for(unsigned i=0;i<eqn_number_range;i++)
911  {
912  int eqn_number=Eqn_number[i];
913  if (eqn_number>=0)
914  {
915  // Note: The spacing around equation number is deliberate.
916  // It allows for searching more easily as Eqn:<space>5<space> would return
917  // a unique dof, whereas Eqn:<space>5 would also return those starting with
918  // 5, such as 500.
919  out<<"Eqn: "<<eqn_number<<" | Value "<<i<<current_string<<std::endl;
920  }
921  }
922 }
923 
924 
925 //================================================================
926 /// Self-test: Have all values been classified as pinned/unpinned?
927 /// Return 0 if OK.
928 
929 //================================================================
930 unsigned Data::self_test()
931 {
932  //Initialise test flag
933  bool passed=true;
934 
935  //Loop over all equation numbers
936  const unsigned eqn_number_range = Nvalue;
937  for(unsigned i=0;i<eqn_number_range;i++)
938  {
939  //If the equation number has not been assigned, issue an error
941  {
942  passed=false;
943  oomph_info
944  << "\n ERROR: Failed Data::self_test() for i=" << i << std::endl;
945  oomph_info
946  << " (Value is not classified as pinned or free)" << std::endl;
947  }
948  }
949 
950  //Return verdict
951  if (passed) {return 0;}
952  else {return 1;}
953 }
954 
955 //================================================================
956 /// Increase the number of data values stored, useful when adding
957 /// additional data at a node, almost always Lagrange multipliers.
958 /// Note if any of the unresized data is copied, then we assume all the
959 /// resized data is copied from the same node as the unresized data.
960 //================================================================
961 void Data::resize(const unsigned &n_value)
962 {
963  //Find current number of values
964  const unsigned n_value_old = nvalue();
965  //Set the desired number of values
966  const unsigned n_value_new = n_value;
967 
968  //If the number of values hasn't changed, do nothing
969  if(n_value_new==n_value_old) {return;}
970 
971  //Put in a little safely check here
972 #ifdef PARANOID
973  if(n_value_new < n_value_old)
974  {
975  std::ostringstream error_stream;
976  error_stream
977  << "Warning : Data cannot be resized to a smaller value!" << std::endl;
978  throw OomphLibError(error_stream.str(),
979  OOMPH_CURRENT_FUNCTION,
980  OOMPH_EXCEPTION_LOCATION);
981  }
982 #endif
983 
984  //Find amount of additional time storage required
985  //N.B. We can't change timesteppers in this process
986  const unsigned t_storage = ntstorage();
987 
988  //Create new sets of pointers of the appropriate (new) size
989  double **value_new_pt = new double*[n_value_new];
990  long *eqn_number_new = new long[n_value_new];
991 
992  //Create new array of values that is contiguous in memory
993  double *values = new double[n_value_new*t_storage];
994 
995  //Copy the old values over into the new storage scheme
996  for(unsigned i=0;i<n_value_old;i++)
997  {
998  //Set pointer for the new values
999  value_new_pt[i] = &values[i*t_storage];
1000  //Copy value
1001  for(unsigned t=0;t<t_storage;t++)
1002  {value_new_pt[i][t] = Value[i][t];}
1003 
1004  //Copy equation number
1005  eqn_number_new[i] = Eqn_number[i];
1006  }
1007 
1008  //Loop over the new entries, set pointers and initialise data
1009  for(unsigned i=n_value_old;i<n_value_new;i++)
1010  {
1011  //Set the pointer
1012  value_new_pt[i] = &values[i*t_storage];
1013  //Initialise the new data values to zero
1014  for(unsigned t=0;t<t_storage;t++) {value_new_pt[i][t] = 0.0;}
1015 
1016  //Initialise the equation number to Is_unclassified
1017  eqn_number_new[i] = Is_unclassified;
1018  }
1019 
1020  //Set the number of new values
1021  Nvalue = n_value_new;
1022 
1023  //Now delete the old storage and set the new pointers
1024  if (n_value_old!=0) delete[] Value[0];
1025  delete[] Value;
1026  Value = value_new_pt;
1027  delete[] Eqn_number;
1028  Eqn_number = eqn_number_new;
1029 
1030  //Now update pointers in any copies of this data
1031  for(unsigned i=0;i<Ncopies;i++)
1032  {
1034  }
1035 }
1036 
1037 //=======================================================================
1038 /// Add pointers to all unpinned and unconstrained data to a map
1039 /// indexed by (global) equation number
1040 //=======================================================================
1041 void Data::add_value_pt_to_map(std::map<unsigned,double*> &map_of_value_pt)
1042 {
1043  //How many values does it have
1044  const unsigned n_value = this->nvalue();
1045  //Find the global equation number
1046  for(unsigned i=0;i<n_value;i++)
1047  {
1048  int global_eqn = this->eqn_number(i);
1049  //If it is a degree of freedom, add it to the map
1050  if(global_eqn >= 0)
1051  {
1052  map_of_value_pt[static_cast<unsigned>(global_eqn)] =
1053  this->value_pt(i);
1054  }
1055  }
1056 }
1057 
1058 #ifdef OOMPH_HAS_MPI
1059 //==================================================================
1060 /// Add all data and time history values to a vector that will be
1061 /// used when communicating data between processors. The function is
1062 /// virtual so that it can be overloaded by Nodes and SolidNodes to
1063 /// add the additional data present in those objects.
1064 //==================================================================
1066 {
1067  //Find the number of stored values
1068  const unsigned n_value = this->nvalue();
1069 
1070 #ifndef PARANOID
1071 
1072  //If no values are stored then return immediately
1073  if(n_value==0) {return;}
1074 
1075 #endif
1076 
1077  //Find the number of stored time data
1078  const unsigned n_tstorage = this->ntstorage();
1079 
1080  //Resize the vector to accommodate the new data
1081  const unsigned n_current_value = vector_of_values.size();
1082 
1083 #ifdef PARANOID
1084  unsigned n_debug=2;
1085 #else
1086  unsigned n_debug=0;
1087 #endif
1088 
1089  vector_of_values.resize(n_current_value + n_tstorage*n_value + n_debug);
1090 
1091  //Now add the data to the vector
1092  unsigned index = n_current_value;
1093 
1094 #ifdef PARANOID
1095  vector_of_values[index++]=n_tstorage;
1096  vector_of_values[index++]=n_value;
1097  // Now return
1098  if(n_value==0) {return;}
1099 #endif
1100 
1101  //Pointer to the first entry in the data array
1102  double* data_pt = Value[0];
1103 
1104  //Loop over values
1105  for(unsigned i=0;i<n_value;i++)
1106  {
1107  //Loop over time histories
1108  for(unsigned t=0;t<n_tstorage;t++)
1109  {
1110  //Add the data to the vector
1111  vector_of_values[index] = *data_pt;
1112 
1113  //Increment the counter and the pointer
1114  ++index;
1115  ++data_pt;
1116  }
1117  }
1118 }
1119 
1120 //==================================================================
1121 /// Read all data and time history values from a vector that will be
1122 /// used when communicating data between processors. The function is
1123 /// virtual so that it can be overloaded by Nodes and SolidNodes to
1124 /// read the additional data present in those objects. The unsigned
1125 /// index is used to indicate the start position for reading in
1126 /// the vector and will be set the end of the data that has been
1127 /// read in on return.
1128 //==================================================================
1129 void Data::read_values_from_vector(const Vector<double> &vector_of_values,
1130  unsigned &index)
1131 {
1132  //Find the number of stored values
1133  unsigned n_value = this->nvalue();
1134 
1135  //Find the number of stored time data
1136  const unsigned n_tstorage = this->ntstorage();
1137 
1138 #ifdef PARANOID
1139  unsigned orig_n_tstorage=unsigned(vector_of_values[index++]);
1140  unsigned orig_n_value=unsigned(vector_of_values[index++]);
1141  if ((orig_n_tstorage!=n_tstorage)||(orig_n_value!=n_value))
1142  {
1143  std::ostringstream error_stream;
1144  error_stream << "Non-matching number of values:\n"
1145  << "sent and local n_tstorage: " << orig_n_tstorage << " "
1146  << n_tstorage << std::endl
1147  << "sent and local n_value: " << orig_n_value << " "
1148  << n_value << std::endl;
1149  throw OomphLibError(error_stream.str(),
1150  OOMPH_CURRENT_FUNCTION,
1151  OOMPH_EXCEPTION_LOCATION);
1152  }
1153 #endif
1154 
1155  //If no values are stored, return immediately
1156  if(n_value==0) {return;}
1157 
1158  //Pointer to the first entry in the data array
1159  double* data_pt = Value[0];
1160 
1161  //Loop over values
1162  for(unsigned i=0;i<n_value;i++)
1163  {
1164  //Loop over time histories
1165  for(unsigned t=0;t<n_tstorage;t++)
1166  {
1167  //Read the data from the vector
1168  *data_pt = vector_of_values[index];
1169  //Increment the counter and the pointer
1170  ++index;
1171  ++data_pt;
1172  }
1173  }
1174 }
1175 
1176 //==================================================================
1177 /// Add all equation numbers to the vector. The function is virtual
1178 /// so that it can be overloaded by SolidNodes to add the additional
1179 /// equation numbers associated with the solid dofs in those objects.
1180 //==================================================================
1181 void Data::add_eqn_numbers_to_vector(Vector<long> &vector_of_eqn_numbers)
1182 {
1183  //Find the number of stored values
1184  const unsigned n_value = this->nvalue();
1185  //If no values are stored then return immediately
1186  if(n_value==0) {return;}
1187 
1188  //Resize the vector to accommodate the new data
1189  const unsigned n_current_value = vector_of_eqn_numbers.size();
1190  vector_of_eqn_numbers.resize(n_current_value + n_value);
1191 
1192  //Now add the data to the vector
1193  unsigned index = n_current_value;
1194  //Pointer to the first entry in the equation number array
1195  long* eqn_number_pt = Eqn_number;
1196  //Loop over values
1197  for(unsigned i=0;i<n_value;i++)
1198  {
1199  //Add the data to the vector
1200  vector_of_eqn_numbers[index] = *eqn_number_pt;
1201  //Increment the counter and the pointer
1202  ++index;
1203  ++eqn_number_pt;
1204  }
1205 }
1206 
1207 //==================================================================
1208 /// Read all equation numbers from the vector. The function is virtual
1209 /// so that it can be overloaded by SolidNodes to add the additional
1210 /// equation numbers associated with the solid dofs in those objects.
1211 /// The unsigned
1212 /// index is used to indicate the start position for reading in
1213 /// the vector and will be set the end of the data that has been
1214 /// read in on return.
1215 //==================================================================
1217  const Vector<long> &vector_of_eqn_numbers, unsigned &index)
1218 {
1219  //Find the number of stored values
1220  const unsigned n_value = this->nvalue();
1221  //If no values are stored then return immediately
1222  if(n_value==0) {return;}
1223 
1224  //Pointer to the first entry in the equation number array
1225  long* eqn_number_pt = Eqn_number;
1226  //Loop over values
1227  for(unsigned i=0;i<n_value;i++)
1228  {
1229  //Read the data from the vector
1230  *eqn_number_pt = vector_of_eqn_numbers[index];
1231  //Increment the counter and the pointer
1232  ++index;
1233  ++eqn_number_pt;
1234  }
1235 }
1236 
1237 #endif
1238 
1239 
1240 //================================================================
1241 /// Reset the pointers to the copied data
1242 //===============================================================
1244 {
1245  //Copy the pointer to the value. This will give the appropriate
1246  //"slice" of the array
1248 
1249  //Copy the pointer to the equation number
1251 }
1252 
1253 
1254 //===============================================================
1255 /// Clear the pointers to the copied data
1256 //===============================================================
1258 {
1259  Copied_data_pt = 0;
1260  Value = 0; Eqn_number = 0;
1261 }
1262 
1263 //================================================================
1264 /// Constructor, creates a HijackedData object with a single value
1265 /// that is copied from another Data object.
1266 /// The ordering of the aguments is used to
1267 /// distinguish this case from that of copying all data values, except one
1268 /// independent value.
1269 //================================================================
1271 HijackedData(const unsigned &copied_index, Data* const &data_pt) :
1272  Data(data_pt->time_stepper_pt(),1,false),
1273  Copied_data_pt(data_pt),
1274  Copied_index(copied_index)
1275 {
1276  //Don't allow copying of a copy
1277  if(data_pt->is_a_copy(copied_index))
1278  {
1279  std::ostringstream error_stream;
1280  error_stream << "The data you are trying to hijack is already a copy"
1281  << std::endl;
1282  error_stream << "Please copy the original data" << std::endl;
1283  error_stream << "In a later version, I might do this for you,"
1284  << " but not today" << std::endl;
1285 
1286  throw OomphLibError(error_stream.str(),
1287  OOMPH_CURRENT_FUNCTION,
1288  OOMPH_EXCEPTION_LOCATION);
1289  }
1290 
1291  //Copy the pointer to the value. This will give the appropriate
1292  //"slice" of the array
1293  Value = &data_pt->Value[copied_index];
1294  //Copy the pointer to the equation number
1295  Eqn_number = &data_pt->Eqn_number[copied_index];
1296  //Inform the original data that it has been copied
1297  data_pt->add_copy(this);
1298 }
1299 
1300 //=================================================================
1301 /// We do not allow Hijacked Data to be resized
1302 //=================================================================
1303 void HijackedData::resize(const unsigned &n_value)
1304 {
1305  throw OomphLibError("HijackedData cannot be resized",
1306  OOMPH_CURRENT_FUNCTION,
1307  OOMPH_EXCEPTION_LOCATION);
1308 }
1309 
1310 
1311 //////////////////////////////////////////////////////////////////////////
1312 //////////////////////////////////////////////////////////////////////////
1313 //Functions for the CopiedData class
1314 //////////////////////////////////////////////////////////////////////////
1315 //////////////////////////////////////////////////////////////////////////
1316 
1317 
1318 //================================================================
1319 /// Reset the pointers to the copied data
1320 //===============================================================
1322 {
1323  //Set the new number of values
1325 
1326  //Copy the pointer to the value. This will give the appropriate
1327  //"slice" of the array
1329 
1330  //Copy the pointer to the equation numbers
1332 }
1333 
1334 
1335 //===============================================================
1336 /// Clear ther pointers to the copied data
1337 //===============================================================
1339 {
1340  Copied_data_pt = 0;
1341  Value = 0; Eqn_number = 0;
1342 }
1343 
1344 //================================================================
1345 /// Constructor, creates a CopiedData object with all values
1346 /// copied from another Data object.
1347 //================================================================
1348 CopiedData::CopiedData(Data* const &data_pt) :
1349  Data(data_pt->time_stepper_pt(),data_pt->nvalue(),false),
1350  Copied_data_pt(data_pt)
1351 {
1352  //Don't allow copying of a copy
1353  if(data_pt->is_a_copy())
1354  {
1355  std::ostringstream error_stream;
1356  error_stream << "The data you are trying to copy is already a copy"
1357  << std::endl;
1358  error_stream << "Please copy the original data" << std::endl;
1359  error_stream << "In a later version, I might do this for you,"
1360  << " but not today" << std::endl;
1361 
1362  throw OomphLibError(error_stream.str(),
1363  OOMPH_CURRENT_FUNCTION,
1364  OOMPH_EXCEPTION_LOCATION);
1365  }
1366 
1367  //Copy the pointer to the value.
1368  Value = data_pt->Value;
1369  //Copy the pointer to the equation number
1370  Eqn_number = data_pt->Eqn_number;
1371  //Inform the original data that it has been copied
1372  data_pt->add_copy(this);
1373 }
1374 
1375 //=================================================================
1376 /// We do not allow Copied Data to be resized
1377 //=================================================================
1378 void CopiedData::resize(const unsigned &n_value)
1379 {
1380  throw OomphLibError("CopiedData cannot be resized",
1381  OOMPH_CURRENT_FUNCTION,
1382  OOMPH_EXCEPTION_LOCATION);
1383 }
1384 
1385 
1386 
1387 //////////////////////////////////////////////////////////////////////////
1388 //////////////////////////////////////////////////////////////////////////
1389 //Functions for the HangInfo class
1390 //////////////////////////////////////////////////////////////////////////
1391 //////////////////////////////////////////////////////////////////////////
1392 
1393 //================================================================
1394 /// Check that the argument is within the range of
1395 /// stored master nodes
1396 //=================================================================
1397 void HangInfo::range_check(const unsigned &i) const
1398 {
1399  //If the argument is negative or greater than the number of stored
1400  //values die
1401  if(i >= Nmaster)
1402  {
1403  std::ostringstream error_message;
1404  error_message << "Range Error: the index " << i
1405  << " is not in the range (0,"
1406  << Nmaster-1 << ")";
1407  throw OomphLibError(error_message.str(),
1408  OOMPH_CURRENT_FUNCTION,
1409  OOMPH_EXCEPTION_LOCATION);
1410  }
1411 }
1412 
1413 
1414 //=====================================================================
1415 /// Set the pointer to the i-th master node and its weight
1416 //=====================================================================
1417 void HangInfo::set_master_node_pt(const unsigned &i,
1418  Node* const &master_node_pt_,
1419  const double &weight)
1420 {
1421 #ifdef RANGE_CHECKING
1422  range_check(i);
1423 #endif
1424  Master_nodes_pt[i] = master_node_pt_;
1425  Master_weights[i] = weight;
1426 }
1427 
1428 //====================================================================
1429 /// Add (pointer to) master node and corresponding weight to
1430 /// the internally stored (pointers to) master nodes and weights.
1431 //====================================================================
1432  void HangInfo::add_master_node_pt(Node* const &master_node_pt_,
1433  const double &weight)
1434 {
1435  //Find the present number of master nodes
1436  const unsigned n_master = Nmaster;
1437  //Make new data
1438  Node* *new_master_nodes_pt = new Node*[n_master+1];
1439  double *new_master_weights = new double[n_master+1];
1440 
1441  //Copy the old values over to the new data
1442  for(unsigned i=0;i<n_master;i++)
1443  {
1444  new_master_nodes_pt[i] = Master_nodes_pt[i];
1445  new_master_weights[i] = Master_weights[i];
1446  }
1447  //Add the new values at the end
1448  new_master_nodes_pt[n_master] = master_node_pt_;
1449  new_master_weights[n_master] = weight;
1450 
1451  //Reset the pointers
1452  delete[] Master_nodes_pt; Master_nodes_pt = new_master_nodes_pt;
1453  delete[] Master_weights; Master_weights = new_master_weights;
1454  //Increase the number of master nodes
1455  ++Nmaster;
1456 }
1457 
1458 //////////////////////////////////////////////////////////////////////////
1459 //////////////////////////////////////////////////////////////////////////
1460 //Functions for the Node class
1461 //////////////////////////////////////////////////////////////////////////
1462 //////////////////////////////////////////////////////////////////////////
1463 
1464 //=================================================================
1465 /// \short Private function to check that the arguments are within
1466 /// the range of the stored coordinates, position types and time history
1467 /// values.
1468 //=================================================================
1469 void Node::x_gen_range_check(const unsigned &t, const unsigned &k,
1470  const unsigned &i) const
1471 {
1472  //Number of stored history values
1473  const unsigned position_ntstorage = Position_time_stepper_pt->ntstorage();
1474  //If any of the coordinates or time values are out of range
1475  if((i >= Ndim) || (k >= Nposition_type) ||
1476  (t >= position_ntstorage))
1477  {
1478  std::ostringstream error_message;
1479  //If it's the dimension
1480  if(i >= Ndim)
1481  {
1482  error_message << "Range Error: X coordinate " << i
1483  << " is not in the range (0,"
1484  << Ndim-1 << ")";
1485  }
1486  //If it's the position type
1487  if(k >= Nposition_type)
1488  {
1489  error_message << "Range Error: Position type " << k
1490  << " is not in the range (0,"
1491  << Nposition_type-1 << ")";
1492  }
1493  //If it's the time
1494  if(t >= position_ntstorage)
1495  {
1496  error_message << "Range Error: Position Time Value " << t
1497  << " is not in the range (0,"
1498  << position_ntstorage - 1 << ")";
1499  }
1500  //Throw the error
1501  throw OomphLibError(error_message.str(),
1502  OOMPH_CURRENT_FUNCTION,
1503  OOMPH_EXCEPTION_LOCATION);
1504  }
1505 }
1506 
1507 //========================================================================
1508 /// Static "Magic number" passed as independent_position when there is
1509 /// no independent position in the periodic node. For example, in a periodic
1510 /// mesh.
1511 //=======================================================================
1512 unsigned Node::No_independent_position=10;
1513 
1514 //========================================================================
1515 /// Default constructor.
1516 //========================================================================
1518  Position_time_stepper_pt(Data::Default_static_time_stepper_pt),
1519  Hanging_pt(0),
1520  Ndim(0), Nposition_type(0), Obsolete(false),
1521  Aux_node_update_fct_pt(0)
1522 {
1523 #ifdef LEAK_CHECK
1525 #endif
1526 }
1527 
1528 //========================================================================
1529 /// Steady Constructor, allocates storage for initial_n_value values
1530 /// at a node of spatial dimension NDim. nposition_type: # of coordinate
1531 /// types needed in the mapping between local and global coordinates
1532 /// (e.g. 1 for Lagrange-type elements; 2 for 1D Hermite elements; 4 for
1533 /// 2D Hermite elements, etc).
1534 //========================================================================
1535 Node::Node(const unsigned &n_dim,
1536  const unsigned &n_position_type,
1537  const unsigned &initial_n_value,
1538  const bool &allocate_x_position) :
1539  Data(initial_n_value),
1540  X_position(0),
1541  Position_time_stepper_pt(Data::Default_static_time_stepper_pt),
1542  Hanging_pt(0),
1543  Ndim(n_dim),
1544  Nposition_type(n_position_type), Obsolete(false), Aux_node_update_fct_pt(0)
1545 {
1546 #ifdef LEAK_CHECK
1548 #endif
1549 
1550  //Determine the total amount of storage required for position variables
1551  const unsigned n_storage = n_dim*n_position_type;
1552 
1553  //If we are in charge of the x coordinates (non-solid node)
1554  //the allocate storage
1555  if(allocate_x_position)
1556  {
1557  //Allocate the pointers to each coordinate and coordinate type
1558  X_position = new double*[n_storage];
1559 
1560  //Create one big array of positions
1561  double *x_positions = new double[n_storage];
1562 
1563  //Set pointers from the contiguous array
1564  for(unsigned j=0;j<n_storage;j++)
1565  {
1566  X_position[j] = &x_positions[j];
1567  //Initialise value to zero
1568  X_position[j][0] = 0.0;
1569  }
1570  }
1571 
1572 }
1573 
1574 //========================================================================
1575 /// Unsteady Constructor for a node of spatial dimension n_dim.
1576 /// Allocates storage
1577 /// for initial_n_value values with history values as required
1578 /// by timestepper. n_position_type: # of coordinate
1579 /// types needed in the mapping between local and global coordinates
1580 /// (e.g. 1 for Lagrange-type elements; 2 for 1D Hermite elements; 4 for
1581 /// 2D Hermite elements)
1582 //========================================================================
1583 Node::Node(TimeStepper* const &time_stepper_pt_,
1584  const unsigned &n_dim,
1585  const unsigned &n_position_type,
1586  const unsigned &initial_n_value,
1587  const bool &allocate_x_position)
1588  : Data(time_stepper_pt_,initial_n_value),
1589  X_position(0),
1590  Position_time_stepper_pt(time_stepper_pt_),
1591  Hanging_pt(0),
1592  Ndim(n_dim), Nposition_type(n_position_type), Obsolete(false), Aux_node_update_fct_pt(0)
1593 {
1594 #ifdef LEAK_CHECK
1596 #endif
1597 
1598  //Determine the total amount of storage required for position variables
1599  const unsigned n_storage = n_dim*n_position_type;
1600 
1601  //If we are allocating the storage (non-solid node)
1602  if(allocate_x_position)
1603  {
1604  //Amount of storage required for history values
1605  const unsigned n_tstorage = Position_time_stepper_pt->ntstorage();
1606 
1607  //Allocate the pointers to each coordinate and coordinate type
1608  X_position = new double*[n_storage];
1609 
1610  //Allocate the positions in one big array
1611  double *x_positions = new double[n_storage*n_tstorage];
1612 
1613  //Set the pointers to the contiguous memory
1614  for(unsigned j=0;j<n_storage;j++)
1615  {
1616  //Set the pointer from the bug array
1617  X_position[j] = &x_positions[j*n_tstorage];
1618  //Initialise all values to zero
1619  for(unsigned t=0;t<n_tstorage;t++) {X_position[j][t] = 0.0;}
1620  }
1621  }
1622 }
1623 
1624 
1625 //========================================================================
1626 /// Destructor to clean up the memory allocated for nodal position
1627 //========================================================================
1629 {
1630 #ifdef LEAK_CHECK
1632 #endif
1633 
1634  //Clean up memory allocated to hanging nodes
1635  if(Hanging_pt!=0)
1636  {
1637  //The number of hanging pointers is the number of values plus one
1638  const unsigned nhang = nvalue() + 1;
1639  for(unsigned ival=1;ival<nhang;ival++)
1640  {
1641  //If the ival-th HangInfo object is not the same as the geometrical
1642  //one, delete it
1643  if(Hanging_pt[ival]!=Hanging_pt[0]) {delete Hanging_pt[ival];}
1644  //Always NULL out the HangInfo pointer
1645  Hanging_pt[ival]=0;
1646  }
1647 
1648  //Delete the Geometrical HangInfo pointer
1649  delete Hanging_pt[0];
1650  Hanging_pt[0]=0;
1651 
1652  //Delete the Hanging_pt
1653  delete[] Hanging_pt;
1654  Hanging_pt=0;
1655  }
1656 
1657  //Free the memory allocated
1658 
1659  //If we did not allocate then the memory must have been freed by the
1660  //destructor of the object that did the allocating and X_position MUST
1661  //have been set back to zero
1662  //Test this and if so, we're done
1663  if(X_position==0) {return;}
1664 
1665  //If we're still here we must free our own memory which was allocated
1666  //in one block
1667  delete[] X_position[0];
1668 
1669  //Now delete the pointer
1670  delete[] X_position; X_position=0;
1671 }
1672 
1673 //================================================================
1674 /// Set a new position TimeStepper be resizing the appropriate storage.
1675 /// The current (zero) values will be unaffected, but all other entries
1676 /// will be set to zero.
1677 //================================================================
1679  const &position_time_stepper_pt,
1680  const bool &preserve_existing_data)
1681 {
1682  //If the timestepper is unchanged do nothing
1683  if(Position_time_stepper_pt==position_time_stepper_pt) {return;}
1684 
1685  //Find the amount of data to be preserved
1686  unsigned n_preserved_tstorage =1;
1687  if(preserve_existing_data)
1688  {n_preserved_tstorage = Position_time_stepper_pt->ntstorage();}
1689 
1690  //Set the new time stepper
1692 
1693  //Determine the total amount of storage required for position variables
1694  const unsigned n_storage = this->ndim()*this->nposition_type();
1695 
1696  //Amount of storage required for history values
1697  const unsigned n_tstorage = Position_time_stepper_pt->ntstorage();
1698 
1699  //Allocate all position data in one big array
1700  double *x_positions = new double[n_storage*n_tstorage];
1701 
1702  //If we have reduced the storage, reduce the size of preserved storage
1703  //to that of the new storage
1704  if(n_tstorage < n_preserved_tstorage) {n_preserved_tstorage = n_tstorage;}
1705 
1706  //Copy the old "preserved" positions into the new storage scheme
1707  for(unsigned j=0;j<n_storage;++j)
1708  {
1709  for(unsigned t=0;t<n_preserved_tstorage;t++)
1710  {
1711  x_positions[j*n_tstorage + t] = this->X_position[j][t];
1712  }
1713  }
1714 
1715  //Now delete the old position storage, which was allocated in one block
1716  delete[] X_position[0];
1717 
1718  //Set the pointers to the contiguous memory
1719  for(unsigned j=0;j<n_storage;j++)
1720  {
1721  //Set the pointer from the bug array
1722  X_position[j] = &x_positions[j*n_tstorage];
1723  //Initialise all new time storgae values to be zero
1724  for(unsigned t=n_preserved_tstorage;t<n_tstorage;t++)
1725  {X_position[j][t] = 0.0;}
1726  }
1727 }
1728 
1729 
1730 
1731 
1732 //================================================================
1733 /// Return the i-th component of nodal velocity: dx/dt
1734 //================================================================
1735 double Node::dx_dt(const unsigned &i) const
1736 {
1737  //Number of timsteps (past & present)
1738  const unsigned n_time = Position_time_stepper_pt->ntstorage();
1739 
1740  double dxdt=0.0;
1741 
1742  //If the timestepper is not steady
1744  {
1745  //Loop over the additional storage and add the appropriate contributions
1746  for(unsigned t=0;t<n_time;t++)
1747  {
1748  dxdt+=Position_time_stepper_pt->weight(1,t)*x(t,i);
1749  }
1750  }
1751 
1752  return dxdt;
1753 }
1754 
1755 //================================================================
1756 /// Return the i-th component of j-th derivative of nodal position:
1757 /// d^jx/dt^j.
1758 //================================================================
1759 double Node::dx_dt(const unsigned &j, const unsigned &i) const
1760 {
1761  // Number of timsteps (past & present)
1762  const unsigned n_time = Position_time_stepper_pt->ntstorage();
1763 
1764  double dxdt=0.0;
1765 
1766  //If the timestepper is not steady
1767  if ((!Position_time_stepper_pt->is_steady()) || (j==0))
1768  {
1769  //Loop over the additional storage and add the appropriate contributions
1770  for(unsigned t=0;t<n_time;t++)
1771  {
1772  dxdt+=Position_time_stepper_pt->weight(j,t)*x(t,i);
1773  }
1774  }
1775 
1776  return dxdt;
1777 }
1778 
1779 //================================================================
1780 /// \short i-th component of time derivative (velocity) of the
1781 /// generalised position, dx(k,i)/dt. `Type': k; Coordinate direction: i.
1782 //================================================================
1783 double Node::dx_gen_dt(const unsigned &k, const unsigned &i) const
1784 {
1785  // Number of timsteps (past & present)
1786  const unsigned n_time = Position_time_stepper_pt->ntstorage();
1787 
1788  double dxdt=0.0;
1789 
1790  //If the timestepper is not steady
1792  {
1793  //Loop over the additional time storage and add the appropriate
1794  //contributions
1795  for(unsigned t=0;t<n_time;t++)
1796  {
1797  dxdt+=Position_time_stepper_pt->weight(1,t)*x_gen(t,k,i);
1798  }
1799  }
1800 
1801  return dxdt;
1802 }
1803 
1804 //================================================================
1805 /// \short i-th component of j-th time derivative (velocity) of the
1806 /// generalised position, d^jx(k,i)/dt^j. `Type': k; Coordinate direction: i.
1807 //================================================================
1808 double Node::dx_gen_dt(const unsigned &j, const unsigned &k,
1809  const unsigned &i) const
1810 {
1811  // Number of timsteps (past & present)
1812  const unsigned n_time = Position_time_stepper_pt->ntstorage();
1813 
1814  double dxdt=0.0;
1815 
1816  //If the timestepper is not steady
1817  if ((!Position_time_stepper_pt->is_steady()) || (j==0))
1818  {
1819  //Loop over the additional storage and add the appropriate contributions
1820  for(unsigned t=0;t<n_time;t++)
1821  {
1822  dxdt+=Position_time_stepper_pt->weight(j,t)*x_gen(t,k,i);
1823  }
1824  }
1825 
1826  return dxdt;
1827 }
1828 
1829 
1830 
1831 //================================================================
1832 /// Copy all nodal data from specified Node object
1833 //================================================================
1834 void Node::copy(Node* orig_node_pt)
1835 {
1836 
1837  // Number of positional values
1838  const unsigned npos_storage = Ndim*Nposition_type;
1839 
1840  // Check # of values:
1841  const unsigned long npos_storage_orig
1842  = orig_node_pt->ndim()*orig_node_pt->nposition_type();
1843  if (npos_storage!=npos_storage_orig)
1844  {
1845  std::ostringstream error_stream;
1846  error_stream << "The allocated positional storage "
1847  << npos_storage << " is not the same as the original Node "
1848  << npos_storage_orig << std::endl;
1849 
1850  throw OomphLibError(error_stream.str(),
1851  OOMPH_CURRENT_FUNCTION,
1852  OOMPH_EXCEPTION_LOCATION);
1853  }
1854 
1855  // Number of time values (incl present)
1856  const unsigned n_time = Position_time_stepper_pt->ntstorage();
1857 
1858  // Check # of values:
1859  const unsigned long n_time_orig =
1860  orig_node_pt->position_time_stepper_pt()->ntstorage();
1861  if (n_time!=n_time_orig)
1862  {
1863  std::ostringstream error_stream;
1864  error_stream << "The number of positional time history values, "
1865  << n_time
1866  << " is not the same of those in the original node "
1867  << n_time_orig << std::endl;
1868 
1869  throw OomphLibError(error_stream.str(),
1870  OOMPH_CURRENT_FUNCTION,
1871  OOMPH_EXCEPTION_LOCATION);
1872  }
1873 
1874  // Copy fixed nodal positions
1875  for(unsigned t=0;t<n_time;t++)
1876  {
1877  for(unsigned j=0;j<npos_storage;j++)
1878  {
1879  X_position[j][t] = orig_node_pt->X_position[j][t];
1880  }
1881  }
1882 
1883  // Read associated data
1884  Data::copy(orig_node_pt);
1885 
1886 }
1887 
1888 
1889 //================================================================
1890 ///Dump nodal positions and associated data to file for restart
1891 //================================================================
1892 void Node::dump(std::ostream& dump_file) const
1893 {
1894  // Number of positional values
1895  const unsigned npos_storage = Ndim*Nposition_type;
1896  dump_file << npos_storage
1897  << " # number of fixed position variables" << std::endl;
1898 
1899  const unsigned Time_steps_range = Position_time_stepper_pt->ntstorage();
1900  dump_file << Time_steps_range
1901  << " # total number of doubles for time history (incl present)"
1902  << std::endl;
1903 
1904  for(unsigned t=0;t<Time_steps_range;t++)
1905  {
1906  for(unsigned j=0;j<npos_storage;j++)
1907  {
1908  dump_file << X_position[j][t] << std::endl;
1909  }
1910  }
1911 
1912  // Dump out data
1913  Data::dump(dump_file);
1914 }
1915 
1916 //================================================================
1917 ///Read nodal positions and associated data from file for restart
1918 //================================================================
1919 void Node::read(std::ifstream& restart_file)
1920 {
1921 
1922  std::string input_string;
1923 
1924  // Number of positional values
1925  const unsigned npos_storage = Ndim*Nposition_type;
1926 
1927  // Read line up to termination sign
1928  getline(restart_file,input_string,'#');
1929  // Ignore rest of line
1930  restart_file.ignore(80,'\n');
1931  // Check # of values:
1932  const unsigned long check_npos_storage=atoi(input_string.c_str());
1933  if (check_npos_storage!=npos_storage)
1934  {
1935  std::ostringstream error_stream;
1936  error_stream << "The allocated positional storage "
1937  << npos_storage <<
1938  " is not the same as that in the input file"
1939  << check_npos_storage << std::endl;
1940 
1941  throw OomphLibError(error_stream.str(),
1942  OOMPH_CURRENT_FUNCTION,
1943  OOMPH_EXCEPTION_LOCATION);
1944 
1945  }
1946 
1947  // Number of time values (incl present)
1948  const unsigned time_steps_range = Position_time_stepper_pt->ntstorage();
1949 
1950  // Read line up to termination sign
1951  getline(restart_file,input_string,'#');
1952  // Ignore rest of line
1953  restart_file.ignore(80,'\n');
1954  // Check # of values:
1955  const unsigned long check_time_steps_range=atoi(input_string.c_str());
1956  if (check_time_steps_range!=time_steps_range)
1957  {
1958  std::ostringstream error_stream;
1959  error_stream
1960  << "Number of positional history values in dump file is less "
1961  << "than the storage allocated in Node object: "
1962  << check_time_steps_range
1963  << " " << time_steps_range << std::endl;
1964 
1965  throw OomphLibError(error_stream.str(),
1966  OOMPH_CURRENT_FUNCTION,
1967  OOMPH_EXCEPTION_LOCATION);
1968  }
1969 
1970  // Read fixed nodal positions
1971  for(unsigned t=0;t<time_steps_range;t++)
1972  {
1973  for(unsigned j=0;j<npos_storage;j++)
1974  {
1975  // Read line
1976  getline(restart_file,input_string);
1977 
1978  // Transform to double
1979  X_position[j][t] = atof(input_string.c_str());
1980  }
1981  }
1982 
1983  // Read associated data
1984  Data::read(restart_file);
1985 }
1986 
1987 //=====================================================================
1988 /// Set the hanging data for the i-th value.
1989 /// If node is already hanging, simply overwrite the appropriate entry.
1990 /// If the node isn't hanging (because it might not be hanging
1991 /// geometrically), create the Vector of hanging pointers
1992 /// and make the other entries point to the node's geometric
1993 /// hanging data. Set hang_pt=0 to make entry explicitly non-hanging.
1994 /// Use Node::set_nonhanging() to unhang everything and clear up
1995 /// storage.
1996 //=====================================================================
1997 void Node::set_hanging_pt(HangInfo* const &hang_pt, const int &i)
1998 {
1999  //The number of hanging values is the number of stored values plus
2000  //one (geometry)
2001  unsigned n_hang = nvalue() + 1;
2002 
2003  //Has the vector of pointers to the HangInfo objects already been created?
2004  //If not create it
2005  if(Hanging_pt==0)
2006  {
2007  Hanging_pt = new HangInfo*[n_hang];
2008  //Initialise all entries to zero
2009  for(unsigned n=0;n<n_hang;n++) {Hanging_pt[n] = 0;}
2010  }
2011 
2012  //Geometric hanging data
2013  if(i==-1)
2014  {
2015  //Setup boolean array to find which pointers match the geometric pointer
2016  std::vector<bool> Same_as_geometric(n_hang,true);
2017 
2018  //Mark up any values that DON'T use the geometric hanging scheme
2019  for(unsigned n=1;n<n_hang;n++)
2020  {if(Hanging_pt[n] != Hanging_pt[0]) {Same_as_geometric[n] = false;}}
2021 
2022  //Remove the old geometric HangInfo
2023  delete Hanging_pt[0];
2024  //Assign the new geometric hanging data
2025  Hanging_pt[0] = hang_pt;
2026 
2027  //Constrain the geometric data (virtual function that is
2028  //overladed in solid nodes)
2029  if (hang_pt!=0)
2030  {
2032  }
2033  else
2034  {
2036  }
2037 
2038  //Loop over the entries again and update all pointers that pointed to
2039  //the geometric data
2040  for(unsigned n=1;n<n_hang;n++)
2041  {
2042  if(Same_as_geometric[n]==true)
2043  {
2044  Hanging_pt[n] = Hanging_pt[0];
2045 
2046  //In addition set the corresponding value to be constrained (hanging)
2047  if (Hanging_pt[n]!=0)
2048  {
2049  constrain(n-1);
2050  }
2051  else
2052  {
2053  unconstrain(n-1);
2054  }
2055  }
2056  }
2057 
2058  }
2059  //Value data
2060  else
2061  {
2062  //If the data is different from geometric, delete it
2063  if(Hanging_pt[i+1] != Hanging_pt[0])
2064  {
2065  delete Hanging_pt[i+1];
2066  Hanging_pt[i+1] = 0;
2067  }
2068 
2069  //Overwrite hanging data for the required value
2070  //Do not need to delete previous value, because it is assigned outside
2071  //the Node class
2072  Hanging_pt[i+1]=hang_pt;
2073 
2074  //In addition set the value to be constrained (hanging)
2075  if (hang_pt!=0)
2076  {
2077  constrain(i);
2078  }
2079  else
2080  {
2081  unconstrain(i);
2082  }
2083  }
2084 }
2085 
2086 //=====================================================================
2087 /// Resize the node to allow it to store n_value unknowns
2088 //=====================================================================
2089 void Node::resize(const unsigned &n_value)
2090 {
2091 
2092  // Old number of values
2093  unsigned old_nvalue=nvalue();
2094 
2095  // Now deal with the hanging Data (if any)
2096  HangInfo** backup_hanging_pt=0;
2097  if (Hanging_pt!=0)
2098  {
2099  // Backup
2100  backup_hanging_pt = new HangInfo*[old_nvalue+1];
2101 
2102  // Copy across existing ones
2103  for(unsigned i=0;i<old_nvalue+1;i++)
2104  {
2105  backup_hanging_pt[i]=Hanging_pt[i];
2106  }
2107 
2108  // Cleanup old one
2109  delete[] Hanging_pt;
2110  Hanging_pt=0;
2111  }
2112 
2113  // Call the resize function for the underlying Data object
2114  Data::resize(n_value);
2115 
2116 
2117  // Now deal with the hanging Data (if any)
2118  if (backup_hanging_pt!=0)
2119  {
2120  //The size of the new hanging point is the number of new values plus one.
2121  const unsigned n_hang = n_value+1;
2122  Hanging_pt = new HangInfo*[n_hang];
2123  //Initialise all entries to zero
2124  for(unsigned i=0;i<n_hang;i++) {Hanging_pt[i] = 0;}
2125 
2126  //Restore the saved data
2127  for(unsigned i=0;i<=old_nvalue;i++)
2128  {
2129  Hanging_pt[i] = backup_hanging_pt[i];
2130  }
2131 
2132  //Loop over the new values and set equal to the geometric hanging data
2133  for(unsigned i=old_nvalue+1;i<n_hang;i++)
2134  {
2135  Hanging_pt[i] = Hanging_pt[0];
2136  //and constrain if necessary
2137  if(Hanging_pt[i]!=0)
2138  {
2139  constrain(i-1);
2140  }
2141  else
2142  {
2143  unconstrain(i-1);
2144  }
2145  }
2146 
2147  //If necessary constrain
2148  //Positions
2149  //if(Hanging_pt[0] != 0) {constrain_positions();}
2150  //Values
2151  //for(unsigned i=0;i<n_value;i++)
2152  // {
2153  // if(Hanging_pt[i+1] != 0) {constrain(i);}
2154  // }
2155 
2156  // Loop over all values and geom hanging data
2157  /*for (int i=-1;i<int(old_nvalue);i++)
2158  {
2159  set_hanging_pt(backup_hanging_pt[i+1],i);
2160  }
2161 
2162  // By default use geometric hanging data for any new entries
2163  for (int i=int(old_nvalue);i<int(n_value);i++)
2164  {
2165  set_hanging_pt(backup_hanging_pt[0],i);
2166  }*/
2167 
2168  delete [] backup_hanging_pt;
2169  }
2170 
2171 }
2172 
2173 
2174 
2175 
2176 //=====================================================================
2177 /// Make the node periodic by copying values from node_pt.
2178 /// Broken virtual (only implemented in BoundaryNodes)
2179 //=====================================================================
2180 void Node::make_periodic(Node* const &node_pt)
2181 {
2182  throw OomphLibError("Only BoundaryNodes can be made periodic",
2183  OOMPH_CURRENT_FUNCTION,
2184  OOMPH_EXCEPTION_LOCATION);
2185 }
2186 
2187 //=====================================================================
2188 /// Make the nodes passed in periodic_nodes_pt periodic by copying values
2189 /// across from this node. At present all the positions will be assumed
2190 /// to be independent.
2191 /// Broken virtual (only implemented in BoundaryNodes)
2192 //=====================================================================
2193 void Node::make_periodic_nodes(const Vector<Node*> &periodic_nodes_pt)
2194 {
2195  throw OomphLibError("Only BoundaryNodes can make periodic nodes",
2196  OOMPH_CURRENT_FUNCTION,
2197  OOMPH_EXCEPTION_LOCATION);
2198 }
2199 
2200 
2201 //====================================================================
2202 /// Label node as non-hanging node by removing all hanging node data.
2203 //====================================================================
2205 {
2206  if(Hanging_pt!=0)
2207  {
2208  //Kill any additional hanging data for values
2209  const unsigned nhang = nvalue() + 1;
2210  for(unsigned ival=1;ival<nhang;ival++)
2211  {
2212  // Only kill it if it's different from the geometric hanging node data
2213  if (Hanging_pt[ival]!= Hanging_pt[0]) {delete Hanging_pt[ival];}
2214  //Always zero the entry
2215  Hanging_pt[ival] = 0;
2216 
2217  //Unconstrain any values that were constrained only because they
2218  //were hanging
2219  unconstrain(ival-1);
2220  }
2221 
2222  //Unconstrain the positions (virtual function that is overloaded for
2223  //solid nodes)
2225 
2226  //Kill the geometric hanging node data
2227  delete Hanging_pt[0];
2228  Hanging_pt[0]=0;
2229 
2230  //Kill the pointer to all hanging data
2231  delete[] Hanging_pt;
2232  Hanging_pt=0;
2233  }
2234 }
2235 
2236 
2237 
2238 //=======================================================================
2239 /// Interface for function to add the node to the mesh boundary b.
2240 /// Broken here in order to report run-time errors. Must be overloaded
2241 /// by all boundary nodes
2242 //=======================================================================
2243 void Node::add_to_boundary(const unsigned &b)
2244 {
2245  std::stringstream ss;
2246  ss << "Cannot add non BoundaryNode<NODE> to boundary " << b <<"\n";
2247  throw OomphLibError(ss.str(),
2248  OOMPH_CURRENT_FUNCTION,
2249  OOMPH_EXCEPTION_LOCATION);
2250 }
2251 
2252 
2253 //=======================================================================
2254 /// Interface for function to remove the node from the mesh boundary b.
2255 /// Broken here in order to report run-time erorrs. Must be overloaded
2256 /// by all boundary nodes
2257 //=======================================================================
2258 void Node::remove_from_boundary(const unsigned &b)
2259 {
2260  throw OomphLibError("Cannot remove non BoundaryNode<NODE> to boundary",
2261  OOMPH_CURRENT_FUNCTION,
2262  OOMPH_EXCEPTION_LOCATION);
2263 }
2264 
2265 
2266 //=========================================================================
2267 /// Interface to get the number of boundary coordinates on mesh boundary b.
2268 /// Broken here in order to provide run-time error reporting. Must
2269 /// be overloaded by all boundary nodes.
2270 //=========================================================================
2271 unsigned Node::ncoordinates_on_boundary(const unsigned &b)
2272 {
2273  throw OomphLibError("Non-boundary Node cannot have boundary coordinates",
2274  OOMPH_CURRENT_FUNCTION,
2275  OOMPH_EXCEPTION_LOCATION);
2276  // dummy return
2277  return 0;
2278 }
2279 
2280 
2281 //=========================================================================
2282 /// Interface for function to get the k-th generalised boundary coordinate
2283 /// of the node on boundary b. Broken here in order to
2284 /// provide run-time error reporting. Must be overloaded by all boundary
2285 /// nodes.
2286 //=========================================================================
2287 void Node::get_coordinates_on_boundary(const unsigned &b, const unsigned& k,
2288  Vector<double> &boundary_zeta)
2289 {
2290  throw OomphLibError("Non-boundary Node cannot have boundary coordinates",
2291  OOMPH_CURRENT_FUNCTION,
2292  OOMPH_EXCEPTION_LOCATION);
2293 }
2294 
2295 
2296 //=========================================================================
2297 /// Interface for function to set the k-th generalised boundary coordinate
2298 /// of the node on boundary b. Broken here to provide
2299 /// run-time error reports. Must be overloaded by all boundary nodes.
2300 //=========================================================================
2301 void Node::set_coordinates_on_boundary(const unsigned &b, const unsigned& k,
2302  const Vector<double> &boundary_zeta)
2303 {
2304  throw OomphLibError("Non-boundary Node cannot have boundary coordinates",
2305  OOMPH_CURRENT_FUNCTION,
2306  OOMPH_EXCEPTION_LOCATION);
2307 }
2308 
2309 
2310 //=================================================================
2311 /// Return i-th value (free or pinned) at this node
2312 /// either directly or via hanging node representation.
2313 //================================================================
2314 double Node::value(const unsigned &i) const
2315 {
2316  //If value is not hanging, just return the underlying value
2317  if(!is_hanging(i)) {return raw_value(i);}
2318  // Hanging node: Use hanging node representation
2319  else
2320  {
2321  // Initialise
2322  double sum=0.0;
2323  // Add contribution from master nodes
2324  const unsigned n_master = hanging_pt(i)->nmaster();
2325  for(unsigned m=0;m<n_master;m++)
2326  {
2327  //A master node cannot be hanging by definition.
2328  //so we get the raw value to avoid an unnecessary it
2329  sum += hanging_pt(i)->master_node_pt(m)->raw_value(i)*
2330  hanging_pt(i)->master_weight(m);
2331  }
2332  return sum;
2333  }
2334 }
2335 
2336 //=================================================================
2337 /// Return i-th value (free or pinned) at this node at time level t
2338 /// either directly or via hanging node representation.
2339 //================================================================
2340 double Node::value(const unsigned &t, const unsigned &i) const
2341 {
2342  //If value is not hanging, just return the raw value
2343  if(!is_hanging(i)) {return raw_value(t,i);}
2344  // Hanging node: Use hanging node representation
2345  else
2346  {
2347  // Initialise
2348  double sum=0.0;
2349 
2350  // Add contribution from master nodes
2351  const unsigned n_master=hanging_pt(i)->nmaster();
2352  for(unsigned m=0;m<n_master;m++)
2353  {
2354  //Get the raw nodal values at each master to avoid un-necessary ifs
2355  sum += hanging_pt(i)->master_node_pt(m)->raw_value(t,i)*
2356  hanging_pt(i)->master_weight(m);
2357  }
2358  return sum;
2359  }
2360 }
2361 
2362 //==================================================================
2363 /// Compute Vector of values (dofs or pinned) at this Data object
2364 /// either directly or via hanging node representation.
2365 //==================================================================
2366 void Node::value(Vector<double>& values) const
2367 {
2368  //Loop over all the values
2369  const unsigned n_value = nvalue();
2370  for(unsigned i=0;i<n_value;i++)
2371  {
2372  //Set the value, using the hanging node representation if necessary
2373  values[i] = value(i);
2374  }
2375 }
2376 
2377 //==================================================================
2378 /// Compute Vector of values (dofs or pinned) at this node
2379 /// at time level t (t=0: present; t>0: previous)
2380 /// either directly or via hanging node representation.
2381 //==================================================================
2382 void Node::value(const unsigned& t, Vector<double>& values) const
2383 {
2384  //Loop over all the values
2385  const unsigned n_value = nvalue();
2386  for(unsigned i=0;i<n_value;i++)
2387  {
2388  //Set the value at the time-level t, using the hanging node representation
2389  //if necessary
2390  values[i] = value(t,i);
2391  }
2392 }
2393 
2394 
2395 //============================================================
2396 /// Compute Vector of nodal positions
2397 /// either directly or via hanging node representation
2398 //===========================================================
2400 {
2401  //Assign all positions using hanging node representation where necessary
2402  const unsigned n_dim = ndim();
2403  for(unsigned i=0;i<n_dim;i++) {pos[i] = position(i);}
2404 }
2405 
2406 //===============================================================
2407 /// Compute Vector of nodal position at timestep t
2408 /// (t=0: current; t>0: previous timestep),
2409 /// either directly or via hanging node representation.
2410 //==============================================================
2411 void Node::position(const unsigned &t, Vector<double>& pos) const
2412 {
2413  //Assign all positions, using hanging node representation where necessary
2414  const unsigned n_dim = ndim();
2415  for(unsigned i=0;i<n_dim;i++) {pos[i] = position(t,i);}
2416 }
2417 
2418 //=======================================================================
2419 /// Return i-th nodal coordinate
2420 /// either directly or via hanging node representation.
2421 //======================================================================
2422 double Node::position(const unsigned &i) const
2423 {
2424  double posn=0.0;
2425 
2426  // Non-hanging node: just return value
2427  if (!is_hanging()) {posn = x(i);}
2428  // Hanging node: Use hanging node representation
2429  else
2430  {
2431  // Initialise
2432  double interpolated_position=0.0;
2433 
2434  // Add contribution from master nodes
2435  const unsigned n_master=hanging_pt()->nmaster();
2436  for (unsigned m=0;m<n_master;m++)
2437  {
2438  interpolated_position += hanging_pt()->master_node_pt(m)->x(i)*
2439  hanging_pt()->master_weight(m);
2440  }
2441  posn=interpolated_position;
2442  }
2443  return posn;
2444 }
2445 
2446 //================================================================
2447 /// Return i-th nodal coordinate at time level t
2448 /// (t=0: current; t>0: previous time level),
2449 /// either directly or via hanging node representation.
2450 //================================================================
2451 double Node::position(const unsigned &t, const unsigned &i) const
2452 {
2453  double posn=0.0;
2454 
2455  // Non-hanging node: just return value
2456  if(!is_hanging()) {posn = x(t,i);}
2457  // Hanging node: Use hanging node representation
2458  else
2459  {
2460  // Initialise
2461  double interpolated_position=0.0;
2462 
2463  // Add contribution from master nodes
2464  const unsigned n_master=hanging_pt()->nmaster();
2465  for (unsigned m=0;m<n_master;m++)
2466  {
2467  interpolated_position+=
2468  hanging_pt()->master_node_pt(m)->x(t,i)*
2469  hanging_pt()->master_weight(m);
2470  }
2471  posn=interpolated_position;
2472  }
2473 
2474  return posn;
2475 }
2476 
2477 //=======================================================================
2478 /// Return generalised nodal coordinate
2479 /// either directly or via hanging node representation.
2480 //======================================================================
2481 double Node::position_gen(const unsigned &k, const unsigned &i) const
2482 {
2483  double posn=0.0;
2484 
2485  // Non-hanging node: just return value
2486  if(!is_hanging()) {posn=x_gen(k,i);}
2487  // Hanging node: Use hanging node representation
2488  else
2489  {
2490  // Initialise
2491  double interpolated_position=0.0;
2492 
2493  // Add contribution from master nodes
2494  const unsigned n_master=hanging_pt()->nmaster();
2495  for (unsigned m=0;m<n_master;m++)
2496  {
2497  interpolated_position+=
2498  hanging_pt()->master_node_pt(m)->x_gen(k,i)*
2499  hanging_pt()->master_weight(m);
2500  }
2501  posn=interpolated_position;
2502  }
2503  return posn;
2504 }
2505 
2506 //================================================================
2507 /// Return generalised nodal coordinate at time level t
2508 /// (t=0: current; t>0: previous time level),
2509 /// either directly or via hanging node representation.
2510 //================================================================
2511 double Node::position_gen(const unsigned &t, const unsigned &k,
2512  const unsigned &i) const
2513 {
2514  double posn=0.0;
2515 
2516  // Non-hanging node: just return value
2517  if(!is_hanging()) {posn=x_gen(t,k,i);}
2518  // Hanging node: Use hanging node representation
2519  else
2520  {
2521  // Initialise
2522  double interpolated_position=0.0;
2523 
2524  // Add contribution from master nodes
2525  const unsigned n_master=hanging_pt()->nmaster();
2526  for (unsigned m=0;m<n_master;m++)
2527  {
2528  interpolated_position+=
2529  hanging_pt()->master_node_pt(m)->x_gen(t,k,i)*
2530  hanging_pt()->master_weight(m);
2531  }
2532  posn=interpolated_position;
2533  }
2534 
2535  return posn;
2536 }
2537 
2538 //================================================================
2539 /// Return the i-th component of nodal velocity: dx/dt
2540 //// Use the hanging node representation if required.
2541 //================================================================
2542 double Node::dposition_dt(const unsigned &i) const
2543 {
2544  //Number of timsteps (past & present)
2545  const unsigned n_time = Position_time_stepper_pt->ntstorage();
2546 
2547  double dxdt=0.0;
2548 
2549  //If the timestepper is not steady
2551  {
2552  //Loop over the additional storage and add the appropriate contributions
2553  for(unsigned t=0;t<n_time;t++)
2554  {
2556  }
2557  }
2558 
2559  return dxdt;
2560 }
2561 
2562 //================================================================
2563 /// Return the i-th component of j-th derivative of nodal position:
2564 /// d^jx/dt^j. Use the hanging node representation.
2565 //================================================================
2566 double Node::dposition_dt(const unsigned &j, const unsigned &i) const
2567 {
2568  // Number of timsteps (past & present)
2569  const unsigned n_time = Position_time_stepper_pt->ntstorage();
2570 
2571  double dxdt=0.0;
2572 
2573  //If the timestepper is not steady
2574  if ((!Position_time_stepper_pt->is_steady()) || (j==0))
2575  {
2576  //Loop over the additional storage and add the appropriate contributions
2577  for(unsigned t=0;t<n_time;t++)
2578  {
2580  }
2581  }
2582 
2583  return dxdt;
2584 }
2585 
2586 //================================================================
2587 /// \short i-th component of time derivative (velocity) of the
2588 /// generalised position, dx(k,i)/dt. `Type': k; Coordinate direction: i.
2589 /// Use the hanging node representation
2590 //================================================================
2591 double Node::dposition_gen_dt(const unsigned &k, const unsigned &i) const
2592 {
2593  // Number of timsteps (past & present)
2594  const unsigned n_time = Position_time_stepper_pt->ntstorage();
2595 
2596  double dxdt=0.0;
2597 
2598  //If the timestepper is not steady
2600  {
2601  //Loop over the additional time storage and add the appropriate
2602  //contributions
2603  for(unsigned t=0;t<n_time;t++)
2604  {
2606  }
2607  }
2608 
2609  return dxdt;
2610 }
2611 
2612 //================================================================
2613 /// \short i-th component of j-th time derivative (velocity) of the
2614 /// generalised position, d^jx(k,i)/dt^j. `Type': k; Coordinate direction: i.
2615 /// Use the hanging node representation.
2616 //================================================================
2617 double Node::dposition_gen_dt(const unsigned &j, const unsigned &k,
2618  const unsigned &i) const
2619 {
2620  // Number of timsteps (past & present)
2621  const unsigned n_time = Position_time_stepper_pt->ntstorage();
2622 
2623  double dxdt=0.0;
2624 
2625  //If the timestepper is not steady
2626  if ((!Position_time_stepper_pt->is_steady()) || (j==0))
2627  {
2628  //Loop over the additional storage and add the appropriate contributions
2629  for(unsigned t=0;t<n_time;t++)
2630  {
2632  }
2633  }
2634 
2635  return dxdt;
2636 }
2637 
2638 
2639 //========================================================================
2640 /// Output nodal coordinates
2641 //========================================================================
2642 void Node::output(std::ostream &outfile)
2643 {
2644  //Loop over the number of dimensions of the node
2645  const unsigned n_dim = this->ndim();
2646  for (unsigned i=0;i<n_dim;i++) {outfile << x(i) << " ";}
2647  outfile << std::endl;
2648 }
2649 
2650 
2651 
2652 #ifdef OOMPH_HAS_MPI
2653 //==================================================================
2654 /// Add position data and time-history values to the vector after the
2655 /// addition of the "standard" data stored at the node
2656 //==================================================================
2658 {
2659  //Firstly add the value data to the vector
2660  Data::add_values_to_vector(vector_of_values);
2661 
2662  //Now add the additional position data to the vector
2663 
2664  //Find the number of stored time data for the position
2665  const unsigned n_tstorage = this->position_time_stepper_pt()->ntstorage();
2666 
2667  //Find the total amount of storage required for the position variables
2668  const unsigned n_storage = this->ndim()*this->nposition_type();
2669 
2670  //Resize the vector to accommodate the new data
2671  const unsigned n_current_value = vector_of_values.size();
2672 
2673 #ifdef PARANOID
2674  unsigned n_debug=2;
2675 #else
2676  unsigned n_debug=0;
2677 #endif
2678 
2679  vector_of_values.resize(n_current_value + n_tstorage*n_storage+n_debug);
2680 
2681  //Now add the data to the vector
2682  unsigned index = n_current_value;
2683 
2684 #ifdef PARANOID
2685  vector_of_values[index++]=n_storage;
2686  vector_of_values[index++]=n_tstorage;
2687 #endif
2688 
2689 
2690  //Pointer to the first entry in the data array
2691  double* data_pt = X_position[0];
2692 
2693  //Loop over values
2694  for(unsigned i=0;i<n_storage;i++)
2695  {
2696  //Loop over time histories
2697  for(unsigned t=0;t<n_tstorage;t++)
2698  {
2699  //Add the position data to the vector
2700  vector_of_values[index] = *data_pt;
2701  //Increment the counter and the pointer
2702  ++index;
2703  ++data_pt;
2704  }
2705  }
2706 }
2707 
2708 //==================================================================
2709 /// Read the position data and its time histories from the vector
2710 /// after reading the "standard" data.
2711 //==================================================================
2712 void Node::read_values_from_vector(const Vector<double> &vector_of_values,
2713  unsigned &index)
2714 {
2715  //Read the standard nodal data
2716  Data::read_values_from_vector(vector_of_values,index);
2717 
2718  //Now read the additional position data to the vector
2719 
2720  //Find the number of stored time data for the position
2721  const unsigned n_tstorage = this->position_time_stepper_pt()->ntstorage();
2722 
2723 //Find the total amount of storage required for the position variables
2724  const unsigned n_storage = this->ndim()*this->nposition_type();
2725 
2726 #ifdef PARANOID
2727  unsigned orig_n_storage=unsigned(vector_of_values[index++]);
2728  unsigned orig_n_tstorage=unsigned(vector_of_values[index++]);
2729  if ((orig_n_tstorage!=n_tstorage)||(orig_n_storage!=n_storage))
2730  {
2731  std::ostringstream error_stream;
2732  error_stream << "Non-matching number of values:\n"
2733  << "sent and local n_tstorage: " << orig_n_tstorage << " "
2734  << n_tstorage << std::endl
2735  << "sent and local n_storage: " << orig_n_storage << " "
2736  << n_storage << std::endl;
2737  throw OomphLibError(error_stream.str(),
2738  OOMPH_CURRENT_FUNCTION,
2739  OOMPH_EXCEPTION_LOCATION);
2740  }
2741 #endif
2742 
2743  //Pointer to the first entry in the data array
2744  double* data_pt = X_position[0];
2745 
2746  //Loop over values
2747  for(unsigned i=0;i<n_storage;i++)
2748  {
2749  //Loop over time histories
2750  for(unsigned t=0;t<n_tstorage;t++)
2751  {
2752  //Read the position data from the vector
2753  *data_pt = vector_of_values[index];
2754  //Increment the counter and the pointer
2755  ++index;
2756  ++data_pt;
2757  }
2758  }
2759 }
2760 
2761 #endif
2762 
2763 
2764 
2765 /////////////////////////////////////////////////////////////////////////
2766 /////////////////////////////////////////////////////////////////////////
2767 //Functions for the BoundaryNodeBase class
2768 /////////////////////////////////////////////////////////////////////////
2769 /////////////////////////////////////////////////////////////////////////
2770 
2771 //==================================================================
2772 /// \short Helper function that is used to turn BoundaryNodes into
2773 /// periodic boundary nodes by setting the data values of the nodes
2774 /// in the vector periodic_copies_pt to be the same as those
2775 /// in copied_node_pt. This function should be used when making doubly periodic
2776 /// sets of nodes.
2777 //==================================================================
2779  Node* const &copied_node_pt,
2780  Vector<Node*> const &periodic_copies_pt)
2781 {
2782  //Don't allow copying if the original or periodic nodes are already
2783  //periodic
2784  bool already_a_copy = false;
2785  already_a_copy |= copied_node_pt->is_a_copy();
2786  const unsigned n_periodic = periodic_copies_pt.size();
2787  for(unsigned n=0;n<n_periodic;n++)
2788  {
2789  already_a_copy |= periodic_copies_pt[n]->is_a_copy();
2790  }
2791 
2792  //If we have a copy bail
2793  if(already_a_copy)
2794  {
2795  std::ostringstream error_stream;
2796  error_stream <<
2797  "The nodes you are trying to make periodic are already periodic\n"
2798  <<
2799  "Or you are trying to make a copy of another already periodic node\n";
2800  error_stream << "Please copy the original data if you can\n";
2801  throw OomphLibError(error_stream.str(),
2802  OOMPH_CURRENT_FUNCTION,
2803  OOMPH_EXCEPTION_LOCATION);
2804 
2805  }
2806 
2807  //Now we simply delete and copy over for each node
2808  for(unsigned n=0;n<n_periodic;n++)
2809  {
2810  //Local cache of the node
2811  Node* const nod_pt = periodic_copies_pt[n];
2812  //Miss out the node itself if it's in the list
2813  if(nod_pt != copied_node_pt)
2814  {
2815  //Delete the storage allocated in the copy
2816  nod_pt->delete_value_storage();
2817  //Now set the Value and Equation number pointers to be the same
2818  nod_pt->Value = copied_node_pt->Value;
2819  nod_pt->Eqn_number = copied_node_pt->Eqn_number;
2820 
2821  //Set the copied node pointer in the copy
2822  BoundaryNodeBase* cast_nod_pt = dynamic_cast<BoundaryNodeBase*>(nod_pt);
2823  cast_nod_pt->Copied_node_pt = copied_node_pt;
2824  //Inform the node that it has been copied
2825  copied_node_pt->add_copy(nod_pt);
2826  }
2827  }
2828 
2829 }
2830 
2831 
2832 //====================================================================
2833 /// Helper function that is used to turn BoundaryNodes into
2834 /// peridic boundary nodes by setting the data values of
2835 /// copy_of_node_pt to those of copied_node_pt.
2836 //=====================================================================
2838  Node* const &copied_node_pt)
2839 {
2840  // Don't allow this node to already be a copy (you should clear it's
2841  // copied status first).
2842  if(node_pt->is_a_copy())
2843  {
2844  std::ostringstream error_stream;
2845  error_stream <<
2846  "The node you are trying to make into a periodic copy is already a copy\n.";
2847  throw OomphLibError(error_stream.str(),
2848  OOMPH_CURRENT_FUNCTION,
2849  OOMPH_EXCEPTION_LOCATION);
2850  }
2851 
2852  // If the node to be copied from is a copy then copy it's "copied_node_pt"
2853  // instead.
2854  if(copied_node_pt->is_a_copy())
2855  {
2856  make_node_periodic(node_pt, copied_node_pt->copied_node_pt());
2857  }
2858 
2859  // Otherwise just do it
2860  else
2861  {
2862  //Set the copied node pointer
2863  Copied_node_pt = copied_node_pt;
2864 
2865  //First copy the data values
2866  //Delete the storage allocated in the copy
2867  node_pt->delete_value_storage();
2868  //Now set the Value and Equation number pointers to be the same
2869  node_pt->Value = copied_node_pt->Value;
2870  node_pt->Eqn_number = copied_node_pt->Eqn_number;
2871 
2872  //Inform the node that it has been copied
2873  copied_node_pt->add_copy(node_pt);
2874  }
2875 }
2876 
2877 //======================================================================
2878 /// Destructor to clean up any memory that might have been allocated.
2879 //=======================================================================
2881 {
2882  //Delete the set of boundaries on which the Node lies
2883  delete Boundaries_pt;
2884  Boundaries_pt=0;
2885 
2886  //If the Boundary coordinates have been set then delete them
2887  if(Boundary_coordinates_pt != 0)
2888  {
2889  //Loop over the boundary coordinate entries and delete them
2890  for(std::map<unsigned,DenseMatrix<double> *>::iterator it
2891  = Boundary_coordinates_pt->begin();
2892  it!=Boundary_coordinates_pt->end();++it)
2893  {
2894  //Delete the vectors that have been allocated for the storage
2895  //of the boundary coordinates
2896  delete it->second;
2897  }
2898 
2899  //Now delete the Boundary coordinates map itself
2900  delete Boundary_coordinates_pt;
2901  //Set the pointer to null to be on the safe side
2903  }
2904 
2905  //Delete the map of face element's first value
2908 }
2909 
2910 
2911 
2912 
2913 
2914 
2915 
2916 //=======================================================================
2917 /// Add the node to the mesh boundary b
2918 //=======================================================================
2919 void BoundaryNodeBase::add_to_boundary(const unsigned &b)
2920 {
2921  //If there is not storage then create storage and set the entry
2922  //to be the boundary b
2923  if(Boundaries_pt==0) {Boundaries_pt = new std::set<unsigned>;}
2924 
2925 // //If the boundary is already stored in the node issue a warning
2926 // if(find(Boundaries_pt->begin(),Boundaries_pt->end(),b) !=
2927 // Boundaries_pt->end())
2928 // {
2929 // // MH: who cares?
2930 // // oomph_info << std::endl << "============================================"
2931 // // << std::endl << "Warning in Node::add_to_boundary() "
2932 // // << std::endl << "Node is already marked as being on boundary " << b
2933 // // << std::endl << "============================================"
2934 // // << std::endl << std::endl;
2935 // }
2936 // else
2937 // {
2938 
2939  Boundaries_pt->insert(b);
2940 
2941 // }
2942 }
2943 
2944 
2945 //=======================================================================
2946 /// Remove the node from the mesh boundary b
2947 //=======================================================================
2949 {
2950 #ifdef PARANOID
2951  if(is_on_boundary(b)==false)
2952  {
2953  std::ostringstream error_stream;
2954  error_stream << "Node is not on boundary " << b << std::endl;
2955 
2956  throw OomphLibError(error_stream.str(),
2957  OOMPH_CURRENT_FUNCTION,
2958  OOMPH_EXCEPTION_LOCATION);
2959  }
2960 #endif
2961 
2962  //Remove the boundary from the set
2963  Boundaries_pt->erase(b);
2964 
2965  //Need to delete the equivalent entry in the Boundary coordinate
2966  //map, if the storage has actually been allocated
2968  {
2969  //Delete the vector storage that has been allocated
2970  delete (*Boundary_coordinates_pt)[b];
2971  //Delete the entry in the map
2972  Boundary_coordinates_pt->erase(b);
2973  }
2974 
2975  //If all entries have been removed, delete the storage
2976  if(Boundaries_pt->size()==0)
2977  {
2978  delete Boundaries_pt;
2979  Boundaries_pt=0;
2980  }
2981 }
2982 
2983 //========================================================================
2984 /// Test whether the node lies on the mesh boundary b
2985 //========================================================================
2986 bool BoundaryNodeBase::is_on_boundary(const unsigned &b) const
2987 {
2988  //If the node lies on any boundary
2989  if(Boundaries_pt!=0)
2990  {
2991  if(find(Boundaries_pt->begin(),Boundaries_pt->end(),b)
2992  != Boundaries_pt->end()) {return true;}
2993  }
2994 
2995  //If we haven't returned yet, then the node does not lie on the boundary
2996  return false;
2997 }
2998 
2999 
3000 //=========================================================================
3001 /// Get the number of boundary coordinates on mesh boundary b
3002 //=========================================================================
3004 {
3005  //Check that the node lies on a boundary
3006 #ifdef PARANOID
3007  if(Boundaries_pt==0)
3008  {
3009  throw OomphLibError("Node does not lie on any boundary",
3010  OOMPH_CURRENT_FUNCTION,
3011  OOMPH_EXCEPTION_LOCATION);
3012  }
3013 
3014 
3015  //Does the node lie on the mesh boundary b
3016  if(!is_on_boundary(b))
3017  {
3018  std::ostringstream error_stream;
3019  error_stream << "Node is not on boundary " << b << std::endl;
3020 
3021  throw OomphLibError(error_stream.str(),
3022  OOMPH_CURRENT_FUNCTION,
3023  OOMPH_EXCEPTION_LOCATION);
3024  }
3025 
3026  //Check that the boundary coordinates have been set
3027  if(Boundary_coordinates_pt == 0)
3028  {
3029  std::ostringstream error_stream;
3030  error_stream << "Boundary coordinates have not been set\n"
3031  << "[Note: In refineable problems, the boundary coordinates\n"
3032  << " will only be interpolated to newly created nodes\n"
3033  << " if Mesh::Boundary_coordinate_exists[...] has been\n"
3034  << " set to true!]\n";
3035  throw OomphLibError(error_stream.str(),
3036  OOMPH_CURRENT_FUNCTION,
3037  OOMPH_EXCEPTION_LOCATION);
3038  }
3039 #endif
3040 
3041  //Find out how may coordinates there are from the map
3042  return (*Boundary_coordinates_pt)[b]->nrow();
3043 
3044 }
3045 
3046 
3047 
3048 //=========================================================================
3049 /// Given the mesh boundary b, return the k-th generalised boundary
3050 /// coordinates of the node in the vector boundary_zeta
3051 //=========================================================================
3052 void BoundaryNodeBase::
3053 get_coordinates_on_boundary(const unsigned &b, const unsigned& k,
3054  Vector<double> &boundary_zeta)
3055 {
3056  //Check that the node lies on a boundary
3057 #ifdef PARANOID
3058  if(Boundaries_pt==0)
3059  {
3060  throw OomphLibError("Node does not lie on any boundary",
3061  OOMPH_CURRENT_FUNCTION,
3062  OOMPH_EXCEPTION_LOCATION);
3063  }
3064 #endif
3065 
3066  //Does the node lie on the mesh boundary b
3067  if(!is_on_boundary(b))
3068  {
3069  std::ostringstream error_stream;
3070  error_stream << "Node is not on boundary " << b << std::endl;
3071 
3072  throw OomphLibError(error_stream.str(),
3073  OOMPH_CURRENT_FUNCTION,
3074  OOMPH_EXCEPTION_LOCATION);
3075  }
3076 
3077 
3078 #ifdef PARANOID
3079  //Check that the boundary coordinates have been set
3080  if(Boundary_coordinates_pt == 0)
3081  {
3082  std::ostringstream error_stream;
3083  error_stream << "Boundary coordinates have not been set\n"
3084  << "[Note: In refineable problems, the boundary coordinates\n"
3085  << " will only be interpolated to newly created nodes\n"
3086  << " if Mesh::Boundary_coordinate_exists[...] has been\n"
3087  << " set to true!]\n";
3088  throw OomphLibError(error_stream.str(),
3089  OOMPH_CURRENT_FUNCTION,
3090  OOMPH_EXCEPTION_LOCATION);
3091  }
3092 #endif
3093 
3094 
3095  //Find out how may coordinates there are from the map
3096  const unsigned nboundary_coord = (*Boundary_coordinates_pt)[b]->nrow();
3097 #ifdef PARANOID
3098  if(nboundary_coord != boundary_zeta.size())
3099  {
3100  std::ostringstream error_stream;
3101  error_stream
3102  << "Wrong number of coordinates in the vector boundary_zeta"
3103  << std::endl << "There are " << nboundary_coord
3104  << " boundary coordinates"
3105  << std::endl << "But bounday_zeta() has size " << boundary_zeta.size()
3106  << std::endl;
3107 
3108  throw OomphLibError(error_stream.str(),
3109  OOMPH_CURRENT_FUNCTION,
3110  OOMPH_EXCEPTION_LOCATION);
3111  }
3112 #endif
3113 
3114  //Loop over and assign the coordinates
3115  for(unsigned i=0;i<nboundary_coord;i++)
3116  {boundary_zeta[i] = (*(*Boundary_coordinates_pt)[b])(i,k);}
3117 }
3118 
3119 
3120 //=========================================================================
3121 /// Given the mesh boundary b, set the k-th generalised boundary
3122 /// coordinates of the node from the vector boundary_zeta
3123 //=========================================================================
3124 void BoundaryNodeBase::
3125 set_coordinates_on_boundary(const unsigned &b, const unsigned& k,
3126  const Vector<double> &boundary_zeta)
3127 {
3128  //Check that the node lies on a boundary
3129 #ifdef PARANOID
3130  if(Boundaries_pt==0)
3131  {
3132  throw OomphLibError("Node does not lie on any boundary",
3133  OOMPH_CURRENT_FUNCTION,
3134  OOMPH_EXCEPTION_LOCATION);
3135  }
3136 #endif
3137 
3138  //Does the node lie on the mesh boundary b
3139  if(!is_on_boundary(b))
3140  {
3141  std::ostringstream error_stream;
3142  error_stream << "Node is not on boundary " << b << std::endl;
3143 
3144  throw OomphLibError(error_stream.str(),
3145  OOMPH_CURRENT_FUNCTION,
3146  OOMPH_EXCEPTION_LOCATION);
3147  }
3148 
3149  //If the storage has not been assigned, then assign it
3150  if(Boundary_coordinates_pt == 0)
3151  {
3152  Boundary_coordinates_pt = new std::map<unsigned,DenseMatrix<double> *>;
3153  }
3154 
3155 
3156  //Find the number of boundary coordinates
3157  const unsigned nboundary_coord = boundary_zeta.size();
3158 
3159  //Allocate the vector for the boundary coordinates, if we haven't already
3160  if((*Boundary_coordinates_pt)[b] == 0)
3161  {
3162  // Need k+1 columns initially
3163  (*Boundary_coordinates_pt)[b] =
3164  new DenseMatrix<double>(nboundary_coord,k+1);
3165  }
3166  //Otherwise resize it, in case the number of boundary coordinates
3167  //or the number of types has changed
3168  else
3169  {
3170  // Adjust number of boundary coordinates -- retain number of types
3171  unsigned ncol=(*Boundary_coordinates_pt)[b]->ncol();
3172  {
3173  (*Boundary_coordinates_pt)[b]->resize(nboundary_coord,ncol);
3174  }
3175 
3176  // Resize number of types if required
3177  if ((k+1)>(*Boundary_coordinates_pt)[b]->ncol())
3178  {
3179  (*Boundary_coordinates_pt)[b]->resize(nboundary_coord,k+1);
3180  }
3181  }
3182 
3183  //Loop over and assign the coordinates
3184  for(unsigned i=0;i<nboundary_coord;i++)
3185  {(*(*Boundary_coordinates_pt)[b])(i,k) = boundary_zeta[i];}
3186 }
3187 
3188 
3189 
3190 
3191 ///////////////////////////////////////////////////////////////////////////
3192 ///////////////////////////////////////////////////////////////////////////
3193 //Functions for the SolidNode class
3194 ///////////////////////////////////////////////////////////////////////////
3195 ///////////////////////////////////////////////////////////////////////////
3196 
3197 
3198 //=================================================================
3199 /// Private function to check that the argument is within the
3200 /// range of stored Lagrangain coordinates and position types.
3201 //=================================================================
3202 void SolidNode::xi_gen_range_check(const unsigned &k, const unsigned &i) const
3203 {
3204  //If either the coordinate or type are out of range
3205  if((i >= Nlagrangian) || (k >= Nlagrangian_type))
3206  {
3207  std::ostringstream error_message;
3208  //If it's the lagrangian coordinate
3209  if(i >= Nlagrangian)
3210  {
3211  error_message << "Range Error: Xi coordinate " << i
3212  << " is not in the range (0,"
3213  << Nlagrangian-1 << ")";
3214  }
3215  //If it's the position type
3216  if(k >= Nlagrangian_type)
3217  {
3218  error_message << "Range Error: Lagrangian type " << k
3219  << " is not in the range (0,"
3220  << Nlagrangian_type-1 << ")";
3221  }
3222 
3223  throw OomphLibError(error_message.str(),
3224  OOMPH_CURRENT_FUNCTION,
3225  OOMPH_EXCEPTION_LOCATION);
3226  }
3227 }
3228 
3229 
3230 
3231 //========================================================================
3232 /// Steady constructor. The node has NLgrangian Lagrangian
3233 /// coordinates of n_lagrangian_type types (1 for Lagrange elements,
3234 /// 2 for 1D Hermite etc.).
3235 /// The Eulerian dimension of the Node is n_dim and we have n_position_type
3236 /// (generalised) Eulerian coordinates. There are
3237 /// initial_n_value values stored at
3238 /// this node and NAdditional_solid additional values associated with the
3239 /// solid equations are stored in a separate Data object at the node.
3240 //========================================================================
3241 SolidNode::SolidNode(const unsigned &n_lagrangian,
3242  const unsigned &n_lagrangian_type,
3243  const unsigned &n_dim,
3244  const unsigned &n_position_type,
3245  const unsigned &initial_n_value)
3246  : Node(n_dim,n_position_type,initial_n_value,false),
3247  Nlagrangian(n_lagrangian), Nlagrangian_type(n_lagrangian_type)
3248 {
3249  //Calculate the total storage required for positions
3250  const unsigned n_storage = Ndim*Nposition_type;
3251 
3252  //Allocate a data object with exactly the coorect number of coordinates
3253  Variable_position_pt = new Data(n_storage);
3254  //Set X_position to point to the data's positions
3256 
3257  //Setup the lagrangian storage
3258  const unsigned n_lagrangian_storage = n_lagrangian*n_lagrangian_type;
3259  Xi_position = new double[n_lagrangian_storage];
3260  //Initialise lagrangian positions to zero
3261  for(unsigned j=0;j<n_lagrangian_storage;j++) {Xi_position[j] = 0.0;}
3262 }
3263 
3264 //========================================================================
3265 /// Unsteady constructor.
3266 /// Allocates storage for initial_n_value nodal values with history values
3267 /// as required by timestepper.
3268 /// The node has NLgrangian Lagrangian coordinates of
3269 /// n_lagrangian_type types (1 for Lagrange elements, 2 for 1D Hermite etc.)/
3270 /// The Eulerian dimension of the Node is n_dim and we have n_position_type
3271 /// generalised Eulerian coordinates.
3272 //========================================================================
3273 SolidNode::SolidNode(TimeStepper* const &time_stepper_pt_,
3274  const unsigned &n_lagrangian,
3275  const unsigned &n_lagrangian_type,
3276  const unsigned &n_dim,
3277  const unsigned &n_position_type,
3278  const unsigned &initial_n_value)
3279  : Node(time_stepper_pt_,n_dim,n_position_type,initial_n_value,false),
3280  Nlagrangian(n_lagrangian), Nlagrangian_type(n_lagrangian_type)
3281 {
3282  //Calculate the total storage required for positions
3283  const unsigned n_storage = n_dim*n_position_type;
3284 
3285  //Allocate a Data value to each element of the Vector
3286  Variable_position_pt = new Data(time_stepper_pt_,n_storage);
3287  //Set the pointer
3289 
3290  //Setup the lagrangian storage
3291  const unsigned n_lagrangian_storage = n_lagrangian*n_lagrangian_type;
3292  Xi_position = new double[n_lagrangian_storage];
3293  //Initialise lagrangian positions to zero
3294  for(unsigned j=0;j<n_lagrangian_storage;j++) {Xi_position[j] = 0.0;}
3295 }
3296 
3297 //========================================================================
3298 ///Destructor to clean up the memory allocated for nodal positions and
3299 ///additional solid variables
3300 //========================================================================
3302 {
3303  //Null out X_position so that the Node destructor doesn't delete it
3304  X_position=0;
3305  //Delete the position data
3307  //Now clean up lagrangian position data
3308  delete[] Xi_position; Xi_position=0;
3309 }
3310 
3311 
3312 //================================================================
3313 /// Copy nodal positions and associated data from specified
3314 /// node object
3315 //================================================================
3316 void SolidNode::copy(SolidNode* orig_node_pt)
3317 {
3318  // Eulerian positions are stored as Data, so copy the data values
3319  // from one data to another
3321 
3322  //Copy the Lagrangian coordinates
3323  const unsigned nlagrangian_storage = Nlagrangian*Nlagrangian_type;
3324 
3325  // Check # of values:
3326  const unsigned long nlagrangian_storage_orig
3327  = orig_node_pt->nlagrangian()*orig_node_pt->nlagrangian_type();
3328  if (nlagrangian_storage!=nlagrangian_storage_orig)
3329  {
3330  std::ostringstream error_stream;
3331  error_stream << "The allocated lagrangian storage "
3332  << nlagrangian_storage
3333  << " is not the same as the original Solid Node "
3334  << nlagrangian_storage_orig << std::endl;
3335 
3336  throw OomphLibError(error_stream.str(),
3337  OOMPH_CURRENT_FUNCTION,
3338  OOMPH_EXCEPTION_LOCATION);
3339  }
3340 
3341  // Copy lagrangian positions
3342  for(unsigned j=0;j<nlagrangian_storage;j++)
3343  {
3344  Xi_position[j] = orig_node_pt->Xi_position[j];
3345  }
3346 
3347  // Copy the associated data
3348  Data::copy(orig_node_pt);
3349 
3350 }
3351 
3352 
3353 //================================================================
3354 ///Dump nodal positions and associated data to file for restart
3355 //================================================================
3356 void SolidNode::dump(std::ostream& dump_file) const
3357 {
3358  //Dump out the Lagrangian coordinates
3359  // Number of lagrangian values
3360  const unsigned nlagrangian_storage = Nlagrangian*Nlagrangian_type;
3361  dump_file << nlagrangian_storage
3362  << " # number of Lagrangian position variables" << std::endl;
3363 
3364  for(unsigned j=0;j<nlagrangian_storage;j++)
3365  {
3366  dump_file << Xi_position[j] << std::endl;
3367  }
3368 
3369  // Dump out Eulerian positions and nodal data
3370  Node::dump(dump_file);
3371 }
3372 
3373 //================================================================
3374 /// Read nodal positions and associated data to file for restart
3375 //================================================================
3376 void SolidNode::read(std::ifstream& restart_file)
3377 {
3378  std::string input_string;
3379 
3380  // Number of lagrangian values
3381  const unsigned nlagrangian_storage = Nlagrangian*Nlagrangian_type;
3382 
3383  // Read line up to termination sign
3384  getline(restart_file,input_string,'#');
3385  // Ignore rest of line
3386  restart_file.ignore(80,'\n');
3387  // Check # of values:
3388  const unsigned long check_nlagrangian_storage=atoi(input_string.c_str());
3389  if(check_nlagrangian_storage!=nlagrangian_storage)
3390  {
3391  std::ostringstream error_stream;
3392  error_stream << "The allocated Lagrangian storage "
3393  << nlagrangian_storage <<
3394  " is not the same as that in the input file"
3395  << check_nlagrangian_storage << std::endl;
3396 
3397  throw OomphLibError(error_stream.str(),
3398  OOMPH_CURRENT_FUNCTION,
3399  OOMPH_EXCEPTION_LOCATION);
3400 
3401  }
3402 
3403  // Read Lagrangian positions
3404  for(unsigned j=0;j<nlagrangian_storage;j++)
3405  {
3406  // Read line
3407  getline(restart_file,input_string);
3408 
3409  // Transform to double
3410  Xi_position[j] = atof(input_string.c_str());
3411  }
3412 
3413  // Read Eulerian positions and nodal data
3414  Node::read(restart_file);
3415 }
3416 
3417 //===================================================================
3418 /// Set the variable position data from an external source.
3419 /// This is mainly used when setting periodic solid problems.
3420 //==================================================================
3422 {
3423  //Wipe the existing value
3424  delete Variable_position_pt;
3425  //Set the new value
3426  Variable_position_pt = new CopiedData(data_pt);
3427  //Set the new value of x
3429 }
3430 
3431 
3432 //================================================================
3433 /// Set a new position TimeStepper be resizing the appropriate storage.
3434 /// The current (zero) values will be unaffected, but all other entries
3435 /// will be set to zero.
3436 //================================================================
3438  const &position_time_stepper_pt,
3439  const bool &preserve_existing_data)
3440 {
3441  //If the timestepper is unchanged do nothing
3442  if(Position_time_stepper_pt==position_time_stepper_pt) {return;}
3443 
3444  //Set the new time stepper
3446 
3447  //Now simply set the time stepper of the variable position data
3448  this->Variable_position_pt->set_time_stepper(position_time_stepper_pt,
3449  preserve_existing_data);
3450  //Need to reset the X_position to point to the variable positions data
3451  //values which have been reassigned
3453 }
3454 
3455 
3456 //====================================================================
3457 /// Check whether the pointer parameter_pt refers to positional data
3458 //====================================================================
3460  double* const &parameter_pt)
3461 {
3462  //Simply pass the test through to the data of the Variabl position pointer
3463  return
3465 }
3466 
3467 
3468 
3469 //=======================================================================
3470 /// Return lagrangian coordinate
3471 /// either directly or via hanging node representation.
3472 //======================================================================
3473 double SolidNode::lagrangian_position(const unsigned &i) const
3474 {
3475  double posn=0.0;
3476 
3477  // Non-hanging node: just return value
3478  if(!is_hanging()) {posn=xi(i);}
3479  // Hanging node: Use hanging node representation
3480  else
3481  {
3482  // Initialise
3483  double lagn_position=0.0;
3484 
3485  // Add contribution from master nodes
3486  const unsigned nmaster=hanging_pt()->nmaster();
3487  for (unsigned imaster=0;imaster<nmaster;imaster++)
3488  {
3489  lagn_position+=
3490  static_cast<SolidNode*>(hanging_pt()->master_node_pt(imaster))->xi(i)*
3491  hanging_pt()->master_weight(imaster);
3492  }
3493  posn=lagn_position;
3494  }
3495  return posn;
3496 }
3497 
3498 //=======================================================================
3499 /// Return generalised lagrangian coordinate
3500 /// either directly or via hanging node representation.
3501 //======================================================================
3502 double SolidNode::lagrangian_position_gen(const unsigned &k,
3503  const unsigned &i) const
3504 {
3505  double posn=0.0;
3506 
3507  // Non-hanging node: just return value
3508  if(!is_hanging()) {posn=xi_gen(k,i);}
3509  // Hanging node: Use hanging node representation
3510  else
3511  {
3512  // Initialise
3513  double lagn_position=0.0;
3514 
3515  // Add contribution from master nodes
3516  const unsigned nmaster=hanging_pt()->nmaster();
3517  for (unsigned imaster=0;imaster<nmaster;imaster++)
3518  {
3519  lagn_position+=
3520  static_cast<SolidNode*>(hanging_pt()->master_node_pt(imaster))
3521  ->xi_gen(k,i)*
3522  hanging_pt()->master_weight(imaster);
3523  }
3524  posn=lagn_position;
3525  }
3526  return posn;
3527 }
3528 
3529 //================================================================
3530 /// Assign (global) equation number, for SolidNodes
3531 //================================================================
3532 void SolidNode::assign_eqn_numbers(unsigned long &global_number,
3533  Vector<double*>& dof_pt)
3534 {
3535  //Let's call position equations first
3536  Variable_position_pt->assign_eqn_numbers(global_number,dof_pt);
3537  //Then call standard Data assign_eqn_numbers
3538  Data::assign_eqn_numbers(global_number,dof_pt);
3539 }
3540 
3541 //================================================================
3542 /// \short Function to describe the dofs of the SolidNode. The ostream
3543 /// specifies the output stream to which the description
3544 /// is written; the string stores the currently
3545 /// assembled output that is ultimately written to the
3546 /// output stream by Data::describe_dofs(...); it is typically
3547 /// built up incrementally as we descend through the
3548 /// call hierarchy of this function when called from
3549 /// Problem::describe_dofs(...)
3550 //================================================================
3551 void SolidNode::describe_dofs(std::ostream& out,
3552  const std::string& current_string) const
3553 {
3554  //Let's call position equations first
3555  {
3556  std::stringstream conversion;
3557  conversion << " of Solid Node Position" << current_string;
3558  std::string in(conversion.str());
3560  }// Let conversion go out of scope.
3561 
3562  //Then call standard Data version
3563  std::stringstream conversion;
3564  conversion << " of Data" << current_string;
3565  std::string in(conversion.str());
3566  Data::describe_dofs(out,in);
3567 
3568 }
3569 
3570 //================================================================
3571 /// Add pointers to all data values (including position data)
3572 /// to a map
3573 //================================================================
3575 (std::map<unsigned,double*> &map_of_value_pt)
3576 {
3577  //Add the position values first
3578  Variable_position_pt->add_value_pt_to_map(map_of_value_pt);
3579  //Then call standard Data values
3580  Data::add_value_pt_to_map(map_of_value_pt);
3581 }
3582 
3583 
3584 #ifdef OOMPH_HAS_MPI
3585 //==================================================================
3586 /// Add Lagrangian coordinates to the vector after positional data
3587 /// and "standard" data
3588 //==================================================================
3590 {
3591  //Firstly add the position and value data to the vector
3592  Node::add_values_to_vector(vector_of_values);
3593 
3594  //Now add the additional Lagrangian data to the vector
3595 
3596  //Find the total amount of storage required for the Lagrangian variables
3597  const unsigned n_lagrangian_storage =
3598  this->nlagrangian()*this->nlagrangian_type();
3599 
3600  //Resize the vector to accommodate the new data
3601  const unsigned n_current_value = vector_of_values.size();
3602 
3603 #ifdef PARANOID
3604  unsigned n_debug=1;
3605 #else
3606  unsigned n_debug=0;
3607 #endif
3608 
3609  vector_of_values.resize(n_current_value + n_lagrangian_storage+n_debug);
3610 
3611  //Now add the data to the vector
3612  unsigned index = n_current_value;
3613 
3614 #ifdef PARANOID
3615  vector_of_values[index++]=n_lagrangian_storage;
3616 #endif
3617 
3618  //Pointer to the first entry in the data array
3619  double* data_pt = Xi_position;
3620 
3621  //Loop over values
3622  for(unsigned i=0;i<n_lagrangian_storage;i++)
3623  {
3624  //Add the position data to the vector
3625  vector_of_values[index] = *data_pt;
3626  //Increment the counter and the pointer
3627  ++index;
3628  ++data_pt;
3629  }
3630 }
3631 
3632 //==================================================================
3633 /// Read the lagrangian coordinates in from the vector after
3634 /// reading in the positional and "standard" data.
3635 //==================================================================
3637  unsigned &index)
3638 {
3639  //Read the standard nodal data and positions
3640  Node::read_values_from_vector(vector_of_values,index);
3641 
3642  //Now read the additional Lagrangian data to the vector
3643 
3644  //Find the total amount of storage required for the Lagrangian variables
3645  const unsigned n_lagrangian_storage =
3646  this->nlagrangian()*this->nlagrangian_type();
3647 
3648 #ifdef PARANOID
3649  unsigned orig_n_lagrangian_storage=unsigned(vector_of_values[index++]);
3650  if (orig_n_lagrangian_storage!=n_lagrangian_storage)
3651  {
3652  std::ostringstream error_stream;
3653  error_stream << "Non-matching number of values:\n"
3654  << "sent and local n_lagrangian_storage: "
3655  << orig_n_lagrangian_storage << " "
3656  << n_lagrangian_storage << std::endl;
3657  throw OomphLibError(error_stream.str(),
3658  OOMPH_CURRENT_FUNCTION,
3659  OOMPH_EXCEPTION_LOCATION);
3660  }
3661 #endif
3662 
3663  //Pointer to the first entry in the data array
3664  double* data_pt = Xi_position;
3665 
3666  //Loop over values
3667  for(unsigned i=0;i<n_lagrangian_storage;i++)
3668  {
3669  //Read the position data from the vector
3670  *data_pt = vector_of_values[index];
3671  //Increment the counter and the pointer
3672  ++index;
3673  ++data_pt;
3674  }
3675 }
3676 
3677 //==================================================================
3678 /// Add equations numbers associated with the node position
3679 /// to the vector after the standard nodal equation numbers
3680 //==================================================================
3682 {
3683  //Firstly add the standard equation numbers
3684  Data::add_eqn_numbers_to_vector(vector_of_eqn_numbers);
3685  //Now add the equation numbers associated with the positional data
3686  Variable_position_pt->add_eqn_numbers_to_vector(vector_of_eqn_numbers);
3687 }
3688 
3689 //=======================================================================
3690 /// Read the equation numbers associated with the node position
3691 /// from the vector after reading in the standard nodal equaiton numbers
3692 //=======================================================================
3694  const Vector<long> &vector_of_eqn_numbers,
3695  unsigned &index)
3696 {
3697  //Read the standard nodal data and positions
3698  Data::read_eqn_numbers_from_vector(vector_of_eqn_numbers,index);
3699  //Now add the equation numbers associated with the positional data
3701  read_eqn_numbers_from_vector(vector_of_eqn_numbers,index);
3702 }
3703 
3704 #endif
3705 
3706 
3707 }
double * Xi_position
Storage for the Lagrangian positions.
Definition: nodes.h:1589
Vector< double > value() const
Return vector of values calculated using value(vector).
Definition: nodes.h:1399
void read_eqn_numbers_from_vector(const Vector< long > &vector_of_eqn_numbers, unsigned &index)
Read all equation numbers from the vector starting from index. On return the index will be set to the...
Definition: nodes.cc:3693
std::map< unsigned, DenseMatrix< double > * > * Boundary_coordinates_pt
Pointer to a map of pointers to intrinsic boundary coordinates of the Node, indexed by the boundary n...
Definition: nodes.h:1864
virtual void dump(std::ostream &dump_file) const
Dump nodal position and associated data to file for restart.
Definition: nodes.cc:1892
void clear_copied_pointers()
Clear the pointers to the copied data.
Definition: nodes.cc:1257
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
Definition: nodes.cc:2301
virtual ~SolidNode()
Destructor that cleans up the additional memory allocated in SolidNodes.
Definition: nodes.cc:3301
TimeStepper * Time_stepper_pt
Pointer to a Timestepper. The inclusion of a Timestepper pointer in the Data class, ensures that time-derivatives can be calculated and storage can be managed at the low (Data) level.
Definition: nodes.h:126
void delete_value_storage()
Delete all storage allocated by the Data object for values and equation numbers.
Definition: nodes.cc:190
void remove_from_boundary(const unsigned &b)
Remove the node from the mesh boundary b.
Definition: nodes.cc:2948
bool is_pinned(const unsigned &i) const
Test whether the i-th variable is pinned (1: true; 0: false).
Definition: nodes.h:403
bool does_pointer_correspond_to_value(double *const &parameter_pt)
Check whether the pointer parameter_pt addresses internal data values.
Definition: nodes.cc:533
void constrain(const unsigned &i)
Constrain the i-th stored variable when making hanging data If the data is already pinned leave it al...
Definition: nodes.h:416
virtual void unconstrain_positions()
Unconstrain the positions when the node is made non-hanging Empty virtual function that is overloaded...
Definition: nodes.h:1271
Data * Copied_data_pt
Pointer to the Data object from which the values are copied.
Definition: nodes.h:605
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:847
void set_nonhanging()
Label node as non-hanging node by removing all hanging node data.
Definition: nodes.cc:2204
long * eqn_number_pt(const unsigned &i)
Return the pointer to the equation number of the i-th stored variable.
Definition: nodes.h:356
void dump(std::ostream &dump_file) const
Dump the data object to a file.
Definition: nodes.cc:608
void add_to_boundary(const unsigned &b)
Add the node to the mesh boundary b.
Definition: nodes.cc:2919
cstr elem_len * i
Definition: cfortran.h:607
double * Master_weights
C-style array of weights for the dofs on the master nodes.
Definition: nodes.h:787
unsigned Nvalue
Number of values stored in the data object.
Definition: nodes.h:141
virtual void add_value_pt_to_map(std::map< unsigned, double * > &map_of_value_pt)
Add pointers to all unpinned and unconstrained data to a map indexed by (global) equation number...
Definition: nodes.cc:1041
virtual void make_periodic(Node *const &node_pt)
Make the node periodic by copying the values from node_pt. Note that the coordinates will always rema...
Definition: nodes.cc:2180
virtual ~BoundaryNodeBase()
Destructor, clean up any allocated storage for the boundaries.
Definition: nodes.cc:2880
virtual void constrain_positions()
Constrain the positions when the node is made hanging Empty virtual function that is overloaded in So...
Definition: nodes.h:1267
void xi_gen_range_check(const unsigned &k, const unsigned &i) const
Private function to check that the arguments to the position functions are in range.
Definition: nodes.cc:3202
virtual Node * copied_node_pt() const
Return pointer to copied node (null if the current node is not a copy – always the case here; it's ov...
Definition: nodes.h:1042
Node()
Default constructor.
Definition: nodes.cc:1517
void read_values_from_vector(const Vector< double > &vector_of_values, unsigned &index)
Read all data and time history values from the vector starting from index. On return the index will b...
Definition: nodes.cc:3636
Data * Variable_position_pt
Pointer to data that will hold variable positions in elastic nodes.
Definition: nodes.h:1586
bool does_pointer_correspond_to_position_data(double *const &parameter_pt)
Overload the check whether the pointer parameter_pt addresses position data values.
Definition: nodes.cc:3459
static TimeStepper * Default_static_time_stepper_pt
Default (static) timestepper used in steady problems.
Definition: nodes.h:169
Node ** Master_nodes_pt
C-style array of pointers to nodes that this hanging node depends on.
Definition: nodes.h:784
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:557
char t
Definition: cfortran.h:572
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
Definition: nodes.h:322
friend class CopiedData
Definition: nodes.h:99
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void add_values_to_vector(Vector< double > &vector_of_values)
Add all data and time history values to the vector. Overloaded to add the position information as wel...
Definition: nodes.cc:2657
unsigned self_test()
Self-test: Have all values been classified as pinned/unpinned? Return 0 if OK.
Definition: nodes.cc:930
static long Is_constrained
Static "Magic number" used in place of the equation number to indicate that the value is constrained ...
Definition: nodes.h:207
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
OomphInfo oomph_info
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:375
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
virtual void assign_eqn_numbers(unsigned long &global_ndof, Vector< double * > &dof_pt)
Assign global equation numbers; increment global number of unknowns, global_ndof; and add any new dof...
Definition: nodes.cc:861
void set_master_node_pt(const unsigned &i, Node *const &master_node_pt, const double &weight)
Set the pointer to the i-th master node and its weight.
Definition: nodes.cc:1417
void assign_eqn_numbers(unsigned long &global_number, Vector< double * > &dof_pt)
Overload the assign equation numbers routine.
Definition: nodes.cc:3532
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
void range_check(const unsigned &i) const
Check that the argument is within the range of stored data values.
Definition: nodes.cc:1397
virtual void get_coordinates_on_boundary(const unsigned &b, const unsigned &k, Vector< double > &boundary_zeta)
Return the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interface provides run-time error checking.
Definition: nodes.cc:2287
virtual unsigned ncoordinates_on_boundary(const unsigned &b)
Get the number of boundary coordinates on mesh boundary b. Broken virtual interface provides run-time...
Definition: nodes.cc:2271
void read(std::ifstream &restart_file)
Read nodal positions (variable and fixed) and associated data from file for restart.
Definition: nodes.cc:3376
double lagrangian_position(const unsigned &i) const
Return lagrangian coordinate either directly or via hanging node representation.
Definition: nodes.cc:3473
HijackedData(const unsigned &copied_value, Data *const &data_pt)
Constructor.
Definition: nodes.cc:1271
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:966
void clear_copied_pointers()
Clear the pointers to the copied data.
Definition: nodes.cc:1338
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
Definition: nodes.h:1654
void get_coordinates_on_boundary(const unsigned &b, Vector< double > &boundary_zeta)
Return the vector of boundary coordinates on mesh boundary b.
Definition: nodes.h:2055
virtual void remove_from_boundary(const unsigned &b)
Broken interface for removing the node from the mesh boundary b Here to provide error reporting...
Definition: nodes.cc:2258
bool is_halo() const
Is this Data a halo?
Definition: nodes.h:490
unsigned Nlagrangian_type
Number of types of Lagrangian coordinates used to interpolate the Lagrangian coordinates within the e...
Definition: nodes.h:1583
void dump(std::ostream &dump_file) const
Dump nodal positions (variable and fixed) and associated data to file for restart.
Definition: nodes.cc:3356
virtual void make_periodic_nodes(const Vector< Node * > &periodic_nodes_pt)
Make the nodes passed in the vector periodic_nodes share the same data as this node.
Definition: nodes.cc:2193
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
double & xi(const unsigned &i)
Reference to i-th Lagrangian position.
Definition: nodes.h:1740
A class that contains the information required by Nodes that are located on Mesh boundaries. A BoundaryNode of a particular type is obtained by combining a given Node with this class. By differentiating between Nodes and BoundaryNodes we avoid a lot of un-necessary storage in the bulk Nodes.
Definition: nodes.h:1854
double & xi_gen(const unsigned &k, const unsigned &i)
Reference to the generalised Lagrangian position. `Type': k; 'Coordinate direction: i...
Definition: nodes.h:1759
void add_values_to_vector(Vector< double > &vector_of_values)
Add all data, position and time history values to the vector Overload to add the Lagrangian coordinat...
Definition: nodes.cc:3589
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
void set_coordinates_on_boundary(const unsigned &b, const Vector< double > &boundary_zeta)
Set the vector of boundary coordinates on mesh boundary b.
Definition: nodes.h:2064
void read(std::ifstream &restart_file)
Read nodal position and associated data from file for restart.
Definition: nodes.cc:1919
virtual bool is_a_copy() const
Return a boolean to indicate whether the Data objact contains any copied values. A base Data object c...
Definition: nodes.h:255
virtual void add_values_to_vector(Vector< double > &vector_of_values)
Add all data and time history values to the vector in the internal storage order. ...
Definition: nodes.cc:1065
unsigned Nposition_type
Number of coordinate types used in the mapping between local and global coordinates (e...
Definition: nodes.h:900
void make_nodes_periodic(Node *const &node_pt, Vector< Node * > const &periodic_copies_pt)
Helper function that is used to turn BoundaryNodes into periodic boundary nodes by setting the data v...
Definition: nodes.cc:2778
bool is_segregated_solve_pinned(const unsigned &i)
Test whether the i-th variable is temporaily pinned for a segregated solve.
Definition: nodes.h:408
unsigned Ncopies
Number of Data that contain copies of this Data object's values.
Definition: nodes.h:136
static long Is_pinned
Static "Magic number" used in place of the equation number to indicate that the value is pinned...
Definition: nodes.h:192
double ** Value
C-style array of pointers to data values and possible history values. The data must be ordered in suc...
Definition: nodes.h:116
void copy(SolidNode *orig_node_pt)
Copy nodal positions and associated data from specified node object.
Definition: nodes.cc:3316
void set_position_time_stepper(TimeStepper *const &position_time_stepper_pt, const bool &preserve_existing_data)
Set a new position timestepper be resizing the appropriate storage Overloaded from the basic implemen...
Definition: nodes.cc:3437
double position_gen(const unsigned &k, const unsigned &i) const
Return generalised nodal coordinate either directly or via hanging node representation.
Definition: nodes.cc:2481
void add_copy(Data *const &data_pt)
Add the pointer data_pt to the array Copy_of_data_pt. This should be used whenever copies are made of...
Definition: nodes.cc:84
long * Eqn_number
C-style array of pointers to the (global) equation numbers of the values.
Definition: nodes.h:120
virtual void describe_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the dofs of the Node. The ostream specifies the output stream to which the descr...
Definition: nodes.cc:905
void copy(Data *orig_data_pt)
Copy Data values from specified Data object.
Definition: nodes.cc:561
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
double raw_value(const unsigned &i) const
Return the i-th value stored at the Node. This interface does NOT take the hanging status of the Node...
Definition: nodes.h:1357
void set_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set a new timestepper by resizing the appropriate storage. If already assigned the equation numbering...
Definition: nodes.cc:392
void make_node_periodic(Node *const &node_pt, Node *const &original_node_pt)
Helper function that is used to turn BoundaryNodes into peridic boundary nodes by setting the data va...
Definition: nodes.cc:2837
unsigned Nmaster
Number of master nodes required by this hanging node.
Definition: nodes.h:790
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:507
static long Is_segregated_solve_pinned
Static "Magic number" used in place of the equation number to indicate that the value is pinned...
Definition: nodes.h:197
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void unconstrain(const unsigned &i)
Unconstrain the i-th stored variable when make the data nonhanging. Only unconstrain if it was actual...
Definition: nodes.h:421
virtual ~Node()
Destructor: Clean up the memory allocated for nodal position.
Definition: nodes.cc:1628
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i).
Definition: nodes.h:1055
virtual void read_values_from_vector(const Vector< double > &vector_of_values, unsigned &index)
Read all data and time history values from the vector starting from index. On return the index will b...
Definition: nodes.cc:1129
void reset_copied_pointers()
Reset the pointers to the copied data.
Definition: nodes.cc:1243
void resize(const unsigned &n_value)
We cannot resize HijackedData, so the resize function throws a warning.
Definition: nodes.cc:1303
void read(std::ifstream &restart_file)
Read data object from a file.
Definition: nodes.cc:635
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
void range_check(const unsigned &t, const unsigned &i) const
Check that the arguments are within the range of the stored data values and timesteps.
Definition: nodes.cc:52
unsigned Copied_index
Index of the value that is copied from within the Data object.
Definition: nodes.h:538
bool is_constrained(const unsigned &i)
Test whether the i-th variable is constrained (1: true; 0: false).
Definition: nodes.h:439
HangInfo ** Hanging_pt
C-style array of pointers to hanging node info. It's set to NULL if the node isn't hanging...
Definition: nodes.h:891
Data ** Copy_of_data_pt
C-style array of any Data objects that contain copies of the current Data object's data values...
Definition: nodes.h:132
Class that contains data for hanging nodes.
Definition: nodes.h:684
unsigned ncoordinates_on_boundary(const unsigned &b)
Get the number of boundary coordinates on mesh boundary b.
Definition: nodes.cc:3003
virtual void set_position_time_stepper(TimeStepper *const &position_time_stepper_pt, const bool &preserve_existing_data)
Set a new position timestepper be resizing the appropriate storage.
Definition: nodes.cc:1678
double lagrangian_position_gen(const unsigned &k, const unsigned &i) const
Return generalised lagrangian coordinate either directly or via hanging node representation.
Definition: nodes.cc:3502
std::map< unsigned, unsigned > * Index_of_first_value_assigned_by_face_element_pt
Pointer to a map, indexed by the face element identifier it returns the position of the first face el...
Definition: nodes.h:1878
void remove_copy(Data *const &data_pt)
Remove the pointer data_pt from the array Copy_of_data_pt. This should be used whenever copies of the...
Definition: nodes.cc:108
void read_values_from_vector(const Vector< double > &vector_of_values, unsigned &index)
Read all data and time history values from the vector starting from index. On return the index will b...
Definition: nodes.cc:2712
void add_master_node_pt(Node *const &master_node_pt, const double &weight)
Add (pointer to) master node and corresponding weight to the internally stored (pointers to) master n...
Definition: nodes.cc:1432
void set_external_variable_position_pt(Data *const &data_pt)
Set the variable position data from an external data object.
Definition: nodes.cc:3421
void x_gen_range_check(const unsigned &t, const unsigned &k, const unsigned &i) const
Private function to check that the arguemnts to the position functions are in range.
Definition: nodes.cc:1469
void copy(Node *orig_node_pt)
Copy all nodal data from specified Node object.
Definition: nodes.cc:1834
void set_hanging_pt(HangInfo *const &hang_pt, const int &i)
Set the hanging data for the i-th value. (hang_pt=0 to make non-hanging)
Definition: nodes.cc:1997
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
unsigned Ndim
Eulerian dimension of the node.
Definition: nodes.h:894
SolidNode()
Default Constructor.
Definition: nodes.h:1594
void resize(const unsigned &n_value)
Resize the number of equations.
Definition: nodes.cc:2089
virtual void clear_copied_pointers()
Helper function that should be overloaded derived classes that contain copies of data. The function must unset (NULL out) the internal pointers to the copied data. This is used when destructing data to ensure that all pointers remain valid. The default implementation throws an error because Data cannot be a copy.
Definition: nodes.cc:179
double dposition_gen_dt(const unsigned &k, const unsigned &i) const
i-th component of time derivative (velocity) of the generalised position, dx(k,i)/dt. `Type': k; Coordinate direction: i. This function uses the hanging node representation if necessary.
Definition: nodes.cc:2591
std::ostream & operator<<(std::ostream &out, const DoubleVector &v)
output operator
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
TimeStepper * Position_time_stepper_pt
Pointer to the timestepper associated with the position data.
Definition: nodes.h:881
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
static long Is_unclassified
Static "Magic number" used in place of the equation number to denote a value that hasn't been class...
Definition: nodes.h:201
Node * Copied_node_pt
If the BoundaryNode is periodic, this pointer is set to the BoundaryNode whose data it shares...
Definition: nodes.h:1882
CopiedData(Data *const &data_pt)
Constructor.
Definition: nodes.cc:1348
static unsigned No_independent_position
Static "Magic number" used to indicate that there is no independent position in a periodic node...
Definition: nodes.h:919
virtual void resize(const unsigned &n_value)
Change (increase) the number of values that may be stored.
Definition: nodes.cc:961
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
unsigned nlagrangian_type() const
Number of types of Lagrangian coordinates used to interpolate the Lagrangian coordinates within the e...
Definition: nodes.h:1737
void reset_copied_pointers()
Reset the pointers to the copied data.
Definition: nodes.cc:1321
unsigned Nlagrangian
Number of Lagrangian coordinates of the node.
Definition: nodes.h:1579
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
bool is_on_boundary() const
Test whether the node lies on a boundary.
Definition: nodes.h:2046
void output(std::ostream &outfile)
Output nodal position.
Definition: nodes.cc:2642
std::set< unsigned > * Boundaries_pt
Pointer to set of mesh boundaries occupied by the Node; NULL if the Node is not on any boundaries...
Definition: nodes.h:1868
double dposition_dt(const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt, either directly or via hanging node representatio...
Definition: nodes.cc:2542
double dx_gen_dt(const unsigned &k, const unsigned &i) const
i-th component of time derivative (velocity) of the generalised position, dx(k,i)/dt. `Type': k; Coordinate direction: i.
Definition: nodes.cc:1783
Vector< double > position() const
Return vector of position of node at current time.
Definition: nodes.h:1422
virtual void add_to_boundary(const unsigned &b)
Broken interface for adding the node to the mesh boundary b Essentially here for error reporting...
Definition: nodes.cc:2243
Data()
Default: Just set pointer to (steady) timestepper. No storage for values is allocated.
Definition: nodes.cc:235
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:291
double dx_dt(const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt.
Definition: nodes.cc:1735
virtual void assign_eqn_numbers(unsigned long &global_ndof, Vector< double * > &dof_pt)
Assign global equation numbers; increment global number of unknowns, global_ndof; and add any new dof...
Definition: nodes.cc:499
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
void add_value_pt_to_map(std::map< unsigned, double * > &map_of_value_pt)
Overload the function add_values_to_map so that it also adds the variable position data...
Definition: nodes.cc:3575
virtual void read_eqn_numbers_from_vector(const Vector< long > &vector_of_eqn_numbers, unsigned &index)
Read all equation numbers from the vector starting from index. On return the index will be set to the...
Definition: nodes.cc:1216
void describe_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the dofs of the Node. The ostream specifies the output stream to which the descr...
Definition: nodes.cc:3551
virtual void add_eqn_numbers_to_vector(Vector< long > &vector_of_eqn_numbers)
Add all equation numbers to the vector in the internal storage order.
Definition: nodes.cc:1181
double ** X_position
Array of pointers to the data holding the Eulerian positions. The storage format must be the same as ...
Definition: nodes.h:878
void add_eqn_numbers_to_vector(Vector< long > &vector_of_eqn_numbers)
Add all equation numbers to the vector in the internal storage order. Overload to add equation number...
Definition: nodes.cc:3681
unsigned nlagrangian() const
Return number of lagrangian coordinates.
Definition: nodes.h:1733
virtual ~Data()
Destructor, deallocates memory assigned for data.
Definition: nodes.cc:454
Data * Copied_data_pt
Pointer to the Data object from which the value is copied.
Definition: nodes.h:535
void resize(const unsigned &n_value)
We cannot resize CopiedData, so the resize function throws a warning.
Definition: nodes.cc:1378
virtual void reset_copied_pointers()
Helper function that should be overloaded in derived classes that can contain copies of Data...
Definition: nodes.cc:165