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