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-2023 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 
38 namespace 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  //=======================================================================
62  unsigned GeneralisedElement::add_internal_data(Data* const& data_pt,
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
144  ++Ninternal_data;
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  //=======================================================================
307  unsigned GeneralisedElement::add_external_data(Data* const& data_pt,
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;
326  oomph_info
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
377  ++Nexternal_data;
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
498  --Nexternal_data;
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  {
536  internal_data_pt(i)->assign_eqn_numbers(global_number, Dof_pt);
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());
559  internal_data_pt(i)->describe_dofs(out, in);
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
786  ElementWithExternalElement* e_el_pt =
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
1032  local_eqn_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
1069  local_eqn_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
1654  unsigned Max_newton_iterations = 10;
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  //=========================================================================
3274  void FiniteElement::d2shape_local_at_knot(const unsigned& ipt,
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  //==========================================================================
3502  double FiniteElement::d2shape_eulerian_at_knot(const unsigned& ipt,
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
3632  local_eqn_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  //==================================================================
4378  FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt,
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 
4445  if (GeneralisedElement::self_test() != 0)
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,
4736  Vector<double>& s,
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() =
5155  face_to_bulk_coordinate_fct_pt(face_index);
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
5774  // and two surface coordinates
5775  // Then we need only to get the derivatives wrt the local coordinates
5776  // in this face element
5777  case 2:
5778  {
5779 #ifdef PARANOID
5780  // Check that we actually have three spatial dimensions
5781  if (spatial_dim != 3)
5782  {
5783  std::ostringstream error_stream;
5784  error_stream << "There are only " << spatial_dim
5785  << "coordinates at the nodes of this 2D FaceElement,\n"
5786  << "which must have come from a 3D Bulk element\n";
5787  throw OomphLibError(error_stream.str(),
5788  OOMPH_CURRENT_FUNCTION,
5789  OOMPH_EXCEPTION_LOCATION);
5790  }
5791 #endif
5792 
5793  // Find the number of nodes in the element
5794  const unsigned n_node = this->nnode();
5795  // Find the number of position types
5796  const unsigned n_position_type = this->nnodal_position_type();
5797 
5798  // Allocate storage for the shape functions and their derivatives wrt
5799  // local coordinates
5800  Shape psi(n_node, n_position_type);
5801  DShape dpsids(n_node, n_position_type, 2);
5802  // Get the value of the shape functions at the given local coordinate
5803  this->dshape_local(s, psi, dpsids);
5804 
5805  // Calculate all derivatives of the spatial coordinates
5806  // wrt local coordinates
5807  Vector<Vector<double>> interpolated_dxds(2, Vector<double>(3, 0));
5808 
5809  // Loop over all nodes
5810  for (unsigned l = 0; l < n_node; l++)
5811  {
5812  // Loop over all position types
5813  for (unsigned k = 0; k < n_position_type; k++)
5814  {
5815  // Loop over derivative directions
5816  for (unsigned j = 0; j < 2; j++)
5817  {
5818  // Loop over coordinate directions
5819  for (unsigned i = 0; i < 3; i++)
5820  {
5821  // Compute the spatial derivative
5822  // Remember that we need to translate the position type
5823  // to its location in the bulk node
5824  interpolated_dxds[j][i] +=
5825  this->nodal_position_gen(l, bulk_position_type(k), i) *
5826  dpsids(l, k, j);
5827  }
5828  }
5829  }
5830  }
5831 
5832  // We now take the cross product of the two normal vectors
5834  interpolated_dxds[0], interpolated_dxds[1], unit_normal);
5835 
5836  // Finally normalise unit normal
5837  double normal_length = VectorHelpers::magnitude(unit_normal);
5838 
5839  for (unsigned i = 0; i < spatial_dim; i++)
5840  {
5841  unit_normal[i] *= Normal_sign / normal_length;
5842  }
5843 
5844  // Next we create the continuous tangent vectors!
5845  if (Tangent_direction_pt != 0)
5846  // There is a general direction that the first tangent vector should
5847  // follow.
5848  {
5849  // Project the Tangent direction vector onto the surface.
5850  // Project Tangent_direction_pt D onto the plane P defined by
5851  // T1 and T2:
5852  // proj_P(D) = proj_T1(D) + proj_T2(D), where D is
5853  // Tangent_direction_pt, recall that proj_u(v) = (u.v)/(u.u) * u
5854 
5855  // Get the direction vector. The vector is NOT copied! :)
5856  Vector<double>& direction_vector = *Tangent_direction_pt;
5857 
5858 #ifdef PARANOID
5859  // Check that the angle between the direction vector and the normal
5860  // is not less than 10 degrees
5861  if (VectorHelpers::angle(direction_vector, unit_normal) <
5862  (10.0 * MathematicalConstants::Pi / 180.0))
5863  {
5864  std::ostringstream err_stream;
5865  err_stream << "The angle between the unit normal and the \n"
5866  << "direction vector is less than 10 degrees.";
5867  throw OomphLibError(err_stream.str(),
5868  OOMPH_CURRENT_FUNCTION,
5869  OOMPH_EXCEPTION_LOCATION);
5870  }
5871 #endif
5872 
5873  // Calculate the two scalings, (u.v) / (u.u)
5874  double t1_scaling =
5875  VectorHelpers::dot(interpolated_dxds[0], direction_vector) /
5876  VectorHelpers::dot(interpolated_dxds[0], interpolated_dxds[0]);
5877 
5878  double t2_scaling =
5879  VectorHelpers::dot(interpolated_dxds[1], direction_vector) /
5880  VectorHelpers::dot(interpolated_dxds[1], interpolated_dxds[1]);
5881 
5882  // Finish off the projection.
5883  tang_vec[0][0] = t1_scaling * interpolated_dxds[0][0] +
5884  t2_scaling * interpolated_dxds[1][0];
5885  tang_vec[0][1] = t1_scaling * interpolated_dxds[0][1] +
5886  t2_scaling * interpolated_dxds[1][1];
5887  tang_vec[0][2] = t1_scaling * interpolated_dxds[0][2] +
5888  t2_scaling * interpolated_dxds[1][2];
5889 
5890  // The second tangent vector is the cross product of
5891  // tang_vec[0] and the normal vector N.
5892  VectorHelpers::cross(tang_vec[0], unit_normal, tang_vec[1]);
5893 
5894  // Normalise the tangent vectors
5895  for (unsigned vec_i = 0; vec_i < 2; vec_i++)
5896  {
5897  // Get the length...
5898  double tang_length = VectorHelpers::magnitude(tang_vec[vec_i]);
5899 
5900  for (unsigned dim_i = 0; dim_i < spatial_dim; dim_i++)
5901  {
5902  tang_vec[vec_i][dim_i] /= tang_length;
5903  }
5904  }
5905  }
5906  else
5907  // A direction for the first tangent vector is not supplied, we do our
5908  // best to keep it continuous. But if the normal vector is aligned with
5909  // the z axis, then the tangent vector may "flip", even within elements.
5910  // This is because we check this by comparing doubles.
5911  // The code is copied from impose_parallel_outflow_element.h,
5912  // NO modification is done except for a few variable name changes.
5913  {
5914  double a, b, c;
5915  a = unit_normal[0];
5916  b = unit_normal[1];
5917  c = unit_normal[2];
5918 
5919  if (a != 0.0 || b != 0.0)
5920  {
5921  double a_sq = a * a;
5922  double b_sq = b * b;
5923  double c_sq = c * c;
5924 
5925  tang_vec[0][0] = -b / sqrt(a_sq + b_sq);
5926  tang_vec[0][1] = a / sqrt(a_sq + b_sq);
5927  tang_vec[0][2] = 0;
5928 
5929  double z = (a_sq + b_sq) / sqrt(a_sq * c_sq + b_sq * c_sq +
5930  (a_sq + b_sq) * (a_sq + b_sq));
5931 
5932  tang_vec[1][0] = -(a * c * z) / (a_sq + b_sq);
5933  tang_vec[1][1] = -(b * c * z) / (a_sq + b_sq);
5934  tang_vec[1][2] = z;
5935  // NB : we didn't use the fact that N is normalized,
5936  // that's why we have these unsimplified formulas
5937  }
5938  else if (c != 0.0)
5939  {
5941  {
5942  std::ostringstream warn_stream;
5943  warn_stream
5944  << "I have detected that your normal vector is [0,0,1].\n"
5945  << "Since this function compares against 0.0, you may\n"
5946  << "get discontinuous tangent vectors.";
5947  OomphLibWarning(warn_stream.str(),
5948  OOMPH_CURRENT_FUNCTION,
5949  OOMPH_EXCEPTION_LOCATION);
5950  }
5951 
5952  tang_vec[0][0] = 1.0;
5953  tang_vec[0][1] = 0.0;
5954  tang_vec[0][2] = 0.0;
5955 
5956  tang_vec[1][0] = 0.0;
5957  tang_vec[1][1] = 1.0;
5958  tang_vec[1][2] = 0.0;
5959  }
5960  else
5961  {
5962  throw OomphLibError("You have a zero normal vector!! ",
5963  OOMPH_CURRENT_FUNCTION,
5964  OOMPH_EXCEPTION_LOCATION);
5965  }
5966  }
5967  }
5968  break;
5969 
5970  default:
5971 
5972  throw OomphLibError(
5973  "Cannot have a FaceElement with dimension higher than 2",
5974  OOMPH_CURRENT_FUNCTION,
5975  OOMPH_EXCEPTION_LOCATION);
5976  break;
5977  }
5978  }
5979 
5980  //=======================================================================
5981  /// Compute the tangent vector(s) and the outer unit normal
5982  /// vector at the ipt-th integration point. This is a wrapper around
5983  /// continuous_tangent_and_outer_unit_normal(...) with the integration points
5984  /// converted into local coordinates.
5985  //=======================================================================
5987  const unsigned& ipt,
5988  Vector<Vector<double>>& tang_vec,
5989  Vector<double>& unit_normal) const
5990  {
5991  // Find the dimension of the element
5992  const unsigned element_dim = dim();
5993  // Find the local coordinates of the ipt-th integration point
5994  Vector<double> s(element_dim);
5995  for (unsigned i = 0; i < element_dim; i++)
5996  {
5997  s[i] = integral_pt()->knot(ipt, i);
5998  }
5999  // Call the outer unit normal function
6000  continuous_tangent_and_outer_unit_normal(s, tang_vec, unit_normal);
6001  }
6002 
6003  //=======================================================================
6004  /// Compute the outer unit normal at the specified local coordinate
6005  //=======================================================================
6007  Vector<double>& unit_normal) const
6008  {
6009  // Find the spatial dimension of the FaceElement
6010  const unsigned element_dim = dim();
6011 
6012  // Find the overall dimension of the problem
6013  //(assume that it's the same for all nodes)
6014  const unsigned spatial_dim = nodal_dimension();
6015 
6016 #ifdef PARANOID
6017  // Check the number of local coordinates passed
6018  if (s.size() != element_dim)
6019  {
6020  std::ostringstream error_stream;
6021  error_stream
6022  << "Local coordinate s passed to outer_unit_normal() has dimension "
6023  << s.size() << std::endl
6024  << "but element dimension is " << element_dim << std::endl;
6025 
6026  throw OomphLibError(
6027  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6028  }
6029 
6030  // Check the dimension of the normal vector
6031  if (unit_normal.size() != spatial_dim)
6032  {
6033  std::ostringstream error_stream;
6034  error_stream << "Unit normal passed to outer_unit_normal() has dimension "
6035  << unit_normal.size() << std::endl
6036  << "but spatial dimension is " << spatial_dim << std::endl;
6037 
6038  throw OomphLibError(
6039  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6040  }
6041 #endif
6042 
6043  /* //The spatial dimension of the bulk element will be element_dim+1
6044  const unsigned bulk_dim = element_dim + 1;
6045 
6046  //Find the number of nodes in the bulk element
6047  const unsigned n_node_bulk = Bulk_element_pt->nnode();
6048  //Find the number of position types in the bulk element
6049  const unsigned n_position_type_bulk =
6050  Bulk_element_pt->nnodal_position_type();
6051 
6052  //Construct the local coordinate in the bulk element
6053  Vector<double> s_bulk(bulk_dim);
6054  //Set the value of the bulk coordinate that is fixed on the face
6055  //s_bulk[s_fixed_index()] = s_fixed_value();
6056 
6057  //Set the other bulk coordinates
6058  //for(unsigned i=0;i<element_dim;i++) {s_bulk[bulk_s_index(i)] = s[i];}
6059 
6060  get_local_coordinate_in_bulk(s,s_bulk);
6061 
6062  //Allocate storage for the shape functions and their derivatives wrt
6063  //local coordinates
6064  Shape psi(n_node_bulk,n_position_type_bulk);
6065  DShape dpsids(n_node_bulk,n_position_type_bulk,bulk_dim);
6066  //Get the value of the shape functions at the given local coordinate
6067  Bulk_element_pt->dshape_local(s_bulk,psi,dpsids);
6068 
6069  //Calculate all derivatives of the spatial coordinates wrt local
6070  coordinates DenseMatrix<double> interpolated_dxds(bulk_dim,spatial_dim);
6071  //Initialise to zero
6072  for(unsigned j=0;j<bulk_dim;j++)
6073  {for(unsigned i=0;i<spatial_dim;i++) {interpolated_dxds(j,i) = 0.0;}}
6074 
6075  //Loop over all parent nodes
6076  for(unsigned l=0;l<n_node_bulk;l++)
6077  {
6078  //Loop over all position types in the bulk
6079  for(unsigned k=0;k<n_position_type_bulk;k++)
6080  {
6081  //Loop over derivative direction
6082  for(unsigned j=0;j<bulk_dim;j++)
6083  {
6084  //Loop over coordinate directions
6085  for(unsigned i=0;i<spatial_dim;i++)
6086  {
6087  //Compute the spatial derivative
6088  interpolated_dxds(j,i) +=
6089  Bulk_element_pt->nodal_position_gen(l,k,i)*dpsids(l,k,j);
6090  }
6091  }
6092  }
6093  }*/
6094 
6095  // Now let's consider the different element dimensions
6096  switch (element_dim)
6097  {
6098  // Point element, derived from a 1D element, in this case
6099  // the normal is merely the tangent to the bulk element
6100  // and there is only one free coordinate in the bulk element
6101  // Hence we will need to calculate the derivatives wrt the
6102  // local coordinates in the BULK element.
6103  case 0:
6104  {
6105  // Find the number of nodes in the bulk element
6106  const unsigned n_node_bulk = Bulk_element_pt->nnode();
6107  // Find the number of position types in the bulk element
6108  const unsigned n_position_type_bulk =
6110 
6111