elements.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// Non-inline member functions for generic elements
27
28#include <float.h>
29
30// oomph-lib includes
31#include "elements.h"
32#include "timesteppers.h"
33#include "integral.h"
34#include "shape.h"
35#include "oomph_definitions.h"
37
38namespace oomph
39{
40 /// Static boolean to suppress warnings about repeated internal
41 /// data. Defaults to false
43 false;
44
45
46 /// Static boolean to suppress warnings about repeated external
47 /// data. Defaults to true
49
50
51 /// ////////////////////////////////////////////////////////////////////////
52 /// ////////////////////////////////////////////////////////////////////////
53 // Functions for generalised elements
54 /// ////////////////////////////////////////////////////////////////////////
55 /// ////////////////////////////////////////////////////////////////////////
56
57 //=======================================================================
58 /// Add a (pointer to an) internal data object to the element and
59 /// return the index required to obtain it from the access
60 /// function \c internal_data_pt()
61 //=======================================================================
63 const bool& fd)
64 {
65 // Local cache of numbers of internal and external data
66 const unsigned n_internal_data = Ninternal_data;
67 const unsigned n_external_data = Nexternal_data;
68
69 // Find out whether the data is already stored in the array
70
71 // Loop over the number of internals
72 // The internal data are stored at the beginning of the array
73 for (unsigned i = 0; i < n_internal_data; i++)
74 {
75 // If the passed pointer is stored in the array
76 if (internal_data_pt(i) == data_pt)
77 {
78#ifdef PARANOID
80 {
81 oomph_info << std::endl << std::endl;
83 << "---------------------------------------------------------------"
84 "--"
85 << std::endl
86 << "Info: Data is already included in this element's internal "
87 "storage."
88 << std::endl
89 << "It's stored as entry " << i << " and I'm not adding it again."
90 << std::endl
91 << std::endl
92 << "Note: You can suppress this message by recompiling without"
93 << "\n PARANOID or setting the boolean \n"
94 << "\n "
95 "GeneralisedElement::Suppress_warning_about_repeated_internal_"
96 "data"
97 << "\n\n to true." << std::endl
98 << "---------------------------------------------------------------"
99 "--"
100 << std::endl
101 << std::endl;
102 }
103#endif
104 // Return the index to the data object
105 return i;
106 }
107 }
108
109 // Allocate new storage for the pointers to data
110 Data** new_data_pt = new Data*[n_internal_data + n_external_data + 1];
111
112 // Copy the old internal values across to the beginning of the array
113 for (unsigned i = 0; i < n_internal_data; i++)
114 {
115 new_data_pt[i] = Data_pt[i];
116 }
117
118 // Now add the new value to the end of the internal data
119 new_data_pt[n_internal_data] = data_pt;
120
121 // Copy the external values across
122 for (unsigned i = 0; i < n_external_data; i++)
123 {
124 new_data_pt[n_internal_data + 1 + i] = Data_pt[n_internal_data + i];
125 }
126
127 // Delete the storage associated with the previous values
128 delete[] Data_pt;
129
130 // Set the pointer to the new storage
131 Data_pt = new_data_pt;
132
133 // Resize the array of boolean flags
134 Data_fd.resize(n_internal_data + n_external_data + 1);
135 // Shuffle the flags for the external data to the end of the array
136 for (unsigned i = n_external_data; i > 0; i--)
137 {
138 Data_fd[n_internal_data + i] = Data_fd[n_internal_data + i - 1];
139 }
140 // Now add the new flag to the end of the internal data
141 Data_fd[n_internal_data] = fd;
142
143 // Increase the number of internals
145
146 // Return the final index to the new internal data
147 return n_internal_data;
148 }
149
150 //=======================================================================
151 /// Add the contents of the queue global_eqn_numbers to
152 /// the local storage for the equation numbers, which represents the
153 /// local-to-global translation scheme. It is essential that the entries
154 /// are added in order, i.e. from the front.
155 //=======================================================================
157 std::deque<unsigned long> const& global_eqn_numbers,
158 std::deque<double*> const& global_dof_pt)
159 {
160 // Find the number of dofs
161 const unsigned n_dof = Ndof;
162 // Find the number of additional dofs
163 const unsigned n_additional_dof = global_eqn_numbers.size();
164 // If there are none, return immediately
165 if (n_additional_dof == 0)
166 {
167 return;
168 }
169
170 // Find the new total number of equation numbers
171 const unsigned new_n_dof = n_dof + n_additional_dof;
172 // Create storage for all equations, initialised to NULL
173 unsigned long* new_eqn_number = new unsigned long[new_n_dof];
174
175 // Copy over the existing values to the start new storage
176 for (unsigned i = 0; i < n_dof; i++)
177 {
178 new_eqn_number[i] = Eqn_number[i];
179 }
180
181 // Set an index to the next position in the new storage
182 unsigned index = n_dof;
183 // Loop over the queue and add it's entries to our new storage
184 for (std::deque<unsigned long>::const_iterator it =
185 global_eqn_numbers.begin();
186 it != global_eqn_numbers.end();
187 ++it)
188 {
189 // Add the value to the storage
190 new_eqn_number[index] = *it;
191 // Increase the array index
192 ++index;
193 }
194
195
196 // If a non-empty dof deque has been passed then do stuff
197 const unsigned n_additional_dof_pt = global_dof_pt.size();
198 if (n_additional_dof_pt > 0)
199 {
200// If it's size is not the same as the equation numbers complain
201#ifdef PARANOID
202 if (n_additional_dof_pt != n_additional_dof)
203 {
204 std::ostringstream error_stream;
205 error_stream
206 << "global_dof_pt is non-empty, yet it does not have the same size\n"
207 << "as global_eqn_numbers.\n"
208 << "There are " << n_additional_dof << " equation numbers,\n"
209 << "but " << n_additional_dof_pt << std::endl;
210
211 throw OomphLibError(
212 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
213 }
214#endif
215
216 // Create storge for all dofs initialised to NULL
217 double** new_dof_pt = new double*[new_n_dof];
218 // Copy over the exisiting values to the start of new storage
219 for (unsigned i = 0; i < n_dof; i++)
220 {
221 new_dof_pt[i] = Dof_pt[i];
222 }
223
224 // Set an index to the next position in the new storage
225 unsigned index = n_dof;
226 // Loop over the queue and add it's entries to our new storage
227 for (std::deque<double*>::const_iterator it = global_dof_pt.begin();
228 it != global_dof_pt.end();
229 ++it)
230 {
231 // Add the value to the storage
232 new_dof_pt[index] = *it;
233 // Increase the array index
234 ++index;
235 }
236
237 // Now delete the old storage
238 delete[] Dof_pt;
239 // Set the pointer to address the new storage
240 Dof_pt = new_dof_pt;
241 }
242
243 // Now delete the old for the equation numbers storage
244 delete[] Eqn_number;
245 // Set the pointer to address the new storage
246 Eqn_number = new_eqn_number;
247 // Finally update the number of degrees of freedom
248 Ndof = new_n_dof;
249 }
250
251 //========================================================================
252 /// Empty dense matrix used as a dummy argument to combined
253 /// residual and jacobian functions in the case when only the residuals
254 /// are being assembled
255 //========================================================================
257
258 //========================================================================
259 /// Static storage used when pointers to the dofs are being assembled by
260 /// add_global_eqn_numbers()
261 //========================================================================
262 std::deque<double*> GeneralisedElement::Dof_pt_deque;
263
264
265 //=========================================================================
266 /// Default value used as the increment for finite difference calculations
267 /// of the jacobian matrices
268 //=========================================================================
270
271 //==========================================================================
272 /// Destructor for generalised elements: Wipe internal data. Pointers
273 /// to external data get NULLed but are not deleted because they
274 /// are (generally) shared by lots of elements.
275 //==========================================================================
276 GeneralisedElement::~GeneralisedElement()
277 {
278 // Delete each of the objects stored as internal data
279 for (unsigned i = Ninternal_data; i > 0; i--)
280 {
281 // The objects are stored at the beginning of the Data_pt array
282 delete Data_pt[i - 1];
283 Data_pt[i - 1] = 0;
284 }
285
286 // Now delete the storage for internal and external data
287 delete[] Data_pt;
288
289 // Now if we have allocated storage for the local equation for
290 // the internal and external data, delete it.
291 if (Data_local_eqn)
292 {
293 delete[] Data_local_eqn[0];
294 delete[] Data_local_eqn;
295 }
296
297 // Delete the storage for the global equation numbers
298 delete[] Eqn_number;
299 }
300
301
302 //=======================================================================
303 /// Add a (pointer to an) external data object to the element and
304 /// return the index required to obtain it from the access
305 /// function \c external_data_pt()
306 //=======================================================================
308 const bool& fd)
309 {
310 // Find the numbers of internal and external data
311 const unsigned n_internal_data = Ninternal_data;
312 const unsigned n_external_data = Nexternal_data;
313 // Find out whether the data is already stored in the array
314
315 // Loop over the number of externals
316 // The external data are stored at the end of the array Data_pt
317 for (unsigned i = 0; i < n_external_data; i++)
318 {
319 // If the passed pointer is stored in the array
320 if (external_data_pt(i) == data_pt)
321 {
322#ifdef PARANOID
324 {
325 oomph_info << std::endl << std::endl;
327 << "---------------------------------------------------------------"
328 "--"
329 << std::endl
330 << "Info: Data is already included in this element's external "
331 "storage."
332 << std::endl
333 << "It's stored as entry " << i << " and I'm not adding it again"
334 << std::endl
335 << std::endl
336 << "Note: You can suppress this message by recompiling without"
337 << "\n PARANOID or setting the boolean \n"
338 << "\n "
339 "GeneralisedElement::Suppress_warning_about_repeated_external_"
340 "data"
341 << "\n\n to true." << std::endl
342 << "---------------------------------------------------------------"
343 "--"
344 << std::endl
345 << std::endl;
346 }
347#endif
348 // Return the index to the data object
349 return i;
350 }
351 }
352
353 // Allocate new storage for the pointers to data
354 Data** new_data_pt = new Data*[n_internal_data + n_external_data + 1];
355
356 // Copy the old internal and external values across to the new array
357 for (unsigned i = 0; i < (n_internal_data + n_external_data); i++)
358 {
359 new_data_pt[i] = Data_pt[i];
360 }
361
362 // Add the new data pointer to the end of the array
363 new_data_pt[n_internal_data + n_external_data] = data_pt;
364
365 // Delete the storage associated with the previous values
366 delete[] Data_pt;
367
368 // Set the pointer to the new storage
369 Data_pt = new_data_pt;
370
371 // Resize the array of boolean flags
372 Data_fd.resize(n_internal_data + n_external_data + 1);
373 // Now add the new flag to the end of the external data
374 Data_fd[n_internal_data + n_external_data] = fd;
375
376 // Increase the number of externals
378
379 // Return the final index to the new external data
380 return n_external_data;
381 }
382
383 //========================================================================
384 /// Flush all the external data, i.e. clear the pointers to external
385 /// data from the internal storage.
386 //========================================================================
388 {
389 // Get the numbers of internal and external data
390 const unsigned n_external_data = Nexternal_data;
391 // If there is external data
392 if (n_external_data > 0)
393 {
394 // Storage for new data, initialised to NULL
395 Data** new_data_pt = 0;
396
397 // Find the number of internal data
398 const unsigned n_internal_data = Ninternal_data;
399 // If there is internal data resize Data_pt and copy over
400 // the pointers
401 if (n_internal_data > 0)
402 {
403 // The new data pointer should only be the size of the internal data
404 new_data_pt = new Data*[n_internal_data];
405 // Copy over the internal data only
406 for (unsigned i = 0; i < n_internal_data; i++)
407 {
408 new_data_pt[i] = Data_pt[i];
409 }
410 }
411
412 // Delete the old storage
413 delete[] Data_pt;
414 // Set the new storage, this will be NULL if there is no internal data
415 Data_pt = new_data_pt;
416 // Set the number of externals to zero
417 Nexternal_data = 0;
418
419 // Resize the array of boolean flags to the number of internals
420 Data_fd.resize(n_internal_data);
421 }
422 }
423
424
425 //=========================================================================
426 /// Remove the object addressed by data_pt from the external data array
427 /// Note that this could mess up the numbering of other external data
428 //========================================================================
430 {
431 // Get the current numbers of external data
432 const unsigned n_external_data = Nexternal_data;
433 // Index of the data to be removed
434 // We initialise this to be n_external_data, and it will be smaller than
435 // n_external_data if the data pointer is found in the array
436 unsigned index = n_external_data;
437 // Loop over the external data and find the argument
438 for (unsigned i = 0; i < n_external_data; i++)
439 {
440 // If we have found the data pointer, set the index and break
441 if (external_data_pt(i) == data_pt)
442 {
443 index = i;
444 break;
445 }
446 }
447
448 // If we have found an index less than Nexternal_data, then we have found
449 // the data in the array
450 if (index < n_external_data)
451 {
452 // Initialise a new array to NULL
453 Data** new_data_pt = 0;
454
455 // Find the number of internals
456 const unsigned n_internal_data = Ninternal_data;
457
458 // Find the new number of total data items (one fewer)
459 const unsigned n_total_data = n_internal_data + n_external_data - 1;
460
461 // Create a new array containing the data items, if we have any
462 if (n_total_data > 0)
463 {
464 new_data_pt = new Data*[n_total_data];
465 }
466
467 // Copy over all the internal data
468 for (unsigned i = 0; i < n_internal_data; i++)
469 {
470 new_data_pt[i] = Data_pt[i];
471 }
472
473 // Now copy over the un-flushed data
474 unsigned counter = 0;
475 for (unsigned i = 0; i < n_external_data; i++)
476 {
477 // If we are not at the deleted index
478 if (i != index)
479 {
480 // Copy the undeleted entry into the next entry in our new
481 // array of pointers to Data
482 new_data_pt[n_internal_data + counter] = Data_pt[i];
483 // Increase the counter
484 ++counter;
485 }
486 }
487
488 // Delete the storage associated with the previous values
489 delete[] Data_pt;
490
491 // Set pointers to the new storage, will be NULL if no data left
492 Data_pt = new_data_pt;
493
494 // Remove the entry from the array of boolean flags
495 Data_fd.erase(Data_fd.begin() + n_internal_data + index);
496
497 // Decrease the number of externals
499
500 // Issue a warning if there will be external data remaining
501 if (Nexternal_data > 1)
502 {
503 std::ostringstream warning_stream;
504 warning_stream
505 << "Data removed from element's external data " << std::endl
506 << "You may have to update the indices for remaining data "
507 << std::endl
508 << "This can be achieved by using add_external_data() "
509 << std::endl;
510 OomphLibWarning(warning_stream.str(),
511 "GeneralisedElement::flush_external_data()",
512 OOMPH_EXCEPTION_LOCATION);
513 }
514 }
515 }
516
517
518 //==========================================================================
519 /// This function loops over the internal data of the element and assigns
520 /// GLOBAL equation numbers to the data objects.
521 ///
522 /// Pass:
523 /// - the current total number of dofs, global_number, which gets
524 /// incremented.
525 /// - Dof_pt, the Vector of pointers to the global dofs
526 /// (to which any internal dofs are added).
527 /// .
528 //==========================================================================
530 unsigned long& global_number, Vector<double*>& Dof_pt)
531 {
532 // Loop over the internal data and assign the equation numbers
533 // The internal data are stored at the beginning of the Data_pt array
534 for (unsigned i = 0; i < Ninternal_data; i++)
535 {
537 }
538 }
539
540
541 //==========================================================================
542 /// Function to describe the dofs of the Element. The ostream
543 /// specifies the output stream to which the description
544 /// is written; the string stores the currently
545 /// assembled output that is ultimately written to the
546 /// output stream by Data::describe_dofs(...); it is typically
547 /// built up incrementally as we descend through the
548 /// call hierarchy of this function when called from
549 /// Problem::describe_dofs(...)
550 //==========================================================================
552 std::ostream& out, const std::string& current_string) const
553 {
554 for (unsigned i = 0; i < Ninternal_data; i++)
555 {
556 std::stringstream conversion;
557 conversion << " of Internal Data " << i << current_string;
558 std::string in(conversion.str());
560 }
561 }
562
563 //==========================================================================
564 /// Function to describe the local dofs of the element. The ostream
565 /// specifies the output stream to which the description
566 /// is written; the string stores the currently
567 /// assembled output that is ultimately written to the
568 /// output stream by Data::describe_dofs(...); it is typically
569 /// built up incrementally as we descend through the
570 /// call hierarchy of this function when called from
571 /// Problem::describe_dofs(...)
572 //==========================================================================
574 std::ostream& out, const std::string& current_string) const
575 {
576 // Find the number of internal and external data
577 const unsigned n_internal_data = Ninternal_data;
578 const unsigned n_external_data = Nexternal_data;
579
580 // Now loop over the internal data and describe local equation numbers
581 for (unsigned i = 0; i < n_internal_data; i++)
582 {
583 // Pointer to the internal data
584 Data* const data_pt = internal_data_pt(i);
585
586 std::stringstream conversion;
587 conversion << " of Internal Data " << i << current_string;
588 std::string in(conversion.str());
589 data_pt->describe_dofs(out, in);
590 }
591
592
593 // Now loop over the external data and assign local equation numbers
594 for (unsigned i = 0; i < n_external_data; i++)
595 {
596 // Pointer to the external data
597 Data* const data_pt = external_data_pt(i);
598
599 std::stringstream conversion;
600 conversion << " of External Data " << i << current_string;
601 std::string in(conversion.str());
602 data_pt->describe_dofs(out, in);
603 }
604 }
605
606 //==========================================================================
607 /// This function loops over the internal data of the element and add
608 /// pointers to their unconstrained values to a map indexed by the global
609 /// equation number.
610 //==========================================================================
612 std::map<unsigned, double*>& map_of_value_pt)
613 {
614 // Loop over the internal data and add their data to the map
615 // The internal data are stored at the beginning of the Data_pt array
616 for (unsigned i = 0; i < Ninternal_data; i++)
617 {
618 internal_data_pt(i)->add_value_pt_to_map(map_of_value_pt);
619 }
620 }
621
622
623#ifdef OOMPH_HAS_MPI
624 //=========================================================================
625 /// Add all internal data and time history values to the vector in
626 /// the internal storage order
627 //=========================================================================
629 Vector<double>& vector_of_values)
630 {
631 for (unsigned i = 0; i < Ninternal_data; i++)
632 {
633 internal_data_pt(i)->add_values_to_vector(vector_of_values);
634 }
635 }
636
637 //========================================================================
638 /// Read all internal data and time history values from the vector
639 /// starting from index. On return the index will be
640 /// set to the value at the end of the data that has been read in
641 //========================================================================
643 const Vector<double>& vector_of_values, unsigned& index)
644 {
645 for (unsigned i = 0; i < Ninternal_data; i++)
646 {
647 internal_data_pt(i)->read_values_from_vector(vector_of_values, index);
648 }
649 }
650
651 //=========================================================================
652 /// Add all equation numbers associated with internal data
653 /// to the vector in the internal storage order
654 //=========================================================================
656 Vector<long>& vector_of_eqn_numbers)
657 {
658 for (unsigned i = 0; i < Ninternal_data; i++)
659 {
660 internal_data_pt(i)->add_eqn_numbers_to_vector(vector_of_eqn_numbers);
661 }
662 }
663
664 //=========================================================================
665 /// Read all equation numbers associated with internal data
666 /// from the vector
667 /// starting from index. On return the index will be
668 /// set to the value at the end of the data that has been read in
669 //=========================================================================
671 const Vector<long>& vector_of_eqn_numbers, unsigned& index)
672 {
673 for (unsigned i = 0; i < Ninternal_data; i++)
674 {
675 internal_data_pt(i)->read_eqn_numbers_from_vector(vector_of_eqn_numbers,
676 index);
677 }
678 }
679
680#endif
681
682
683 //====================================================================
684 /// Setup the arrays of local equation numbers for the element.
685 /// In general, this function should not need to be overloaded. Instead
686 /// the function assign_all_generic_local_eqn_numbers() will be overloaded
687 /// by different types of Element.
688 /// That said, the function is virtual so that it
689 /// may be overloaded by the user if they *really* know what they are doing.
690 //==========================================================================
692 const bool& store_local_dof_pt)
693 {
695 assign_all_generic_local_eqn_numbers(store_local_dof_pt);
697
698 // Check that no global equation numbers are repeated
699#ifdef PARANOID
700
701 std::ostringstream error_stream;
702
703 // Loop over the array of equation numbers and add to set to assess
704 // uniqueness
705 std::map<unsigned, bool> is_repeated;
706 std::set<unsigned long> set_of_global_eqn_numbers;
707 unsigned old_ndof = 0;
708 for (unsigned n = 0; n < Ndof; ++n)
709 {
710 set_of_global_eqn_numbers.insert(Eqn_number[n]);
711 if (set_of_global_eqn_numbers.size() == old_ndof)
712 {
713 error_stream << "Repeated global eqn: " << Eqn_number[n] << std::endl;
714 is_repeated[Eqn_number[n]] = true;
715 }
716 old_ndof = set_of_global_eqn_numbers.size();
717 }
718
719 // If the sizes do not match we have a repeat, throw an error
720 if (set_of_global_eqn_numbers.size() != Ndof)
721 {
722#ifdef OOMPH_HAS_MPI
723 error_stream << "Element is ";
724 if (!is_halo()) error_stream << "not ";
725 error_stream << "a halo element\n\n";
726#endif
727 error_stream << "\nLocal/lobal equation numbers: " << std::endl;
728 for (unsigned n = 0; n < Ndof; ++n)
729 {
730 error_stream << " " << n << " " << Eqn_number[n] << std::endl;
731 }
732
733 // It's helpful for debugging purposes to output more about this element
734 error_stream << std::endl << " element address is " << this << std::endl;
735
736 // Check if the repeated dofs are among the internal Data values
737 unsigned nint = this->ninternal_data();
738 error_stream << "Number of internal data " << nint << std::endl;
739 for (unsigned i = 0; i < nint; i++)
740 {
741 Data* data_pt = this->internal_data_pt(i);
742 unsigned nval = data_pt->nvalue();
743 for (unsigned j = 0; j < nval; j++)
744 {
745 int eqn_no = data_pt->eqn_number(j);
746 error_stream << "Internal dof: " << eqn_no << std::endl;
747 if (is_repeated[unsigned(eqn_no)])
748 {
749 error_stream << "Repeated internal dof: " << eqn_no << std::endl;
750 }
751 }
752 }
753
754 // Check if the repeated dofs are among the external Data values
755 unsigned next = this->nexternal_data();
756 error_stream << "Number of external data " << next << std::endl;
757 for (unsigned i = 0; i < next; i++)
758 {
759 Data* data_pt = this->external_data_pt(i);
760 unsigned nval = data_pt->nvalue();
761 for (unsigned j = 0; j < nval; j++)
762 {
763 int eqn_no = data_pt->eqn_number(j);
764 error_stream << "External dof: " << eqn_no << std::endl;
765 if (is_repeated[unsigned(eqn_no)])
766 {
767 error_stream << "Repeated external dof: " << eqn_no;
768 Node* nod_pt = dynamic_cast<Node*>(data_pt);
769 if (nod_pt != 0)
770 {
771 error_stream << " (is a node at: ";
772 unsigned ndim = nod_pt->ndim();
773 for (unsigned ii = 0; ii < ndim; ii++)
774 {
775 error_stream << nod_pt->x(ii) << " ";
776 }
777 }
778 error_stream << ")\n";
779 }
780 }
781 }
782
783
784 // If it's an element with external element check the associated
785 // Data
787 dynamic_cast<ElementWithExternalElement*>(this);
788 if (e_el_pt != 0)
789 {
790 // Check if the repeated dofs are among the external Data values
791 {
792 unsigned next = e_el_pt->nexternal_interaction_field_data();
793 error_stream << "Number of external element data " << next
794 << std::endl;
796 for (unsigned i = 0; i < next; i++)
797 {
798 unsigned nval = data_pt[i]->nvalue();
799 for (unsigned j = 0; j < nval; j++)
800 {
801 int eqn_no = data_pt[i]->eqn_number(j);
802 error_stream << "External element dof: " << eqn_no << std::endl;
803 if (is_repeated[unsigned(eqn_no)])
804 {
805 error_stream << "Repeated external element dof: " << eqn_no;
806 Node* nod_pt = dynamic_cast<Node*>(data_pt[i]);
807 if (nod_pt != 0)
808 {
809 error_stream << " (is a node at: ";
810 unsigned ndim = nod_pt->ndim();
811 for (unsigned ii = 0; ii < ndim; ii++)
812 {
813 error_stream << nod_pt->x(ii) << " ";
814 }
815 }
816 error_stream << ")\n";
817 }
818 }
819 }
820 }
821
822
823 // Check if the repeated dofs are among the external geom Data values
824 {
825 unsigned next = e_el_pt->nexternal_interaction_geometric_data();
826 error_stream << "Number of external element geom data " << next
827 << std::endl;
828 Vector<Data*> data_pt(
830 for (unsigned i = 0; i < next; i++)
831 {
832 unsigned nval = data_pt[i]->nvalue();
833 for (unsigned j = 0; j < nval; j++)
834 {
835 int eqn_no = data_pt[i]->eqn_number(j);
836 error_stream << "External element geometric dof: " << eqn_no
837 << std::endl;
838 if (is_repeated[unsigned(eqn_no)])
839 {
840 error_stream
841 << "Repeated external element geometric dof: " << eqn_no
842 << " " << data_pt[i]->value(j);
843 Node* nod_pt = dynamic_cast<Node*>(data_pt[i]);
844 if (nod_pt != 0)
845 {
846 error_stream << " (is a node at: ";
847 unsigned ndim = nod_pt->ndim();
848 for (unsigned ii = 0; ii < ndim; ii++)
849 {
850 error_stream << nod_pt->x(i) << " ";
851 }
852 error_stream << ")";
853 }
854 error_stream << "\n";
855 }
856 }
857 }
858 }
859 }
860
861 // If it's a FiniteElement then output its nodes
862 FiniteElement* f_el_pt = dynamic_cast<FiniteElement*>(this);
863 if (f_el_pt != 0)
864 {
865 unsigned n_node = f_el_pt->nnode();
866 for (unsigned n = 0; n < n_node; n++)
867 {
868 Node* nod_pt = f_el_pt->node_pt(n);
869 unsigned nval = nod_pt->nvalue();
870 for (unsigned j = 0; j < nval; j++)
871 {
872 int eqn_no = nod_pt->eqn_number(j);
873 error_stream << "Node " << n << ": Nodal dof: " << eqn_no;
874 if (eqn_no >= 0)
875 {
876 if (is_repeated[unsigned(eqn_no)])
877 {
878 error_stream << "Node " << n
879 << ": Repeated nodal dof: " << eqn_no;
880 if (j >= f_el_pt->required_nvalue(n))
881 {
882 error_stream << " (resized)";
883 }
884 error_stream << std::endl;
885 }
886 }
887 }
888 SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
889 if (solid_nod_pt != 0)
890 {
891 Data* data_pt = solid_nod_pt->variable_position_pt();
892 unsigned nval = data_pt->nvalue();
893 for (unsigned j = 0; j < nval; j++)
894 {
895 int eqn_no = data_pt->eqn_number(j);
896 error_stream << "Node " << n << ": Positional dof: " << eqn_no
897 << std::endl;
898 if (is_repeated[unsigned(eqn_no)])
899 {
900 error_stream << "Repeated positional dof: " << eqn_no << " "
901 << data_pt->value(j) << std::endl;
902 }
903 }
904 }
905 }
906
907 // Output nodal coordinates of element
908 n_node = f_el_pt->nnode();
909 for (unsigned n = 0; n < n_node; n++)
910 {
911 Node* nod_pt = f_el_pt->node_pt(n);
912 unsigned n_dim = nod_pt->ndim();
913 error_stream << "Node " << n << " at ( ";
914 for (unsigned i = 0; i < n_dim; i++)
915 {
916 error_stream << nod_pt->x(i) << " ";
917 }
918 error_stream << ")" << std::endl;
919 }
920 }
921
922
923 throw OomphLibError(
924 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
925 }
926#endif
927 }
928
929
930 //==========================================================================
931 /// This function loops over the internal and external data of the element,
932 /// adds the GLOBAL equation numbers to the local-to-global look-up scheme and
933 /// fills in the look-up schemes for the local equation
934 /// numbers.
935 /// If the boolean argument is true then pointers to the dofs will be
936 /// stored in Dof_pt
937 //==========================================================================
939 const bool& store_local_dof_pt)
940 {
941 // Find the number of internal and external data
942 const unsigned n_internal_data = Ninternal_data;
943 const unsigned n_external_data = Nexternal_data;
944 // Find the total number of data
945 const unsigned n_total_data = n_internal_data + n_external_data;
946
947 // If there is data
948 if (n_total_data > 0)
949 {
950 // Find the number of local equations assigned so far
951 unsigned local_eqn_number = ndof();
952
953 // We need to find the total number of values stored in all the
954 // internal and external data
955 // Initialise to the number stored in the first data object
956 unsigned n_total_values = Data_pt[0]->nvalue();
957 // Loop over the other data and add the number of values stored
958 for (unsigned i = 1; i < n_total_data; ++i)
959 {
960 n_total_values += Data_pt[i]->nvalue();
961 }
962
963 // If allocated delete the old storage
964 if (Data_local_eqn)
965 {
966 delete[] Data_local_eqn[0];
967 delete[] Data_local_eqn;
968 }
969
970 // If there are no values then we are done, null out the storage and
971 // return
972 if (n_total_values == 0)
973 {
974 Data_local_eqn = 0;
975 return;
976 }
977
978 // Allocate the storage for the local equation numbers
979 // The idea is that we store internal equation numbers followed by
980 // external equation numbers
981
982 // Firstly allocate pointers to the rows for each data item
983 Data_local_eqn = new int*[n_total_data];
984 // Now allocate storage for all the equation numbers
985 Data_local_eqn[0] = new int[n_total_values];
986 // Set all values to be unclassified
987 for (unsigned i = 0; i < n_total_values; ++i)
988 {
990 }
991
992 // Loop over the remaining rows and set their pointers
993 for (unsigned i = 1; i < n_total_data; ++i)
994 {
995 // Initially set the pointer to the i-th row to the pointer
996 // to the i-1th row
998 // Now increase the row pointer by the number of values
999 // stored at the i-1th data object
1000 Data_local_eqn[i] += Data_pt[i - 1]->nvalue();
1001 }
1002
1003 // A local queue to store the global equation numbers
1004 std::deque<unsigned long> global_eqn_number_queue;
1005
1006 // Now loop over the internal data and assign local equation numbers
1007 for (unsigned i = 0; i < n_internal_data; i++)
1008 {
1009 // Pointer to the internal data
1010 Data* const data_pt = internal_data_pt(i);
1011 // Find the number of values stored at the internal data
1012 unsigned n_value = data_pt->nvalue();
1013
1014 // Loop over the number of values
1015 for (unsigned j = 0; j < n_value; j++)
1016 {
1017 // Get the GLOBAL equation number
1018 long eqn_number = data_pt->eqn_number(j);
1019 // If the GLOBAL equation number is positive (a free variable)
1020 if (eqn_number >= 0)
1021 {
1022 // Add the GLOBAL equation number to the queue
1023 global_eqn_number_queue.push_back(eqn_number);
1024 // Add pointer to the dof to the queue if required
1025 if (store_local_dof_pt)
1026 {
1027 GeneralisedElement::Dof_pt_deque.push_back(data_pt->value_pt(j));
1028 }
1029 // Add the local equation number to the storage scheme
1031 // Increase the local number
1033 }
1034 else
1035 {
1036 // Set the local scheme to be pinned
1038 }
1039 }
1040 } // End of loop over internal data
1041
1042
1043 // Now loop over the external data and assign local equation numbers
1044 for (unsigned i = 0; i < n_external_data; i++)
1045 {
1046 // Pointer to the external data
1047 Data* const data_pt = external_data_pt(i);
1048 // Find the number of values stored at the external data
1049 unsigned n_value = data_pt->nvalue();
1050
1051 // Loop over the number of values
1052 for (unsigned j = 0; j < n_value; j++)
1053 {
1054 // Get the GLOBAL equation number
1055 long eqn_number = data_pt->eqn_number(j);
1056 // If the GLOBAL equation number is positive (a free variable)
1057 if (eqn_number >= 0)
1058 {
1059 // Add the GLOBAL equation number to the queue
1060 global_eqn_number_queue.push_back(eqn_number);
1061 // Add pointer to the dof to the queue if required
1062 if (store_local_dof_pt)
1063 {
1064 GeneralisedElement::Dof_pt_deque.push_back(data_pt->value_pt(j));
1065 }
1066 // Add the local equation number to the local scheme
1067 Data_local_eqn[n_internal_data + i][j] = local_eqn_number;
1068 // Increase the local number
1070 }
1071 else
1072 {
1073 // Set the local scheme to be pinned
1074 Data_local_eqn[n_internal_data + i][j] = Data::Is_pinned;
1075 }
1076 }
1077 }
1078
1079 // Now add our global equations numbers to the internal element storage
1080 add_global_eqn_numbers(global_eqn_number_queue,
1082 // Clear the memory used in the deque
1083 if (store_local_dof_pt)
1084 {
1085 std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);
1086 }
1087 }
1088 }
1089
1090
1091 //============================================================================
1092 /// This function calculates the entries of Jacobian matrix, used in
1093 /// the Newton method, associated with the internal degrees of freedom.
1094 /// It does this using finite differences,
1095 /// rather than an analytical formulation, so can be done in total generality.
1096 /// If the boolean argument is true, the finite
1097 /// differencing will be performed for all internal data, irrespective of
1098 /// the information in Data_fd. The default value (false)
1099 /// uses the information in Data_fd to selectively difference only
1100 /// certain data.
1101 //==========================================================================
1103 Vector<double>& residuals,
1104 DenseMatrix<double>& jacobian,
1105 const bool& fd_all_data)
1106 {
1107 // Locally cache the number of internal data
1108 const unsigned n_internal_data = Ninternal_data;
1109
1110 // If there aren't any internal data, then return straight away
1111 if (n_internal_data == 0)
1112 {
1113 return;
1114 }
1115
1116 // Call the update function to ensure that the element is in
1117 // a consistent state before finite differencing starts
1119
1120 // Find the number of dofs in the element
1121 const unsigned n_dof = ndof();
1122
1123 // Create newres vector
1124 Vector<double> newres(n_dof);
1125
1126 // Integer storage for local unknown
1127 int local_unknown = 0;
1128
1129 // Use the default finite difference step
1130 const double fd_step = Default_fd_jacobian_step;
1131
1132 // Loop over the internal data
1133 for (unsigned i = 0; i < n_internal_data; i++)
1134 {
1135 // If we are doing all finite differences or
1136 // the variable is included in the finite difference loop, do it
1137 if (fd_all_data || internal_data_fd(i))
1138 {
1139 // Get the number of value at the internal data
1140 const unsigned n_value = internal_data_pt(i)->nvalue();
1141
1142 // Loop over the number of values
1143 for (unsigned j = 0; j < n_value; j++)
1144 {
1145 // Get the local equation number
1146 local_unknown = internal_local_eqn(i, j);
1147 // If it's not pinned
1148 if (local_unknown >= 0)
1149 {
1150 // Get a pointer to the value of the internal data
1151 double* const value_pt = internal_data_pt(i)->value_pt(j);
1152
1153 // Save the old value of the Internal data
1154 const double old_var = *value_pt;
1155
1156 // Increment the value of the Internal data
1157 *value_pt += fd_step;
1158
1159 // Now update any dependent variables
1161
1162 // Calculate the new residuals
1163 get_residuals(newres);
1164
1165 // Do finite differences
1166 for (unsigned m = 0; m < n_dof; m++)
1167 {
1168 double sum = (newres[m] - residuals[m]) / fd_step;
1169 // Stick the entry into the Jacobian matrix
1170 jacobian(m, local_unknown) = sum;
1171 }
1172
1173 // Reset the Internal data
1174 *value_pt = old_var;
1175
1176 // Reset any dependent variables
1178 }
1179 }
1180 } // End of finite-differencing for datum (if block)
1181 }
1182
1183 // End of finite difference loop
1184 // Final reset of any dependent data
1186 }
1187
1188 //============================================================================
1189 /// This function calculates the entries of Jacobian matrix, used in
1190 /// the Newton method, associated with the external degrees of freedom.
1191 /// It does this using finite differences,
1192 /// rather than an analytical formulation, so can be done in total generality.
1193 /// If the boolean argument is true, the finite
1194 /// differencing will be performed for all internal data, irrespective of
1195 /// the information in Data_fd. The default value (false)
1196 /// uses the information in Data_fd to selectively difference only
1197 /// certain data.
1198 //==========================================================================
1200 Vector<double>& residuals,
1201 DenseMatrix<double>& jacobian,
1202 const bool& fd_all_data)
1203 {
1204 // Locally cache the number of external data
1205 const unsigned n_external_data = Nexternal_data;
1206 // If there aren't any external data, then return straight away
1207 if (n_external_data == 0)
1208 {
1209 return;
1210 }
1211
1212 // Call the update function to ensure that the element is in
1213 // a consistent state before finite differencing starts
1215
1216 // Find the number of dofs in the element
1217 const unsigned n_dof = ndof();
1218
1219 // Create newres vector
1220 Vector<double> newres(n_dof);
1221
1222 // Integer storage for local unknown
1223 int local_unknown = 0;
1224
1225 // Use the default finite difference step
1226 const double fd_step = Default_fd_jacobian_step;
1227
1228 // Loop over the external data
1229 for (unsigned i = 0; i < n_external_data; i++)
1230 {
1231 // If we are doing all finite differences or
1232 // the variable is included in the finite difference loop, do it
1233 if (fd_all_data || external_data_fd(i))
1234 {
1235 // Get the number of value at the external data
1236 const unsigned n_value = external_data_pt(i)->nvalue();
1237
1238 // Loop over the number of values
1239 for (unsigned j = 0; j < n_value; j++)
1240 {
1241 // Get the local equation number
1242 local_unknown = external_local_eqn(i, j);
1243 // If it's not pinned
1244 if (local_unknown >= 0)
1245 {
1246 // Get a pointer to the External data value
1247 double* const value_pt = external_data_pt(i)->value_pt(j);
1248
1249 // Save the old value of the External data
1250 const double old_var = *value_pt;
1251
1252 // Increment the value of the External data
1253 *value_pt += fd_step;
1254
1255 // Now update any dependent variables
1257
1258 // Calculate the new residuals
1259 get_residuals(newres);
1260
1261 // Do finite differences
1262 for (unsigned m = 0; m < n_dof; m++)
1263 {
1264 double sum = (newres[m] - residuals[m]) / fd_step;
1265 // Stick the entry into the Jacobian matrix
1266 jacobian(m, local_unknown) = sum;
1267 }
1268
1269 // Reset the External data
1270 *value_pt = old_var;
1271
1272 // Reset any dependent variables
1274 }
1275 }
1276 } // End of finite differencing for datum (if block)
1277 }
1278
1279 // End of finite difference loop
1280 // Final reset of any dependent data
1282 }
1283
1284 //=====================================================================
1285 /// Add the elemental contribution to the mass matrix
1286 /// and the residuals vector. Note that
1287 /// this function will NOT initialise the residuals vector or the mass
1288 /// matrix. It must be called after the residuals vector and
1289 /// jacobian matrix have been initialised to zero. The default
1290 /// is deliberately broken.
1291 //=====================================================================
1293 Vector<double>& residuals, DenseMatrix<double>& mass_matrix)
1294 {
1295 std::string error_message =
1296 "Empty fill_in_contribution_to_mass_matrix() has been called.\n";
1297 error_message +=
1298 "This function is called from the default implementation of\n";
1299 error_message += "get_mass_matrix();\n";
1300 error_message += "and must calculate the residuals vector and mass matrix ";
1301 error_message += "without initialising any of their entries.\n\n";
1302
1303 error_message +=
1304 "If you do not wish to use these defaults, you must overload\n";
1305 error_message += "get_mass_matrix(), which must initialise the entries\n";
1306 error_message += "of the residuals vector and mass matrix to zero.\n";
1307
1308 throw OomphLibError(
1309 error_message,
1310 "GeneralisedElement::fill_in_contribution_to_mass_matrix()",
1311 OOMPH_EXCEPTION_LOCATION);
1312 }
1313
1314 //=====================================================================
1315 /// Add the elemental contribution to the jacobian matrix,
1316 /// mass matrix and the residuals vector. Note that
1317 /// this function should NOT initialise any entries.
1318 /// It must be called after the residuals vector and
1319 /// matrices have been initialised to zero. The default is deliberately
1320 /// broken.
1321 //=====================================================================
1323 Vector<double>& residuals,
1324 DenseMatrix<double>& jacobian,
1325 DenseMatrix<double>& mass_matrix)
1326 {
1327 std::string error_message =
1328 "Empty fill_in_contribution_to_jacobian_and_mass_matrix() has been ";
1329 error_message += "called.\n";
1330 error_message +=
1331 "This function is called from the default implementation of\n";
1332 error_message += "get_jacobian_and_mass_matrix();\n";
1333 error_message +=
1334 "and must calculate the residuals vector and mass and jacobian matrices ";
1335 error_message += "without initialising any of their entries.\n\n";
1336
1337 error_message +=
1338 "If you do not wish to use these defaults, you must overload\n";
1339 error_message +=
1340 "get_jacobian_and_mass_matrix(), which must initialise the entries\n";
1341 error_message +=
1342 "of the residuals vector, jacobian and mass matrix to zero.\n";
1343
1344 throw OomphLibError(
1345 error_message,
1346 "GeneralisedElement::fill_in_contribution_to_jacobian_and_mass_matrix()",
1347 OOMPH_EXCEPTION_LOCATION);
1348 }
1349
1350
1351 //=====================================================================
1352 /// Add the elemental contribution to the derivatives of
1353 /// the residuals with respect to a parameter. This function should
1354 /// NOT initialise any entries and must be called after the entries
1355 /// have been initialised to zero
1356 /// The default implementation is deliberately broken
1357 //========================================================================
1359 double* const& parameter_pt, Vector<double>& dres_dparam)
1360 {
1361 std::string error_message =
1362 "Empty fill_in_contribution_to_dresiduals_dparameter() has been ";
1363 error_message += "called.\n";
1364 error_message +=
1365 "This function is called from the default implementation of\n";
1366 error_message += "get_dresiduals_dparameter();\n";
1367 error_message += "and must calculate the derivatives of the residuals "
1368 "vector with respect\n";
1369 error_message += "to the passed parameter ";
1370 error_message += "without initialising any values.\n\n";
1371
1372 error_message +=
1373 "If you do not wish to use these defaults, you must overload\n";
1374 error_message +=
1375 "get_dresiduals_dparameter(), which must initialise the entries\n";
1376 error_message += "of the returned vector to zero.\n";
1377
1378 error_message +=
1379 "This function is intended for use instead of the default (global) \n";
1380 error_message +=
1381 "finite-difference routine when analytic expressions are to be used\n";
1382 error_message += "in continuation and bifurcation tracking problems.\n\n";
1383 error_message += "This function is only called when the function\n";
1384 error_message +=
1385 "Problem::set_analytic_dparameter() has been called in the driver code\n";
1386
1387 throw OomphLibError(
1388 error_message,
1389 "GeneralisedElement::fill_in_contribution_to_dresiduals_dparameter()",
1390 OOMPH_EXCEPTION_LOCATION);
1391 }
1392
1393 //======================================================================
1394 /// Add the elemental contribution to the derivatives of
1395 /// the elemental Jacobian matrix
1396 /// and residuals with respect to a parameter. This function should
1397 /// NOT initialise any entries and must be called after the entries
1398 /// have been initialised to zero
1399 /// The default implementation is to use finite differences to calculate
1400 /// the derivatives.
1401 //========================================================================
1403 double* const& parameter_pt,
1404 Vector<double>& dres_dparam,
1405 DenseMatrix<double>& djac_dparam)
1406 {
1407 std::string error_message =
1408 "Empty fill_in_contribution_to_djacobian_dparameter() has been ";
1409 error_message += "called.\n";
1410 error_message +=
1411 "This function is called from the default implementation of\n";
1412 error_message += "get_djacobian_dparameter();\n";
1413 error_message +=
1414 "and must calculate the derivatives of residuals vector and jacobian ";
1415 error_message += "matrix\n";
1416 error_message += "with respect to the passed parameter ";
1417 error_message += "without initialising any values.\n\n";
1418
1419 error_message +=
1420 "If you do not wish to use these defaults, you must overload\n";
1421 error_message +=
1422 "get_djacobian_dparameter(), which must initialise the entries\n";
1423 error_message += "of the returned vector and matrix to zero.\n\n";
1424
1425 error_message +=
1426 "This function is intended for use instead of the default (global) \n";
1427 error_message +=
1428 "finite-difference routine when analytic expressions are to be used\n";
1429 error_message += "in continuation and bifurcation tracking problems.\n\n";
1430 error_message += "This function is only called when the function\n";
1431 error_message +=
1432 "Problem::set_analytic_dparameter() has been called in the driver code\n";
1433
1434
1435 throw OomphLibError(
1436 error_message,
1437 "GeneralisedElement::fill_in_contribution_to_djacobian_dparameter()",
1438 OOMPH_EXCEPTION_LOCATION);
1439 }
1440
1441
1442 //=====================================================================
1443 /// Add the elemental contribution to the derivative of the
1444 /// jacobian matrix,
1445 /// mass matrix and the residuals vector with respect to a parameter.
1446 /// Note that
1447 /// this function should NOT initialise any entries.
1448 /// It must be called after the residuals vector and
1449 /// matrices have been initialised to zero. The default is deliberately
1450 /// broken.
1451 //=====================================================================
1454 double* const& parameter_pt,
1455 Vector<double>& dres_dparam,
1456 DenseMatrix<double>& djac_dparam,
1457 DenseMatrix<double>& dmass_matrix_dparam)
1458 {
1459 std::string error_message =
1460 "Empty fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter() "
1461 "has";
1462 error_message += " been called.\n";
1463 error_message +=
1464 "This function is called from the default implementation of\n";
1465 error_message += "get_djacobian_and_dmass_matrix_dparameter();\n";
1466 error_message +=
1467 "and must calculate the residuals vector and mass and jacobian matrices ";
1468 error_message += "without initialising any of their entries.\n\n";
1469
1470 error_message +=
1471 "If you do not wish to use these defaults, you must overload\n";
1472 error_message += "get_djacobian_and_dmass_matrix_dparameter(), which must "
1473 "initialise the\n";
1474 error_message += "entries of the returned vector and matrices to zero.\n";
1475
1476
1477 error_message +=
1478 "This function is intended for use instead of the default (global) \n";
1479 error_message +=
1480 "finite-difference routine when analytic expressions are to be used\n";
1481 error_message += "in continuation and bifurcation tracking problems.\n\n";
1482 error_message += "This function is only called when the function\n";
1483 error_message +=
1484 "Problem::set_analytic_dparameter() has been called in the driver code\n";
1485
1486
1487 throw OomphLibError(error_message,
1488 "GeneralisedElement::fill_in_contribution_to_djacobian_"
1489 "and_dmass_matrix_dparameter()",
1490 OOMPH_EXCEPTION_LOCATION);
1491 }
1492
1493
1494 //========================================================================
1495 /// Fill in contribution to the product of the Hessian
1496 /// (derivative of Jacobian with
1497 /// respect to all variables) an eigenvector, Y, and
1498 /// other specified vectors, C
1499 /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}
1500 /// At the moment the dof pointer is passed in to enable
1501 /// easy calculation of finite difference default
1502 //==========================================================================
1504 Vector<double> const& Y,
1505 DenseMatrix<double> const& C,
1506 DenseMatrix<double>& product)
1507 {
1508 std::string error_message =
1509 "Empty fill_in_contribution_to_hessian_vector_products() has been ";
1510 error_message += "called.\n";
1511 error_message +=
1512 "This function is called from the default implementation of\n";
1513 error_message += "get_hessian_vector_products(); ";
1514 error_message += " and must calculate the products \n";
1515 error_message += "of the hessian matrix with the (global) ";
1516 error_message += "vectors Y and C\n";
1517 error_message += "without initialising any values.\n\n";
1518
1519 error_message +=
1520 "If you do not wish to use these defaults, you must overload\n";
1521 error_message +=
1522 "get_hessian_vector_products(), which must initialise the entries\n";
1523 error_message += "of the returned vector to zero.\n\n";
1524
1525 error_message +=
1526 "This function is intended for use instead of the default (global) \n";
1527 error_message +=
1528 "finite-difference routine when analytic expressions are to be used.\n";
1529 error_message += "This function is only called when the function\n";
1530 error_message += "Problem::set_analytic_hessian_products() has been called "
1531 "in the driver code\n";
1532
1533 throw OomphLibError(
1534 error_message,
1535 "GeneralisedElement::fill_in_contribution_to_hessian_vector_product()",
1536 OOMPH_EXCEPTION_LOCATION);
1537 }
1538
1539 //========================================================================
1540 /// Fill in the contribution to the inner products between given
1541 /// pairs of history values
1542 //==========================================================================
1544 Vector<std::pair<unsigned, unsigned>> const& history_index,
1545 Vector<double>& inner_product)
1546 {
1547 std::string error_message =
1548 "Empty fill_in_contribution_to_inner_products() has been called.\n";
1549 error_message +=
1550 "This function is called from the default implementation of\n";
1551 error_message += "get_inner_products();\n ";
1552 error_message += "It must calculate the inner products over the element\n";
1553 error_message += "of given pairs of history values\n";
1554 error_message += "without initialision of the output vector.\n\n";
1555
1556 error_message +=
1557 "If you do not wish to use these defaults, you must overload\n";
1558 error_message +=
1559 "get_inner_products(), which must initialise the entries\n";
1560 error_message += "of the returned vector to zero.\n\n";
1561
1562 throw OomphLibError(
1563 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1564 }
1565
1566 //========================================================================
1567 /// Fill in the contributions to the vectors that when taken
1568 /// as dot product with other history values give the inner product
1569 /// over the element
1570 //==========================================================================
1572 Vector<unsigned> const& history_index,
1573 Vector<Vector<double>>& inner_product_vector)
1574 {
1575 std::string error_message =
1576 "Empty fill_in_contribution_to_inner_product_vectors() has been "
1577 "called.\n";
1578 error_message +=
1579 "This function is called from the default implementation of\n";
1580 error_message += "get_inner_product_vectors(); ";
1581 error_message += " and must calculate the vectors \n";
1582 error_message += "corresponding to the input history values\n ";
1583 error_message +=
1584 "that when multiplied by other vectors of history values\n";
1585 error_message += "return the appropriate dot products.\n\n";
1586 error_message += "The output vectors must not be initialised.\n\n";
1587
1588 error_message +=
1589 "If you do not wish to use these defaults, you must overload\n";
1590 error_message +=
1591 "get_inner_products(), which must initialise the entries\n";
1592 error_message += "of the returned vector to zero.\n\n";
1593
1594 throw OomphLibError(
1595 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1596 }
1597
1598
1599 //==========================================================================
1600 /// Self-test: Have all internal values been classified as
1601 /// pinned/unpinned? Return 0 if OK.
1602 //==========================================================================
1604 {
1605 // Initialise
1606 bool passed = true;
1607
1608 // Loop over internal Data
1609 for (unsigned i = 0; i < Ninternal_data; i++)
1610 {
1611 if (internal_data_pt(i)->self_test() != 0)
1612 {
1613 passed = false;
1614 oomph_info << "\n ERROR: Failed GeneralisedElement::self_test()!"
1615 << std::endl;
1616 oomph_info << "for internal data object number: " << i << std::endl;
1617 }
1618 }
1619
1620 // Loop over external Data
1621 for (unsigned i = 0; i < Nexternal_data; i++)
1622 {
1623 if (external_data_pt(i)->self_test() != 0)
1624 {
1625 passed = false;
1626 oomph_info << "\n ERROR: Failed GeneralisedElement::self_test()!"
1627 << std::endl;
1628 oomph_info << "for external data object number: " << i << std::endl;
1629 }
1630 }
1631
1632 // Return verdict
1633 if (passed)
1634 {
1635 return 0;
1636 }
1637 else
1638 {
1639 return 1;
1640 }
1641 }
1642
1643
1644 //======================================================================
1645 /// Helper namespace for tolerances, number of iterations, etc
1646 /// used in the locate_zeta function in FiniteElement
1647 //======================================================================
1648 namespace Locate_zeta_helpers
1649 {
1650 /// Convergence tolerance for the newton solver
1651 double Newton_tolerance = 1.0e-7;
1652
1653 /// Maximum number of newton iterations
1655
1656 /// Multiplier for (zeta-based) outer radius of element used for
1657 /// deciding that point is outside element. Set to negative value
1658 /// to suppress test.
1660
1661 /// Number of points along one dimension of each element used
1662 /// to create coordinates within the element in order to see
1663 /// which has the smallest initial residual (and is therefore used
1664 /// as the initial guess in the Newton method when locating coordinate)
1665 unsigned N_local_points = 5;
1666 } // namespace Locate_zeta_helpers
1667
1668
1669 /// ////////////////////////////////////////////////////////////////////////
1670 /// ////////////////////////////////////////////////////////////////////////
1671 // Functions for finite elements
1672 /// ////////////////////////////////////////////////////////////////////////
1673 /// ////////////////////////////////////////////////////////////////////////
1674
1675 //======================================================================
1676 /// Return the i-th coordinate at local node n. Do not use
1677 /// the hanging node representation.
1678 /// NOTE: Moved to cc file because of a possible compiler bug
1679 /// in gcc (yes, really!). The move to the cc file avoids inlining
1680 /// which appears to cause problems (only) when compiled with gcc
1681 /// and -O3; offensive "illegal read" is in optimised-out section
1682 /// of code and data that is allegedly illegal is readily readable
1683 /// (by other means) just before this function is called so I can't
1684 /// really see how we could possibly be responsible for this...
1685 //======================================================================
1686 double FiniteElement::raw_nodal_position(const unsigned& n,
1687 const unsigned& i) const
1688 {
1689 /* oomph_info << "n: "<< n << std::endl; */
1690 /* oomph_info << "i: "<< i << std::endl; */
1691 /* oomph_info << "node_pt(n): "<< node_pt(n) << std::endl; */
1692 /* oomph_info << "node_pt(n)->x(i): "<< node_pt(n)->x(i) << std::endl; */
1693 // double tmp=node_pt(n)->x(i);
1694 // oomph_info << "tmp: "<< tmp << std::endl;
1695 return node_pt(n)->x(i);
1696 }
1697
1698
1699 //======================================================================
1700 /// Function to describe the local dofs of the element. The ostream
1701 /// specifies the output stream to which the description
1702 /// is written; the string stores the currently
1703 /// assembled output that is ultimately written to the
1704 /// output stream by Data::describe_dofs(...); it is typically
1705 /// built up incrementally as we descend through the
1706 /// call hierarchy of this function when called from
1707 /// Problem::describe_dofs(...)
1708 //======================================================================
1710 std::ostream& out, const std::string& current_string) const
1711 {
1712 // Call the standard finite element classification.
1713 GeneralisedElement::describe_local_dofs(out, current_string);
1714 describe_nodal_local_dofs(out, current_string);
1715 }
1716
1717 //======================================================================
1718 // Function to describe the local dofs of the element. The ostream
1719 /// specifies the output stream to which the description
1720 /// is written; the string stores the currently
1721 /// assembled output that is ultimately written to the
1722 /// output stream by Data::describe_dofs(...); it is typically
1723 /// built up incrementally as we descend through the
1724 /// call hierarchy of this function when called from
1725 /// Problem::describe_dofs(...)
1726 //======================================================================
1728 std::ostream& out, const std::string& current_string) const
1729 {
1730 // Find the number of nodes
1731 const unsigned n_node = nnode();
1732 for (unsigned n = 0; n < n_node; n++)
1733 {
1734 // Pointer to node
1735 Node* const nod_pt = node_pt(n);
1736
1737 std::stringstream conversion;
1738 conversion << " of Node " << n << current_string;
1739 std::string in(conversion.str());
1740 nod_pt->describe_dofs(out, in);
1741 } // End if for n_node
1742 } // End describe_nodal_local_dofs
1743
1744 //========================================================================
1745 /// Internal function used to check for singular or negative values
1746 /// of the determinant of the Jacobian of the mapping between local and
1747 /// global or lagrangian coordinates. Negative jacobians are allowed if the
1748 /// Accept_negative_jacobian flag is set to true.
1749 //========================================================================
1750 void FiniteElement::check_jacobian(const double& jacobian) const
1751 {
1752 // First check for a zero jacobian
1753 if (std::fabs(jacobian) < Tolerance_for_singular_jacobian)
1754 {
1756 {
1757 throw OomphLibQuietException();
1758 }
1759 else
1760 {
1761 std::ostringstream error_stream;
1762 error_stream
1763 << "Determinant of Jacobian matrix is zero --- "
1764 << "singular mapping!\nThe determinant of the "
1765 << "jacobian is " << std::fabs(jacobian)
1766 << " which is smaller than the treshold "
1768 << "You can change this treshold, by specifying "
1769 << "FiniteElement::Tolerance_for_singular_jacobian \n"
1770 << "Here are the nodal coordinates of the inverted element\n"
1771 << "in the form \n\n x,y[,z], hang_status\n\n"
1772 << "where hang_status = 1 or 2 for non-hanging or hanging\n"
1773 << "nodes respectively (useful for automatic sizing of\n"
1774 << "tecplot markers to identify the hanging nodes). \n\n";
1775 unsigned n_dim_nod = node_pt(0)->ndim();
1776 unsigned n_nod = nnode();
1777 unsigned hang_count = 0;
1778 for (unsigned j = 0; j < n_nod; j++)
1779 {
1780 for (unsigned i = 0; i < n_dim_nod; i++)
1781 {
1782 error_stream << node_pt(j)->x(i) << " ";
1783 }
1784 if (node_pt(j)->is_hanging())
1785 {
1786 error_stream << " 2";
1787 hang_count++;
1788 }
1789 else
1790 {
1791 error_stream << " 1";
1792 }
1793 error_stream << std::endl;
1794 }
1795 error_stream << std::endl << std::endl;
1796 if ((Macro_elem_pt != 0) && (0 != hang_count))
1797 {
1798 error_stream
1799 << "NOTE: Offending element is associated with a MacroElement\n"
1800 << " AND the element has hanging nodes! \n"
1801 << " If an element is thin and highly curved, the \n"
1802 << " constraints imposed by\n \n"
1803 << " (1) inter-element continuity (imposed by the hanging\n"
1804 << " node constraints) and \n\n"
1805 << " (2) the need to respect curvilinear domain boundaries\n"
1806 << " during mesh refinement (imposed by the element's\n"
1807 << " macro element mapping)\n\n"
1808 << " may be irreconcilable! \n \n"
1809 << " You may have to re-design your base mesh to avoid \n"
1810 << " the creation of thin, highly curved elements during\n"
1811 << " the refinement process.\n"
1812 << std::endl;
1813 }
1814 throw OomphLibError(
1815 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1816 }
1817 }
1818
1819
1820 // Now check for negative jacobians, if we're not allowing them (default)
1821 if ((Accept_negative_jacobian == false) && (jacobian < 0.0))
1822 {
1824 {
1825 throw OomphLibQuietException();
1826 }
1827 else
1828 {
1829 std::ostringstream error_stream;
1830 error_stream << "Negative Jacobian in transform from "
1831 << "local to global coordinates" << std::endl
1832 << " You have an inverted coordinate system!"
1833 << std::endl
1834 << std::endl;
1835 error_stream
1836 << "Here are the nodal coordinates of the inverted element\n"
1837 << "in the form \n\n x,y[,z], hang_status\n\n"
1838 << "where hang_status = 1 or 2 for non-hanging or hanging\n"
1839 << "nodes respectively (useful for automatic sizing of\n"
1840 << "tecplot markers to identify the hanging nodes). \n\n";
1841 unsigned n_dim_nod = node_pt(0)->ndim();
1842 unsigned n_nod = nnode();
1843 unsigned hang_count = 0;
1844 for (unsigned j = 0; j < n_nod; j++)
1845 {
1846 for (unsigned i = 0; i < n_dim_nod; i++)
1847 {
1848 error_stream << node_pt(j)->x(i) << " ";
1849 }
1850 if (node_pt(j)->is_hanging())
1851 {
1852 error_stream << " 2";
1853 hang_count++;
1854 }
1855 else
1856 {
1857 error_stream << " 1";
1858 }
1859 error_stream << std::endl;
1860 }
1861 error_stream << std::endl << std::endl;
1862 if ((Macro_elem_pt != 0) && (0 != hang_count))
1863 {
1864 error_stream
1865 << "NOTE: The inverted element is associated with a MacroElement\n"
1866 << " AND the element has hanging nodes! \n"
1867 << " If an element is thin and highly curved, the \n"
1868 << " constraints imposed by\n \n"
1869 << " (1) inter-element continuity (imposed by the hanging\n"
1870 << " node constraints) and \n\n"
1871 << " (2) the need to respect curvilinear domain boundaries\n"
1872 << " during mesh refinement (imposed by the element's\n"
1873 << " macro element mapping)\n\n"
1874 << " may be irreconcilable! \n \n"
1875 << " You may have to re-design your base mesh to avoid \n"
1876 << " the creation of thin, highly curved elements during\n"
1877 << " the refinement process.\n"
1878 << std::endl;
1879 }
1880
1881 error_stream
1882 << std::endl
1883 << std::endl
1884 << "If you believe that inverted elements do not cause any\n"
1885 << "problems in your specific application you can \n "
1886 << "suppress this test by: " << std::endl
1887 << " i) setting the (static) flag "
1888 << "FiniteElement::Accept_negative_jacobian to be true" << std::endl;
1889 error_stream << " ii) switching OFF the PARANOID flag" << std::endl
1890 << std::endl;
1891
1892 throw OomphLibError(
1893 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1894 }
1895 }
1896 }
1897
1898 //=========================================================================
1899 /// Internal function that is used to assemble the jacobian of the mapping
1900 /// from local coordinates (s) to the eulerian coordinates (x), given the
1901 /// derivatives of the shape functions. The entire jacobian matrix is
1902 /// constructed and this function will only work if there are the same number
1903 /// of local coordinates as global coordinates (i.e. for "bulk" elements).
1904 //=========================================================================
1906 const DShape& dpsids, DenseMatrix<double>& jacobian) const
1907 {
1908 // Locally cache the elemental dimension
1909 const unsigned el_dim = dim();
1910 // The number of shape functions must be equal to the number
1911 // of nodes (by definition)
1912 const unsigned n_shape = nnode();
1913 // The number of shape function types must be equal to the number
1914 // of nodal position types (by definition)
1915 const unsigned n_shape_type = nnodal_position_type();
1916
1917#ifdef PARANOID
1918 // Check for dimensional compatibility
1920 {
1921 std::ostringstream error_message;
1922 error_message << "Dimension mismatch" << std::endl;
1923 error_message << "The elemental dimension: " << Elemental_dimension
1924 << " must equal the nodal dimension: " << Nodal_dimension
1925 << " for the jacobian of the mapping to be well-defined"
1926 << std::endl;
1927 throw OomphLibError(
1928 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1929 }
1930#endif
1931
1932
1933 // Loop over the rows of the jacobian
1934 for (unsigned i = 0; i < el_dim; i++)
1935 {
1936 // Loop over the columns of the jacobian
1937 for (unsigned j = 0; j < el_dim; j++)
1938 {
1939 // Zero the entry
1940 jacobian(i, j) = 0.0;
1941 // Loop over the shape functions
1942 for (unsigned l = 0; l < n_shape; l++)
1943 {
1944 for (unsigned k = 0; k < n_shape_type; k++)
1945 {
1946 // Jacobian is dx_j/ds_i, which is represented by the sum
1947 // over the dpsi/ds_i of the nodal points X j
1948 // Call the Non-hanging version of positions
1949 // This is overloaded in refineable elements
1950 jacobian(i, j) += raw_nodal_position_gen(l, k, j) * dpsids(l, k, i);
1951 }
1952 }
1953 }
1954 }
1955 }
1956
1957 //=========================================================================
1958 /// Internal function that is used to assemble the jacobian of second
1959 /// derivatives of the mapping from local coordinates (s) to the
1960 /// eulerian coordinates (x), given the second derivatives of the
1961 /// shape functions.
1962 //=========================================================================
1964 const DShape& d2psids, DenseMatrix<double>& jacobian2) const
1965 {
1966 // Find the dimension of the element
1967 const unsigned el_dim = dim();
1968 // Find the number of shape functions and shape functions types
1969 // Must be equal to the number of nodes and their position types
1970 // by the definition of the shape function.
1971 const unsigned n_shape = nnode();
1972 const unsigned n_shape_type = nnodal_position_type();
1973 // Find the number of second derivatives
1974 const unsigned n_row = N2deriv[el_dim];
1975
1976 // Assemble the "jacobian" (d^2 x_j/ds_i^2) for second derivatives of
1977 // shape functions
1978 // Loop over the rows (number of second derivatives)
1979 for (unsigned i = 0; i < n_row; i++)
1980 {
1981 // Loop over the columns (element dimension
1982 for (unsigned j = 0; j < el_dim; j++)
1983 {
1984 // Zero the entry
1985 jacobian2(i, j) = 0.0;
1986 // Loop over the shape functions
1987 for (unsigned l = 0; l < n_shape; l++)
1988 {
1989 // Loop over the shape function types
1990 for (unsigned k = 0; k < n_shape_type; k++)
1991 {
1992 // Add the terms to the jacobian entry
1993 // Call the Non-hanging version of positions
1994 // This is overloaded in refineable elements
1995 jacobian2(i, j) +=
1996 raw_nodal_position_gen(l, k, j) * d2psids(l, k, i);
1997 }
1998 }
1999 }
2000 }
2001 }
2002
2003 //=====================================================================
2004 /// Assemble the covariant Eulerian base vectors and return them in the
2005 /// matrix interpolated_G. The derivatives of the shape functions with
2006 /// respect to the local coordinate should already have been calculated
2007 /// before calling this function
2008 //=====================================================================
2010 const DShape& dpsids, DenseMatrix<double>& interpolated_G) const
2011 {
2012 // Find the number of nodes and position types
2013 const unsigned n_node = nnode();
2014 const unsigned n_position_type = nnodal_position_type();
2015 // Find the dimension of the node and element
2016 const unsigned n_dim_node = nodal_dimension();
2017 const unsigned n_dim_element = dim();
2018
2019 // Loop over the dimensions and compute the entries of the
2020 // base vector matrix
2021 for (unsigned i = 0; i < n_dim_element; i++)
2022 {
2023 for (unsigned j = 0; j < n_dim_node; j++)
2024 {
2025 // Initialise the j-th component of the i-th base vector to zero
2026 interpolated_G(i, j) = 0.0;
2027 for (unsigned l = 0; l < n_node; l++)
2028 {
2029 for (unsigned k = 0; k < n_position_type; k++)
2030 {
2031 interpolated_G(i, j) +=
2032 raw_nodal_position_gen(l, k, j) * dpsids(l, k, i);
2033 }
2034 }
2035 }
2036 }
2037 }
2038
2039
2040 //============================================================================
2041 /// Zero-d specialisation of function to calculate inverse of jacobian mapping
2042 //============================================================================
2043 template<>
2044 double FiniteElement::invert_jacobian<0>(
2045 const DenseMatrix<double>& jacobian,
2046 DenseMatrix<double>& inverse_jacobian) const
2047 {
2048 // Issue a warning
2049 oomph_info << "\nWarning: You are trying to invert the jacobian for "
2050 << "a 'point' element" << std::endl
2051 << "This makes no sense and is almost certainly an error"
2052 << std::endl
2053 << std::endl;
2054
2055 // Dummy return
2056 return (1.0);
2057 }
2058
2059
2060 //===========================================================================
2061 /// One-d specialisation of function to calculate inverse of jacobian mapping
2062 //===========================================================================
2063 template<>
2064 double FiniteElement::invert_jacobian<1>(
2065 const DenseMatrix<double>& jacobian,
2066 DenseMatrix<double>& inverse_jacobian) const
2067 {
2068 // Calculate the determinant of the matrix
2069 const double det = jacobian(0, 0);
2070
2071// Report if Matrix is singular or negative
2072#ifdef PARANOID
2073 check_jacobian(det);
2074#endif
2075
2076 // Calculate the inverse --- trivial in 1D
2077 inverse_jacobian(0, 0) = 1.0 / jacobian(0, 0);
2078
2079 // Return the determinant
2080 return (det);
2081 }
2082
2083 //===========================================================================
2084 /// Two-d specialisation of function to calculate inverse of jacobian mapping
2085 //===========================================================================
2086 template<>
2087 double FiniteElement::invert_jacobian<2>(
2088 const DenseMatrix<double>& jacobian,
2089 DenseMatrix<double>& inverse_jacobian) const
2090 {
2091 // Calculate the determinant of the matrix
2092 const double det =
2093 jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0);
2094
2095// Report if Matrix is singular or negative
2096#ifdef PARANOID
2097 check_jacobian(det);
2098#endif
2099
2100 // Calculate the inverset of the 2x2 matrix
2101 inverse_jacobian(0, 0) = jacobian(1, 1) / det;
2102 inverse_jacobian(0, 1) = -jacobian(0, 1) / det;
2103 inverse_jacobian(1, 0) = -jacobian(1, 0) / det;
2104 inverse_jacobian(1, 1) = jacobian(0, 0) / det;
2105
2106 // Return the jacobian
2107 return (det);
2108 }
2109
2110 //=============================================================================
2111 /// Three-d specialisation of function to calculate inverse of jacobian
2112 /// mapping
2113 //=============================================================================
2114 template<>
2115 double FiniteElement::invert_jacobian<3>(
2116 const DenseMatrix<double>& jacobian,
2117 DenseMatrix<double>& inverse_jacobian) const
2118 {
2119 // Calculate the determinant of the matrix
2120 const double det = jacobian(0, 0) * jacobian(1, 1) * jacobian(2, 2) +
2121 jacobian(0, 1) * jacobian(1, 2) * jacobian(2, 0) +
2122 jacobian(0, 2) * jacobian(1, 0) * jacobian(2, 1) -
2123 jacobian(0, 0) * jacobian(1, 2) * jacobian(2, 1) -
2124 jacobian(0, 1) * jacobian(1, 0) * jacobian(2, 2) -
2125 jacobian(0, 2) * jacobian(1, 1) * jacobian(2, 0);
2126
2127 // Report if Matrix is singular or negative
2128#ifdef PARANOID
2129 check_jacobian(det);
2130#endif
2131
2132 // Calculate the inverse of the 3x3 matrix
2133 inverse_jacobian(0, 0) =
2134 (jacobian(1, 1) * jacobian(2, 2) - jacobian(1, 2) * jacobian(2, 1)) / det;
2135 inverse_jacobian(0, 1) =
2136 -(jacobian(0, 1) * jacobian(2, 2) - jacobian(0, 2) * jacobian(2, 1)) /
2137 det;
2138 inverse_jacobian(0, 2) =
2139 (jacobian(0, 1) * jacobian(1, 2) - jacobian(0, 2) * jacobian(1, 1)) / det;
2140 inverse_jacobian(1, 0) =
2141 -(jacobian(1, 0) * jacobian(2, 2) - jacobian(1, 2) * jacobian(2, 0)) /
2142 det;
2143 inverse_jacobian(1, 1) =
2144 (jacobian(0, 0) * jacobian(2, 2) - jacobian(0, 2) * jacobian(2, 0)) / det;
2145 inverse_jacobian(1, 2) =
2146 -(jacobian(0, 0) * jacobian(1, 2) - jacobian(0, 2) * jacobian(1, 0)) /
2147 det;
2148 inverse_jacobian(2, 0) =
2149 (jacobian(1, 0) * jacobian(2, 1) - jacobian(1, 1) * jacobian(2, 0)) / det;
2150 inverse_jacobian(2, 1) =
2151 -(jacobian(0, 0) * jacobian(2, 1) - jacobian(0, 1) * jacobian(2, 0)) /
2152 det;
2153 inverse_jacobian(2, 2) =
2154 (jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0)) / det;
2155
2156 // Return the determinant
2157 return (det);
2158 }
2159
2160 //========================================================================
2161 /// Template-free interface for inversion of the jacobian of a mapping.
2162 /// This is slightly inefficient, given that it uses a switch statement.
2163 /// It can always be overloaded in specific geometric elements,
2164 /// for efficiency reasons.
2165 //========================================================================
2167 const DenseMatrix<double>& jacobian,
2168 DenseMatrix<double>& inverse_jacobian) const
2169 {
2170 // Find the spatial dimension of the element
2171 const unsigned el_dim = dim();
2172 // Call the appropriate templated function, depending on the
2173 // element dimension
2174 switch (el_dim)
2175 {
2176 case 0:
2177 return invert_jacobian<0>(jacobian, inverse_jacobian);
2178 break;
2179 case 1:
2180 return invert_jacobian<1>(jacobian, inverse_jacobian);
2181 break;
2182 case 2:
2183 return invert_jacobian<2>(jacobian, inverse_jacobian);
2184 break;
2185 case 3:
2186 return invert_jacobian<3>(jacobian, inverse_jacobian);
2187 break;
2188 // Catch-all default case: issue warning and die
2189 default:
2190 std::ostringstream error_stream;
2191 error_stream << "Dimension of the element must be 0,1,2 or 3, not "
2192 << el_dim << std::endl;
2193
2194 throw OomphLibError(
2195 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2196 }
2197 // Dummy return for Intel compiler
2198 return 1.0;
2199 }
2200
2201 //============================================================================
2202 /// Zero-d specialisation of function to calculate the derivative of the
2203 /// jacobian of a mapping with respect to the nodal coordinates X_ij.
2204 //============================================================================
2205 template<>
2206 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<0>(
2207 const DenseMatrix<double>& jacobian,
2208 const DShape& dpsids,
2209 DenseMatrix<double>& djacobian_dX) const
2210 {
2211 // Issue a warning
2212 oomph_info << "\nWarning: You are trying to calculate derivatives of "
2213 << "a jacobian w.r.t. nodal coordinates for a 'point' "
2214 << "element." << std::endl
2215 << "This makes no sense and is almost certainly an error."
2216 << std::endl
2217 << std::endl;
2218 }
2219
2220 //===========================================================================
2221 /// One-d specialisation of function to calculate the derivative of the
2222 /// jacobian of a mapping with respect to the nodal coordinates X_ij.
2223 //===========================================================================
2224 template<>
2225 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<1>(
2226 const DenseMatrix<double>& jacobian,
2227 const DShape& dpsids,
2228 DenseMatrix<double>& djacobian_dX) const
2229 {
2230 // Determine the number of nodes in the element
2231 const unsigned n_node = nnode();
2232
2233 // Loop over nodes
2234 for (unsigned j = 0; j < n_node; j++)
2235 {
2236 djacobian_dX(0, j) = dpsids(j, 0);
2237 }
2238 }
2239
2240 //===========================================================================
2241 /// Two-d specialisation of function to calculate the derivative of the
2242 /// jacobian of a mapping with respect to the nodal coordinates X_ij.
2243 //===========================================================================
2244 template<>
2245 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<2>(
2246 const DenseMatrix<double>& jacobian,
2247 const DShape& dpsids,
2248 DenseMatrix<double>& djacobian_dX) const
2249 {
2250 // Determine the number of nodes in the element
2251 const unsigned n_node = nnode();
2252
2253 // Loop over nodes
2254 for (unsigned j = 0; j < n_node; j++)
2255 {
2256 // i=0
2257 djacobian_dX(0, j) =
2258 dpsids(j, 0) * jacobian(1, 1) - dpsids(j, 1) * jacobian(0, 1);
2259
2260 // i=1
2261 djacobian_dX(1, j) =
2262 dpsids(j, 1) * jacobian(0, 0) - dpsids(j, 0) * jacobian(1, 0);
2263 }
2264 }
2265
2266 //=============================================================================
2267 /// Three-d specialisation of function to calculate the derivative of the
2268 /// jacobian of a mapping with respect to the nodal coordinates X_ij.
2269 //=============================================================================
2270 template<>
2271 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<3>(
2272 const DenseMatrix<double>& jacobian,
2273 const DShape& dpsids,
2274 DenseMatrix<double>& djacobian_dX) const
2275 {
2276 // Determine the number of nodes in the element
2277 const unsigned n_node = nnode();
2278
2279 // Loop over nodes
2280 for (unsigned j = 0; j < n_node; j++)
2281 {
2282 // i=0
2283 djacobian_dX(0, j) =
2284 dpsids(j, 0) *
2285 (jacobian(1, 1) * jacobian(2, 2) - jacobian(1, 2) * jacobian(2, 1)) +
2286 dpsids(j, 1) *
2287 (jacobian(0, 2) * jacobian(2, 1) - jacobian(0, 1) * jacobian(2, 2)) +
2288 dpsids(j, 2) *
2289 (jacobian(0, 1) * jacobian(1, 2) - jacobian(0, 2) * jacobian(1, 1));
2290
2291 // i=1
2292 djacobian_dX(1, j) =
2293 dpsids(j, 0) *
2294 (jacobian(1, 2) * jacobian(2, 0) - jacobian(1, 0) * jacobian(2, 2)) +
2295 dpsids(j, 1) *
2296 (jacobian(0, 0) * jacobian(2, 2) - jacobian(0, 2) * jacobian(2, 0)) +
2297 dpsids(j, 2) *
2298 (jacobian(0, 2) * jacobian(1, 0) - jacobian(0, 0) * jacobian(1, 2));
2299
2300 // i=2
2301 djacobian_dX(2, j) =
2302 dpsids(j, 0) *
2303 (jacobian(1, 0) * jacobian(2, 1) - jacobian(1, 1) * jacobian(2, 0)) +
2304 dpsids(j, 1) *
2305 (jacobian(0, 1) * jacobian(2, 0) - jacobian(0, 0) * jacobian(2, 1)) +
2306 dpsids(j, 2) *
2307 (jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0));
2308 }
2309 }
2310
2311 //============================================================================
2312 /// Zero-d specialisation of function to calculate the derivative w.r.t. the
2313 /// nodal coordinates \f$ X_{pq} \f$ of the derivative of the shape functions
2314 /// w.r.t. the global eulerian coordinates \f$ x_i \f$.
2315 //============================================================================
2316 template<>
2317 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<0>(
2318 const double& det_jacobian,
2319 const DenseMatrix<double>& jacobian,
2320 const DenseMatrix<double>& djacobian_dX,
2321 const DenseMatrix<double>& inverse_jacobian,
2322 const DShape& dpsids,
2323 RankFourTensor<double>& d_dpsidx_dX) const
2324 {
2325 // Issue a warning
2326 oomph_info << "\nWarning: You are trying to calculate derivatives of "
2327 << "eulerian derivatives of shape functions w.r.t. nodal "
2328 << "coordinates for a 'point' element." << std::endl
2329 << "This makes no sense and is almost certainly an error."
2330 << std::endl
2331 << std::endl;
2332 }
2333
2334 //===========================================================================
2335 /// One-d specialisation of function to calculate the derivative w.r.t. the
2336 /// nodal coordinates \f$ X_{pq} \f$ of the derivative of the shape functions
2337 /// w.r.t. the global eulerian coordinates \f$ x_i \f$.
2338 //===========================================================================
2339 template<>
2340 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<1>(
2341 const double& det_jacobian,
2342 const DenseMatrix<double>& jacobian,
2343 const DenseMatrix<double>& djacobian_dX,
2344 const DenseMatrix<double>& inverse_jacobian,
2345 const DShape& dpsids,
2346 RankFourTensor<double>& d_dpsidx_dX) const
2347 {
2348 // Find inverse of determinant of jacobian of mapping
2349 const double inv_det_jac = 1.0 / det_jacobian;
2350
2351 // Determine the number of nodes in the element
2352 const unsigned n_node = nnode();
2353
2354 // Loop over the shape functions
2355 for (unsigned q = 0; q < n_node; q++)
2356 {
2357 // Loop over the shape functions
2358 for (unsigned j = 0; j < n_node; j++)
2359 {
2360 d_dpsidx_dX(0, q, j, 0) =
2361 -djacobian_dX(0, q) * dpsids(j, 0) * inv_det_jac * inv_det_jac;
2362 }
2363 }
2364 }
2365
2366 //===========================================================================
2367 /// Two-d specialisation of function to calculate the derivative w.r.t. the
2368 /// nodal coordinates \f$ X_{pq} \f$ of the derivative of the shape functions
2369 /// w.r.t. the global eulerian coordinates \f$ x_i \f$.
2370 //===========================================================================
2371 template<>
2372 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<2>(
2373 const double& det_jacobian,
2374 const DenseMatrix<double>& jacobian,
2375 const DenseMatrix<double>& djacobian_dX,
2376 const DenseMatrix<double>& inverse_jacobian,
2377 const DShape& dpsids,
2378 RankFourTensor<double>& d_dpsidx_dX) const
2379 {
2380 // Find inverse of determinant of jacobian of mapping
2381 const double inv_det_jac = 1.0 / det_jacobian;
2382
2383 // Determine the number of nodes in the element
2384 const unsigned n_node = nnode();
2385
2386 // Loop over the spatial dimension (this must be 2)
2387 for (unsigned p = 0; p < 2; p++)
2388 {
2389 // Loop over the shape functions
2390 for (unsigned q = 0; q < n_node; q++)
2391 {
2392 // Loop over the shape functions
2393 for (unsigned j = 0; j < n_node; j++)
2394 {
2395 // i=0
2396 d_dpsidx_dX(p, q, j, 0) =
2397 -djacobian_dX(p, q) * (inverse_jacobian(0, 0) * dpsids(j, 0) +
2398 inverse_jacobian(0, 1) * dpsids(j, 1));
2399
2400 if (p == 1)
2401 {
2402 d_dpsidx_dX(p, q, j, 0) +=
2403 dpsids(j, 0) * dpsids(q, 1) - dpsids(j, 1) * dpsids(q, 0);
2404 }
2405 d_dpsidx_dX(p, q, j, 0) *= inv_det_jac;
2406
2407 // i=1
2408 d_dpsidx_dX(p, q, j, 1) =
2409 -djacobian_dX(p, q) * (inverse_jacobian(1, 1) * dpsids(j, 1) +
2410 inverse_jacobian(1, 0) * dpsids(j, 0));
2411
2412 if (p == 0)
2413 {
2414 d_dpsidx_dX(p, q, j, 1) +=
2415 dpsids(j, 1) * dpsids(q, 0) - dpsids(j, 0) * dpsids(q, 1);
2416 }
2417 d_dpsidx_dX(p, q, j, 1) *= inv_det_jac;
2418 }
2419 }
2420 }
2421 }
2422
2423 //=============================================================================
2424 /// Three-d specialisation of function to calculate the derivative w.r.t. the
2425 /// nodal coordinates \f$ X_{pq} \f$ of the derivative of the shape functions
2426 /// w.r.t. the global eulerian coordinates \f$ x_i \f$.
2427 //=============================================================================
2428 template<>
2429 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<3>(
2430 const double& det_jacobian,
2431 const DenseMatrix<double>& jacobian,
2432 const DenseMatrix<double>& djacobian_dX,
2433 const DenseMatrix<double>& inverse_jacobian,
2434 const DShape& dpsids,
2435 RankFourTensor<double>& d_dpsidx_dX) const
2436 {
2437 // Find inverse of determinant of jacobian of mapping
2438 const double inv_det_jac = 1.0 / det_jacobian;
2439
2440 // Determine the number of nodes in the element
2441 const unsigned n_node = nnode();
2442
2443 // Loop over the spatial dimension (this must be 3)
2444 for (unsigned p = 0; p < 3; p++)
2445 {
2446 // Loop over the shape functions
2447 for (unsigned q = 0; q < n_node; q++)
2448 {
2449 // Loop over the shape functions
2450 for (unsigned j = 0; j < n_node; j++)
2451 {
2452 // Terms not multiplied by delta function
2453 for (unsigned i = 0; i < 3; i++)
2454 {
2455 d_dpsidx_dX(p, q, j, i) =
2456 -djacobian_dX(p, q) * (inverse_jacobian(i, 0) * dpsids(j, 0) +
2457 inverse_jacobian(i, 1) * dpsids(j, 1) +
2458 inverse_jacobian(i, 2) * dpsids(j, 2));
2459 }
2460
2461 // Delta function terms
2462 switch (p)
2463 {
2464 case 0:
2465 d_dpsidx_dX(p, q, j, 1) += ((dpsids(q, 2) * jacobian(1, 2) -
2466 dpsids(q, 1) * jacobian(2, 2)) *
2467 dpsids(j, 0) +
2468 (dpsids(q, 0) * jacobian(2, 2) -
2469 dpsids(q, 2) * jacobian(0, 2)) *
2470 dpsids(j, 1) +
2471 (dpsids(q, 1) * jacobian(0, 2) -
2472 dpsids(q, 0) * jacobian(1, 2)) *
2473 dpsids(j, 2));
2474
2475 d_dpsidx_dX(p, q, j, 2) += ((dpsids(q, 1) * jacobian(2, 1) -
2476 dpsids(q, 2) * jacobian(1, 1)) *
2477 dpsids(j, 0) +
2478 (dpsids(q, 2) * jacobian(0, 1) -
2479 dpsids(q, 0) * jacobian(2, 1)) *
2480 dpsids(j, 1) +
2481 (dpsids(q, 0) * jacobian(1, 1) -
2482 dpsids(q, 1) * jacobian(0, 1)) *
2483 dpsids(j, 2));
2484 break;
2485
2486 case 1:
2487
2488 d_dpsidx_dX(p, q, j, 0) += ((dpsids(q, 1) * jacobian(2, 2) -
2489 dpsids(q, 2) * jacobian(1, 2)) *
2490 dpsids(j, 0) +
2491 (dpsids(q, 2) * jacobian(0, 2) -
2492 dpsids(q, 0) * jacobian(2, 2)) *
2493 dpsids(j, 1) +
2494 (dpsids(q, 0) * jacobian(1, 2) -
2495 dpsids(q, 1) * jacobian(0, 2)) *
2496 dpsids(j, 2));
2497
2498 d_dpsidx_dX(p, q, j, 2) += ((dpsids(q, 2) * jacobian(1, 0) -
2499 dpsids(q, 1) * jacobian(2, 0)) *
2500 dpsids(j, 0) +
2501 (dpsids(q, 0) * jacobian(2, 0) -
2502 dpsids(q, 2) * jacobian(0, 0)) *
2503 dpsids(j, 1) +
2504 (dpsids(q, 1) * jacobian(0, 0) -
2505 dpsids(q, 0) * jacobian(1, 0)) *
2506 dpsids(j, 2));
2507 break;
2508
2509 case 2:
2510
2511 d_dpsidx_dX(p, q, j, 0) += ((dpsids(q, 2) * jacobian(1, 1) -
2512 dpsids(q, 1) * jacobian(2, 1)) *
2513 dpsids(j, 0) +
2514 (dpsids(q, 0) * jacobian(2, 1) -
2515 dpsids(q, 2) * jacobian(0, 1)) *
2516 dpsids(j, 1) +
2517 (dpsids(q, 1) * jacobian(0, 1) -
2518 dpsids(q, 0) * jacobian(1, 1)) *
2519 dpsids(j, 2));
2520
2521 d_dpsidx_dX(p, q, j, 1) += ((dpsids(q, 1) * jacobian(2, 0) -
2522 dpsids(q, 2) * jacobian(1, 0)) *
2523 dpsids(j, 0) +
2524 (dpsids(q, 2) * jacobian(0, 0) -
2525 dpsids(q, 0) * jacobian(2, 0)) *
2526 dpsids(j, 1) +
2527 (dpsids(q, 0) * jacobian(1, 0) -
2528 dpsids(q, 1) * jacobian(0, 0)) *
2529 dpsids(j, 2));
2530 break;
2531 }
2532
2533 // Divide through by the determinant of the Jacobian mapping
2534 for (unsigned i = 0; i < 3; i++)
2535 {
2536 d_dpsidx_dX(p, q, j, i) *= inv_det_jac;
2537 }
2538 }
2539 }
2540 }
2541 }
2542
2543 //=======================================================================
2544 /// Default value for the number of values at a node
2545 //=======================================================================
2546 const unsigned FiniteElement::Default_Initial_Nvalue = 0;
2547
2548 //======================================================================
2549 /// Default value that is used for the tolerance required when
2550 /// locating nodes via local coordinates
2551 const double FiniteElement::Node_location_tolerance = 1.0e-14;
2552
2553 //=======================================================================
2554 /// Set the default tolerance for a singular jacobian
2555 //=======================================================================
2557
2558 //======================================================================
2559 /// Set the default value of the Accept_negative_jacobian flag to be
2560 /// false
2561 //======================================================================
2563
2564
2565 //======================================================================
2566 /// Set default for static boolean to suppress output while checking
2567 /// for inverted elements
2568 //======================================================================
2570 false;
2571
2572 //========================================================================
2573 /// Static array that holds the number of rows in the second derivative
2574 /// matrix as a function of spatial dimension. In one-dimension, there is
2575 /// only one possible second derivative. In two-dimensions, there are three,
2576 /// the two second derivatives and the mixed derivatives. In three
2577 /// dimensions there are six.
2578 //=========================================================================
2579 const unsigned FiniteElement::N2deriv[4] = {0, 1, 3, 6};
2580
2581 //==========================================================================
2582 /// Calculate the mapping from local to eulerian coordinates
2583 /// assuming that the coordinates are aligned in the direction of the local
2584 /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
2585 /// The local derivatives are passed as dpsids and the jacobian and
2586 /// inverse jacobian are returned.
2587 //==========================================================================
2589 const DShape& dpsids,
2590 DenseMatrix<double>& jacobian,
2591 DenseMatrix<double>& inverse_jacobian) const
2592 {
2593 // Find the dimension of the element
2594 const unsigned el_dim = dim();
2595 // Find the number of shape functions and shape functions types
2596 // Equal to the number of nodes and their position types by definition
2597 const unsigned n_shape = nnode();
2598 const unsigned n_shape_type = nnodal_position_type();
2599
2600#ifdef PARANOID
2601 // Check for dimension compatibility
2603 {
2604 std::ostringstream error_message;
2605 error_message << "Dimension mismatch" << std::endl;
2606 error_message << "The elemental dimension: " << Elemental_dimension
2607 << " must equal the nodal dimension: " << Nodal_dimension
2608 << " for the jacobian of the mapping to be well-defined"
2609 << std::endl;
2610 throw OomphLibError(
2611 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2612 }
2613#endif
2614
2615 // In this case we assume that there are no cross terms, that is
2616 // global coordinate 0 is always in the direction of local coordinate 0
2617
2618 // Loop over the coordinates
2619 for (unsigned i = 0; i < el_dim; i++)
2620 {
2621 // Zero the jacobian and inverse jacobian entries
2622 for (unsigned j = 0; j < el_dim; j++)
2623 {
2624 jacobian(i, j) = 0.0;
2625 inverse_jacobian(i, j) = 0.0;
2626 }
2627
2628 // Loop over the shape functions
2629 for (unsigned l = 0; l < n_shape; l++)
2630 {
2631 // Loop over the types of dof
2632 for (unsigned k = 0; k < n_shape_type; k++)
2633 {
2634 // Derivatives are always dx_{i}/ds_{i}
2635 jacobian(i, i) += raw_nodal_position_gen(l, k, i) * dpsids(l, k, i);
2636 }
2637 }
2638 }
2639
2640 // Now calculate the determinant of the matrix
2641 double det = 1.0;
2642 for (unsigned i = 0; i < el_dim; i++)
2643 {
2644 det *= jacobian(i, i);
2645 }
2646
2647// Report if Matrix is singular, or negative
2648#ifdef PARANOID
2649 check_jacobian(det);
2650#endif
2651
2652 // Calculate the inverse mapping (trivial in this case)
2653 for (unsigned i = 0; i < el_dim; i++)
2654 {
2655 inverse_jacobian(i, i) = 1.0 / jacobian(i, i);
2656 }
2657
2658 // Return the value of the Jacobian
2659 return (det);
2660 }
2661
2662 //========================================================================
2663 /// Template-free interface calculating the derivative of the jacobian
2664 /// of a mapping with respect to the nodal coordinates X_ij. This is
2665 /// slightly inefficient, given that it uses a switch statement. It can
2666 /// always be overloaded in specific geometric elements, for efficiency
2667 /// reasons.
2668 //========================================================================
2670 const DenseMatrix<double>& jacobian,
2671 const DShape& dpsids,
2672 DenseMatrix<double>& djacobian_dX) const
2673 {
2674 // Determine the spatial dimension of the element
2675 const unsigned el_dim = dim();
2676
2677#ifdef PARANOID
2678 // Determine the number of nodes in the element
2679 const unsigned n_node = nnode();
2680
2681 // Check that djacobian_dX has the correct number of rows (= el_dim)
2682 if (djacobian_dX.nrow() != el_dim)
2683 {
2684 std::ostringstream error_message;
2685 error_message << "djacobian_dX must have the same number of rows as the"
2686 << "\nspatial dimension of the element.";
2687 throw OomphLibError(
2688 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2689 }
2690 // Check that djacobian_dX has the correct number of columns (= n_node)
2691 if (djacobian_dX.ncol() != n_node)
2692 {
2693 std::ostringstream error_message;
2694 error_message
2695 << "djacobian_dX must have the same number of columns as the"
2696 << "\nnumber of nodes in the element.";
2697 throw OomphLibError(
2698 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2699 }
2700#endif
2701
2702 // Call the appropriate templated function, depending on the
2703 // element dimension
2704 switch (el_dim)
2705 {
2706 case 0:
2707 dJ_eulerian_dnodal_coordinates_templated_helper<0>(
2708 jacobian, dpsids, djacobian_dX);
2709 break;
2710 case 1:
2711 dJ_eulerian_dnodal_coordinates_templated_helper<1>(
2712 jacobian, dpsids, djacobian_dX);
2713 break;
2714 case 2:
2715 dJ_eulerian_dnodal_coordinates_templated_helper<2>(
2716 jacobian, dpsids, djacobian_dX);
2717 break;
2718 case 3:
2719 dJ_eulerian_dnodal_coordinates_templated_helper<3>(
2720 jacobian, dpsids, djacobian_dX);
2721 break;
2722 // Catch-all default case: issue warning and die
2723 default:
2724 std::ostringstream error_stream;
2725 error_stream << "Dimension of the element must be 0,1,2 or 3, not "
2726 << el_dim << std::endl;
2727
2728 throw OomphLibError(
2729 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2730 }
2731 }
2732
2733 //========================================================================
2734 /// Template-free interface calculating the derivative w.r.t. the nodal
2735 /// coordinates \f$ X_{pq} \f$ of the derivative of the shape functions
2736 /// \f$ \psi_j \f$ w.r.t. the global eulerian coordinates \f$ x_i \f$.
2737 /// I.e. this function calculates
2738 /// \f[ \frac{\partial}{\partial X_{pq}} \left( \frac{\partial \psi_j}{\partial x_i} \right). \f]
2739 /// To do this it requires the determinant of the jacobian mapping, its
2740 /// derivative w.r.t. the nodal coordinates \f$ X_{pq} \f$, the inverse
2741 /// jacobian and the derivatives of the shape functions w.r.t. the local
2742 /// coordinates. The result is returned as a tensor of rank four.
2743 /// Numbering:
2744 /// d_dpsidx_dX(p,q,j,i) = \f$ \frac{\partial}{\partial X_{pq}} \left( \frac{\partial \psi_j}{\partial x_i} \right) \f$
2745 /// This function is slightly inefficient, given that it uses a switch
2746 /// statement. It can always be overloaded in specific geometric elements,
2747 /// for efficiency reasons.
2748 //========================================================================
2750 const double& det_jacobian,
2751 const DenseMatrix<double>& jacobian,
2752 const DenseMatrix<double>& djacobian_dX,
2753 const DenseMatrix<double>& inverse_jacobian,
2754 const DShape& dpsids,
2755 RankFourTensor<double>& d_dpsidx_dX) const
2756 {
2757 // Determine the spatial dimension of the element
2758 const unsigned el_dim = dim();
2759
2760#ifdef PARANOID
2761 // Determine the number of nodes in the element
2762 const unsigned n_node = nnode();
2763
2764 // Check that d_dpsidx_dX is of the correct size
2765 if (d_dpsidx_dX.nindex1() != el_dim || d_dpsidx_dX.nindex2() != n_node ||
2766 d_dpsidx_dX.nindex3() != n_node || d_dpsidx_dX.nindex4() != el_dim)
2767 {
2768 std::ostringstream error_message;
2769 error_message << "d_dpsidx_dX must be of the following dimensions:"
2770 << "\nd_dpsidx_dX(el_dim,n_node,n_node,el_dim)";
2771 throw OomphLibError(
2772 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2773 }
2774#endif
2775
2776 // Call the appropriate templated function, depending on the
2777 // element dimension
2778 switch (el_dim)
2779 {
2780 case 0:
2781 d_dshape_eulerian_dnodal_coordinates_templated_helper<0>(
2782 det_jacobian,
2783 jacobian,
2784 djacobian_dX,
2785 inverse_jacobian,
2786 dpsids,
2787 d_dpsidx_dX);
2788 break;
2789 case 1:
2790 d_dshape_eulerian_dnodal_coordinates_templated_helper<1>(
2791 det_jacobian,
2792 jacobian,
2793 djacobian_dX,
2794 inverse_jacobian,
2795 dpsids,
2796 d_dpsidx_dX);
2797 break;
2798 case 2:
2799 d_dshape_eulerian_dnodal_coordinates_templated_helper<2>(
2800 det_jacobian,
2801 jacobian,
2802 djacobian_dX,
2803 inverse_jacobian,
2804 dpsids,
2805 d_dpsidx_dX);
2806 break;
2807 case 3:
2808 d_dshape_eulerian_dnodal_coordinates_templated_helper<3>(
2809 det_jacobian,
2810 jacobian,
2811 djacobian_dX,
2812 inverse_jacobian,
2813 dpsids,
2814 d_dpsidx_dX);
2815 break;
2816 // Catch-all default case: issue warning and die
2817 default:
2818 std::ostringstream error_stream;
2819 error_stream << "Dimension of the element must be 0,1,2 or 3, not "
2820 << el_dim << std::endl;
2821
2822 throw OomphLibError(
2823 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2824 }
2825 }
2826
2827 //=====================================================================
2828 /// Convert derivatives w.r.t local coordinates to derivatives w.r.t the
2829 /// coordinates used to assemble the inverse jacobian mapping passed as
2830 /// inverse_jacobian. The derivatives passed in dbasis will be
2831 /// modified in this function from dbasisds to dbasisdX.
2832 //======================================================================
2834 const DenseMatrix<double>& inverse_jacobian, DShape& dbasis) const
2835 {
2836 // Find the number of basis functions and basis function types
2837 const unsigned n_basis = dbasis.nindex1();
2838 const unsigned n_basis_type = dbasis.nindex2();
2839 // Find the dimension of the array (Must be the elemental dimension)
2840 const unsigned n_dim = dim();
2841
2842 // Allocate temporary (stack) storage of the dimension of the element
2843 double new_derivatives[n_dim];
2844
2845 // Loop over the number of basis functions
2846 for (unsigned l = 0; l < n_basis; l++)
2847 {
2848 // Loop over type of basis functions
2849 for (unsigned k = 0; k < n_basis_type; k++)
2850 {
2851 // Loop over the coordinates
2852 for (unsigned j = 0; j < n_dim; j++)
2853 {
2854 // Zero the new transformed derivatives
2855 new_derivatives[j] = 0.0;
2856 // Do premultiplication by inverse jacobian
2857 for (unsigned i = 0; i < n_dim; i++)
2858 {
2859 new_derivatives[j] += inverse_jacobian(j, i) * dbasis(l, k, i);
2860 }
2861 }
2862 // We now copy the new derivatives into the shape functions
2863 for (unsigned j = 0; j < n_dim; j++)
2864 {
2865 dbasis(l, k, j) = new_derivatives[j];
2866 }
2867 }
2868 }
2869 }
2870
2871 //=====================================================================
2872 /// Convert derivatives w.r.t local coordinates to derivatives w.r.t the
2873 /// coordinates used to assemble the inverse jacobian mapping passed as
2874 /// inverse_jacobian, assuming that the mapping is diagonal. This merely
2875 /// saves a few loops, but is probably worth it.
2876 //======================================================================
2878 const DenseMatrix<double>& inverse_jacobian, DShape& dbasis) const
2879 {
2880 // Find the number of basis functions and basis function types
2881 const unsigned n_basis = dbasis.nindex1();
2882 const unsigned n_basis_type = dbasis.nindex2();
2883 // Find the dimension of the array (must be the elemental dimension)
2884 const unsigned n_dim = dim();
2885
2886 // Loop over the number of basis functions
2887 for (unsigned l = 0; l < n_basis; l++)
2888 {
2889 // Loop over type of basis functions
2890 for (unsigned k = 0; k < n_basis_type; k++)
2891 {
2892 // Loop over the coordinates
2893 for (unsigned j = 0; j < n_dim; j++)
2894 {
2895 dbasis(l, k, j) *= inverse_jacobian(j, j);
2896 }
2897 }
2898 }
2899 }
2900
2901 //=======================================================================
2902 /// Convert derivatives and second derivatives w.r.t local coordinates to
2903 /// derivatives w.r.t. the coordinates used to assemble the jacobian,
2904 /// inverse_jacobian and jacobian 2 passed. This must be specialised
2905 /// for each dimension, otherwise it gets very ugly
2906 /// Specialisation to one dimension.
2907 //=======================================================================
2908 template<>
2909 void FiniteElement::transform_second_derivatives_template<1>(
2910 const DenseMatrix<double>& jacobian,
2911 const DenseMatrix<double>& inverse_jacobian,
2912 const DenseMatrix<double>& jacobian2,
2913 DShape& dbasis,
2914 DShape& d2basis) const
2915 {
2916 // Find the number of basis functions and basis function types
2917 const unsigned n_basis = dbasis.nindex1();
2918 const unsigned n_basis_type = dbasis.nindex2();
2919
2920 // The second derivatives are easy, because there can be no mixed terms
2921 // Loop over number of basis functions
2922 for (unsigned l = 0; l < n_basis; l++)
2923 {
2924 // Loop over number of basis function types
2925 for (unsigned k = 0; k < n_basis_type; k++)
2926 {
2927 d2basis(l, k, 0) =
2928 d2basis(l, k, 0) / (jacobian(0, 0) * jacobian(0, 0))
2929 // Second term comes from taking d/dx of (dpsi/ds / dx/ds)
2930 - dbasis(l, k, 0) * jacobian2(0, 0) /
2931 (jacobian(0, 0) * jacobian(0, 0) * jacobian(0, 0));
2932 }
2933 }
2934
2935 // Assemble the first derivatives (do this last so that we don't
2936 // overwrite the dphids before we use it in the above)
2937 transform_derivatives(inverse_jacobian, dbasis);
2938 }
2939
2940 //=======================================================================
2941 /// Convert derivatives and second derivatives w.r.t local coordinates to
2942 /// derivatives w.r.t. the coordinates used to assemble the jacobian,
2943 /// inverse_jacobian and jacobian 2 passed. This must be specialised
2944 /// for each dimension, otherwise it gets very ugly.
2945 /// Specialisation to two spatial dimensions
2946 //=======================================================================
2947 template<>
2948 void FiniteElement::transform_second_derivatives_template<2>(
2949 const DenseMatrix<double>& jacobian,
2950 const DenseMatrix<double>& inverse_jacobian,
2951 const DenseMatrix<double>& jacobian2,
2952 DShape& dbasis,
2953 DShape& d2basis) const
2954 {
2955 // Find the number of basis functions and basis function types
2956 const unsigned n_basis = dbasis.nindex1();
2957 const unsigned n_basis_type = dbasis.nindex2();
2958
2959 // Calculate the determinant
2960 const double det =
2961 jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0);
2962
2963 // Second derivatives ... the approach taken here is to construct
2964 // dphidX/ds which can then be used to calculate the second derivatives
2965 // using the relationship d/dX = inverse_jacobian*d/ds
2966
2967 double ddetds[2];
2968
2969 ddetds[0] =
2970 jacobian2(0, 0) * jacobian(1, 1) + jacobian(0, 0) * jacobian2(2, 1) -
2971 jacobian2(0, 1) * jacobian(1, 0) - jacobian(0, 1) * jacobian2(2, 0);
2972 ddetds[1] =
2973 jacobian2(2, 0) * jacobian(1, 1) + jacobian(0, 0) * jacobian2(1, 1) -
2974 jacobian2(2, 1) * jacobian(1, 0) - jacobian(0, 1) * jacobian2(1, 0);
2975
2976 // Calculate the derivatives of the inverse jacobian wrt the local
2977 // coordinates
2978 double dinverse_jacobiands[2][2][2];
2979
2980 dinverse_jacobiands[0][0][0] =
2981 jacobian2(2, 1) / det - jacobian(1, 1) * ddetds[0] / (det * det);
2982 dinverse_jacobiands[0][1][0] =
2983 -jacobian2(0, 1) / det + jacobian(0, 1) * ddetds[0] / (det * det);
2984 dinverse_jacobiands[1][0][0] =
2985 -jacobian2(2, 0) / det + jacobian(1, 0) * ddetds[0] / (det * det);
2986 dinverse_jacobiands[1][1][0] =
2987 jacobian2(0, 0) / det - jacobian(0, 0) * ddetds[0] / (det * det);
2988
2989 dinverse_jacobiands[0][0][1] =
2990 jacobian2(1, 1) / det - jacobian(1, 1) * ddetds[1] / (det * det);
2991 dinverse_jacobiands[0][1][1] =
2992 -jacobian2(2, 1) / det + jacobian(0, 1) * ddetds[1] / (det * det);
2993 dinverse_jacobiands[1][0][1] =
2994 -jacobian2(1, 0) / det + jacobian(1, 0) * ddetds[1] / (det * det);
2995 dinverse_jacobiands[1][1][1] =
2996 jacobian2(2, 0) / det - jacobian(0, 0) * ddetds[1] / (det * det);
2997
2998 // Set up derivatives of dpsidx wrt local coordinates
2999 DShape dphidXds0(n_basis, n_basis_type, 2),
3000 dphidXds1(n_basis, n_basis_type, 2);
3001
3002 for (unsigned l = 0; l < n_basis; l++)
3003 {
3004 for (unsigned k = 0; k < n_basis_type; k++)
3005 {
3006 for (unsigned j = 0; j < 2; j++)
3007 {
3008 // Note that we can't have an inner loop because of the
3009 // convention I've chosen for the mixed derivatives!
3010 dphidXds0(l, k, j) = dinverse_jacobiands[j][0][0] * dbasis(l, k, 0) +
3011 dinverse_jacobiands[j][1][0] * dbasis(l, k, 1) +
3012 inverse_jacobian(j, 0) * d2basis(l, k, 0) +
3013 inverse_jacobian(j, 1) * d2basis(l, k, 2);
3014
3015 dphidXds1(l, k, j) = dinverse_jacobiands[j][0][1] * dbasis(l, k, 0) +
3016 dinverse_jacobiands[j][1][1] * dbasis(l, k, 1) +
3017 inverse_jacobian(j, 0) * d2basis(l, k, 2) +
3018 inverse_jacobian(j, 1) * d2basis(l, k, 1);
3019 }
3020 }
3021 }
3022
3023 // Now calculate the DShape d2phidX
3024 for (unsigned l = 0; l < n_basis; l++)
3025 {
3026 for (unsigned k = 0; k < n_basis_type; k++)
3027 {
3028 // Zero dpsidx
3029 for (unsigned j = 0; j < 3; j++)
3030 {
3031 d2basis(l, k, j) = 0.0;
3032 }
3033
3034 // Do premultiplication by inverse jacobian
3035 for (unsigned i = 0; i < 2; i++)
3036 {
3037 d2basis(l, k, i) = inverse_jacobian(i, 0) * dphidXds0(l, k, i) +
3038 inverse_jacobian(i, 1) * dphidXds1(l, k, i);
3039 }
3040 // Calculate mixed derivative term
3041 d2basis(l, k, 2) += inverse_jacobian(0, 0) * dphidXds0(l, k, 1) +
3042 inverse_jacobian(0, 1) * dphidXds1(l, k, 1);
3043 }
3044 }
3045
3046 // Assemble the first derivatives second, so that we don't
3047 // overwrite dphids
3048 transform_derivatives(inverse_jacobian, dbasis);
3049 }
3050
3051
3052 //=======================================================================
3053 /// Convert derivatives and second derivatives w.r.t local coordinates to
3054 /// derivatives w.r.t. the coordinates used to assemble the jacobian,
3055 /// inverse_jacobian and jacobian 2 passed. This must be specialised
3056 /// for each dimension, otherwise it gets very ugly
3057 /// Specialisation to one dimension.
3058 //=======================================================================
3059 template<>
3060 void FiniteElement::transform_second_derivatives_diagonal<1>(
3061 const DenseMatrix<double>& jacobian,
3062 const DenseMatrix<double>& inverse_jacobian,
3063 const DenseMatrix<double>& jacobian2,
3064 DShape& dbasis,
3065 DShape& d2basis) const
3066 {
3067 FiniteElement::transform_second_derivatives_template<1>(
3068 jacobian, inverse_jacobian, jacobian2, dbasis, d2basis);
3069 }
3070
3071
3072 //=========================================================================
3073 /// Convert second derivatives w.r.t. local coordinates to
3074 /// second derivatives w.r.t. the coordinates passed in the tensor
3075 /// coordinate. Specialised to two spatial dimension
3076 //=========================================================================
3077 template<>
3078 void FiniteElement::transform_second_derivatives_diagonal<2>(
3079 const DenseMatrix<double>& jacobian,
3080 const DenseMatrix<double>& inverse_jacobian,
3081 const DenseMatrix<double>& jacobian2,
3082 DShape& dbasis,
3083 DShape& d2basis) const
3084 {
3085 // Find the number of basis functions and basis function types
3086 const unsigned n_basis = dbasis.nindex1();
3087 const unsigned n_basis_type = dbasis.nindex2();
3088
3089 // Again we assume that there are no cross terms and that coordinate
3090 // i depends only upon local coordinate i
3091
3092 // Now calculate the DShape d2phidx
3093 for (unsigned l = 0; l < n_basis; l++)
3094 {
3095 for (unsigned k = 0; k < n_basis_type; k++)
3096 {
3097 // Second derivatives
3098 d2basis(l, k, 0) =
3099 d2basis(l, k, 0) / (jacobian(0, 0) * jacobian(0, 0)) -
3100 dbasis(l, k, 0) * jacobian2(0, 0) /
3101 (jacobian(0, 0) * jacobian(0, 0) * jacobian(0, 0));
3102
3103 d2basis(l, k, 1) =
3104 d2basis(l, k, 1) / (jacobian(1, 1) * jacobian(1, 1)) -
3105 dbasis(l, k, 1) * jacobian2(1, 1) /
3106 (jacobian(1, 1) * jacobian(1, 1) * jacobian(1, 1));
3107
3108 d2basis(l, k, 2) = d2basis(l, k, 2) / (jacobian(0, 0) * jacobian(1, 1));
3109 }
3110 }
3111
3112
3113 // Assemble the first derivatives
3114 transform_derivatives_diagonal(inverse_jacobian, dbasis);
3115 }
3116
3117
3118 //=============================================================================
3119 /// Convert derivatives and second derivatives w.r.t. local coordiantes
3120 /// to derivatives and second derivatives w.r.t. the coordinates used to
3121 /// assemble the jacobian, inverse jacobian and jacobian2 passed to the
3122 /// function. This is a template-free general interface, that should be
3123 /// overloaded for efficiency
3124 //============================================================================
3126 const DenseMatrix<double>& jacobian,
3127 const DenseMatrix<double>& inverse_jacobian,
3128 const DenseMatrix<double>& jacobian2,
3129 DShape& dbasis,
3130 DShape& d2basis) const
3131 {
3132 // Find the dimension of the element
3133 const unsigned el_dim = dim();
3134 // Choose the appropriate function based on the dimension of the element
3135 switch (el_dim)
3136 {
3137 case 1:
3138 transform_second_derivatives_template<1>(
3139 jacobian, inverse_jacobian, jacobian2, dbasis, d2basis);
3140 break;
3141 case 2:
3142 transform_second_derivatives_template<2>(
3143 jacobian, inverse_jacobian, jacobian2, dbasis, d2basis);
3144 break;
3145
3146 case 3:
3147 throw OomphLibError("Not implemented yet ... maybe one day",
3148 OOMPH_CURRENT_FUNCTION,
3149 OOMPH_EXCEPTION_LOCATION);
3150
3151 // transform_second_derivatives_template<3>(dphids,d2phids,jacobian,
3152 // inverse_jacobian,jacobian2,
3153 // dphidX,d2phidX);
3154 break;
3155 // Catch-all default case: issue warning and die
3156 default:
3157 std::ostringstream error_stream;
3158 error_stream << "Dimension of the element must be 1,2 or 3, not "
3159 << el_dim << std::endl;
3160
3161 throw OomphLibError(
3162 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3163 }
3164 }
3165
3166
3167 //======================================================================
3168 /// The destructor cleans up the memory allocated
3169 /// for storage of pointers to nodes. Internal and external data get
3170 /// wiped by the GeneralisedElement destructor; nodes get
3171 /// killed in mesh destructor.
3172 //=======================================================================
3174 {
3175 // Delete the storage allocated for the pointers to the loca nodes
3176 delete[] Node_pt;
3177
3178 // Delete the storage allocated for the nodal numbering schemes
3179 if (Nodal_local_eqn)
3180 {
3181 delete[] Nodal_local_eqn[0];
3182 delete[] Nodal_local_eqn;
3183 }
3184 }
3185
3186
3187 //==============================================================
3188 /// Get the local fraction of the node j in the element;
3189 /// vector sets its own size
3190 //==============================================================
3192 Vector<double>& s_fraction)
3193 {
3194 // Default implementation is rather dumb
3195 // Get the local coordinate and scale by local coordinate range
3196 local_coordinate_of_node(j, s_fraction);
3197 unsigned n_coordinates = s_fraction.size();
3198 for (unsigned i = 0; i < n_coordinates; i++)
3199 {
3200 s_fraction[i] = (s_fraction[i] - s_min()) / (s_max() - s_min());
3201 }
3202 }
3203
3204 //=======================================================================
3205 /// Set the spatial integration scheme and also calculate the values of the
3206 /// shape functions and their derivatives w.r.t. the local coordinates,
3207 /// placing the values into storage so that they may be re-used,
3208 /// without recalculation
3209 //=======================================================================
3211 {
3212 // Assign the integration scheme
3214 }
3215
3216 //=========================================================================
3217 /// Return the shape function stored at the ipt-th integration
3218 /// point.
3219 //=========================================================================
3220 void FiniteElement::shape_at_knot(const unsigned& ipt, Shape& psi) const
3221 {
3222 // Find the dimension of the element
3223 const unsigned el_dim = dim();
3224 // Storage for the local coordinates of the integration point
3225 Vector<double> s(el_dim);
3226 // Set the local coordinate
3227 for (unsigned i = 0; i < el_dim; i++)
3228 {
3229 s[i] = integral_pt()->knot(ipt, i);
3230 }
3231 // Get the shape function
3232 shape(s, psi);
3233 }
3234
3235 //=========================================================================
3236 /// Return the shape function and its derivatives w.r.t. the local
3237 /// coordinates at the ipt-th integration point.
3238 //=========================================================================
3239 void FiniteElement::dshape_local_at_knot(const unsigned& ipt,
3240 Shape& psi,
3241 DShape& dpsids) const
3242 {
3243 // Find the dimension of the element
3244 const unsigned el_dim = dim();
3245 // Storage for the local coordinates of the integration point
3246 Vector<double> s(el_dim);
3247 // Set the local coordinate
3248 for (unsigned i = 0; i < el_dim; i++)
3249 {
3250 s[i] = integral_pt()->knot(ipt, i);
3251 }
3252 // Get the shape function and derivatives
3253 dshape_local(s, psi, dpsids);
3254 }
3255
3256 //=========================================================================
3257 /// Calculate the shape function and its first and second derivatives
3258 /// w.r.t. local coordinates at the ipt-th integration point.
3259 /// Numbering:
3260 /// \b 1D:
3261 /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
3262 /// \b 2D:
3263 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
3264 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
3265 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
3266 /// \b 3D:
3267 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
3268 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
3269 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
3270 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
3271 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
3272 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
3273 //=========================================================================
3275 Shape& psi,
3276 DShape& dpsids,
3277 DShape& d2psids) const
3278 {
3279 // Find the dimension of the element
3280 const unsigned el_dim = dim();
3281 // Storage for the local coordinates of the integration point
3282 Vector<double> s(el_dim);
3283 // Set the local coordinate
3284 for (unsigned i = 0; i < el_dim; i++)
3285 {
3286 s[i] = integral_pt()->knot(ipt, i);
3287 }
3288 // Get the shape function and first and second derivatives
3289 d2shape_local(s, psi, dpsids, d2psids);
3290 }
3291
3292 //=========================================================================
3293 /// Compute the geometric shape functions and also
3294 /// first derivatives w.r.t. global coordinates at local coordinate s;
3295 /// Returns Jacobian of mapping from global to local coordinates.
3296 /// Most general form of the function, but may be over-loaded, if desired
3297 //=========================================================================
3299 Shape& psi,
3300 DShape& dpsi) const
3301 {
3302 // Find the element dimension
3303 const unsigned el_dim = dim();
3304
3305 // Get the values of the shape functions and their local derivatives
3306 // Temporarily stored in dpsi
3307 dshape_local(s, psi, dpsi);
3308
3309 // Allocate memory for the inverse jacobian
3310 DenseMatrix<double> inverse_jacobian(el_dim);
3311 // Now calculate the inverse jacobian
3312 const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
3313
3314 // Now set the values of the derivatives to be dpsidx
3315 transform_derivatives(inverse_jacobian, dpsi);
3316 // Return the determinant of the jacobian
3317 return det;
3318 }
3319
3320 //========================================================================
3321 /// Compute the geometric shape functions and also first
3322 /// derivatives w.r.t. global coordinates at integration point ipt.
3323 /// Most general form of function, but may be over-loaded if desired
3324 //========================================================================
3325 double FiniteElement::dshape_eulerian_at_knot(const unsigned& ipt,
3326 Shape& psi,
3327 DShape& dpsi) const
3328 {
3329 // Find the element dimension
3330 const unsigned el_dim = dim();
3331
3332 // Get the values of the shape function and local derivatives
3333 // Temporarily store it in dpsi
3334 dshape_local_at_knot(ipt, psi, dpsi);
3335
3336 // Allocate memory for the inverse jacobian
3337 DenseMatrix<double> inverse_jacobian(el_dim);
3338 // Now calculate the inverse jacobian
3339 const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
3340
3341 // Now set the values of the derivatives to dpsidx
3342 transform_derivatives(inverse_jacobian, dpsi);
3343 // Return the determinant of the jacobian
3344 return det;
3345 }
3346
3347
3348 //========================================================================
3349 /// Compute the geometric shape functions (psi) at integration point
3350 /// ipt. Return the determinant of the jacobian of the mapping (detJ).
3351 /// Additionally calculate the derivatives of "detJ" w.r.t. the
3352 /// nodal coordinates.
3353 //========================================================================
3355 const unsigned& ipt, Shape& psi, DenseMatrix<double>& djacobian_dX) const
3356 {
3357 // Find the element dimension
3358 const unsigned el_dim = dim();
3359
3360 // Get the values of the shape function and local derivatives
3361 unsigned nnod = nnode();
3362 DShape dpsi(nnod, el_dim);
3363 dshape_local_at_knot(ipt, psi, dpsi);
3364
3365 // Allocate memory for the jacobian and the inverse of the jacobian
3366 DenseMatrix<double> jacobian(el_dim), inverse_jacobian(el_dim);
3367
3368 // Now calculate the inverse jacobian
3369 const double det =
3370 local_to_eulerian_mapping(dpsi, jacobian, inverse_jacobian);
3371
3372 // Calculate the derivative of the jacobian w.r.t. nodal coordinates
3373 // Note: must call this before "transform_derivatives(...)" since this
3374 // function requires dpsids rather than dpsidx
3375 dJ_eulerian_dnodal_coordinates(jacobian, dpsi, djacobian_dX);
3376
3377 // Return the determinant of the jacobian
3378 return det;
3379 }
3380
3381
3382 //========================================================================
3383 /// Compute the geometric shape functions (psi) and first
3384 /// derivatives w.r.t. global coordinates (dpsidx) at integration point
3385 /// ipt. Return the determinant of the jacobian of the mapping (detJ).
3386 /// Additionally calculate the derivatives of both "detJ" and "dpsidx"
3387 /// w.r.t. the nodal coordinates.
3388 /// Most general form of function, but may be over-loaded if desired.
3389 //========================================================================
3391 const unsigned& ipt,
3392 Shape& psi,
3393 DShape& dpsi,
3394 DenseMatrix<double>& djacobian_dX,
3395 RankFourTensor<double>& d_dpsidx_dX) const
3396 {
3397 // Find the element dimension
3398 const unsigned el_dim = dim();
3399
3400 // Get the values of the shape function and local derivatives
3401 // Temporarily store in dpsi
3402 dshape_local_at_knot(ipt, psi, dpsi);
3403
3404 // Allocate memory for the jacobian and the inverse of the jacobian
3405 DenseMatrix<double> jacobian(el_dim), inverse_jacobian(el_dim);
3406
3407 // Now calculate the inverse jacobian
3408 const double det =
3409 local_to_eulerian_mapping(dpsi, jacobian, inverse_jacobian);
3410
3411 // Calculate the derivative of the jacobian w.r.t. nodal coordinates
3412 // Note: must call this before "transform_derivatives(...)" since this
3413 // function requires dpsids rather than dpsidx
3414 dJ_eulerian_dnodal_coordinates(jacobian, dpsi, djacobian_dX);
3415
3416 // Calculate the derivative of dpsidx w.r.t. nodal coordinates
3417 // Note: this function also requires dpsids rather than dpsidx
3419 det, jacobian, djacobian_dX, inverse_jacobian, dpsi, d_dpsidx_dX);
3420
3421 // Now set the values of the derivatives to dpsidx
3422 transform_derivatives(inverse_jacobian, dpsi);
3423
3424 // Return the determinant of the jacobian
3425 return det;
3426 }
3427
3428
3429 //===========================================================================
3430 /// Compute the geometric shape functions and also first
3431 /// and second derivatives w.r.t. global coordinates at local coordinate s;
3432 /// Also returns Jacobian of mapping from global to local coordinates.
3433 /// Numbering:
3434 /// \b 1D:
3435 /// d2psidx(i,0) = \f$ d^2 \psi_j / d x^2 \f$
3436 /// \b 2D:
3437 /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$
3438 /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$
3439 /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$
3440 /// \b 3D:
3441 /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$
3442 /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$
3443 /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_2^2 \f$
3444 /// d2psidx(i,3) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$
3445 /// d2psidx(i,4) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_2 \f$
3446 /// d2psidx(i,5) = \f$ \partial^2 \psi_j / \partial x_1 \partial x_2 \f$
3447 //===========================================================================
3449 Shape& psi,
3450 DShape& dpsi,
3451 DShape& d2psi) const
3452 {
3453 // Find the values of the indices of the shape functions
3454 // Locally cached.
3455 // Find the element dimension
3456 const unsigned el_dim = dim();
3457 // Find the number of second derivatives required
3458 const unsigned n_deriv = N2deriv[el_dim];
3459
3460 // Get the values of the shape function and local derivatives
3461 d2shape_local(s, psi, dpsi, d2psi);
3462
3463 // Allocate memory for the jacobian and inverse jacobian
3464 DenseMatrix<double> jacobian(el_dim), inverse_jacobian(el_dim);
3465 // Calculate the jacobian and inverse jacobian
3466 const double det =
3467 local_to_eulerian_mapping(dpsi, jacobian, inverse_jacobian);
3468
3469 // Allocate memory for the jacobian of second derivatives
3470 DenseMatrix<double> jacobian2(n_deriv, el_dim);
3471 // Assemble the jacobian of second derivatives
3472 assemble_local_to_eulerian_jacobian2(d2psi, jacobian2);
3473
3474 // Now set the value of the derivatives
3476 jacobian, inverse_jacobian, jacobian2, dpsi, d2psi);
3477 // Return the determinant of the mapping
3478 return det;
3479 }
3480
3481 //===========================================================================
3482 /// Compute the geometric shape functions and also first
3483 /// and second derivatives w.r.t. global coordinates at ipt-th integration
3484 /// point
3485 /// Returns Jacobian of mapping from global to local coordinates.
3486 /// This is the most general version, may be overloaded, if desired.
3487 /// Numbering:
3488 /// \b 1D:
3489 /// d2psidx(i,0) = \f$ d^2 \psi_j / d x^2 \f$
3490 /// \b 2D:
3491 /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$
3492 /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$
3493 /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$
3494 /// \b 3D:
3495 /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$
3496 /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$
3497 /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_2^2 \f$
3498 /// d2psidx(i,3) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$
3499 /// d2psidx(i,4) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_2 \f$
3500 /// d2psidx(i,5) = \f$ \partial^2 \psi_j / \partial x_1 \partial x_2 \f$
3501 //==========================================================================
3503 Shape& psi,
3504 DShape& dpsi,
3505 DShape& d2psi) const
3506 {
3507 // Find the values of the indices of the shape functions
3508 // Locally cached
3509 // Find the element dimension
3510 const unsigned el_dim = dim();
3511 // Find the number of second derivatives required
3512 const unsigned n_deriv = N2deriv[el_dim];
3513
3514 // Get the values of the shape function and local derivatives
3515 d2shape_local_at_knot(ipt, psi, dpsi, d2psi);
3516
3517 // Allocate memory for the jacobian and inverse jacobian
3518 DenseMatrix<double> jacobian(el_dim), inverse_jacobian(el_dim);
3519 // Calculate the jacobian and inverse jacobian
3520 const double det =
3521 local_to_eulerian_mapping(dpsi, jacobian, inverse_jacobian);
3522
3523 // Allocate memory for the jacobian of second derivatives
3524 DenseMatrix<double> jacobian2(n_deriv, el_dim);
3525 // Assemble the jacobian of second derivatives
3526 assemble_local_to_eulerian_jacobian2(d2psi, jacobian2);
3527
3528 // Now set the value of the derivatives
3530 jacobian, inverse_jacobian, jacobian2, dpsi, d2psi);
3531 // Return the determinant of the mapping
3532 return det;
3533 }
3534
3535
3536 //==========================================================================
3537 /// This function loops over the nodal data of the element, adds the
3538 /// GLOBAL equation numbers to the local-to-global look-up scheme and
3539 /// fills in the Nodal_local_eqn look-up scheme for the local equation
3540 /// numbers
3541 /// If the boolean argument is true then pointers to the dofs will be
3542 /// stored in Dof_pt
3543 //==========================================================================
3545 const bool& store_local_dof_pt)
3546 {
3547 // Find the number of nodes
3548 const unsigned n_node = nnode();
3549 // If there are nodes
3550 if (n_node > 0)
3551 {
3552 // Find the number of local equations assigned so far
3553 unsigned local_eqn_number = ndof();
3554
3555 // We need to find the total number of values stored at the node
3556 // Initialise to the number of values stored at the first node
3557 unsigned n_total_values = node_pt(0)->nvalue();
3558 // Loop over the other nodes and add the values stored
3559 for (unsigned n = 1; n < n_node; n++)
3560 {
3561 n_total_values += node_pt(n)->nvalue();
3562 }
3563
3564 // If allocated delete the old storage
3565 if (Nodal_local_eqn)
3566 {
3567 delete[] Nodal_local_eqn[0];
3568 delete[] Nodal_local_eqn;
3569 }
3570
3571 // If there are no values, we are done, null out the storage and
3572 // return
3573 if (n_total_values == 0)
3574 {
3575 Nodal_local_eqn = 0;
3576 return;
3577 }
3578
3579 // Resize the storage for the nodal local equation numbers
3580 // Firstly allocate pointers to rows for each node
3581 Nodal_local_eqn = new int*[n_node];
3582 // Now allocate storage for the equation numbers
3583 Nodal_local_eqn[0] = new int[n_total_values];
3584 // initially all local equations are unclassified
3585 for (unsigned i = 0; i < n_total_values; i++)
3586 {
3588 }
3589
3590 // Loop over the remaining rows and set their pointers
3591 for (unsigned n = 1; n < n_node; ++n)
3592 {
3593 // Initially set the pointer to the i-th row to the pointer
3594 // to the i-1th row
3595 Nodal_local_eqn[n] = Nodal_local_eqn[n - 1];
3596 // Now increase the row pointer by the number of values
3597 // stored at the i-1th node
3598 Nodal_local_eqn[n] += Node_pt[n - 1]->nvalue();
3599 }
3600
3601
3602 // A local queue to store the global equation numbers
3603 std::deque<unsigned long> global_eqn_number_queue;
3604
3605 // Now loop over the nodes again and assign local equation numbers
3606 for (unsigned n = 0; n < n_node; n++)
3607 {
3608 // Pointer to node
3609 Node* const nod_pt = node_pt(n);
3610
3611 // Find the number of values stored at the node
3612 unsigned n_value = nod_pt->nvalue();
3613
3614 // Loop over the number of values
3615 for (unsigned j = 0; j < n_value; j++)
3616 {
3617 // Get the GLOBAL equation number
3618 long eqn_number = nod_pt->eqn_number(j);
3619 // If the GLOBAL equation number is positive (a free variable)
3620 if (eqn_number >= 0)
3621 {
3622 // Add the GLOBAL equation number to the queue
3623 global_eqn_number_queue.push_back(eqn_number);
3624 // Add pointer to the dof to the queue if required
3625 if (store_local_dof_pt)
3626 {
3627 GeneralisedElement::Dof_pt_deque.push_back(nod_pt->value_pt(j));
3628 }
3629 // Add the local equation number to the local scheme
3631 // Increase the local number
3633 }
3634 else
3635 {
3636 // Set the local scheme to be pinned
3638 }
3639 }
3640 }
3641
3642 // Now add our global equations numbers to the internal element storage
3643 add_global_eqn_numbers(global_eqn_number_queue,
3645 // Clear the memory used in the deque
3646 if (store_local_dof_pt)
3647 {
3648 std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);
3649 }
3650 }
3651 }
3652
3653
3654 //============================================================================
3655 /// This function calculates the entries of Jacobian matrix, used in
3656 /// the Newton method, associated with the nodal degrees of freedom.
3657 /// It does this using finite differences,
3658 /// rather than an analytical formulation, so can be done in total generality.
3659 //==========================================================================
3661 Vector<double>& residuals, DenseMatrix<double>& jacobian)
3662 {
3663 // Find the number of nodes
3664 const unsigned n_node = nnode();
3665 // If there aren't any nodes, then return straight awayy
3666 if (n_node == 0)
3667 {
3668 return;
3669 }
3670
3671 // Call the update function to ensure that the element is in
3672 // a consistent state before finite differencing starts
3674
3675 // Find the number of dofs in the element
3676 const unsigned n_dof = ndof();
3677 // Create newres vector
3678 Vector<double> newres(n_dof);
3679
3680 // Integer storage for local unknown
3681 int local_unknown = 0;
3682
3683 // Use the default finite difference step
3684 const double fd_step = Default_fd_jacobian_step;
3685
3686 // Loop over the nodes
3687 for (unsigned n = 0; n < n_node; n++)
3688 {
3689 // Get the number of values stored at the node
3690 const unsigned n_value = node_pt(n)->nvalue();
3691
3692 // Loop over the number of values
3693 for (unsigned i = 0; i < n_value; i++)
3694 {
3695 // Get the local equation number
3696 local_unknown = nodal_local_eqn(n, i);
3697 // If it's not pinned
3698 if (local_unknown >= 0)
3699 {
3700 // Store a pointer to the nodal data value
3701 double* const value_pt = node_pt(n)->value_pt(i);
3702
3703 // Save the old value of the Nodal data
3704 const double old_var = *value_pt;
3705
3706 // Increment the value of the Nodal data
3707 *value_pt += fd_step;
3708
3709 // Now update any dependent variables
3711
3712 // Calculate the new residuals
3713 get_residuals(newres);
3714
3715 // Do finite differences
3716 for (unsigned m = 0; m < n_dof; m++)
3717 {
3718 double sum = (newres[m] - residuals[m]) / fd_step;
3719 // Stick the entry into the Jacobian matrix
3720 jacobian(m, local_unknown) = sum;
3721 }
3722
3723 // Reset the Nodal data
3724 *value_pt = old_var;
3725
3726 // Reset any dependent variables
3728 }
3729 }
3730 }
3731
3732 // End of finite difference loop
3733 // Final reset of any dependent data
3735 }
3736
3737
3738 //=======================================================================
3739 /// Compute derivatives of elemental residual vector with respect
3740 /// to nodal coordinates. Default implementation by FD can be overwritten
3741 /// for specific elements.
3742 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
3743 /// /=======================================================================
3745 RankThreeTensor<double>& dresidual_dnodal_coordinates)
3746 {
3747 // Number of nodes
3748 unsigned n_nod = nnode();
3749
3750 // If the element has no nodes (why??!!) return straightaway
3751 if (n_nod == 0) return;
3752
3753 // Get dimension from first node
3754 unsigned dim_nod = node_pt(0)->ndim();
3755
3756 // Number of dofs
3757 unsigned n_dof = ndof();
3758
3759 // Get reference residual
3760 Vector<double> res(n_dof);
3761 Vector<double> res_pls(n_dof);
3762 get_residuals(res);
3763
3764 // FD step
3766
3767 // Do FD loop
3768 for (unsigned j = 0; j < n_nod; j++)
3769 {
3770 // Get node
3771 Node* nod_pt = node_pt(j);
3772
3773 // Loop over coordinate directions
3774 for (unsigned i = 0; i < dim_nod; i++)
3775 {
3776 // Make backup
3777 double backup = nod_pt->x(i);
3778
3779 // Do FD step. No node update required as we're
3780 // attacking the coordinate directly...
3781 nod_pt->x(i) += eps_fd;
3782
3783 // Perform auxiliary node update function
3785
3786 // Get advanced residual
3787 get_residuals(res_pls);
3788
3789 // Fill in FD entries [Loop order is "wrong" here as l is the
3790 // slow index but this is in a function that's costly anyway
3791 // and gives us the fastest loop outside where these tensor
3792 // is actually used.]
3793 for (unsigned l = 0; l < n_dof; l++)
3794 {
3795 dresidual_dnodal_coordinates(l, i, j) =
3796 (res_pls[l] - res[l]) / eps_fd;
3797 }
3798
3799 // Reset coordinate. No node update required as we're
3800 // attacking the coordinate directly...
3801 nod_pt->x(i) = backup;
3802
3803 // Perform auxiliary node update function
3805 }
3806 }
3807 }
3808
3809
3810 //===============================================================
3811 /// Return the number of the node located at *node_pt
3812 /// if this node is in the element, else return \f$ -1 \f$
3813 //===============================================================
3814 int FiniteElement::get_node_number(Node* const& global_node_pt) const
3815 {
3816 // Initialise the number to -1
3817 int number = -1;
3818 // Find the number of nodes
3819 unsigned n_node = nnode();
3820#ifdef PARANOID
3821 {
3822 // Error check that node does not appear in element more than once
3823 unsigned count = 0;
3824 // Storage for the local node numbers of the element
3825 std::vector<int> local_node_number;
3826 // Loop over the nodes
3827 for (unsigned i = 0; i < n_node; i++)
3828 {
3829 // If the node is present increase the counter
3830 // and store the local node number
3831 if (node_pt(i) == global_node_pt)
3832 {
3833 ++count;
3834 local_node_number.push_back(i);
3835 }
3836 }
3837
3838 // If the node appears more than once, complain
3839 if (count > 1)
3840 {
3841 std::ostringstream error_stream;
3842 error_stream << "Node " << global_node_pt << " appears " << count
3843 << " times in an element." << std::endl
3844 << "In positions: ";
3845 for (std::vector<int>::iterator it = local_node_number.begin();
3846 it != local_node_number.end();
3847 ++it)
3848 {
3849 error_stream << *it << " ";
3850 }
3851 error_stream << std::endl << "That seems very odd." << std::endl;
3852
3853 throw OomphLibError(
3854 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3855 }
3856 }
3857#endif
3858
3859 // Loop over the nodes
3860 for (unsigned i = 0; i < n_node; i++)
3861 {
3862 // If the passed node pointer is present in the element
3863 // set number to be its local node number
3864 if (node_pt(i) == global_node_pt)
3865 {
3866 number = i;
3867 break;
3868 }
3869 }
3870
3871 // Return the node number
3872 return number;
3873 }
3874
3875
3876 //==========================================================================
3877 /// If there is a node at the local coordinate, s, return the pointer
3878 /// to the node. If not return 0. Note that this is a default, brute
3879 /// force implementation, can almost certainly be made more efficient for
3880 /// specific elements.
3881 //==========================================================================
3883 const Vector<double>& s) const
3884 {
3885 // Locally cache the tolerance
3886 const double tol = Node_location_tolerance;
3887 Vector<double> s_node;
3888 // Locally cache the member data
3889 const unsigned el_dim = Elemental_dimension;
3890 const unsigned n_node = Nnode;
3891 // Loop over the nodes
3892 for (unsigned n = 0; n < n_node; n++)
3893 {
3894 bool Match = true;
3895 // Find the local coordinate of the node
3896 local_coordinate_of_node(n, s_node);
3897 for (unsigned i = 0; i < el_dim; i++)
3898 {
3899 // Calculate the difference between coordinates
3900 // and if it's bigger than our tolerance
3901 // break out of the (inner)loop
3902 if (std::fabs(s[i] - s_node[i]) > tol)
3903 {
3904 Match = false;
3905 break;
3906 }
3907 }
3908 // If we haven't complained then we have a match
3909 if (Match)
3910 {
3911 return node_pt(n);
3912 }
3913 }
3914 // If we get here, we have no match
3915 return 0;
3916 }
3917
3918 //======================================================================
3919 /// Compute centre of gravity of all nodes and radius of node that
3920 /// is furthest from it. Used to assess approximately if a point
3921 /// is likely to be contained with an element in locate_zeta-like
3922 /// operations. NOTE: All computed in terms of zeta!
3923 //======================================================================
3925 Vector<double>& cog, double& max_radius) const
3926 {
3927 // Initialise
3928 cog.resize(Elemental_dimension);
3929 max_radius = 0.0;
3930
3931 // Get cog
3932 unsigned nnod = nnode();
3933 for (unsigned j = 0; j < nnod; j++)
3934 {
3935 for (unsigned i = 0; i < Elemental_dimension; i++)
3936 {
3937 cog[i] += zeta_nodal(j, 0, i);
3938 }
3939 }
3940 for (unsigned i = 0; i < Elemental_dimension; i++)
3941 {
3942 cog[i] /= double(nnod);
3943 }
3944
3945 // Get max distance
3946 for (unsigned j = 0; j < nnod; j++)
3947 {
3948 double dist_squared = 0.0;
3949 for (unsigned i = 0; i < Elemental_dimension; i++)
3950 {
3951 dist_squared +=
3952 (cog[i] - zeta_nodal(j, 0, i)) * (cog[i] - zeta_nodal(j, 0, i));
3953 }
3954 if (dist_squared > max_radius) max_radius = dist_squared;
3955 }
3956 max_radius = sqrt(max_radius);
3957 }
3958
3959 //======================================================================
3960 /// Return FE interpolated coordinate x[i] at local coordinate s
3961 //======================================================================
3963 const unsigned& i) const
3964 {
3965 // Find the number of nodes
3966 const unsigned n_node = nnode();
3967 // Find the number of positional types
3968 const unsigned n_position_type = nnodal_position_type();
3969 // Assign storage for the local shape function
3970 Shape psi(n_node, n_position_type);
3971 // Find the values of shape function
3972 shape(s, psi);
3973
3974 // Initialise value of x
3975 double interpolated_x = 0.0;
3976 // Loop over the local nodes
3977 for (unsigned l = 0; l < n_node; l++)
3978 {
3979 // Loop over the number of dofs
3980 for (unsigned k = 0; k < n_position_type; k++)
3981 {
3982 interpolated_x += nodal_position_gen(l, k, i) * psi(l, k);
3983 }
3984 }
3985
3986 return (interpolated_x);
3987 }
3988
3989 //=========================================================================
3990 /// Return FE interpolated coordinate x[i] at local coordinate s
3991 /// at previous timestep t (t=0: present; t>0: previous timestep)
3992 //========================================================================
3993 double FiniteElement::interpolated_x(const unsigned& t,
3994 const Vector<double>& s,
3995 const unsigned& i) const
3996 {
3997 // Find the number of nodes
3998 const unsigned n_node = nnode();
3999 // Find the number of positional types
4000 const unsigned n_position_type = nnodal_position_type();
4001
4002 // Assign storage for the local shape function
4003 Shape psi(n_node, n_position_type);
4004 // Find the values of shape function
4005 shape(s, psi);
4006
4007 // Initialise value of x
4008 double interpolated_x = 0.0;
4009 // Loop over the local nodes
4010 for (unsigned l = 0; l < n_node; l++)
4011 {
4012 // Loop over the number of dofs
4013 for (unsigned k = 0; k < n_position_type; k++)
4014 {
4015 interpolated_x += nodal_position_gen(t, l, k, i) * psi(l, k);
4016 }
4017 }
4018
4019 return (interpolated_x);
4020 }
4021
4022 //=======================================================================
4023 /// Return FE interpolated position x[] at local coordinate s as Vector
4024 //=======================================================================
4026 Vector<double>& x) const
4027 {
4028 // Find the number of nodes
4029 const unsigned n_node = nnode();
4030 // Find the number of positional types
4031 const unsigned n_position_type = nnodal_position_type();
4032 // Find the dimension stored in the node
4033 const unsigned nodal_dim = nodal_dimension();
4034
4035 // Assign storage for the local shape function
4036 Shape psi(n_node, n_position_type);
4037 // Find the values of shape function
4038 shape(s, psi);
4039
4040 // Loop over the dimensions
4041 for (unsigned i = 0; i < nodal_dim; i++)
4042 {
4043 // Initilialise value of x[i] to zero
4044 x[i] = 0.0;
4045 // Loop over the local nodes
4046 for (unsigned l = 0; l < n_node; l++)
4047 {
4048 // Loop over the number of dofs
4049 for (unsigned k = 0; k < n_position_type; k++)
4050 {
4051 x[i] += nodal_position_gen(l, k, i) * psi(l, k);
4052 }
4053 }
4054 }
4055 }
4056
4057 //==========================================================================
4058 /// Return FE interpolated position x[] at local coordinate s
4059 /// at previous timestep t as Vector (t=0: present; t>0: previous timestep)
4060 //==========================================================================
4061 void FiniteElement::interpolated_x(const unsigned& t,
4062 const Vector<double>& s,
4063 Vector<double>& x) const
4064 {
4065 // Find the number of nodes
4066 const unsigned n_node = nnode();
4067 // Find the number of positional types
4068 const unsigned n_position_type = nnodal_position_type();
4069 // Find the dimensions of the nodes
4070 const unsigned nodal_dim = nodal_dimension();
4071
4072 // Assign storage for the local shape function
4073 Shape psi(n_node, n_position_type);
4074 // Find the values of shape function
4075 shape(s, psi);
4076
4077 // Loop over the dimensions
4078 for (unsigned i = 0; i < nodal_dim; i++)
4079 {
4080 // Initilialise value of x[i] to zero
4081 x[i] = 0.0;
4082 // Loop over the local nodes
4083 for (unsigned l = 0; l < n_node; l++)
4084 {
4085 // Loop over the number of dofs
4086 for (unsigned k = 0; k < n_position_type; k++)
4087 {
4088 x[i] += nodal_position_gen(t, l, k, i) * psi(l, k);
4089 }
4090 }
4091 }
4092 }
4093
4094 //========================================================================
4095 /// Calculate the determinant of the
4096 /// Jacobian of the mapping between local and global
4097 /// coordinates at the position. Works directly from the base vectors
4098 /// without assuming that coordinates match spatial dimension. Will
4099 /// be overloaded in FaceElements, in which the elemental dimension does
4100 /// not match the spatial dimension. WARNING: this is always positive
4101 /// and cannot be used to check if the element is inverted, say!
4102 //========================================================================
4104 {
4105 // Find the number of nodes and position types
4106 const unsigned n_node = nnode();
4107 const unsigned n_position_type = nnodal_position_type();
4108 // Find the dimension of the node and element
4109 const unsigned n_dim_node = nodal_dimension();
4110 const unsigned n_dim_element = dim();
4111
4112 // Set up dummy memory for the shape functions
4113 Shape psi(n_node, n_position_type);
4114 DShape dpsids(n_node, n_position_type, n_dim_element);
4115 // Get the shape functions and local derivatives
4116 dshape_local(s, psi, dpsids);
4117
4118 // Right calculate the base vectors
4119 DenseMatrix<double> interpolated_G(n_dim_element, n_dim_node);
4120 assemble_eulerian_base_vectors(dpsids, interpolated_G);
4121
4122 // Calculate the metric tensor of the element
4123 DenseMatrix<double> G(n_dim_element, n_dim_element, 0.0);
4124 for (unsigned i = 0; i < n_dim_element; i++)
4125 {
4126 for (unsigned j = 0; j < n_dim_element; j++)
4127 {
4128 for (unsigned k = 0; k < n_dim_node; k++)
4129 {
4130 G(i, j) += interpolated_G(i, k) * interpolated_G(j, k);
4131 }
4132 }
4133 }
4134
4135 // Calculate the determinant of the metric tensor
4136 double det = 0.0;
4137 switch (n_dim_element)
4138 {
4139 case 0:
4140 throw OomphLibError("Cannot calculate J_eulerian() for point element\n",
4141 OOMPH_CURRENT_FUNCTION,
4142 OOMPH_EXCEPTION_LOCATION);
4143 break;
4144 case 1:
4145 det = G(0, 0);
4146 break;
4147 case 2:
4148 det = G(0, 0) * G(1, 1) - G(0, 1) * G(1, 0);
4149 break;
4150 case 3:
4151 det = G(0, 0) * G(1, 1) * G(2, 2) + G(0, 1) * G(1, 2) * G(2, 0) +
4152 G(0, 2) * G(1, 0) * G(2, 1) - G(0, 0) * G(1, 2) * G(2, 1) -
4153 G(0, 1) * G(1, 0) * G(2, 2) - G(0, 2) * G(1, 1) * G(2, 0);
4154 break;
4155 default:
4156 oomph_info << "More than 3 dimensions in J_eulerian()" << std::endl;
4157 break;
4158 }
4159
4160 // Return the Jacobian (square-root of the determinant of the metric tensor)
4161 return sqrt(det);
4162 }
4163
4164 //========================================================================
4165 /// Compute the Jacobian of the mapping between the local and global
4166 /// coordinates at the ipt-th integration point
4167 //========================================================================
4168 double FiniteElement::J_eulerian_at_knot(const unsigned& ipt) const
4169 {
4170 // Find the number of nodes
4171 const unsigned n_node = nnode();
4172 // Find the number of position types
4173 const unsigned n_position_type = nnodal_position_type();
4174 // Find the dimension of the node and element
4175 const unsigned n_dim_node = nodal_dimension();
4176 const unsigned n_dim_element = dim();
4177
4178 // Set up dummy memory for the shape functions
4179 Shape psi(n_node, n_position_type);
4180 DShape dpsids(n_node, n_position_type, n_dim_element);
4181 // Get the shape functions and local derivatives at the knot
4182 // This call may use the stored versions, which is why this entire
4183 // function doesn't just call J_eulerian(s), after reading out s from
4184 // the knots.
4185 dshape_local_at_knot(ipt, psi, dpsids);
4186
4187 // Right calculate the base vectors
4188 DenseMatrix<double> interpolated_G(n_dim_element, n_dim_node);
4189 assemble_eulerian_base_vectors(dpsids, interpolated_G);
4190
4191 // Calculate the metric tensor of the element
4192 DenseMatrix<double> G(n_dim_element, n_dim_element, 0.0);
4193 for (unsigned i = 0; i < n_dim_element; i++)
4194 {
4195 for (unsigned j = 0; j < n_dim_element; j++)
4196 {
4197 for (unsigned k = 0; k < n_dim_node; k++)
4198 {
4199 G(i, j) += interpolated_G(i, k) * interpolated_G(j, k);
4200 }
4201 }
4202 }
4203
4204 // Calculate the determinant of the metric tensor
4205 double det = 0.0;
4206 switch (n_dim_element)
4207 {
4208 case 0:
4209 throw OomphLibError("Cannot calculate J_eulerian() for point element\n",
4210 OOMPH_CURRENT_FUNCTION,
4211 OOMPH_EXCEPTION_LOCATION);
4212 break;
4213 case 1:
4214 det = G(0, 0);
4215 break;
4216 case 2:
4217 det = G(0, 0) * G(1, 1) - G(0, 1) * G(1, 0);
4218 break;
4219 case 3:
4220 det = G(0, 0) * G(1, 1) * G(2, 2) + G(0, 1) * G(1, 2) * G(2, 0) +
4221 G(0, 2) * G(1, 0) * G(2, 1) - G(0, 0) * G(1, 2) * G(2, 1) -
4222 G(0, 1) * G(1, 0) * G(2, 2) - G(0, 2) * G(1, 1) * G(2, 0);
4223 break;
4224 default:
4225 oomph_info << "More than 3 dimensions in J_eulerian()" << std::endl;
4226 break;
4227 }
4228
4229 // Return the Jacobian (square-root of the determinant of the metric tensor)
4230 return sqrt(det);
4231 }
4232
4233 //========================================================================
4234 /// Check that Jacobian of mapping between local and Eulerian
4235 /// coordinates at all integration points is positive.
4236 //========================================================================
4238 {
4239 // Bypass check deep down in the guts...
4242
4243 // Let's be optimistic...
4244 passed = true;
4245
4246 // Find the number of nodes
4247 const unsigned n_node = nnode();
4248
4249 // Find the number of position types
4250 const unsigned n_position_type = nnodal_position_type();
4251
4252 // DRAIG: Unused variable
4253 // const unsigned n_dim_node = nodal_dimension();
4254
4255 // Find the dimension of the node and element
4256 const unsigned n_dim_element = dim();
4257
4258 // Set up dummy memory for the shape functions
4259 Shape psi(n_node, n_position_type);
4260 DShape dpsi(n_node, n_dim_element);
4261
4262 unsigned nintpt = integral_pt()->nweight();
4263 for (unsigned ipt = 0; ipt < nintpt; ipt++)
4264 {
4265 double jac = dshape_eulerian_at_knot(ipt, psi, dpsi);
4266
4267 // Are we dead yet?
4268 if (jac <= 0.0)
4269 {
4270 passed = false;
4271
4272 // Reset
4274
4275 return;
4276 }
4277 }
4278
4279 // Reset
4281 }
4282
4283 //====================================================================
4284 /// Calculate the size of the element (in Eulerian computational
4285 /// coordinates. Use suitably overloaded compute_physical_size()
4286 /// function to compute the actual size (taking into account
4287 /// factors such as 2pi or radii the integrand). Such function
4288 /// can only be implemented on an equation-by-equation basis.
4289 //====================================================================
4290 double FiniteElement::size() const
4291 {
4292 // Initialise the sum to zero
4293 double sum = 0.0;
4294
4295 // Loop over the integration points
4296 const unsigned n_intpt = integral_pt()->nweight();
4297 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
4298 {
4299 // Get the integral weight
4300 double w = integral_pt()->weight(ipt);
4301 // Get the value of the Jacobian of the mapping to global coordinates
4302 double J = J_eulerian_at_knot(ipt);
4303
4304 // Add the product to the sum
4305 sum += w * J;
4306 }
4307
4308 // Return the answer
4309 return (sum);
4310 }
4311
4312 //==================================================================
4313 /// Integrate Vector-valued time-dep function over element
4314 //==================================================================
4317 const double& time,
4318 Vector<double>& integral)
4319 {
4320 // Initialise all components of integral Vector and setup integrand vector
4321 const unsigned ncomponents = integral.size();
4322 Vector<double> integrand(ncomponents);
4323 for (unsigned i = 0; i < ncomponents; i++)
4324 {
4325 integral[i] = 0.0;
4326 }
4327
4328 // Figure out the global (Eulerian) spatial dimension of the
4329 // element
4330 const unsigned n_dim_eulerian = nodal_dimension();
4331
4332 // Allocate Vector of global Eulerian coordinates
4333 Vector<double> x(n_dim_eulerian);
4334
4335 // Set the value of n_intpt
4336 const unsigned n_intpt = integral_pt()->nweight();
4337
4338 // Vector of local coordinates
4339 const unsigned n_dim = dim();
4340 Vector<double> s(n_dim);
4341
4342 // Loop over the integration points
4343 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
4344 {
4345 // Assign the values of s
4346 for (unsigned i = 0; i < n_dim; i++)
4347 {
4348 s[i] = integral_pt()->knot(ipt, i);
4349 }
4350
4351 // Assign the values of the global Eulerian coordinates
4352 for (unsigned i = 0; i < n_dim_eulerian; i++)
4353 {
4354 x[i] = interpolated_x(s, i);
4355 }
4356
4357 // Get the integral weight
4358 double w = integral_pt()->weight(ipt);
4359
4360 // Get Jacobian of mapping
4361 double J = J_eulerian(s);
4362
4363 // Evaluate the integrand at the knot points
4364 integrand_fct_pt(time, x, integrand);
4365
4366 // Add to components of integral Vector
4367 for (unsigned i = 0; i < ncomponents; i++)
4368 {
4369 integral[i] += integrand[i] * w * J;
4370 }
4371 }
4372 }
4373
4374 //==================================================================
4375 /// Integrate Vector-valued function over element
4376 //==================================================================
4379 Vector<double>& integral)
4380 {
4381 // Initialise all components of integral Vector
4382 const unsigned ncomponents = integral.size();
4383 Vector<double> integrand(ncomponents);
4384 for (unsigned i = 0; i < ncomponents; i++)
4385 {
4386 integral[i] = 0.0;
4387 }
4388
4389 // Figure out the global (Eulerian) spatial dimension of the
4390 // element by checking the Eulerian dimension of the nodes
4391 const unsigned n_dim_eulerian = nodal_dimension();
4392
4393 // Allocate Vector of global Eulerian coordinates
4394 Vector<double> x(n_dim_eulerian);
4395
4396 // Set the value of n_intpt
4397 const unsigned n_intpt = integral_pt()->nweight();
4398
4399 // Vector of local coordinates
4400 const unsigned n_dim = dim();
4401 Vector<double> s(n_dim);
4402
4403 // Loop over the integration points
4404 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
4405 {
4406 // Assign the values of s
4407 for (unsigned i = 0; i < n_dim; i++)
4408 {
4409 s[i] = integral_pt()->knot(ipt, i);
4410 }
4411
4412 // Assign the values of the global Eulerian coordinates
4413 for (unsigned i = 0; i < n_dim_eulerian; i++)
4414 {
4415 x[i] = interpolated_x(s, i);
4416 }
4417
4418 // Get the integral weight
4419 double w = integral_pt()->weight(ipt);
4420
4421 // Get Jacobian of mapping
4422 double J = J_eulerian(s);
4423
4424 // Evaluate the integrand at the knot points
4425 integrand_fct_pt(x, integrand);
4426
4427 // Add to components of integral Vector
4428 for (unsigned i = 0; i < ncomponents; i++)
4429 {
4430 integral[i] += integrand[i] * w * J;
4431 }
4432 }
4433 }
4434
4435 //==========================================================================
4436 /// Self-test: Have all internal values been classified as
4437 /// pinned/unpinned? Has pointer to spatial integration scheme
4438 /// been set? Return 0 if OK.
4439 //==========================================================================
4441 {
4442 // Initialise
4443 bool passed = true;
4444
4446 {
4447 passed = false;
4448 }
4449
4450 // Check that pointer to spatial integration scheme has been set
4451 if (integral_pt() == 0)
4452 {
4453 passed = false;
4454
4455 OomphLibWarning("Pointer to spatial integration scheme has not been set.",
4456 "FiniteElement::self_test()",
4457 OOMPH_EXCEPTION_LOCATION);
4458 }
4459
4460 // If the dimension of the element is zero (point element), there
4461 // is not jacobian
4462 const unsigned dim_el = dim();
4463
4464 if (dim_el > 0)
4465 {
4466 // Loop over integration points to check sign of Jacobian
4467 //-------------------------------------------------------
4468
4469 // Set the value of n_intpt
4470 const unsigned n_intpt = integral_pt()->nweight();
4471
4472 // Set the Vector to hold local coordinates
4473 Vector<double> s(dim_el);
4474
4475 // Find the number of local nodes
4476 const unsigned n_node = nnode();
4477 const unsigned n_position_type = nnodal_position_type();
4478 // Set up memory for the shape and test functions
4479 Shape psi(n_node, n_position_type);
4480 DShape dpsidx(n_node, dim_el);
4481
4482 // Jacobian
4483 double jacobian;
4484
4485
4486 // Two ways of testing for negative Jacobian for non-FaceElements
4487 unsigned ntest = 1;
4488
4489 // For FaceElements checking the Jacobian via dpsidx doesn't
4490 // make sense
4491 FiniteElement* tmp_pt = const_cast<FiniteElement*>(this);
4492 FaceElement* face_el_pt = dynamic_cast<FaceElement*>(tmp_pt);
4493 if (face_el_pt == 0)
4494 {
4495 ntest = 2;
4496 }
4497
4498 // For now overwrite -- the stuff above fails for Bretherton.
4499 // Not sure why.
4500 ntest = 1;
4501
4502 // Loop over the integration points
4503 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
4504 {
4505 // Assign values of s
4506 for (unsigned i = 0; i < dim_el; i++)
4507 {
4508 s[i] = integral_pt()->knot(ipt, i);
4509 }
4510
4511
4512 // Do tests
4513 for (unsigned test = 0; test < ntest; test++)
4514 {
4515 switch (test)
4516 {
4517 case 0:
4518
4519 // Get the jacobian from the mapping between local and Eulerian
4520 // coordinates
4521 jacobian = J_eulerian(s);
4522
4523 break;
4524
4525 case 1:
4526
4527 // Call the geometrical shape functions and derivatives
4528 // This also computes the Jacobian by a slightly different
4529 // method
4530 jacobian = dshape_eulerian_at_knot(ipt, psi, dpsidx);
4531
4532 break;
4533
4534 default:
4535
4536 throw OomphLibError("Never get here",
4537 OOMPH_CURRENT_FUNCTION,
4538 OOMPH_EXCEPTION_LOCATION);
4539 }
4540
4541
4542 // Check for a singular jacobian
4543 if (std::fabs(jacobian) < 1.0e-16)
4544 {
4545 std::ostringstream warning_stream;
4546 warning_stream << "Determinant of Jacobian matrix is zero at ipt "
4547 << ipt << std::endl;
4548 OomphLibWarning(warning_stream.str(),
4549 "FiniteElement::self_test()",
4550 OOMPH_EXCEPTION_LOCATION);
4551 passed = false;
4552 // Skip the next test
4553 continue;
4554 }
4555
4556 // Check sign of Jacobian
4557 if ((Accept_negative_jacobian == false) && (jacobian < 0.0))
4558 {
4559 std::ostringstream warning_stream;
4560 warning_stream << "Jacobian negative at integration point ipt="
4561 << ipt << std::endl;
4562 warning_stream
4563 << "If you think that this is what you want you may: "
4564 << std::endl
4565 << "set the (static) flag "
4566 << "FiniteElement::Accept_negative_jacobian to be true"
4567 << std::endl;
4568
4569 OomphLibWarning(warning_stream.str(),
4570 "FiniteElement::self_test()",
4571 OOMPH_EXCEPTION_LOCATION);
4572 passed = false;
4573 }
4574
4575 } // end of loop over two tests
4576 }
4577 } // End of non-zero dimension check
4578
4579 // Return verdict
4580 if (passed)
4581 {
4582 return 0;
4583 }
4584 else
4585 {
4586 return 1;
4587 }
4588 }
4589
4590
4591 //=======================================================================
4592 /// Return the t-th time-derivative of the
4593 /// i-th FE-interpolated Eulerian coordinate at
4594 /// local coordinate s.
4595 //=======================================================================
4597 const unsigned& i,
4598 const unsigned& t_deriv)
4599 {
4600 // Find the number of nodes and positions (locally cached)
4601 const unsigned n_node = nnode();
4602 const unsigned n_position_type = nnodal_position_type();
4603 // Get shape functions: specify # of nodes, # of positional dofs
4604 Shape psi(n_node, n_position_type);
4605 shape(s, psi);
4606
4607 // Initialise
4608 double drdt = 0.0;
4609
4610 // Assemble time derivative
4611 // Loop over nodes
4612 for (unsigned l = 0; l < n_node; l++)
4613 {
4614 // Loop over types of dof
4615 for (unsigned k = 0; k < n_position_type; k++)
4616 {
4617 drdt += dnodal_position_gen_dt(t_deriv, l, k, i) * psi(l, k);
4618 }
4619 }
4620 return drdt;
4621 }
4622
4623
4624 //=======================================================================
4625 /// Compute t-th time-derivative of the
4626 /// FE-interpolated Eulerian coordinate vector at
4627 /// local coordinate s.
4628 //=======================================================================
4630 const unsigned& t_deriv,
4631 Vector<double>& dxdt)
4632 {
4633 // Find the number of nodes and positions (locally cached)
4634 const unsigned n_node = nnode();
4635 const unsigned n_position_type = nnodal_position_type();
4636 const unsigned nodal_dim = nodal_dimension();
4637
4638 // Get shape functions: specify # of nodes, # of positional dofs
4639 Shape psi(n_node, n_position_type);
4640 shape(s, psi);
4641
4642 // Loop over directions
4643 for (unsigned i = 0; i < nodal_dim; i++)
4644 {
4645 // Initialise
4646 dxdt[i] = 0.0;
4647
4648 // Assemble time derivative
4649 // Loop over nodes
4650 for (unsigned l = 0; l < n_node; l++)
4651 {
4652 // Loop over types of dof
4653 for (unsigned k = 0; k < n_position_type; k++)
4654 {
4655 dxdt[i] += dnodal_position_gen_dt(t_deriv, l, k, i) * psi(l, k);
4656 }
4657 }
4658 }
4659 }
4660
4661 //============================================================================
4662 /// Calculate the interpolated value of zeta, the intrinsic coordinate
4663 /// of the element when viewed as a compound geometric object within a Mesh
4664 /// as a function of the local coordinate of the element, s.
4665 /// The default
4666 /// assumption is the zeta is interpolated using the shape functions of
4667 /// the element with the values given by zeta_nodal().
4668 /// A MacroElement representation of the intrinsic coordinate parametrised
4669 /// by the local coordinate s is used if available.
4670 /// Choosing the MacroElement representation of zeta (Eulerian x by default)
4671 /// allows a correspondence to be established between elements on different
4672 /// Meshes covering the same curvilinear domain in cases where one element
4673 /// is much coarser than the other.
4674 //==========================================================================
4676 Vector<double>& zeta) const
4677 {
4678 // If there is a macro element use it
4679 if (Macro_elem_pt != 0)
4680 {
4681 this->get_x_from_macro_element(s, zeta);
4682 }
4683 // Otherwise interpolate zeta_nodal using the shape functions
4684 else
4685 {
4686 // Find the number of nodes
4687 const unsigned n_node = this->nnode();
4688 // Find the number of positional types
4689 const unsigned n_position_type = this->nnodal_position_type();
4690 // Storage for the shape functions
4691 Shape psi(n_node, n_position_type);
4692 // Get the values of the shape functions at the local coordinate s
4693 this->shape(s, psi);
4694
4695 // Find the number of coordinates
4696 const unsigned ncoord = this->dim();
4697 // Initialise the value of zeta to zero
4698 for (unsigned i = 0; i < ncoord; i++)
4699 {
4700 zeta[i] = 0.0;
4701 }
4702
4703 // Add the contributions from each nodal dof to the interpolated value
4704 // of zeta.
4705 for (unsigned l = 0; l < n_node; l++)
4706 {
4707 for (unsigned k = 0; k < n_position_type; k++)
4708 {
4709 // Locally cache the value of the shape function
4710 const double psi_ = psi(l, k);
4711 for (unsigned i = 0; i < ncoord; i++)
4712 {
4713 zeta[i] += this->zeta_nodal(l, k, i) * psi_;
4714 }
4715 }
4716 }
4717 }
4718 }
4719
4720 //==========================================================================
4721 /// For a given value of zeta, the "global" intrinsic coordinate of
4722 /// a mesh of FiniteElements represented as a compound geometric object,
4723 /// find the local coordinate in this element that corresponds to the
4724 /// requested value of zeta.
4725 /// This is achieved in generality by using Newton's method to find the value
4726 /// of the local coordinate, s, such that
4727 /// interpolated_zeta(s) is equal to the requested value of zeta.
4728 /// If zeta cannot be located in this element, geom_object_pt is set
4729 /// to NULL. If zeta is located in this element, we return its "this"
4730 /// pointer.
4731 /// Setting the optional bool argument to true means that the coordinate
4732 /// argument "s" is used as the initial guess. (Default is false).
4733 //=========================================================================
4735 GeomObject*& geom_object_pt,
4737 const bool& use_coordinate_as_initial_guess)
4738 {
4739 // Find the number of coordinates, the dimension of the element
4740 // This must be the same for the local and intrinsic global coordinate
4741 unsigned ncoord = this->dim();
4742
4743 // Fast return based on centre of gravity and max. radius of any nodal
4744 // point
4746 0.0)
4747 {
4748 Vector<double> cog(ncoord);
4749 double max_radius = 0.0;
4751
4752 // Get radius
4753 double radius = 0.0;
4754 for (unsigned i = 0; i < ncoord; i++)
4755 {
4756 radius += (cog[i] - zeta[i]) * (cog[i] - zeta[i]);
4757 }
4758 radius = sqrt(radius);
4759 if (radius > Locate_zeta_helpers::
4761 max_radius)
4762 {
4763 geom_object_pt = 0;
4764 return;
4765 }
4766 }
4767
4768 // Assign storage for the vector and matrix used in Newton's method
4769 Vector<double> dx(ncoord, 0.0);
4770 DenseDoubleMatrix jacobian(ncoord, ncoord, 0.0);
4771
4772 // // Make a list of (equally-spaced) local coordinates inside the element
4773 // unsigned n_list=Locate_zeta_helpers::N_local_points;
4774
4775 // double list_space=(1.0/(double(n_list)-1.0))*(s_max()-s_min());
4776
4777 // Possible initial guesses for Newton's method
4778 Vector<Vector<double>> s_list;
4779
4780 // If the boolean argument use_coordinate_as_initial_guess was set
4781 // to true then we don't need to initialise s
4782 if (!use_coordinate_as_initial_guess)
4783 {
4784 // Vector of local coordinates
4785 Vector<double> s_c(ncoord);
4786
4787 // Loop over plot points
4788 unsigned num_plot_points =
4790 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
4791 {
4792 // Get local coordinates of plot point
4794 s_list.push_back(s_c);
4795 }
4796 }
4797
4798 // Counter for the number of Newton steps
4799 unsigned count = 0;
4800
4801 // Control flag for the Newton loop
4802 bool keep_going = true;
4803
4804 // Storage for the interpolated value of x
4805 Vector<double> inter_x(ncoord);
4806
4807 // Ifwe have specified the coordinate already
4808 if (use_coordinate_as_initial_guess)
4809 {
4810 // Get the value of x at the initial guess
4811 this->interpolated_zeta(s, inter_x);
4812
4813 // Set up the residuals
4814 for (unsigned i = 0; i < ncoord; i++)
4815 {
4816 dx[i] = zeta[i] - inter_x[i];
4817 }
4818 }
4819 else
4820 {
4821 // Find the smallest residual from the list of coordinates made earlier
4822 double my_min_resid = DBL_MAX;
4823 Vector<double> s_local(ncoord);
4824 Vector<double> work_x(ncoord);
4825 Vector<double> work_dx(ncoord);
4826
4827 unsigned n_list_coord = s_list.size();
4828
4829 for (unsigned i_coord = 0; i_coord < n_list_coord; i_coord++)
4830 {
4831 for (unsigned i = 0; i < ncoord; i++)
4832 {
4833 s_local[i] = s_list[i_coord][i];
4834 }
4835 // get_x for this coordinate
4836 this->interpolated_zeta(s_local, work_x);
4837
4838 // calculate residuals
4839 for (unsigned i = 0; i < ncoord; i++)
4840 {
4841 work_dx[i] = zeta[i] - work_x[i];
4842 }
4843
4844 double maxres = std::fabs(
4845 *std::max_element(work_dx.begin(), work_dx.end(), AbsCmp<double>()));
4846
4847 // test against previous residuals
4848 if (maxres < my_min_resid)
4849 {
4850 my_min_resid = maxres;
4851 dx = work_dx;
4852 inter_x = work_x;
4853 s = s_local;
4854 }
4855 }
4856 }
4857
4858 // Main Newton Loop
4859 do // start of do while loop
4860 {
4861 // Increase loop counter
4862 count++;
4863
4864 // Bail out if necessary (without an error for now...)
4866 {
4867 keep_going = false;
4868 continue;
4869 }
4870
4871 // If it's the first time round the loop, check the initial residuals
4872 if (count == 1)
4873 {
4874 double maxres =
4875 std::fabs(*std::max_element(dx.begin(), dx.end(), AbsCmp<double>()));
4876
4877 // If it's small enough exit
4879 {
4880 keep_going = false;
4881 continue;
4882 }
4883 }
4884
4885 // Is there a macro element? If so, assemble the Jacobian by FD-ing
4886 if (macro_elem_pt() != 0)
4887 {
4888 // Assemble jacobian on the fly by finite differencing
4889 Vector<double> work_s = s;
4890 Vector<double> r = inter_x; // i.e. the result of previous call to get_x
4891
4892 // Finite difference step
4894
4895 // Storage for calculated r from incremented s
4896 Vector<double> work_r(ncoord, 0.0);
4897
4898 // Loop over s coordinates
4899 for (unsigned i = 0; i < ncoord; i++)
4900 {
4901 // Increment work_s by a small amount
4902 work_s[i] += fd_step;
4903
4904 // Calculate work_r from macro element
4905 this->interpolated_zeta(work_s, work_r);
4906
4907 // Loop over r to fill Jacobian
4908 for (unsigned j = 0; j < ncoord; j++)
4909 {
4910 jacobian(j, i) = -(work_r[j] - r[j]) / fd_step;
4911 }
4912
4913 // Reset work_s
4914 work_s[i] = s[i];
4915 }
4916 }
4917 else // no macro element, so compute Jacobian with shape functions etc.
4918 {
4919 // Compute the entries of the Jacobian matrix
4920 unsigned n_node = this->nnode();
4921 unsigned n_position_type = this->nnodal_position_type();
4922 Shape psi(n_node, n_position_type);
4923 DShape dpsids(n_node, n_position_type, ncoord);
4924
4925 // Get the local shape functions and their derivatives
4926 dshape_local(s, psi, dpsids);
4927
4928 // Calculate the values of dxds
4929 DenseMatrix<double> interpolated_dxds(ncoord, ncoord, 0.0);
4930
4931 // MH: No longer needed
4932 // //This implementation will only work for n_position_type=1
4933 // //since the function nodal_position_gen does not yet exist
4934 // #ifdef PARANOID
4935 // if (n_position_type!=1)
4936 // {
4937 // std::ostringstream error_stream;
4938 // error_stream << "This implementation does not exist
4939 // yet;\n"
4940 // << "it currently uses
4941 // raw_nodal_position_gen\n"
4942 // << "which does not take hangingness into
4943 // account\n"
4944 // << "It will work if n_position_type=1\n";
4945 // throw OomphLibError(error_stream.str(),
4946 // OOMPH_CURRENT_FUNCTION,
4947 // OOMPH_EXCEPTION_LOCATION);
4948 // }
4949 // #endif
4950
4951 // Loop over the nodes
4952 for (unsigned l = 0; l < n_node; l++)
4953 {
4954 // Loop over position type even though it should be 1; the
4955 // functionality for n_position_type>1 will exist in the future
4956 for (unsigned k = 0; k < n_position_type; k++)
4957 {
4958 // Add the contribution from the nodal coordinates to the matrix
4959 for (unsigned i = 0; i < ncoord; i++)
4960 {
4961 for (unsigned j = 0; j < ncoord; j++)
4962 {
4963 interpolated_dxds(i, j) +=
4964 this->zeta_nodal(l, k, i) * dpsids(l, k, j);
4965 }
4966 }
4967 }
4968 }
4969
4970 // The entries of the Jacobian matrix are merely dresiduals/ds
4971 // i.e. \f$ -dx/ds \f$
4972 for (unsigned i = 0; i < ncoord; i++)
4973 {
4974 for (unsigned j = 0; j < ncoord; j++)
4975 {
4976 jacobian(i, j) = -interpolated_dxds(i, j);
4977 }
4978 }
4979 }
4980
4981 // Now solve the damn thing
4982 try
4983 {
4984 jacobian.solve(dx);
4985 }
4986 catch (OomphLibError& error)
4987 {
4988 // I've caught the error so shut up!
4989 error.disable_error_message();
4990#ifdef PARANOID
4991 oomph_info << "Error in linear solve for "
4992 << "FiniteElement::locate_zeta()" << std::endl;
4993 oomph_info << "Should not affect the result!" << std::endl;
4994#endif
4995 }
4996
4997 // Add the correction to the local coordinates
4998 for (unsigned i = 0; i < ncoord; i++)
4999 {
5000 s[i] -= dx[i];
5001 }
5002
5003 // Get the new residuals
5004 this->interpolated_zeta(s, inter_x);
5005 for (unsigned i = 0; i < ncoord; i++)
5006 {
5007 dx[i] = zeta[i] - inter_x[i];
5008 }
5009
5010 // Get the maximum residuals
5011 double maxres =
5012 std::fabs(*std::max_element(dx.begin(), dx.end(), AbsCmp<double>()));
5013
5014 // If we have converged jump straight to the test at the end of the loop
5016 {
5017 keep_going = false;
5018 continue;
5019 }
5020 } while (keep_going);
5021
5022 // Test whether the local coordinates are valid or not
5023 bool valid = local_coord_is_valid(s);
5024
5025 // If not valid, experimentally push back into element
5026 // and see if the result is still valid (within the Newton tolerance)
5027 if (!valid)
5028 {
5030
5031 // Check residuals again
5032 this->interpolated_zeta(s, inter_x);
5033 for (unsigned i = 0; i < ncoord; i++)
5034 {
5035 dx[i] = zeta[i] - inter_x[i];
5036 }
5037
5038 // Get the maximum residuals
5039 double maxres =
5040 std::fabs(*std::max_element(dx.begin(), dx.end(), AbsCmp<double>()));
5041
5042 // Are we still OK?
5044 {
5045 // oomph_info
5046 // << "Pushing back inside has violated the Newton tolerance: max_res =
5047 // "
5048 // << maxres << std::endl;
5049 geom_object_pt = 0;
5050 return;
5051 }
5052 }
5053
5054 // It is also possible now that it may not have converged "correctly",
5055 // i.e. count is greater than Max_newton_iterations
5057 {
5058 // Don't trust the current answer, return null
5059 geom_object_pt = 0;
5060 return;
5061 }
5062
5063 // Otherwise the required point is located in "this" element:
5064 geom_object_pt = this;
5065 }
5066
5067
5068 //=======================================================================
5069 /// Loop over all nodes in the element and update their positions
5070 /// using each node's (algebraic) update function
5071 //=======================================================================
5073 {
5074 const unsigned n_node = nnode();
5075 for (unsigned n = 0; n < n_node; n++)
5076 {
5077 node_pt(n)->node_update();
5078 }
5079 }
5080
5081 //======================================================================
5082 /// The purpose of this function is to identify all possible
5083 /// Data that can affect the fields interpolated by the FiniteElement.
5084 /// The information will typically be used in interaction problems in
5085 /// which the FiniteElement provides a forcing term for an
5086 /// ElementWithExternalElement. The Data must be provided as
5087 /// \c paired_load data containing (a) the pointer to a Data object
5088 /// and (b) the index of the value in that Data object.
5089 /// The generic implementation (should be overloaded in more specific
5090 /// applications) is to include all nodal and internal Data stored in
5091 /// the FiniteElement. Note that the geometric data,
5092 /// which includes the positions
5093 /// of SolidNodes, is treated separately by the function
5094 /// \c identify_geometric_data()
5095 //======================================================================
5097 std::set<std::pair<Data*, unsigned>>& paired_field_data)
5098 {
5099 // Loop over all internal data
5100 const unsigned n_internal = this->ninternal_data();
5101 for (unsigned n = 0; n < n_internal; n++)
5102 {
5103 // Cache the data pointer
5104 Data* const dat_pt = this->internal_data_pt(n);
5105 // Find the number of data values stored in the data object
5106 const unsigned n_value = dat_pt->nvalue();
5107 // Add the index of each data value and the pointer to the set
5108 // of pairs
5109 for (unsigned i = 0; i < n_value; i++)
5110 {
5111 paired_field_data.insert(std::make_pair(dat_pt, i));
5112 }
5113 }
5114
5115 // Loop over all the nodes
5116 const unsigned n_node = this->nnode();
5117 for (unsigned n = 0; n < n_node; n++)
5118 {
5119 // Cache the node pointer
5120 Node* const nod_pt = this->node_pt(n);
5121 // Find the number of values stored at the node
5122 const unsigned n_value = nod_pt->nvalue();
5123 // Add the index of each data value and the pointer to the set
5124 // of pairs
5125 for (unsigned i = 0; i < n_value; i++)
5126 {
5127 paired_field_data.insert(std::make_pair(nod_pt, i));
5128 }
5129 }
5130 }
5131
5132 void FiniteElement::build_face_element(const int& face_index,
5133 FaceElement* face_element_pt)
5134 {
5135 // Set the nodal dimension
5136 face_element_pt->set_nodal_dimension(nodal_dimension());
5137
5138 // Set the pointer to the orginal "bulk" element
5139 face_element_pt->bulk_element_pt() = this;
5140
5141#ifdef OOMPH_HAS_MPI
5142 // Pass on non-halo proc ID
5143 face_element_pt->set_halo(Non_halo_proc_ID);
5144#endif
5145
5146 // Set the face index
5147 face_element_pt->face_index() = face_index;
5148
5149 // Get number of bulk nodes on a face of this element
5150 const unsigned nnode_face = nnode_on_face();
5151
5152 // Set the function pointer for face coordinate to bulk coordinate
5153 // mapping
5154 face_element_pt->face_to_bulk_coordinate_fct_pt() =
5156
5157 // Set the function pointer for the derivative of the face coordinate to
5158 // bulk coordinate mapping
5159 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
5161
5162 // Resize storage for the number of values originally stored each of the
5163 // face element's nodes.
5164 face_element_pt->nbulk_value_resize(nnode_face);
5165
5166 // Resize storage for the bulk node numbers corresponding to the face
5167 // element's nodes.
5168 face_element_pt->bulk_node_number_resize(nnode_face);
5169
5170 // Copy bulk_node_numbers and nbulk_values
5171 for (unsigned i = 0; i < nnode_face; i++)
5172 {
5173 // Find the corresponding bulk node's number
5174 unsigned bulk_number = get_bulk_node_number(face_index, i);
5175
5176 // Assign the pointer and number into the face element
5177 face_element_pt->node_pt(i) = node_pt(bulk_number);
5178 face_element_pt->bulk_node_number(i) = bulk_number;
5179
5180 // Set the number of values originally stored at this node
5181 face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
5182 }
5183
5184 // Set the outer unit normal sign
5185 face_element_pt->normal_sign() = face_outer_unit_normal_sign(face_index);
5186 }
5187
5188 /// ///////////////////////////////////////////////////////////////////////
5189 /// ///////////////////////////////////////////////////////////////////////
5190 /// ///////////////////////////////////////////////////////////////////////
5191
5192 // Initialise the static variable.
5193 // Do not ignore warning for discontinuous tangent vectors.
5195
5196
5197 //========================================================================
5198 /// Output boundary coordinate zeta
5199 //========================================================================
5200 void FaceElement::output_zeta(std::ostream& outfile, const unsigned& nplot)
5201 {
5202 // Vector of local coordinates
5203 unsigned n_dim = dim();
5204 Vector<double> s(n_dim);
5205
5206 // Tecplot header info
5207 outfile << tecplot_zone_string(nplot);
5208
5209 // Loop over plot points
5210 unsigned num_plot_points = nplot_points(nplot);
5211 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
5212 {
5213 // Get local coordinates of plot point
5214 get_s_plot(iplot, nplot, s);
5215
5216 // Spatial coordinates are one higher
5217 for (unsigned i = 0; i < n_dim + 1; i++)
5218 {
5219 outfile << interpolated_x(s, i) << " ";
5220 }
5221
5222 // Boundary coordinate
5223 Vector<double> zeta(n_dim);
5224 interpolated_zeta(s, zeta);
5225 for (unsigned i = 0; i < n_dim; i++)
5226 {
5227 outfile << zeta[i] << " ";
5228 }
5229 outfile << std::endl;
5230 }
5231
5232 // Write tecplot footer (e.g. FE connectivity lists)
5233 write_tecplot_zone_footer(outfile, nplot);
5234 }
5235
5236
5237 //========================================================================
5238 /// Calculate the determinant of the
5239 /// Jacobian of the mapping between local and global
5240 /// coordinates at the position s. Overloaded from FiniteElement.
5241 //========================================================================
5243 {
5244 // Find out the sptial dimension of the element
5245 unsigned n_dim_el = this->dim();
5246
5247 // Bail out if we're in a point element -- not sure what
5248 // J_eulerian actually is, but this is harmless
5249 if (n_dim_el == 0) return 1.0;
5250
5251 // Find out how many nodes there are
5252 unsigned n_node = nnode();
5253
5254 // Find out how many positional dofs there are
5255 unsigned n_position_type = this->nnodal_position_type();
5256
5257 // Find out the dimension of the node
5258 unsigned n_dim = this->nodal_dimension();
5259
5260 // Set up memory for the shape functions
5261 Shape psi(n_node, n_position_type);
5262 DShape dpsids(n_node, n_position_type, n_dim_el);
5263
5264 // Only need to call the local derivatives
5265 dshape_local(s, psi, dpsids);
5266
5267 // Also calculate the surface Vectors (derivatives wrt local coordinates)
5268 DenseMatrix<double> interpolated_A(n_dim_el, n_dim, 0.0);
5269
5270 // Calculate positions and derivatives
5271 for (unsigned l = 0; l < n_node; l++)
5272 {
5273 // Loop over positional dofs
5274 for (unsigned k = 0; k < n_position_type; k++)
5275 {
5276 // Loop over coordinates
5277 for (unsigned i = 0; i < n_dim; i++)
5278 {
5279 // Loop over LOCAL derivative directions, to calculate the tangent(s)
5280 for (unsigned j = 0; j < n_dim_el; j++)
5281 {
5282 interpolated_A(j, i) +=
5283 nodal_position_gen(l, bulk_position_type(k), i) * dpsids(l, k, j);
5284 }
5285 }
5286 }
5287 }
5288 // Now find the local deformed metric tensor from the tangent Vectors
5289 DenseMatrix<double> A(n_dim_el, n_dim_el, 0.0);
5290 for (unsigned i = 0; i < n_dim_el; i++)
5291 {
5292 for (unsigned j = 0; j < n_dim_el; j++)
5293 {
5294 // Take the dot product
5295 for (unsigned k = 0; k < n_dim; k++)
5296 {
5297 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
5298 }
5299 }
5300 }
5301
5302 // Find the determinant of the metric tensor
5303 double Adet = 0.0;
5304 switch (n_dim_el)
5305 {
5306 case 1:
5307 Adet = A(0, 0);
5308 break;
5309 case 2:
5310 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
5311 break;
5312 default:
5313 throw OomphLibError("Wrong dimension in FaceElement",
5314 "FaceElement::J_eulerian()",
5315 OOMPH_EXCEPTION_LOCATION);
5316 }
5317
5318 // Return
5319 return sqrt(Adet);
5320 }
5321
5322
5323 //========================================================================
5324 /// Compute the Jacobian of the mapping between the local and global
5325 /// coordinates at the ipt-th integration point. Overloaded from
5326 /// FiniteElement.
5327 //========================================================================
5328 double FaceElement::J_eulerian_at_knot(const unsigned& ipt) const
5329 {
5330 // Find the number of nodes
5331 const unsigned n_node = nnode();
5332
5333 // Find the number of position types
5334 const unsigned n_position_type = nnodal_position_type();
5335
5336 // Find the dimension of the node and element
5337 const unsigned n_dim = nodal_dimension();
5338 const unsigned n_dim_el = dim();
5339
5340 // Set up dummy memory for the shape functions
5341 Shape psi(n_node, n_position_type);
5342 DShape dpsids(n_node, n_position_type, n_dim_el);
5343
5344 // Get the shape functions and local derivatives at the knot
5345 // This call may use the stored versions, which is why this entire
5346 // function doesn't just call J_eulerian(s), after reading out s from
5347 // the knots.
5348 dshape_local_at_knot(ipt, psi, dpsids);
5349
5350 // Also calculate the surface Vectors (derivatives wrt local coordinates)
5351 DenseMatrix<double> interpolated_A(n_dim_el, n_dim, 0.0);
5352
5353 // Calculate positions and derivatives
5354 for (unsigned l = 0; l < n_node; l++)
5355 {
5356 // Loop over positional dofs
5357 for (unsigned k = 0; k < n_position_type; k++)
5358 {
5359 // Loop over coordinates
5360 for (unsigned i = 0; i < n_dim; i++)
5361 {
5362 // Loop over LOCAL derivative directions, to calculate the tangent(s)
5363 for (unsigned j = 0; j < n_dim_el; j++)
5364 {
5365 interpolated_A(j, i) +=
5366 nodal_position_gen(l, bulk_position_type(k), i) * dpsids(l, k, j);
5367 }
5368 }
5369 }
5370 }
5371
5372 // Now find the local deformed metric tensor from the tangent Vectors
5373 DenseMatrix<double> A(n_dim_el, n_dim_el, 0.0);
5374 for (unsigned i = 0; i < n_dim_el; i++)
5375 {
5376 for (unsigned j = 0; j < n_dim_el; j++)
5377 {
5378 // Take the dot product
5379 for (unsigned k = 0; k < n_dim; k++)
5380 {
5381 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
5382 }
5383 }
5384 }
5385
5386 // Find the determinant of the metric tensor
5387 double Adet = 0.0;
5388 switch (n_dim_el)
5389 {
5390 case 1:
5391 Adet = A(0, 0);
5392 break;
5393 case 2:
5394 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
5395 break;
5396 default:
5397 throw OomphLibError("Wrong dimension in FaceElement",
5398 "FaceElement::J_eulerian_at_knot()",
5399 OOMPH_EXCEPTION_LOCATION);
5400 }
5401
5402 // Return
5403 return sqrt(Adet);
5404 }
5405
5406 //========================================================================
5407 /// Check that Jacobian of mapping between local and Eulerian
5408 /// coordinates at all integration points is positive.
5409 //========================================================================
5411 {
5412 // Let's be optimistic...
5413 passed = true;
5414
5415 // Find the number of nodes
5416 const unsigned n_node = nnode();
5417
5418 // Find the number of position types
5419 const unsigned n_position_type = nnodal_position_type();
5420
5421 // Find the dimension of the node and element
5422 const unsigned n_dim = nodal_dimension();
5423 const unsigned n_dim_el = dim();
5424
5425 // Set up dummy memory for the shape functions
5426 Shape psi(n_node, n_position_type);
5427 DShape dpsids(n_node, n_position_type, n_dim_el);
5428
5429
5430 unsigned nintpt = integral_pt()->nweight();
5431 for (unsigned ipt = 0; ipt < nintpt; ipt++)
5432 {
5433 // Get the shape functions and local derivatives at the knot
5434 // This call may use the stored versions, which is why this entire
5435 // function doesn't just call J_eulerian(s), after reading out s from
5436 // the knots.
5437 dshape_local_at_knot(ipt, psi, dpsids);
5438
5439 // Also calculate the surface Vectors (derivatives wrt local coordinates)
5440 DenseMatrix<double> interpolated_A(n_dim_el, n_dim, 0.0);
5441
5442 // Calculate positions and derivatives
5443 for (unsigned l = 0; l < n_node; l++)
5444 {
5445 // Loop over positional dofs
5446 for (unsigned k = 0; k < n_position_type; k++)
5447 {
5448 // Loop over coordinates
5449 for (unsigned i = 0; i < n_dim; i++)
5450 {
5451 // Loop over LOCAL derivative directions, to calculate the
5452 // tangent(s)
5453 for (unsigned j = 0; j < n_dim_el; j++)
5454 {
5455 interpolated_A(j, i) +=
5457 dpsids(l, k, j);
5458 }
5459 }
5460 }
5461 }
5462
5463 // Now find the local deformed metric tensor from the tangent Vectors
5464 DenseMatrix<double> A(n_dim_el, n_dim_el, 0.0);
5465 for (unsigned i = 0; i < n_dim_el; i++)
5466 {
5467 for (unsigned j = 0; j < n_dim_el; j++)
5468 {
5469 // Take the dot product
5470 for (unsigned k = 0; k < n_dim; k++)
5471 {
5472 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
5473 }
5474 }
5475 }
5476
5477 // Find the determinant of the metric tensor
5478 double Adet = 0.0;
5479 switch (n_dim_el)
5480 {
5481 case 1:
5482 Adet = A(0, 0);
5483 break;
5484 case 2:
5485 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
5486 break;
5487 default:
5488 throw OomphLibError("Wrong dimension in FaceElement",
5489 "FaceElement::J_eulerian_at_knot()",
5490 OOMPH_EXCEPTION_LOCATION);
5491 }
5492
5493
5494 // Are we dead yet?
5495 if (Adet <= 0.0)
5496 {
5497 passed = false;
5498 return;
5499 }
5500 }
5501 }
5502
5503 //=======================================================================
5504 /// Compute the tangent vector(s) and the outer unit normal
5505 /// vector at the specified local coordinate.
5506 /// In two spatial dimensions, a "tangent direction" is not required.
5507 /// In three spatial dimensions, a tangent direction is required
5508 /// (set via set_tangent_direction(...)), and we project the tanent direction
5509 /// on to the surface. The second tangent vector is taken to be the cross
5510 /// product of the projection and the unit normal.
5511 //=======================================================================
5513 const Vector<double>& s,
5514 Vector<Vector<double>>& tang_vec,
5515 Vector<double>& unit_normal) const
5516 {
5517 // Find the spatial dimension of the FaceElement
5518 const unsigned element_dim = dim();
5519
5520 // Find the overall dimension of the problem
5521 //(assume that it's the same for all nodes)
5522 const unsigned spatial_dim = nodal_dimension();
5523
5524#ifdef PARANOID
5525 // Check the number of local coordinates passed
5526 if (s.size() != element_dim)
5527 {
5528 std::ostringstream error_stream;
5529 error_stream
5530 << "Local coordinate s passed to outer_unit_normal() has dimension "
5531 << s.size() << std::endl
5532 << "but element dimension is " << element_dim << std::endl;
5533
5534 throw OomphLibError(
5535 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5536 }
5537
5538 // Check that if the Tangent_direction_pt is set, then
5539 // it is the same length as the spatial dimension.
5540 if (Tangent_direction_pt != 0 &&
5541 Tangent_direction_pt->size() != spatial_dim)
5542 {
5543 std::ostringstream error_stream;
5544 error_stream << "Tangent direction vector has dimension "
5545 << Tangent_direction_pt->size() << std::endl
5546 << "but spatial dimension is " << spatial_dim << std::endl;
5547
5548 throw OomphLibError(
5549 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5550 }
5551
5552 // Check the dimension of the normal vector
5553 if (unit_normal.size() != spatial_dim)
5554 {
5555 std::ostringstream error_stream;
5556 error_stream << "Unit normal passed to outer_unit_normal() has dimension "
5557 << unit_normal.size() << std::endl
5558 << "but spatial dimension is " << spatial_dim << std::endl;
5559
5560 throw OomphLibError(
5561 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5562 }
5563
5564
5565 // The number of tangent vectors.
5566 unsigned ntangent_vec = tang_vec.size();
5567
5568 // For the tangent vector,
5569 // if element_dim = 2, tang_vec is a 2D Vector of length 3.
5570 // if element_dim = 1, tang_vec is a 1D Vector of length 2.
5571 // if element_dim = 0, tang_vec is a 1D Vector of length 1.
5572 switch (element_dim)
5573 {
5574 // Point element, derived from a 1D element.
5575 case 0:
5576 {
5577 // Check that tang_vec is a 1D vector.
5578 if (ntangent_vec != 1)
5579 {
5580 std::ostringstream error_stream;
5581 error_stream
5582 << "This is a 0D FaceElement, we need one tangent vector.\n"
5583 << "You have given me " << tang_vec.size() << " vector(s).\n";
5584 throw OomphLibError(error_stream.str(),
5585 OOMPH_CURRENT_FUNCTION,
5586 OOMPH_EXCEPTION_LOCATION);
5587 }
5588 }
5589 break;
5590
5591 // Line element, derived from a 2D element.
5592 case 1:
5593 {
5594 // Check that tang_vec is a 1D vector.
5595 if (ntangent_vec != 1)
5596 {
5597 std::ostringstream error_stream;
5598 error_stream
5599 << "This is a 1D FaceElement, we need one tangent vector.\n"
5600 << "You have given me " << tang_vec.size() << " vector(s).\n";
5601 throw OomphLibError(error_stream.str(),
5602 OOMPH_CURRENT_FUNCTION,
5603 OOMPH_EXCEPTION_LOCATION);
5604 }
5605 }
5606 break;
5607
5608 // Plane element, derived from 3D element.
5609 case 2:
5610 {
5611 // Check that tang_vec is a 2D vector.
5612 if (ntangent_vec != 2)
5613 {
5614 std::ostringstream error_stream;
5615 error_stream
5616 << "This is a 2D FaceElement, we need two tangent vectors.\n"
5617 << "You have given me " << tang_vec.size() << " vector(s).\n";
5618 throw OomphLibError(error_stream.str(),
5619 OOMPH_CURRENT_FUNCTION,
5620 OOMPH_EXCEPTION_LOCATION);
5621 }
5622 }
5623 break;
5624 // There are no other elements!
5625 default:
5626
5627 throw OomphLibError(
5628 "Cannot have a FaceElement with dimension higher than 2",
5629 OOMPH_CURRENT_FUNCTION,
5630 OOMPH_EXCEPTION_LOCATION);
5631 break;
5632 }
5633
5634 // Check the lengths of each sub vectors.
5635 for (unsigned vec_i = 0; vec_i < ntangent_vec; vec_i++)
5636 {
5637 if (tang_vec[vec_i].size() != spatial_dim)
5638 {
5639 std::ostringstream error_stream;
5640 error_stream << "This problem has " << spatial_dim
5641 << " spatial dimensions.\n"
5642 << "But the " << vec_i << " tangent vector has size "
5643 << tang_vec[vec_i].size() << std::endl;
5644 throw OomphLibError(
5645 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5646 }
5647 }
5648
5649#endif
5650
5651
5652 // Now let's consider the different element dimensions
5653 switch (element_dim)
5654 {
5655 // Point element, derived from a 1D element.
5656 case 0:
5657 {
5658 std::ostringstream error_stream;
5659 error_stream
5660 << "It is unclear how to calculate a tangent and normal vectors for "
5661 << "a point element.\n";
5662 throw OomphLibError(
5663 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5664 }
5665 break;
5666
5667 // Line element, derived from a 2D element, in this case
5668 // the normal is a mess of cross products
5669 // We need an interior direction, so we must find the local
5670 // derivatives in the BULK element
5671 case 1:
5672 {
5673 // Find the number of nodes in the bulk element
5674 const unsigned n_node_bulk = Bulk_element_pt->nnode();
5675 // Find the number of position types in the bulk element
5676 const unsigned n_position_type_bulk =
5678
5679 // Construct the local coordinate in the bulk element
5680 Vector<double> s_bulk(2);
5681 // Get the local coordinates in the bulk element
5683
5684 // Allocate storage for the shape functions and their derivatives wrt
5685 // local coordinates
5686 Shape psi(n_node_bulk, n_position_type_bulk);
5687 DShape dpsids(n_node_bulk, n_position_type_bulk, 2);
5688 // Get the value of the shape functions at the given local coordinate
5689 Bulk_element_pt->dshape_local(s_bulk, psi, dpsids);
5690
5691 // Calculate all derivatives of the spatial coordinates
5692 // wrt local coordinates
5693 DenseMatrix<double> interpolated_dxds(2, spatial_dim, 0.0);
5694
5695 // Loop over all parent nodes
5696 for (unsigned l = 0; l < n_node_bulk; l++)
5697 {
5698 // Loop over all position types in the bulk
5699 for (unsigned k = 0; k < n_position_type_bulk; k++)
5700 {
5701 // Loop over derivative direction
5702 for (unsigned j = 0; j < 2; j++)
5703 {
5704 // Loop over coordinate directions
5705 for (unsigned i = 0; i < spatial_dim; i++)
5706 {
5707 // Compute the spatial derivative
5708 interpolated_dxds(j, i) +=
5710 dpsids(l, k, j);
5711 }
5712 }
5713 }
5714 }
5715
5716 // Initialise the tangent, interior tangent and normal vectors to zero
5717 // The idea is that even if the element is in a two-dimensional space,
5718 // the normal cannot be calculated without embedding the element in
5719 // three dimensions, in which case, the tangent and interior tangent
5720 // will have zero z-components.
5721 Vector<double> tangent(3, 0.0), interior_tangent(3, 0.0),
5722 normal(3, 0.0);
5723
5724 // We must get the relationship between the coordinate along the face
5725 // and the local coordinates in the bulk element
5726 // We must also find an interior direction
5727 DenseMatrix<double> dsbulk_dsface(2, 1, 0.0);
5728 unsigned interior_direction = 0;
5729 get_ds_bulk_ds_face(s, dsbulk_dsface, interior_direction);
5730 // Load in the values for the tangents
5731 for (unsigned i = 0; i < spatial_dim; i++)
5732 {
5733 // Tangent to the face is the derivative wrt to the face coordinate
5734 // which is calculated using dsbulk_dsface and the chain rule
5735 tangent[i] = interpolated_dxds(0, i) * dsbulk_dsface(0, 0) +
5736 interpolated_dxds(1, i) * dsbulk_dsface(1, 0);
5737 // Interior tangent to the face is the derivative in the interior
5738 // direction
5739 interior_tangent[i] = interpolated_dxds(interior_direction, i);
5740 }
5741
5742 // Now the (3D) normal to the element is the interior tangent
5743 // crossed with the tangent
5744 VectorHelpers::cross(interior_tangent, tangent, normal);
5745
5746 // We find the line normal by crossing the element normal with the
5747 // tangent
5748 Vector<double> full_normal(3);
5749 VectorHelpers::cross(normal, tangent, full_normal);
5750
5751 // Copy the appropriate entries into the unit normal
5752 // Two or Three depending upon the spatial dimension of the system
5753 for (unsigned i = 0; i < spatial_dim; i++)
5754 {
5755 unit_normal[i] = full_normal[i];
5756 }
5757
5758 // Finally normalise unit normal and multiply by the Normal_sign
5759 double length = VectorHelpers::magnitude(unit_normal);
5760 for (unsigned i = 0; i < spatial_dim; i++)
5761 {
5762 unit_normal[i] *= Normal_sign / length;
5763 }
5764
5765 // Create the tangent vector
5766 tang_vec[0][0] = -unit_normal[1];
5767 tang_vec[0][1] = unit_normal[0];
5768 }
5769 break;
5770
5771 // Plane element, derived from 3D element, in this case the normal
5772 // is just the cross product of the two surface tangents
5773 // We assume, therefore, that we have three spatial coordinates