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