elements.h
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 //Header file for classes that define element objects
31 
32 //Include guard to prevent multiple inclusions of the header
33 #ifndef OOMPH_GENERIC_ELEMENTS_HEADER
34 #define OOMPH_GENERIC_ELEMENTS_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38  #include <oomph-lib-config.h>
39 #endif
40 
41 #include <map>
42 #include <deque>
43 #include<string>
44 #include <list>
45 
46 //oomph-lib includes
47 #include "integral.h"
48 #include "nodes.h"
49 #include "geom_objects.h"
50 
51 namespace oomph
52 {
53  //Forward definition for time
54  class Time;
55 
56 
57 //========================================================================
58 /// \short A Generalised Element class.
59 ///
60 /// The main components of a GeneralisedElement are:
61 /// - pointers to its internal and external Data, which are stored together
62 /// in a single array so that we need only store one pointer.
63 /// The internal Data are placed at the beginning
64 /// of the array and the external Data are stored at the end.
65 /// - a pointer to a global Time
66 /// - a lookup table that establishes the relation between local
67 /// and global equation numbers.
68 ///
69 /// We also provide interfaces for functions that compute the
70 /// element's Jacobian matrix and/or the Vector of residuals.
71 /// In addition, an interface that returns a mass matrix --- the matrix
72 /// of terms that multiply any time derivatives in the problem --- is
73 /// also provided to permit explicit time-stepping and the solution
74 /// of the generalised eigenproblems.
75 //========================================================================
77 {
78  private:
79 
80  /// \short Storage for the global equation numbers represented by the
81  /// degrees of freedom in the element.
82  unsigned long *Eqn_number;
83 
84  /// \short Storage for array of pointers to degrees of freedom ordered
85  /// by local equation number. This information is not needed, except in
86  /// continuation, bifurcation tracking and periodic orbit calculations.
87  /// It is not set up unless the control flag
88  /// Problem::Store_local_dof_pts_in_elements = true
89  double **Dof_pt;
90 
91  /// \short Storage for pointers to internal and external data.
92  /// The data is of the same type and so can be addressed by
93  /// a single pointer. The idea is that the array is of
94  /// total size Ninternal_data + Nexternal_data.
95  /// The internal data are stored at the beginning of the array
96  /// and the external data are stored at the end of the array.
98 
99  /// \short Pointer to array storage for local equation numbers associated
100  /// with internal and external data. Again, we save storage by using
101  /// a single pointer to access this information. The first index of the
102  /// array is of dimension Nineternal_data + Nexternal_data and the second
103  /// index varies with the number of values stored at the data object.
104  /// In the most general case, however, the scheme is as memory efficient
105  /// as possible.
107 
108  /// Number of degrees of freedom
109  unsigned Ndof;
110 
111  /// Number of internal data
112  unsigned Ninternal_data;
113 
114  /// Number of external data
115  unsigned Nexternal_data;
116 
117  /// \short Storage for boolean flags of size Ninternal_data + Nexternal_data
118  /// that correspond to the data used in the element. The flags
119  /// indicate whether the particular
120  /// internal or external data should be included in the general
121  /// finite-difference loop in fill_in_jacobian_from_internal_by_fd() or
122  /// fill_in_jacobian_from_external_by_fd(). The default is that all
123  /// data WILL be included in the finite-difference loop, but in many
124  /// circumstances it is possible to treat certain (external) data
125  /// analytically. The use of an STL vector is optimal for memory
126  /// use because the booleans are represented as single-bits.
127  std::vector<bool> Data_fd;
128 
129  protected:
130 
131 #ifdef OOMPH_HAS_MPI
132 
133  /// \short Non-halo processor ID for Data; -1 if it's not a halo.
135 
136  /// Does this element need to be kept as a halo element during a distribute?
138 
139 #endif
140 
141  /// \short Add a (pointer to an) internal data object to the element and
142  /// return the index required to obtain it from the access
143  /// function \c internal_data_pt(). The boolean indicates whether
144  /// the datum should be included in the general finite-difference loop
145  /// when calculating the jacobian. The default value is true, i.e.
146  /// the data will be included in the finite differencing.
147  unsigned add_internal_data(Data* const &data_pt, const bool &fd=true);
148 
149 
150  /// \short Return the status of the boolean flag indicating whether
151  /// the internal data is included in the finite difference loop
152  inline bool internal_data_fd(const unsigned &i) const
153  {
154 #ifdef RANGE_CHECKING
155  if(i >= Ninternal_data)
156  {
157  std::ostringstream error_message;
158  error_message << "Range Error: Internal data " << i
159  << " is not in the range (0,"
160  << Ninternal_data-1 << ")";
161  throw OomphLibError(error_message.str(),
162  OOMPH_CURRENT_FUNCTION,
163  OOMPH_EXCEPTION_LOCATION);
164  }
165 #endif
166  //Return the value
167  return Data_fd[i];
168  }
169 
170 
171  /// \short Set the boolean flag to exclude the internal datum from
172  /// the finite difference loop when computing the jacobian matrix
173  inline void exclude_internal_data_fd(const unsigned &i)
174  {
175 #ifdef RANGE_CHECKING
176  if(i >= Ninternal_data)
177  {
178  std::ostringstream error_message;
179  error_message << "Range Error: Internal data " << i
180  << " is not in the range (0,"
181  << Ninternal_data-1 << ")";
182  throw OomphLibError(error_message.str(),
183  OOMPH_CURRENT_FUNCTION,
184  OOMPH_EXCEPTION_LOCATION);
185  }
186 #endif
187  //Set the value
188  Data_fd[i] = false;
189  }
190 
191  /// \short Set the boolean flag to include the internal datum in
192  /// the finite difference loop when computing the jacobian matrix
193  inline void include_internal_data_fd(const unsigned &i)
194  {
195 #ifdef RANGE_CHECKING
196  if(i >= Ninternal_data)
197  {
198  std::ostringstream error_message;
199  error_message << "Range Error: Internal data " << i
200  << " is not in the range (0,"
201  << Ninternal_data-1 << ")";
202  throw OomphLibError(error_message.str(),
203  OOMPH_CURRENT_FUNCTION,
204  OOMPH_EXCEPTION_LOCATION);
205  }
206 #endif
207  //Set the value
208  Data_fd[i] = true;
209  }
210 
211  /// \short Clear the storage for the global equation numbers
212  /// and pointers to dofs (if stored)
214  {delete[] Eqn_number; Eqn_number=0; delete[] Dof_pt; Dof_pt=0; Ndof=0;}
215 
216  /// \short Add the contents of the queue global_eqn_numbers
217  /// to the local storage for the local-to-global translation scheme.
218  /// It is essential that the entries in the queue are added IN ORDER
219  /// i.e. from the front.
220  void add_global_eqn_numbers(std::deque<unsigned long>
221  const &global_eqn_numbers,
222  std::deque<double*>
223  const &global_dof_pt);
224 
225  /// \short Empty dense matrix used as a dummy argument to combined
226  /// residual and jacobian functions in the case when only the residuals
227  /// are being assembled
229 
230  /// \short Static storage for deque used to add_global_equation_numbers
231  /// when pointers to the dofs in each element are not required
232  static std::deque<double*> Dof_pt_deque;
233 
234  /// \short Assign the local equation numbers for the internal
235  /// and external Data
236  /// This must be called after the global equation numbers have all been
237  /// assigned. It is virtual so that it can be overloaded by
238  /// ElementWithExternalElements so that any external data from the
239  /// external elements in included in the numbering scheme.
240  /// If the boolean argument is true then pointers to the dofs will be
241  /// stored in Dof_pt
243  const bool &store_local_dof_pt);
244 
245  /// \short Assign all the local equation numbering schemes that can
246  /// be applied generically for the element. In most cases, this is the
247  /// function that will be overloaded by inherited classes. It is required
248  /// to ensure that assign_additional_local_eqn_numbers() can always be
249  /// called after ALL other local equation numbering has been performed.
250  /// The default for the GeneralisedElement is simply to call internal
251  /// and external local equation numbering.
252  /// If the boolean argument is true then pointers to the dofs will be stored
253  /// in Dof_pt
255  const bool &store_local_dof_pt)
256  {
257  this->assign_internal_and_external_local_eqn_numbers(store_local_dof_pt);
258  }
259 
260  /// \short Setup any additional look-up schemes for local equation numbers.
261  /// Examples of use include using local storage to refer to explicit
262  /// degrees of freedom. The additional memory cost of such storage may
263  /// or may not be offset by fast local access.
265 
266  /// \short Return the local equation number corresponding to the j-th
267  /// value stored at the i-th internal data
268  inline int internal_local_eqn(const unsigned &i, const unsigned &j) const
269  {
270 #ifdef RANGE_CHECKING
271  if(i >= Ninternal_data)
272  {
273  std::ostringstream error_message;
274  error_message << "Range Error: Internal data " << i
275  << " is not in the range (0,"
276  << Ninternal_data-1 << ")";
277  throw OomphLibError(error_message.str(),
278  OOMPH_CURRENT_FUNCTION,
279  OOMPH_EXCEPTION_LOCATION);
280  }
281  else
282  {
283  unsigned n_value = internal_data_pt(i)->nvalue();
284  if(j >= n_value)
285  {
286  std::ostringstream error_message;
287  error_message << "Range Error: value " << j << " at internal data "
288  << i
289  << " is not in the range (0,"
290  << n_value -1 << ")";
291  throw OomphLibError(error_message.str(),
292  OOMPH_CURRENT_FUNCTION,
293  OOMPH_EXCEPTION_LOCATION);
294  }
295  }
296 #endif
297  //Check that the data has been allocated
298 #ifdef PARANOID
299  if(Data_local_eqn==0)
300  {
301  throw OomphLibError(
302  "Internal local equation numbers have not been allocated",
303  OOMPH_CURRENT_FUNCTION,
304  OOMPH_EXCEPTION_LOCATION);
305  }
306 #endif
307  //Return the local equation number as stored in the Data_local_eqn array
308  return Data_local_eqn[i][j];
309  }
310 
311  /// \short Return the local equation number corresponding to the j-th
312  /// value stored at the i-th external data
313  inline int external_local_eqn(const unsigned &i, const unsigned &j)
314  {
315  //Check that the data has been allocated
316 #ifdef RANGE_CHECKING
317  if(i >= Nexternal_data)
318  {
319  std::ostringstream error_message;
320  error_message << "Range Error: External data " << i
321  << " is not in the range (0,"
322  << Nexternal_data-1 << ")";
323  throw OomphLibError(error_message.str(),
324  OOMPH_CURRENT_FUNCTION,
325  OOMPH_EXCEPTION_LOCATION);
326  }
327  else
328  {
329  unsigned n_value = external_data_pt(i)->nvalue();
330  if(j >= n_value)
331  {
332  std::ostringstream error_message;
333  error_message << "Range Error: value " << j << " at internal data "
334  << i
335  << " is not in the range (0,"
336  << n_value -1 << ")";
337  throw OomphLibError(error_message.str(),
338  OOMPH_CURRENT_FUNCTION,
339  OOMPH_EXCEPTION_LOCATION);
340  }
341  }
342 #endif
343 #ifdef PARANOID
344  if(Data_local_eqn==0)
345  {
346  throw OomphLibError(
347  "External local equation numbers have not been allocated",
348  OOMPH_CURRENT_FUNCTION,
349  OOMPH_EXCEPTION_LOCATION);
350  }
351 #endif
352  //Return the local equation number stored as the j-th value of the
353  //i-th external data object.
354  return Data_local_eqn[Ninternal_data +i][j];
355  }
356 
357  /// \short Add the elemental contribution to the residuals vector. Note that
358  /// this function will NOT initialise the residuals vector. It must be
359  /// called after the residuals vector has been initialised to zero.
361  {
362  std::string error_message =
363  "Empty fill_in_contribution_to_residuals() has been called.\n";
364  error_message +=
365  "This function is called from the default implementations of\n";
366  error_message +=
367  "get_residuals() and get_jacobian();\n";
368  error_message +=
369  "and must calculate the residuals vector without initialising any of ";
370  error_message += "its entries.\n\n";
371 
372  error_message +=
373  "If you do not wish to use these defaults, you must overload both\n";
374  error_message +=
375  "get_residuals() and get_jacobian(), which must initialise the entries\n";
376  error_message +=
377  "of the residuals vector and jacobian matrix to zero.\n";
378 
379  error_message +=
380  "N.B. the default get_jacobian() function employs finite differencing\n";
381 
382  throw OomphLibError(error_message,
383  OOMPH_CURRENT_FUNCTION,
384  OOMPH_EXCEPTION_LOCATION);
385  }
386 
387  /// \short Calculate the contributions to the jacobian from the internal
388  /// degrees of freedom using finite differences.
389  /// This version of the function assumes that the residuals vector has
390  /// already been calculated. If the boolean argument is true, the finite
391  /// differencing will be performed for all internal data, irrespective of
392  /// the information in Data_fd. The default value (false)
393  /// uses the information in Data_fd to selectively difference only
394  /// certain data.
396  DenseMatrix<double> &jacobian,
397  const bool &fd_all_data=false);
398 
399  /// \short Calculate the contributions to the jacobian from the internal
400  /// degrees of freedom using finite differences. This version computes
401  /// the residuals vector before calculating the jacobian terms.
402  /// If the boolean argument is true, the finite
403  /// differencing will be performed for all internal data, irrespective of
404  /// the information in Data_fd. The default value (false)
405  /// uses the information in Data_fd to selectively difference only
406  /// certain data.
408  const bool &fd_all_data=false)
409  {
410  //Allocate storage for the residuals
411  Vector<double> residuals(Ndof,0.0);
412  //Get the residuals for the entire element
413  get_residuals(residuals);
414  //Call the jacobian calculation
415  fill_in_jacobian_from_internal_by_fd(residuals,jacobian,fd_all_data);
416  }
417 
418 
419  /// \short Calculate the contributions to the jacobian from the external
420  /// degrees of freedom using finite differences.
421  /// This version of the function assumes that the residuals vector has
422  /// already been calculated.
423  /// If the boolean argument is true, the finite
424  /// differencing will be performed for all external data, irrespective of
425  /// the information in Data_fd. The default value (false)
426  /// uses the information in Data_fd to selectively difference only
427  /// certain data.
429  DenseMatrix<double> &jacobian,
430  const bool &fd_all_data=false);
431 
432 
433  /// \short Calculate the contributions to the jacobian from the external
434  /// degrees of freedom using finite differences. This version computes
435  /// the residuals vector before calculating the jacobian terms.
436  /// If the boolean argument is true, the finite
437  /// differencing will be performed for all internal data, irrespective of
438  /// the information in Data_fd. The default value (false)
439  /// uses the information in Data_fd to selectively difference only
440  /// certain data.
442  const bool &fd_all_data=false)
443  {
444  //Allocate storage for a residuals vector
445  Vector<double> residuals(Ndof);
446  //Get the residuals for the entire element
447  get_residuals(residuals);
448  //Call the jacobian calculation
449  fill_in_jacobian_from_external_by_fd(residuals,jacobian,fd_all_data);
450  }
451 
452  /// \short Function that is called before the finite differencing of
453  /// any internal data. This may be overloaded to update any slaved
454  /// data before finite differencing takes place.
455  virtual inline void update_before_internal_fd() { }
456 
457  /// \short Function that is call after the finite differencing of
458  /// the internal data. This may be overloaded to reset any slaved
459  /// variables that may have changed during the finite differencing.
460  virtual inline void reset_after_internal_fd() { }
461 
462  /// \short Function called within the finite difference loop for
463  /// internal data after a change in any values in the i-th
464  /// internal data object.
465  virtual inline void update_in_internal_fd(const unsigned &i) { }
466 
467  /// \short Function called within the finite difference loop for
468  /// internal data after the values in the i-th external data object
469  /// are reset. The default behaviour is to call the update function.
470  virtual inline void reset_in_internal_fd(const unsigned &i)
472 
473 
474  /// \short Function that is called before the finite differencing of
475  /// any external data. This may be overloaded to update any slaved
476  /// data before finite differencing takes place.
477  virtual inline void update_before_external_fd() { }
478 
479  /// \short Function that is call after the finite differencing of
480  /// the external data. This may be overloaded to reset any slaved
481  /// variables that may have changed during the finite differencing.
482  virtual inline void reset_after_external_fd() { }
483 
484  /// \short Function called within the finite difference loop for
485  /// external data after a change in any values in the i-th
486  /// external data object
487  virtual inline void update_in_external_fd(const unsigned &i) { }
488 
489  /// \short Function called within the finite difference loop for
490  /// external data after the values in the i-th external data object
491  /// are reset. The default behaviour is to call the update function.
492  virtual inline void reset_in_external_fd(const unsigned &i)
494 
495  /// \short Add the elemental contribution to the jacobian matrix.
496  /// and the residuals vector. Note that
497  /// this function will NOT initialise the residuals vector or the jacobian
498  /// matrix. It must be called after the residuals vector and
499  /// jacobian matrix have been initialised to zero. The default
500  /// is to use finite differences to calculate the jacobian
502  DenseMatrix<double> &jacobian)
503  {
504  //Add the contribution to the residuals
506  //Allocate storage for the full residuals (residuals of entire element)
507  Vector<double> full_residuals(Ndof);
508  //Get the residuals for the entire element
509  get_residuals(full_residuals);
510  //Calculate the contributions from the internal dofs
511  //(finite-difference the lot by default)
512  fill_in_jacobian_from_internal_by_fd(full_residuals,jacobian,true);
513  //Calculate the contributions from the external dofs
514  //(finite-difference the lot by default)
515  fill_in_jacobian_from_external_by_fd(full_residuals,jacobian,true);
516  }
517 
518  /// \short Add the elemental contribution to the mass matrix matrix.
519  /// and the residuals vector. Note that
520  /// this function should NOT initialise the residuals vector or the mass
521  /// matrix. It must be called after the residuals vector and
522  /// jacobian matrix have been initialised to zero. The default
523  /// is deliberately broken
525  Vector<double> &residuals, DenseMatrix<double> &mass_matrix);
526 
527  /// \short Add the elemental contribution to the jacobian matrix,
528  /// mass matrix and the residuals vector. Note that
529  /// this function should NOT initialise any entries.
530  /// It must be called after the residuals vector and
531  /// matrices have been initialised to zero.
533  Vector<double> &residuals,
534  DenseMatrix<double> &jacobian, DenseMatrix<double> &mass_matrix);
535 
536  /// \short Add the elemental contribution to the derivatives of
537  /// the residuals with respect to a parameter. This function should
538  /// NOT initialise any entries and must be called after the entries
539  /// have been initialised to zero
540  /// The default implementation is to use finite differences to calculate
541  /// the derivatives.
543  double* const &parameter_pt,
544  Vector<double> &dres_dparam);
545 
546 
547  /// \short Add the elemental contribution to the derivatives of
548  /// the elemental Jacobian matrix
549  /// and residuals with respect to a parameter. This function should
550  /// NOT initialise any entries and must be called after the entries
551  /// have been initialised to zero
552  /// The default implementation is to use finite differences to calculate
553  /// the derivatives.
555  double* const &parameter_pt,
556  Vector<double> &dres_dparam,
557  DenseMatrix<double> &djac_dparam);
558 
559  /// \short Add the elemental contribution to the derivative of the
560  /// jacobian matrix,
561  /// mass matrix and the residuals vector with respect to the
562  /// passed parameter. Note that
563  /// this function should NOT initialise any entries.
564  /// It must be called after the residuals vector and
565  /// matrices have been initialised to zero.
567  double* const &parameter_pt,
568  Vector<double> &dres_dparam,
569  DenseMatrix<double> &djac_dparam,
570  DenseMatrix<double> &dmass_matrix_dparam);
571 
572 
573  /// \short Fill in contribution to the product of the Hessian
574  /// (derivative of Jacobian with
575  /// respect to all variables) an eigenvector, Y, and
576  /// other specified vectors, C
577  /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}
579  Vector<double> const &Y,
580  DenseMatrix<double> const &C,
581  DenseMatrix<double> &product);
582 
583  /// \short Fill in the contribution to the inner products between given
584  /// pairs of history values
586  Vector<std::pair<unsigned,unsigned> > const &history_index,
587  Vector<double> &inner_product);
588 
589  /// \short Fill in the contributions to the vectors that when taken
590  /// as dot product with other history values give the inner product
591  /// over the element
593  Vector<unsigned> const &history_index,
594  Vector<Vector<double> > &inner_product_vector);
595 
596 public:
597 
598  /// \short Constructor: Initialise all pointers and all values to zero
599  GeneralisedElement() :
600  Eqn_number(0), Dof_pt(0), Data_pt(0),
602 #ifdef OOMPH_HAS_MPI
604 #endif
605  {}
606 
607  /// Virtual destructor to clean up any memory allocated by the object.
608  virtual ~GeneralisedElement();
609 
610  /// Broken copy constructor
612  {
613  BrokenCopy::broken_copy("GeneralisedElement");
614  }
615 
616  /// Broken assignment operator
618  {
619  BrokenCopy::broken_assign("GeneralisedElement");
620  }
621 
622  /// Return a pointer to i-th internal data object.
623  Data*& internal_data_pt(const unsigned &i)
624  {
625 #ifdef RANGE_CHECKING
626  if(i >= Ninternal_data)
627  {
628  std::ostringstream error_message;
629  error_message << "Range Error: Internal data " << i
630  << " is not in the range (0,"
631  << Ninternal_data - 1 << ")";
632  throw OomphLibError(error_message.str(),
633  OOMPH_CURRENT_FUNCTION,
634  OOMPH_EXCEPTION_LOCATION);
635 
636  }
637 #endif
638  return Data_pt[i];
639  }
640 
641  /// Return a pointer to i-th internal data object (const version)
642  Data* const &internal_data_pt(const unsigned &i) const
643  {
644 #ifdef RANGE_CHECKING
645  if(i >= Ninternal_data)
646  {
647  std::ostringstream error_message;
648  error_message << "Range Error: Internal data " << i
649  << " is not in the range (0,"
650  << Ninternal_data - 1 << ")";
651  throw OomphLibError(error_message.str(),
652  OOMPH_CURRENT_FUNCTION,
653  OOMPH_EXCEPTION_LOCATION);
654 
655  }
656 #endif
657  return Data_pt[i];
658  }
659 
660 
661  /// Return a pointer to i-th external data object.
662  Data*& external_data_pt(const unsigned &i)
663  {
664 #ifdef RANGE_CHECKING
665  if(i >= Nexternal_data)
666  {
667  std::ostringstream error_message;
668  error_message << "Range Error: External data " << i
669  << " is not in the range (0,"
670  << Nexternal_data - 1 << ")";
671  throw OomphLibError(error_message.str(),
672  OOMPH_CURRENT_FUNCTION,
673  OOMPH_EXCEPTION_LOCATION);
674 
675  }
676 #endif
677  return Data_pt[Ninternal_data + i];
678  }
679 
680  /// Return a pointer to i-th external data object (const version)
681  Data* const &external_data_pt(const unsigned &i) const
682  {
683 #ifdef RANGE_CHECKING
684  if(i >= Nexternal_data)
685  {
686  std::ostringstream error_message;
687  error_message << "Range Error: External data " << i
688  << " is not in the range (0,"
689  << Nexternal_data - 1 << ")";
690  throw OomphLibError(error_message.str(),
691  OOMPH_CURRENT_FUNCTION,
692  OOMPH_EXCEPTION_LOCATION);
693 
694  }
695 #endif
696  return Data_pt[Ninternal_data + i];
697  }
698 
699  /// \short Static boolean to suppress warnings about repeated internal
700  /// data. Defaults to false.
702 
703  /// \short Static boolean to suppress warnings about repeated external
704  /// data. Defaults to true.
706 
707  /// \short Return the global equation number corresponding to the
708  /// ieqn_local-th local equation number
709  unsigned long eqn_number(const unsigned &ieqn_local) const
710  {
711 #ifdef RANGE_CHECKING
712  if(ieqn_local >= Ndof)
713  {
714  std::ostringstream error_message;
715  error_message << "Range Error: Equation number " << ieqn_local
716  << " is not in the range (0,"
717  << Ndof - 1 << ")";
718  throw OomphLibError(error_message.str(),
719  OOMPH_CURRENT_FUNCTION,
720  OOMPH_EXCEPTION_LOCATION);
721  }
722 #endif
723  return Eqn_number[ieqn_local];
724  }
725 
726 
727  ///\short Return the local equation number corresponding to the ieqn_global-th
728  ///global equation number. Returns minus one (-1) if there is no
729  ///local degree of freedom corresponding to the chosen global equation
730  ///number
731  int local_eqn_number(const unsigned long &ieqn_global) const
732  {
733  //Cache the number of degrees of freedom in the element
734  const unsigned n_dof = this->Ndof;
735  //Loop over the local equation numbers
736  for(unsigned n=0;n<n_dof;n++)
737  {
738  //If the global equation numbers match return
739  if(ieqn_global==Eqn_number[n]) {return n;}
740  }
741 
742  //If we've got all the way to the end the number has not been found
743  //return minus one.
744  return -1;
745  }
746 
747 
748  /// Add a (pointer to an) external data object to the element and return its
749  /// index (i.e. the index required to obtain it from
750  /// the access function \c external_data_pt(...). The optional boolean
751  /// flag indicates whether the data should be included in the
752  /// general finite-difference loop when calculating the jacobian.
753  /// The default value is true, i.e. the data will be included in the
754  /// finite-differencing
755  unsigned add_external_data(Data* const &data_pt, const bool &fd=true);
756 
757  /// \short Return the status of the boolean flag indicating whether
758  /// the external data is included in the finite difference loop
759  inline bool external_data_fd(const unsigned &i) const
760  {
761 #ifdef RANGE_CHECKING
762  if(i >= Nexternal_data)
763  {
764  std::ostringstream error_message;
765  error_message << "Range Error: Internal data " << i
766  << " is not in the range (0,"
767  << Nexternal_data-1 << ")";
768  throw OomphLibError(error_message.str(),
769  OOMPH_CURRENT_FUNCTION,
770  OOMPH_EXCEPTION_LOCATION);
771  }
772 #endif
773  //Return the value
774  return Data_fd[Ninternal_data + i];
775  }
776 
777 
778 
779  /// \short Set the boolean flag to exclude the external datum from the
780  /// the finite difference loop when computing the jacobian matrix
781  inline void exclude_external_data_fd(const unsigned &i)
782  {
783 #ifdef RANGE_CHECKING
784  if(i >= Nexternal_data)
785  {
786  std::ostringstream error_message;
787  error_message << "Range Error: External data " << i
788  << " is not in the range (0,"
789  << Nexternal_data - 1 << ")";
790  throw OomphLibError(error_message.str(),
791  OOMPH_CURRENT_FUNCTION,
792  OOMPH_EXCEPTION_LOCATION);
793 
794  }
795 #endif
796  //Set the value
797  Data_fd[Ninternal_data + i] = false;
798  }
799 
800  /// \short Set the boolean flag to include the external datum in the
801  /// the finite difference loop when computing the jacobian matrix
802  inline void include_external_data_fd(const unsigned &i)
803  {
804 #ifdef RANGE_CHECKING
805  if(i >= Nexternal_data)
806  {
807  std::ostringstream error_message;
808  error_message << "Range Error: External data " << i
809  << " is not in the range (0,"
810  << Nexternal_data - 1 << ")";
811  throw OomphLibError(error_message.str(),
812  OOMPH_CURRENT_FUNCTION,
813  OOMPH_EXCEPTION_LOCATION);
814 
815  }
816 #endif
817  //Set the value
818  Data_fd[Ninternal_data + i] = true;
819  }
820 
821  /// Flush all external data
822  void flush_external_data();
823 
824  /// Flush the object addressed by data_pt from the external data array
825  void flush_external_data(Data* const &data_pt);
826 
827  /// Return the number of internal data objects.
828  unsigned ninternal_data() const {return Ninternal_data;}
829 
830  /// Return the number of external data objects.
831  unsigned nexternal_data() const {return Nexternal_data;}
832 
833  /// Return the number of equations/dofs in the element.
834  unsigned ndof() const {return Ndof;}
835 
836  ///\short Return the vector of dof values at time level t
837  void dof_vector(const unsigned &t, Vector<double> &dof)
838  {
839  //Check that the internal storage has been set up
840 #ifdef PARANOID
841  if(Dof_pt==0)
842  {
843  std::stringstream error_stream;
844  error_stream << "Internal dof array not set up in element.\n"
845  << "In order to set it up you must call\n"
846  << " Problem::enable_store_local_dof_in_elements()\n"
847  << "before the call to Problem::assign_eqn_numbers()\n";
848  throw OomphLibError(error_stream.str(),
849  OOMPH_CURRENT_FUNCTION,
850  OOMPH_EXCEPTION_LOCATION);
851  }
852 #endif
853  //Resize the vector
854  dof.resize(this->Ndof);
855  //Loop over the vector and fill in the desired values
856  for(unsigned i=0;i<this->Ndof;++i)
857  {
858  dof[i] = Dof_pt[i][t];
859  }
860  }
861 
862  ///\short Return the vector of pointers to dof values
864  {
865  //Check that the internal storage has been set up
866 #ifdef PARANOID
867  if(Dof_pt==0)
868  {
869  std::stringstream error_stream;
870  error_stream << "Internal dof array not set up in element.\n"
871  << "In order to set it up you must call\n"
872  << " Problem::enable_store_local_dof_in_elements()\n"
873  << "before the call to Problem::assign_eqn_numbers()\n";
874  throw OomphLibError(error_stream.str(),
875  OOMPH_CURRENT_FUNCTION,
876  OOMPH_EXCEPTION_LOCATION);
877  }
878 #endif
879  //Resize the vector
880  dof_pt.resize(this->Ndof);
881  //Loop over the vector and fill in the desired values
882  for(unsigned i=0;i<this->Ndof;++i)
883  {
884  dof_pt[i] = Dof_pt[i];
885  }
886  }
887 
888 
889 
890  /// \short Set the timestepper associated with the i-th internal data
891  /// object
892  void set_internal_data_time_stepper(const unsigned &i,
893  TimeStepper* const &time_stepper_pt,
894  const bool &preserve_existing_data)
895  {this->internal_data_pt(i)->set_time_stepper(time_stepper_pt,
896  preserve_existing_data);}
897 
898  /// \short Assign the global equation numbers to the internal Data.
899  /// The arguments are the current highest global equation number
900  /// (which will be incremented) and a Vector of pointers to the
901  /// global variables (to which any unpinned values in the internal Data are
902  /// added).
903  void assign_internal_eqn_numbers(unsigned long &global_number,
905 
906  /// \short Function to describe the dofs of the element. The ostream
907  /// specifies the output stream to which the description
908  /// is written; the string stores the currently
909  /// assembled output that is ultimately written to the
910  /// output stream by Data::describe_dofs(...); it is typically
911  /// built up incrementally as we descend through the
912  /// call hierarchy of this function when called from
913  /// Problem::describe_dofs(...)
914  void describe_dofs(std::ostream& out,const std::string& current_string) const;
915 
916  /// \short Function to describe the local dofs of the element. The ostream
917  /// specifies the output stream to which the description
918  /// is written; the string stores the currently
919  /// assembled output that is ultimately written to the
920  /// output stream by Data::describe_dofs(...); it is typically
921  /// built up incrementally as we descend through the
922  /// call hierarchy of this function when called from
923  /// Problem::describe_dofs(...)
924  virtual void describe_local_dofs(std::ostream& out,
925  const std::string& current_string) const;
926 
927  /// \short Add pointers to the internal data values to map indexed
928  /// by the global equation number.
929  void add_internal_value_pt_to_map(std::map<unsigned,double*> &
930  map_of_value_pt);
931 
932 #ifdef OOMPH_HAS_MPI
933  /// \short Add all internal data and time history values to the vector in
934  /// the internal storage order
935  void add_internal_data_values_to_vector(Vector<double> &vector_of_values);
936 
937  /// \short Read all internal data and time history values from the vector
938  /// starting from index. On return the index will be
939  /// set to the value at the end of the data that has been read in
941  const Vector<double> & vector_of_values, unsigned &index);
942 
943  /// \short Add all equation numbers associated with internal data
944  /// to the vector in the internal storage order
945  void add_internal_eqn_numbers_to_vector(Vector<long> &vector_of_eqn_numbers);
946 
947  /// \short Read all equation numbers associated with internal data
948  /// from the vector
949  /// starting from index. On return the index will be
950  /// set to the value at the end of the data that has been read in
952  const Vector<long> & vector_of_eqn_numbers, unsigned &index);
953 
954 #endif
955 
956  /// \short Setup the arrays of local equation numbers for the element.
957  /// If the optional boolean argument is true, then pointers to the associated
958  /// degrees of freedom are stored locally in the array Dof_pt
959  virtual void assign_local_eqn_numbers(const bool &store_local_dof_pt);
960 
961  /// \short Complete the setup of any additional dependencies
962  /// that the element may have. Empty virtual function that may be
963  /// overloaded for specific derived elements. Used, e.g., for elements
964  /// with algebraic node update functions to determine the "geometric
965  /// Data", i.e. the Data that affects the element's shape.
966  /// This function is called (for all elements) at the very beginning of the
967  /// equation numbering procedure to ensure that all dependencies
968  /// are accounted for.
970 
971  /// \short Calculate the vector of residuals of the equations in the element.
972  /// By default initialise the vector to zero and then call the
973  /// fill_in_contribution_to_residuals() function. Note that this entire
974  /// function can be overloaded if desired.
975  virtual void get_residuals(Vector<double> &residuals)
976  {
977  //Zero the residuals vector
978  residuals.initialise(0.0);
979  //Add the elemental contribution to the residuals vector
981  }
982 
983  /// \short Calculate the elemental Jacobian matrix "d equation / d variable".
984  virtual void get_jacobian(Vector<double> &residuals,
985  DenseMatrix<double> &jacobian)
986  {
987  //Zero the residuals vector
988  residuals.initialise(0.0);
989  //Zero the jacobian matrix
990  jacobian.initialise(0.0);
991  //Add the elemental contribution to the residuals vector
992  fill_in_contribution_to_jacobian(residuals,jacobian);
993  }
994 
995  /// \short Calculate the residuals and the elemental "mass" matrix, the
996  /// matrix that multiplies the time derivative terms in a problem.
997  virtual void get_mass_matrix(Vector<double> &residuals,
998  DenseMatrix<double> &mass_matrix)
999  {
1000  //Zero the residuals vector
1001  residuals.initialise(0.0);
1002  //Zero the mass matrix
1003  mass_matrix.initialise(0.0);
1004  //Add the elemental contribution to the vector and matrix
1005  fill_in_contribution_to_mass_matrix(residuals,mass_matrix);
1006  }
1007 
1008  /// \short Calculate the residuals and jacobian and elemental "mass" matrix,
1009  /// the matrix that multiplies the time derivative terms.
1011  DenseMatrix<double> &jacobian,
1012  DenseMatrix<double> &mass_matrix)
1013  {
1014  //Zero the residuals vector
1015  residuals.initialise(0.0);
1016  //Zero the jacobian matrix
1017  jacobian.initialise(0.0);
1018  //Zero the mass matrix
1019  mass_matrix.initialise(0.0);
1020  //Add the elemental contribution to the vector and matrix
1022  mass_matrix);
1023  }
1024 
1025 
1026  /// \short Calculate the derivatives of the residuals with respect to
1027  /// a parameter
1028  virtual void get_dresiduals_dparameter(double* const &parameter_pt,
1029  Vector<double> &dres_dparam)
1030  {
1031  //Zero the dres_dparam vector
1032  dres_dparam.initialise(0.0);
1033  //Add the elemental contribution to the residuals vector
1035  dres_dparam);
1036  }
1037 
1038  /// \short Calculate the derivatives of the elemental Jacobian matrix
1039  /// and residuals with respect to a parameter
1040  virtual void get_djacobian_dparameter(double* const &parameter_pt,
1041  Vector<double> &dres_dparam,
1042  DenseMatrix<double> &djac_dparam)
1043  {
1044  //Zero the residuals vector
1045  dres_dparam.initialise(0.0);
1046  //Zero the jacobian matrix
1047  djac_dparam.initialise(0.0);
1048  //Add the elemental contribution to the residuals vector
1049  this->fill_in_contribution_to_djacobian_dparameter(parameter_pt,dres_dparam,
1050  djac_dparam);
1051  }
1052 
1053  /// \short Calculate the derivatives of the elemental Jacobian matrix
1054  /// mass matrix and residuals with respect to a parameter
1056  double* const &parameter_pt,
1057  Vector<double> &dres_dparam,
1058  DenseMatrix<double> &djac_dparam,
1059  DenseMatrix<double> &dmass_matrix_dparam)
1060  {
1061  //Zero the residuals derivative vector
1062  dres_dparam.initialise(0.0);
1063  //Zero the jacobian derivative matrix
1064  djac_dparam.initialise(0.0);
1065  //Zero the mass matrix derivative
1066  dmass_matrix_dparam.initialise(0.0);
1067  //Add the elemental contribution to the residuals vector and matrices
1069  parameter_pt,dres_dparam,djac_dparam,dmass_matrix_dparam);
1070  }
1071 
1072 
1073  /// \short Calculate the product of the Hessian (derivative of Jacobian with
1074  /// respect to all variables) an eigenvector, Y, and
1075  /// other specified vectors, C
1076  /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}
1078  DenseMatrix<double> const &C,
1079  DenseMatrix<double> &product)
1080  {
1081  //Initialise the product to zero
1082  product.initialise(0.0);
1083  ///Add the elemental contribution to the product
1085  }
1086 
1087  /// \short Return the vector of inner product of the given pairs of
1088  /// history values
1089  virtual void get_inner_products(Vector<std::pair<unsigned,unsigned> >
1090  const &history_index,
1091  Vector<double> &inner_product)
1092  {
1093  inner_product.initialise(0.0);
1094  this->fill_in_contribution_to_inner_products(history_index,
1095  inner_product);
1096  }
1097 
1098  /// \short Compute the vectors that when taken as a dot product with
1099  /// other history values give the inner product over the element.
1100  virtual void get_inner_product_vectors(Vector<unsigned> const &history_index,
1101  Vector<Vector<double> > &inner_product_vector)
1102  {
1103  const unsigned n_inner_product = inner_product_vector.size();
1104  for(unsigned i=0;i<n_inner_product;++i)
1105  {inner_product_vector[i].initialise(0.0);}
1107  inner_product_vector);
1108  }
1109 
1110 
1111  /// \short Self-test: Have all internal values been classified as
1112  /// pinned/unpinned? Return 0 if OK.
1113  virtual unsigned self_test();
1114 
1115 
1116  /// \short Compute norm of solution -- broken virtual can be overloaded
1117  /// by element writer to implement whatever norm is desired for
1118  /// the specific element
1119  virtual void compute_norm(double& norm)
1120  {
1121  std::string error_message = "compute_norm undefined for this element \n";
1122  throw OomphLibError(error_message,
1123  OOMPH_CURRENT_FUNCTION,
1124  OOMPH_EXCEPTION_LOCATION);
1125 
1126  }
1127 
1128 #ifdef OOMPH_HAS_MPI
1129 
1130  /// \short Label the element as halo and specify processor that holds
1131  /// non-halo counterpart
1132  void set_halo(const unsigned& non_halo_proc_ID)
1133  {
1135  }
1136 
1137  /// \short Label the element as not being a halo
1139 
1140  /// \short Is this element a halo?
1141  bool is_halo() const {return (Non_halo_proc_ID!=-1);}
1142 
1143  /// \short ID of processor ID that holds non-halo counterpart
1144  /// of halo element; negative if not a halo.
1146  {
1147  return Non_halo_proc_ID;
1148  }
1149 
1150  /// Insist that this element be kept as a halo element during a distribute?
1152 
1153  /// \short Do not insist that this element be kept as a halo element during
1154  /// distribution
1156 
1157  /// \short Test whether the element must be kept as a halo element
1159 
1160 #endif
1161 
1162  /// \short Double used for the default finite difference step in elemental
1163  /// jacobian calculations
1165 
1166  /// \short The number of types of degrees of freedom in this element
1167  /// are sub-divided into
1168  virtual unsigned ndof_types() const
1169  {
1170  // error message stream
1171  std::ostringstream error_message;
1172  error_message << "ndof_types() const has not been implemented for this \n"
1173  << "element\n" << std::endl;
1174  // throw error
1175  throw OomphLibError(error_message.str(),
1176  OOMPH_CURRENT_FUNCTION,
1177  OOMPH_EXCEPTION_LOCATION);
1178  }
1179 
1180  /// \short Create a list of pairs for the unknowns that this element
1181  /// is "in charge of" -- ignore any unknowns associated with external
1182  /// \c Data. The first entry in each pair must contain the global equation
1183  /// number of the unknown, while the second one contains the number
1184  /// of the DOF type that this unknown is associated with.
1185  /// (The function can obviously only be called if the equation numbering
1186  /// scheme has been set up.)
1188  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
1189  {
1190  // error message stream
1191  std::ostringstream error_message;
1192  error_message << "get_dof_numbers_for_unknowns() const has not been \n"
1193  << " implemented for this element\n" << std::endl;
1194  // throw error
1195  throw OomphLibError(error_message.str(),
1196  OOMPH_CURRENT_FUNCTION,
1197  OOMPH_EXCEPTION_LOCATION);
1198  }
1199 
1200 };
1201 
1202  /// Enumeration a finite element's geometry "type". Either "Q" (square,
1203  /// cubeoid like) or "T" (triangle, tetrahedron).
1205  {
1207  }
1208 
1209 //Forward class definitions that are used in FiniteElements
1210 class FaceElement;
1211 class MacroElement;
1212 class TimeStepper;
1213 class Shape;
1214 class DShape;
1215 class Integral;
1216 
1217 //===========================================================================
1218 /// \short Helper namespace for tolerances and iterations within the Newton
1219 /// method used in the locate_zeta function in FiniteElements.
1220 //===========================================================================
1221 namespace Locate_zeta_helpers
1222 {
1223  /// Convergence tolerance for the Newton solver
1224  extern double Newton_tolerance;
1225 
1226  /// Maximum number of Newton iterations
1227  extern unsigned Max_newton_iterations;
1228 
1229  /// Rounding tolerance for whether coordinate is in element or not
1230  extern double Rounding_tolerance;
1231 
1232  /// Number of points along one dimension of each element used
1233  /// to create coordinates within the element in order to see
1234  /// which has the smallest initial residual (and is therefore used
1235  /// as the initial guess in the Newton method for locate_zeta)
1236  extern unsigned N_local_points;
1237 }
1238 
1239 
1240  /// \short Typedef for the function that translates the face coordinate
1241  /// to the coordinate in the bulk element
1242  typedef void (*CoordinateMappingFctPt)(const Vector<double> &s,
1243  Vector<double> &s_bulk);
1244 
1245  /// \short Typedef for the function that returns the partial derivative
1246  /// of the local coordinates in the bulk element
1247  /// with respect to the coordinates along the face.
1248  /// In addition this function returns an index of one of the
1249  /// bulk local coordinates that varies away from the edge
1251  DenseMatrix<double> &ds_bulk_dsface,
1252  unsigned &interior_direction);
1253 
1254 
1255 //========================================================================
1256 /// \short A general Finite Element class.
1257 ///
1258 /// The main components of a FiniteElement are:
1259 /// - pointers to its Nodes
1260 /// - pointers to its internal Data (inherited from GeneralisedElement)
1261 /// - pointers to its external Data (inherited from GeneralisedElement)
1262 /// - a pointer to a spatial integration scheme
1263 /// - a pointer to the global Time object (inherited from GeneralisedElement)
1264 /// - a lookup table which establishes the relation between local
1265 /// and global equation numbers (inherited from GeneralisedElement)
1266 ///
1267 /// We also provide interfaces for functions that compute the
1268 /// element's Jacobian matrix and/or the Vector of residuals
1269 /// (inherited from GeneralisedElement) plus various output routines.
1270 //========================================================================
1271 class FiniteElement : public virtual GeneralisedElement, public GeomObject
1272 {
1273  private:
1274 
1275  /// Pointer to the spatial integration scheme
1277 
1278  /// Storage for pointers to the nodes in the element
1280 
1281  /// \short Storage for the local equation numbers associated with
1282  /// the values stored at the nodes
1284 
1285  /// Number of nodes in the element
1286  unsigned Nnode;
1287 
1288  /// \short The spatial dimension of the element, i.e. the number
1289  /// of local coordinates used to parametrize it.
1291 
1292  /// \short The spatial dimension of the nodes in the element.
1293  /// We assume that nodes have the same spatial dimension, because
1294  /// we cannot think of any "real" problems for which that would not
1295  /// be the case.
1297 
1298  /// \short The number of coordinate types required to interpolate
1299  /// the element's geometry between the nodes. For Lagrange elements
1300  /// it is 1 (the default). It must be over-ridden by using
1301  /// the set_nposition_type() function in the constructors of elements
1302  /// that use generalised coordinate, e.g. for 1D Hermite elements
1303  /// Nnodal_position_types =2.
1305 
1306  protected:
1307 
1308  /// \short Assemble the jacobian matrix for the mapping from local
1309  /// to Eulerian coordinates, given the derivatives of the shape function
1310  /// w.r.t the local coordinates.
1312  const DShape &dpsids, DenseMatrix<double> &jacobian) const;
1313 
1314  /// \short Assemble the the "jacobian" matrix of second derivatives of the
1315  /// mapping from local to Eulerian coordinates, given
1316  /// the second derivatives of the shape functions w.r.t. local coordinates.
1318  const DShape &d2psids, DenseMatrix<double> &jacobian2) const;
1319 
1320  /// \short Assemble the covariant Eulerian base vectors, assuming that
1321  /// the derivatives of the shape functions with respect to the local
1322  /// coordinates have already been constructed.
1323  virtual void assemble_eulerian_base_vectors(
1324  const DShape &dpsids, DenseMatrix<double> &interpolated_G) const;
1325 
1326  /// \short Default return value for required_nvalue(n) which gives the number
1327  /// of "data" values required by the element at node n; for example, solving
1328  /// a Poisson equation would required only one "data" value at each node. The
1329  /// defaults is set to zero, because a general element is problem-less.
1330  static const unsigned Default_Initial_Nvalue;
1331 
1332  /// \short Default value for the tolerance to be used when locating nodes
1333  /// via local coordinates
1334  static const double Node_location_tolerance;
1335 
1336 public:
1337 
1338  /// \short Set the dimension of the element and initially set
1339  /// the dimension of the nodes to be the same as the dimension of the
1340  /// element.
1341  void set_dimension(const unsigned &dim)
1343 
1344  /// \short Set the dimension of the nodes in the element. This will
1345  /// typically only be required when constructing FaceElements or
1346  /// in beam and shell type elements where a lower dimensional surface
1347  /// is embedded in a higher dimensional space.
1348  void set_nodal_dimension(const unsigned &nodal_dim)
1349  {Nodal_dimension = nodal_dim;}
1350 
1351  /// \short Set the number of types required to interpolate the coordinate
1352  void set_nnodal_position_type(const unsigned &nposition_type)
1353  {Nnodal_position_type = nposition_type;}
1354 
1355 
1356  /// \short Set the number of nodes in the element to n, by resizing
1357  /// the storage for pointers to the Node objects.
1358  void set_n_node(const unsigned &n)
1359  {
1360  //This should only be done once, in a Node constructor
1361 //#ifdef PARANOID
1362  //if(Node_pt)
1363  // {
1364  // OomphLibWarning(
1365  // "You are resetting the number of nodes in an element -- Be careful",
1366  // "FiniteElement::set_n_node()",
1367  // OOMPH_EXCEPTION_LOCATION);
1368  // }
1369 //#endif
1370  //Delete any previous storage to avoid memory leaks
1371  //This will only happen in very special cases
1372  delete[] Node_pt;
1373  //Set the number of nodes
1374  Nnode = n;
1375  //Allocate the storage
1376  Node_pt = new Node*[n];
1377  //Initialise all the pointers to NULL
1378  for(unsigned i=0;i<n;i++) {Node_pt[i] = 0;}
1379  }
1380 
1381  /// \short Return the local equation number corresponding to the i-th
1382  /// value at the n-th local node.
1383  inline int nodal_local_eqn(const unsigned &n, const unsigned &i) const
1384  {
1385 #ifdef RANGE_CHECKING
1386  if(n >= Nnode)
1387  {
1388  std::ostringstream error_message;
1389  error_message << "Range Error: Node number " << n
1390  << " is not in the range (0,"
1391  << Nnode-1 << ")";
1392  throw OomphLibError(error_message.str(),
1393  OOMPH_CURRENT_FUNCTION,
1394  OOMPH_EXCEPTION_LOCATION);
1395  }
1396  else
1397  {
1398  unsigned n_value = node_pt(n)->nvalue();
1399  if(i >= n_value)
1400  {
1401  std::ostringstream error_message;
1402  error_message << "Range Error: value " << i << " at node " << n
1403  << " is not in the range (0,"
1404  << n_value -1 << ")";
1405  throw OomphLibError(error_message.str(),
1406  OOMPH_CURRENT_FUNCTION,
1407  OOMPH_EXCEPTION_LOCATION);
1408  }
1409  }
1410 #endif
1411 #ifdef PARANOID
1412  //Check that the equations have been allocated
1413  if(Nodal_local_eqn==0)
1414  {
1415  throw OomphLibError(
1416  "Nodal local equation numbers have not been allocated",
1417  OOMPH_CURRENT_FUNCTION,
1418  OOMPH_EXCEPTION_LOCATION);
1419  }
1420 #endif
1421  return Nodal_local_eqn[n][i];
1422  }
1423 
1424 
1425 
1426  /// \short Compute the geometric shape functions (psi) at integration point
1427  /// ipt. Return the determinant of the jacobian of the mapping (detJ).
1428  /// Additionally calculate the derivatives of "detJ" w.r.t. the
1429  /// nodal coordinates.
1430  double dJ_eulerian_at_knot(const unsigned &ipt,
1431  Shape &psi,
1432  DenseMatrix<double> &djacobian_dX) const;
1433  protected:
1434 
1435  /// \short Static array that holds the number of second derivatives
1436  /// as a function of the dimension of the element
1437  static const unsigned N2deriv[];
1438 
1439  /// \short Take the matrix passed as jacobian and return its inverse in
1440  /// inverse_jacobian. This function is templated by the dimension of the
1441  /// element because matrix inversion cannot be written efficiently in a
1442  /// generic manner.
1443  template<unsigned DIM>
1444  double invert_jacobian(const DenseMatrix<double> &jacobian,
1445  DenseMatrix<double> &inverse_jacobian) const;
1446 
1447 
1448  /// \short A template-free interface that takes the matrix passed as jacobian
1449  /// and return its inverse in inverse_jacobian. By default the function will
1450  /// use the dimension of the element to call the correct invert_jacobian(..)
1451  /// function. This should be overloaded for efficiency (removal of a switch
1452  /// statement) in specific elements.
1453  virtual double invert_jacobian_mapping(const DenseMatrix<double> &jacobian,
1454  DenseMatrix<double> &inverse_jacobian)
1455  const;
1456 
1457 
1458  /// \short Calculate the mapping from local to Eulerian coordinates,
1459  /// given the derivatives of the shape functions w.r.t. local coordinates.
1460  /// Returns the determinant of the jacobian, the jacobian and inverse jacobian
1461  virtual double local_to_eulerian_mapping(const DShape &dpsids,
1462  DenseMatrix<double> &jacobian,
1463  DenseMatrix<double> &inverse_jacobian) const
1464  {
1465  //Assemble the jacobian
1466  assemble_local_to_eulerian_jacobian(dpsids,jacobian);
1467  //Invert the jacobian (use the template-free interface)
1468  return invert_jacobian_mapping(jacobian,inverse_jacobian);
1469  }
1470 
1471 
1472  /// \short Calculate the mapping from local to Eulerian coordinates,
1473  /// given the derivatives of the shape functions w.r.t. local coordinates,
1474  /// Return only the determinant of the jacobian and the inverse of the
1475  /// mapping (ds/dx).
1476  double local_to_eulerian_mapping(const DShape &dpsids,
1477  DenseMatrix<double> &inverse_jacobian) const
1478  {
1479  //Find the dimension of the element
1480  const unsigned el_dim = dim();
1481  //Assign memory for the jacobian
1482  DenseMatrix<double> jacobian(el_dim);
1483  //Calculate the jacobian and inverse
1484  return local_to_eulerian_mapping(dpsids,jacobian,inverse_jacobian);
1485  }
1486 
1487  /// \short Calculate the mapping from local to Eulerian coordinates given
1488  /// the derivatives of the shape functions w.r.t the local coordinates.
1489  /// assuming that the coordinates are aligned in the direction of the local
1490  /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
1491  /// This function returns the determinant of the jacobian, the jacobian
1492  /// and the inverse jacobian.
1493  virtual double local_to_eulerian_mapping_diagonal(
1494  const DShape &dpsids,DenseMatrix<double> &jacobian,
1495  DenseMatrix<double> &inverse_jacobian) const;
1496 
1497  /// \short A template-free interface that calculates the derivative of the
1498  /// jacobian of a mapping with respect to the nodal coordinates X_ij.
1499  /// To do this it requires the jacobian matrix and the derivatives of the
1500  /// shape functions w.r.t. the local coordinates. By default the function
1501  /// will use the dimension of the element to call the correct
1502  /// dJ_eulerian_dnodal_coordinates_templated_helper(..) function. This
1503  /// should be overloaded for efficiency (removal of a switch statement)
1504  /// in specific elements.
1505  virtual void dJ_eulerian_dnodal_coordinates(
1506  const DenseMatrix<double> &jacobian,const DShape &dpsids,
1507  DenseMatrix<double> &djacobian_dX) const;
1508 
1509  /// \short Calculate the derivative of the jacobian of a mapping with
1510  /// respect to the nodal coordinates X_ij using the jacobian matrix
1511  /// and the derivatives of the shape functions w.r.t. the local
1512  /// coordinates. This function is templated by the dimension of the
1513  /// element.
1514  template<unsigned DIM>
1516  const DenseMatrix<double> &jacobian,const DShape &dpsids,
1517  DenseMatrix<double> &djacobian_dX) const;
1518 
1519  /// \short A template-free interface that calculates the derivative w.r.t.
1520  /// the nodal coordinates \f$ X_{pq} \f$ of the derivative of the shape
1521  /// functions \f$ \psi_j \f$ w.r.t. the global eulerian coordinates
1522  /// \f$ x_i \f$. I.e. this function calculates
1523  /// \f[
1524  /// \frac{\partial}{\partial X_{pq}}
1525  /// \left( \frac{\partial \psi_j}{\partial x_i} \right).
1526  /// \f]
1527  /// To do this it requires the determinant of the jacobian mapping, its
1528  /// derivative w.r.t. the nodal coordinates \f$ X_{pq} \f$, the inverse
1529  /// jacobian and the derivatives of the shape functions w.r.t. the local
1530  /// coordinates. The result is returned as a tensor of rank four.
1531  /// Numbering:
1532  /// d_dpsidx_dX(p,q,j,i) = \f$ \frac{\partial}{\partial X_{pq}}
1533  /// \left( \frac{\partial \psi_j}{\partial x_i} \right) \f$
1534  /// By default the function will use the dimension of the element to call
1535  /// the correct d_dshape_eulerian_dnodal_coordinates_templated_helper(..)
1536  /// function. This should be overloaded for efficiency (removal of a
1537  /// switch statement) in specific elements.
1539  const double &det_jacobian,
1540  const DenseMatrix<double> &jacobian,
1541  const DenseMatrix<double> &djacobian_dX,
1542  const DenseMatrix<double> &inverse_jacobian,
1543  const DShape &dpsids,
1544  RankFourTensor<double> &d_dpsidx_dX) const;
1545 
1546  /// \short Calculate the derivative w.r.t. the nodal coordinates
1547  /// \f$ X_{pq} \f$ of the derivative of the shape functions w.r.t. the
1548  /// global eulerian coordinates \f$ x_i \f$, using the determinant of
1549  /// the jacobian mapping, its derivative w.r.t. the nodal coordinates
1550  /// \f$ X_{pq} \f$, the inverse jacobian and the derivatives of the
1551  /// shape functions w.r.t. the local coordinates. The result is returned
1552  /// as a tensor of rank four.
1553  /// Numbering:
1554  /// d_dpsidx_dX(p,q,j,i) = \f$ \frac{\partial}{\partial X_{pq}}
1555  /// \left( \frac{\partial \psi_j}{\partial x_i} \right) \f$
1556  /// This function is templated by the dimension of the element.
1557  template<unsigned DIM>
1559  const double &det_jacobian,
1560  const DenseMatrix<double> &jacobian,
1561  const DenseMatrix<double> &djacobian_dX,
1562  const DenseMatrix<double> &inverse_jacobian,
1563  const DShape &dpsids,
1564  RankFourTensor<double> &d_dpsidx_dX) const;
1565 
1566  /// \short Convert derivative w.r.t.local coordinates to derivatives
1567  /// w.r.t the coordinates used to assemble the inverse_jacobian passed
1568  /// in the mapping. On entry, dbasis must contain the basis function
1569  /// derivatives w.r.t. the local coordinates; it will contain the
1570  /// derivatives w.r.t. the new coordinates on exit.
1571  /// This is virtual so that it may be overloaded if desired
1572  /// for efficiency reasons.
1573  virtual void transform_derivatives(const DenseMatrix<double>
1574  &inverse_jacobian, DShape &dbasis) const;
1575 
1576  /// \short Convert derivative w.r.t local coordinates to derivatives
1577  /// w.r.t the coordinates used to assemble the inverse jacobian passed
1578  /// in the mapping, assuming that the coordinates are aligned in the
1579  /// direction of the local coordinates. On entry dbasis must contain
1580  /// the derivatives of the basis functions w.r.t. the local coordinates;
1581  /// it will contain the derivatives w.r.t. the new coordinates.
1582  /// are converted into the new using the mapping inverse_jacobian.
1584  &inverse_jacobian, DShape &dbasis) const;
1585 
1586  /// \short Convert derivatives and second derivatives w.r.t. local coordiantes
1587  /// to derivatives and second derivatives w.r.t. the coordinates used to
1588  /// assemble the jacobian, inverse jacobian and jacobian2 passed to the
1589  /// function. By default this function will call
1590  /// transform_second_derivatives_template<>(...) using the dimension of
1591  /// the element as the template parameter. It is virtual so that it can
1592  /// be overloaded by a specific element to save using a switch statement.
1593  /// Optionally, the element writer may wish to use the
1594  /// transform_second_derivatives_diagonal<>(...) function
1595  /// On entry dbasis and d2basis must contain the derivatives w.r.t. the
1596  /// local coordinates; on exit they will be the derivatives w.r.t. the
1597  /// transformed coordinates.
1599  &jacobian,
1600  const DenseMatrix<double>
1601  &inverse_jacobian,
1602  const DenseMatrix<double>
1603  &jacobian2,
1604  DShape &dbasis,
1605  DShape &d2basis) const;
1606 
1607  /// \short Convert derivatives and second derivatives w.r.t. local coordinates
1608  /// to derivatives and second derivatives w.r.t. the coordinates used to
1609  /// asssmble the jacobian, inverse jacobian and jacobian2 passed in the
1610  /// mapping. This is templated by dimension because the method of calculation
1611  /// varies significantly with the dimension.
1612  /// On entry dbasis and d2basis must contain the derivatives w.r.t. the
1613  /// local coordinates; on exit they will be the derivatives w.r.t. the
1614  /// transformed coordinates.
1615  template<unsigned DIM>
1617  &jacobian,
1618  const DenseMatrix<double>
1619  &inverse_jacobian,
1620  const DenseMatrix<double>
1621  &jacobian2,
1622  DShape &dbasis,
1623  DShape &d2basis) const;
1624 
1625  /// \short Convert derivatives and second derivatives w.r.t. local coordinates
1626  /// to derivatives and second derivatives w.r.t. the coordinates used to
1627  /// asssmble the jacobian, inverse jacobian and jacobian2 passed in the
1628  /// mapping. This version of the function assumes that the local coordinates
1629  /// are aligned with the global coordinates, i.e. the jacobians are diagonal
1630  /// On entry dbasis and d2basis must contain the derivatives w.r.t. the
1631  /// local coordinates; on exit they will be the derivatives w.r.t. the
1632  /// transformed coordinates.
1633  template<unsigned DIM>
1635  &jacobian,
1636  const DenseMatrix<double>
1637  &inverse_jacobian,
1638  const DenseMatrix<double>
1639  &jacobian2,
1640  DShape &dbasis,
1641  DShape &d2basis) const;
1642 
1643  /// Pointer to the element's macro element (NULL by default)
1645 
1646  /// \short Calculate the contributions to the jacobian from the nodal
1647  /// degrees of freedom using finite differences.
1648  /// This version of the function assumes that the residuals vector has
1649  /// already been calculated
1650  virtual void fill_in_jacobian_from_nodal_by_fd(Vector<double> &residuals,
1651  DenseMatrix<double> &jacobian);
1652 
1653  /// \short Calculate the contributions to the jacobian from the nodal
1654  /// degrees of freedom using finite differences. This version computes
1655  /// the residuals vector before calculating the jacobian terms.
1657  {
1658  //Allocate storage for a residuals vector and initialise to zero
1659  unsigned n_dof = ndof();
1660  Vector<double> residuals(n_dof,0.0);
1661  //Get the residuals for the entire element
1662  get_residuals(residuals);
1663  //Call the jacobian calculation
1664  fill_in_jacobian_from_nodal_by_fd(residuals,jacobian);
1665  }
1666 
1667  /// \short Function that is called before the finite differencing of
1668  /// any nodal data. This may be overloaded to update any slaved
1669  /// data before finite differencing takes place.
1670  virtual inline void update_before_nodal_fd() { }
1671 
1672  /// \short Function that is call after the finite differencing of
1673  /// the nodal data. This may be overloaded to reset any slaved
1674  /// variables that may have changed during the finite differencing.
1675  virtual inline void reset_after_nodal_fd() { }
1676 
1677  /// \short Function called within the finite difference loop for
1678  /// nodal data after a change in the i-th nodal value.
1679  virtual inline void update_in_nodal_fd(const unsigned &i) { }
1680 
1681  /// \short Function called within the finite difference loop for
1682  /// nodal data after the i-th nodal values is reset.
1683  /// The default behaviour is to call the update function.
1684  virtual inline void reset_in_nodal_fd(const unsigned &i)
1685  {update_in_nodal_fd(i);}
1686 
1687 
1688  /// \short Add the elemental contribution to the jacobian matrix.
1689  /// and the residuals vector. Note that
1690  /// this function will NOT initialise the residuals vector or the jacobian
1691  /// matrix. It must be called after the residuals vector and
1692  /// jacobian matrix have been initialised to zero. The default
1693  /// is to use finite differences to calculate the jacobian
1695  DenseMatrix<double> &jacobian)
1696  {
1697  //Add the contribution to the residuals
1699  //Allocate storage for the full residuals (residuals of entire element)
1700  unsigned n_dof = ndof();
1701  Vector<double> full_residuals(n_dof);
1702  //Get the residuals for the entire element
1703  get_residuals(full_residuals);
1704  //Calculate the contributions from the internal dofs
1705  //(finite-difference the lot by default)
1706  fill_in_jacobian_from_internal_by_fd(full_residuals,jacobian,true);
1707  //Calculate the contributions from the external dofs
1708  //(finite-difference the lot by default)
1709  fill_in_jacobian_from_external_by_fd(full_residuals,jacobian,true);
1710  //Calculate the contributions from the nodal dofs
1711  fill_in_jacobian_from_nodal_by_fd(full_residuals,jacobian);
1712  }
1713 
1714 public:
1715 
1716 
1717  /// \short Function pointer for function that computes vector-valued
1718  /// steady "exact solution" \f$ {\bf f}({\bf x}) \f$
1719  /// as \f$ \mbox{\tt fct}({\bf x}, {\bf f}) \f$.
1721  Vector<double>&);
1722 
1723  /// \short Function pointer for function that computes Vector-valued
1724  /// time-dependent function \f$ {\bf f}(t,{\bf x}) \f$
1725  /// as \f$ \mbox{\tt fct}(t, {\bf x}, {\bf f}) \f$.
1726  typedef void (*UnsteadyExactSolutionFctPt)(const double&,
1727  const Vector<double>&,
1728  Vector<double>&);
1729 
1730  /// \short Tolerance below which the jacobian is considered singular
1732 
1733  /// \short Boolean that if set to true allows a negative jacobian
1734  /// in the transform between global and local coordinates (negative surface
1735  /// area = left-handed coordinate system).
1737 
1738  /// \short Static boolean to suppress output while checking
1739  /// for inverted elements
1741 
1742  /// Constructor
1744  Node_pt(0), Nodal_local_eqn(0), Nnode(0),
1746  Macro_elem_pt(0) {}
1747 
1748  /// \short The destructor cleans up the static memory allocated
1749  /// for shape function storage. Internal and external data get
1750  /// wiped by the GeneralisedElement destructor; nodes get
1751  /// killed in mesh destructor.
1752  virtual ~FiniteElement();
1753 
1754  /// Broken copy constructor
1756  {
1757  BrokenCopy::broken_copy("FiniteElement");
1758  }
1759 
1760  /// Broken assignment operator
1761 //Commented out broken assignment operator because this can lead to a conflict warning
1762 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
1763 //realise that two separate implementations of the broken function are the same and so,
1764 //quite rightly, it shouts.
1765  /*void operator=(const FiniteElement&)
1766  {
1767  BrokenCopy::broken_assign("FiniteElement");
1768  }*/
1769 
1770  ///Check whether the local coordinate are valid or not
1772  {
1773  throw OomphLibError(
1774  "local_coord_is_valid is not implemented for this element\n",
1775  OOMPH_CURRENT_FUNCTION,
1776  OOMPH_EXCEPTION_LOCATION);
1777  return true;
1778  }
1779 
1780  ///\short Check whether the local coordinate are valid or not, allowing for
1781  /// a rounding tolerance. If the point is outside the element by less
1782  /// than the tolerance, we move it back into the element.
1784  const double &rounding_tolerance)
1785  {
1786  throw OomphLibError(
1787  "local_coord_is_valid is not implemented for this element\n",
1788  OOMPH_CURRENT_FUNCTION,
1789  OOMPH_EXCEPTION_LOCATION);
1790  return true;
1791  }
1792 
1793  /// \short Get local coordinates of node j in the element; vector
1794  /// sets its own size (broken virtual)
1795  virtual void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
1796  {
1797  throw OomphLibError(
1798  "local_coordinate_of_node(...) is not implemented for this element\n",
1799  OOMPH_CURRENT_FUNCTION,
1800  OOMPH_EXCEPTION_LOCATION);
1801  }
1802 
1803  /// \short Get the local fraction of the node j in the element
1804  /// A dumb, but correct default implementation is provided.
1805  virtual void local_fraction_of_node(const unsigned &j,
1806  Vector<double> &s_fraction);
1807 
1808  /// \short Get the local fraction of any node in the n-th position
1809  /// in a one dimensional expansion along the i-th local coordinate
1810  virtual double local_one_d_fraction_of_node(const unsigned &n1d,
1811  const unsigned &i)
1812  {
1813  std::string error_message =
1814  "local one_d_fraction_of_node is not implemented for this element\n";
1815  error_message +=
1816  "It only makes sense for elements that use tensor-product expansions\n";
1817 
1818  throw OomphLibError(error_message,
1819  OOMPH_CURRENT_FUNCTION,
1820  OOMPH_EXCEPTION_LOCATION);
1821  }
1822 
1823  /// \short Set pointer to macro element -- can be overloaded in derived
1824  /// elements to perform additional tasks
1827 
1828  ///Access function to pointer to macro element
1830 
1831  /// \short Global coordinates as function of local coordinates.
1832  /// Either via FE representation or via macro-element (if Macro_elem_pt!=0)
1833  void get_x(const Vector<double> &s, Vector<double> &x) const
1834  {
1835  //If there is no macro element then return interpolated x
1836  if(Macro_elem_pt==0) {interpolated_x(s,x);}
1837  //Otherwise call the macro element representation
1838  else {get_x_from_macro_element(s,x);}
1839  }
1840 
1841 
1842  /// \short Global coordinates as function of local coordinates
1843  /// at previous time "level" t (t=0: present; t>0: previous).
1844  /// Either via FE representation of QElement or
1845  /// via macro-element (if Macro_elem_pt!=0).
1846  void get_x(const unsigned &t, const Vector<double> &s,
1847  Vector<double> &x)
1848  {
1849  // Get timestepper from first node
1851 
1852  // Number of previous values
1853  const unsigned nprev=time_stepper_pt->nprev_values();
1854 
1855  // If t > nprev_values(), we're not dealing with a previous value
1856  // but a generalised history value -- this cannot be recovered from
1857  // macro element but must be determined by finite element interpolation
1858 
1859  //If there is no macro element, or we're dealing with a generalised
1860  //history value then use the FE representation
1861  if((Macro_elem_pt==0)||(t>nprev)) {interpolated_x(t,s,x);}
1862  //Otherwise use the macro element representation
1863  else {get_x_from_macro_element(t,s,x);}
1864  }
1865 
1866 
1867  /// \short Global coordinates as function of local coordinates
1868  /// using macro element representation.
1869  /// (Broken virtual --- this must be overloaded in specific geometric
1870  /// element classes)
1872  Vector<double>& x) const
1873  {
1874  throw OomphLibError(
1875  "get_x_from_macro_element(...) is not implemented for this element\n",
1876  OOMPH_CURRENT_FUNCTION,
1877  OOMPH_EXCEPTION_LOCATION);
1878  }
1879 
1880  /// \short Global coordinates as function of local coordinates
1881  /// at previous time "level" t (t=0: present; t>0: previous).
1882  /// using macro element representation
1883  /// (Broken virtual -- overload in specific geometric element class
1884  /// if you want to use this functionality.)
1885  virtual void get_x_from_macro_element(const unsigned& t,
1886  const Vector<double>& s,
1887  Vector<double>& x)
1888  {
1889  throw OomphLibError(
1890  "get_x_from_macro_element(...) is not implemented for this element\n",
1891  OOMPH_CURRENT_FUNCTION,
1892  OOMPH_EXCEPTION_LOCATION);
1893  }
1894 
1895 
1896  ///Set the spatial integration scheme
1897  virtual void set_integration_scheme(Integral* const &integral_pt);
1898 
1899  /// Return the pointer to the integration scheme (const version)
1900  Integral* const &integral_pt() const {return Integral_pt;}
1901 
1902  /// \short Calculate the geometric shape functions
1903  /// at local coordinate s. This function must be overloaded for each specific
1904  /// geometric element.
1905  virtual void shape(const Vector<double> &s, Shape &psi) const=0;
1906 
1907  /// \short Return the geometric shape function at the ipt-th integration point
1908  virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const;
1909 
1910  /// \short Function to compute the geometric shape functions and
1911  /// derivatives w.r.t. local coordinates at local coordinate s.
1912  /// This function must be overloaded for each specific geometric element.
1913  /// (Broken virtual function --- specifies the interface)
1914  virtual void dshape_local(const Vector<double> &s, Shape &psi,
1915  DShape &dpsids) const
1916  {
1917  throw OomphLibError(
1918  "dshape_local() is not implemented for this element\n",
1919  OOMPH_CURRENT_FUNCTION,
1920  OOMPH_EXCEPTION_LOCATION);
1921  }
1922 
1923  /// \short Return the geometric shape function and its derivative w.r.t.
1924  /// the local coordinates at the ipt-th integration point.
1925  virtual void dshape_local_at_knot(const unsigned &ipt, Shape &psi,
1926  DShape &dpsids) const;
1927 
1928  /// \short Function to compute the geometric shape functions and also
1929  /// first and second derivatives w.r.t.
1930  /// local coordinates at local coordinate s. This function must
1931  /// be overloaded for each specific geometric element (if required).
1932  /// (Broken virtual function --- specifies the interface).
1933  /// Numbering:
1934  /// \b 1D:
1935  /// d2psids(i,0) = \f$ d^2 \psi_j / ds^2 \f$
1936  /// \b 2D:
1937  /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
1938  /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
1939  /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
1940  /// \b 3D:
1941  /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
1942  /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
1943  /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
1944  /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
1945  /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
1946  /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
1947  virtual void d2shape_local(const Vector<double> &s,
1948  Shape &psi,
1949  DShape &dpsids,
1950  DShape &d2psids) const
1951  {
1952  throw OomphLibError(
1953  "d2shape_local() is not implemented for this element\n",
1954  OOMPH_CURRENT_FUNCTION,
1955  OOMPH_EXCEPTION_LOCATION);
1956  }
1957 
1958  /// \short Return the geometric shape function and its first and
1959  /// second derivatives w.r.t.
1960  /// the local coordinates at the ipt-th integration point.
1961  /// Numbering:
1962  /// \b 1D:
1963  /// d2psids(i,0) = \f$ d^2 \psi_j / ds^2 \f$
1964  /// \b 2D:
1965  /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
1966  /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
1967  /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
1968  /// \b 3D:
1969  /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
1970  /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
1971  /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
1972  /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
1973  /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
1974  /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
1975  virtual void d2shape_local_at_knot(const unsigned &ipt, Shape &psi,
1976  DShape &dpsids, DShape &d2psids) const;
1977 
1978  /// \short Return the Jacobian of mapping from local to global
1979  /// coordinates at local position s.
1980  virtual double J_eulerian(const Vector<double> &s) const;
1981 
1982  /// \short Return the Jacobian of the mapping from local to global
1983  /// coordinates at the ipt-th integration point
1984  virtual double J_eulerian_at_knot(const unsigned &ipt) const;
1985 
1986  /// \short Helper function used to check for singular or negative
1987  /// Jacobians in the transform from local to global or Lagrangian
1988  /// coordinates.
1989  void check_jacobian(const double &jacobian) const;
1990 
1991  /// \short Compute the geometric shape functions and also
1992  /// first derivatives w.r.t. global coordinates at local coordinate s;
1993  /// Returns Jacobian of mapping from global to local coordinates.
1994  double dshape_eulerian(const Vector<double> &s, Shape &psi,
1995  DShape &dpsidx) const;
1996 
1997  /// \short Return the geometric shape functions and also first
1998  /// derivatives w.r.t. global coordinates at the ipt-th integration point.
1999  virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi,
2000  DShape &dpsidx) const;
2001 
2002  /// \short Compute the geometric shape functions (psi) and first derivatives
2003  /// w.r.t. global coordinates (dpsidx) at the ipt-th integration point.
2004  /// Return the determinant of the jacobian of the mapping (detJ).
2005  /// Additionally calculate the derivatives of both "detJ" and "dpsidx"
2006  /// w.r.t. the nodal coordinates.
2007  virtual double dshape_eulerian_at_knot(
2008  const unsigned &ipt,Shape &psi,DShape &dpsi,
2009  DenseMatrix<double> &djacobian_dX,RankFourTensor<double> &d_dpsidx_dX) const;
2010 
2011  /// \short Compute the geometric shape functions and also first
2012  /// and second derivatives w.r.t. global coordinates at local coordinate s;
2013  /// Returns Jacobian of mapping from global to local coordinates.
2014  /// Numbering:
2015  /// \b 1D:
2016  /// d2psidx(i,0) = \f$ d^2 \psi_j / d x^2 \f$
2017  /// \b 2D:
2018  /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$
2019  /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$
2020  /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$
2021  /// \b 3D:
2022  /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$
2023  /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$
2024  /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_2^2 \f$
2025  /// d2psidx(i,3) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$
2026  /// d2psidx(i,4) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_2 \f$
2027  /// d2psidx(i,5) = \f$ \partial^2 \psi_j / \partial x_1 \partial x_2 \f$
2028  double d2shape_eulerian(const Vector<double> &s, Shape &psi,
2029  DShape &dpsidx, DShape &d2psidx) const;
2030 
2031 
2032  /// \short Return the geometric shape functions and also first
2033  /// and second derivatives w.r.t. global coordinates at ipt-th integration
2034  /// point.
2035  /// Numbering:
2036  /// \b 1D:
2037  /// d2psidx(i,0) = \f$ d^2 \psi_j / d s^2 \f$
2038  /// \b 2D:
2039  /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$
2040  /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$
2041  /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$
2042  /// \b 3D:
2043  /// d2psidx(i,0) = \f$ \partial^2 \psi_j / \partial x_0^2 \f$
2044  /// d2psidx(i,1) = \f$ \partial^2 \psi_j / \partial x_1^2 \f$
2045  /// d2psidx(i,2) = \f$ \partial^2 \psi_j / \partial x_2^2 \f$
2046  /// d2psidx(i,3) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_1 \f$
2047  /// d2psidx(i,4) = \f$ \partial^2 \psi_j / \partial x_0 \partial x_2 \f$
2048  /// d2psidx(i,5) = \f$ \partial^2 \psi_j / \partial x_1 \partial x_2 \f$
2049  virtual double d2shape_eulerian_at_knot(const unsigned &ipt, Shape &psi,
2050  DShape &dpsidx,
2051  DShape &d2psidx) const;
2052 
2053  /// \short Assign the local equation numbers for Data stored at the nodes
2054  /// Virtual so that it can be overloaded by RefineableFiniteElements.
2055  /// If the boolean is true then the pointers to the degrees of freedom
2056  /// associated with each equation number are stored in Dof_pt
2057  virtual void assign_nodal_local_eqn_numbers(
2058  const bool &store_local_dof_pt);
2059 
2060  /// \short Function to describe the local dofs of the element[s]. The ostream
2061  /// specifies the output stream to which the description
2062  /// is written; the string stores the currently
2063  /// assembled output that is ultimately written to the
2064  /// output stream by Data::describe_dofs(...); it is typically
2065  /// built up incrementally as we descend through the
2066  /// call hierarchy of this function when called from
2067  /// Problem::describe_dofs(...)
2068  virtual void describe_local_dofs(std::ostream& out,
2069  const std::string& current_string) const;
2070 
2071  /// \short Function to describe the local dofs of the element[s]. The ostream
2072  /// specifies the output stream to which the description
2073  /// is written; the string stores the currently
2074  /// assembled output that is ultimately written to the
2075  /// output stream by Data::describe_dofs(...); it is typically
2076  /// built up incrementally as we descend through the
2077  /// call hierarchy of this function when called from
2078  /// Problem::describe_dofs(...)
2079  virtual void describe_nodal_local_dofs(
2080  std::ostream& out,
2081  const std::string& current_string) const;
2082 
2083  /// \short Overloaded version of the calculation of the local equation
2084  /// numbers. If the boolean argument is true then pointers to the degrees
2085  /// of freedom associated with each equation number are stored locally
2086  /// in the array Dof_pt.
2088  const bool &store_local_dof_pt)
2089  {
2090  //GeneralisedElement's version assigns internal and external data
2092  //Need simply to add the nodal numbering scheme
2093  assign_nodal_local_eqn_numbers(store_local_dof_pt);
2094  }
2095 
2096  /// Return a pointer to the local node n
2097  Node* &node_pt(const unsigned &n)
2098  {
2099 #ifdef RANGE_CHECKING
2100  if(n >= Nnode)
2101  {
2102  std::ostringstream error_message;
2103  error_message << "Range Error: " << n
2104  << " is not in the range (0,"
2105  << Nnode-1 << ")";
2106  throw OomphLibError(error_message.str(),
2107  OOMPH_CURRENT_FUNCTION,
2108  OOMPH_EXCEPTION_LOCATION);
2109  }
2110 #endif
2111  return Node_pt[n];
2112  }
2113 
2114  /// Return a pointer to the local node n (const version)
2115  Node* const &node_pt(const unsigned &n) const
2116  {
2117 #ifdef RANGE_CHECKING
2118  if(n >= Nnode)
2119  {
2120  std::ostringstream error_message;
2121  error_message << "Range Error: " << n
2122  << " is not in the range (0,"
2123  << Nnode-1 << ")";
2124  throw OomphLibError(error_message.str(),
2125  OOMPH_CURRENT_FUNCTION,
2126  OOMPH_EXCEPTION_LOCATION);
2127  }
2128 #endif
2129 
2130  return Node_pt[n];
2131  }
2132 
2133  /// Return the number of nodes
2134  unsigned nnode() const {return Nnode;}
2135 
2136  /// \short Return the number of nodes along one edge of the element
2137  /// Default is to return zero --- must be overloaded by geometric
2138  /// elements
2139  virtual unsigned nnode_1d() const {return 0;}
2140 
2141  /// \short Return the i-th coordinate at local node n. Do not use
2142  /// the hanging node representation.
2143  /// NOTE: Moved to cc file because of a possible compiler bug
2144  /// in gcc (yes, really!). The move to the cc file avoids inlining
2145  /// which appears to cause problems (only) when compiled with gcc
2146  /// and -O3; offensive "illegal read" is in optimised-out section
2147  /// of code and data that is allegedly illegal is readily readable
2148  /// (by other means) just before this function is called so I can't
2149  /// really see how we could possibly be responsible for this...
2150  double raw_nodal_position(const unsigned &n, const unsigned &i) const;
2151  /* { */
2152  /* /\* oomph_info << "n: "<< n << std::endl; *\/ */
2153  /* /\* oomph_info << "i: "<< i << std::endl; *\/ */
2154  /* /\* oomph_info << "node_pt(n): "<< node_pt(n) << std::endl; *\/ */
2155  /* /\* oomph_info << "node_pt(n)->x(i): "<< node_pt(n)->x(i) << std::endl; *\/ */
2156  /* double tmp=node_pt(n)->x(i); */
2157  /* //oomph_info << "tmp: "<< tmp << std::endl; */
2158  /* return tmp; // node_pt(n)->x(i); */
2159  /* } */
2160 
2161  /// \short Return the i-th coordinate at local node n, at time level t
2162  /// (t=0: present; t>0: previous time level).
2163  /// Do not use the hanging node representation.
2164  double raw_nodal_position(const unsigned &t, const unsigned &n,
2165  const unsigned &i) const
2166  {return node_pt(n)->x(t,i);}
2167 
2168  /// \short Return the i-th component of nodal velocity: dx/dt at local node n.
2169  /// Do not use the hanging node representation.
2170  double raw_dnodal_position_dt(const unsigned &n, const unsigned &i) const
2171  {return node_pt(n)->dx_dt(i);}
2172 
2173  /// \short Return the i-th component of j-th derivative of nodal position:
2174  /// d^jx/dt^j at node n. Do not use the hanging node representation.
2175  double raw_dnodal_position_dt(const unsigned &n, const unsigned &j,
2176  const unsigned &i) const
2177  {return node_pt(n)->dx_dt(j,i);}
2178 
2179  /// \short Return the value of the k-th type of the i-th positional variable
2180  /// at the local node n. Do not use the hanging node representation.
2181  double raw_nodal_position_gen(const unsigned &n, const unsigned &k,
2182  const unsigned &i) const
2183  {return node_pt(n)->x_gen(k,i);}
2184 
2185  /// \short Return the generalised nodal position (type k, i-th variable)
2186  /// at previous timesteps at local node n. Do not use the hanging node
2187  /// representation.
2188  double raw_nodal_position_gen(const unsigned &t, const unsigned &n,
2189  const unsigned &k, const unsigned &i) const
2190  {return node_pt(n)->x_gen(t,k,i);}
2191 
2192  /// \short i-th component of time derivative (velocity) of the
2193  /// generalised position, dx(k,i)/dt at local node n.
2194  /// `Type': k; Coordinate direction: i. Do not use the hanging node
2195  /// representation.
2196  double raw_dnodal_position_gen_dt(const unsigned &n,
2197  const unsigned &k,
2198  const unsigned &i) const
2199  {return node_pt(n)->dx_gen_dt(k,i);}
2200 
2201  /// \short i-th component of j-th time derivative of the
2202  /// generalised position, dx(k,i)/dt at local node n.
2203  /// `Type': k; Coordinate direction: i. Do not use the hanging node
2204  /// representation.
2205  double raw_dnodal_position_gen_dt(const unsigned &j,
2206  const unsigned &n,
2207  const unsigned &k,
2208  const unsigned &i) const
2209  {return node_pt(n)->dx_gen_dt(j,k,i);}
2210 
2211 
2212  /// \short Return the i-th coordinate at local node n. If the
2213  /// node is hanging, the appropriate interpolation is handled by the
2214  /// position function in the Node class.
2215  double nodal_position(const unsigned &n, const unsigned &i) const
2216  {return node_pt(n)->position(i);}
2217 
2218  /// \short Return the i-th coordinate at local node n, at time level t
2219  /// (t=0: present; t>0: previous time level)
2220  /// Returns suitably interpolated version for hanging nodes.
2221  double nodal_position(const unsigned &t, const unsigned &n,
2222  const unsigned &i) const
2223  {return node_pt(n)->position(t,i);}
2224 
2225  /// Return the i-th component of nodal velocity: dx/dt at local node n.
2226  double dnodal_position_dt(const unsigned &n, const unsigned &i) const
2227  {return node_pt(n)->dposition_dt(i);}
2228 
2229  /// \short Return the i-th component of j-th derivative of nodal position:
2230  /// d^jx/dt^j at node n.
2231  double dnodal_position_dt(const unsigned &n, const unsigned &j,
2232  const unsigned &i) const
2233  {return node_pt(n)->dposition_dt(j,i);}
2234 
2235  /// \short Return the value of the k-th type of the i-th positional variable
2236  /// at the local node n.
2237  double nodal_position_gen(const unsigned &n, const unsigned &k,
2238  const unsigned &i) const
2239  {return node_pt(n)->position_gen(k,i);}
2240 
2241  /// \short Return the generalised nodal position (type k, i-th variable)
2242  /// at previous timesteps at local node n
2243  double nodal_position_gen(const unsigned &t, const unsigned &n,
2244  const unsigned &k, const unsigned &i) const
2245  {return node_pt(n)->position_gen(t,k,i);}
2246 
2247  /// \short i-th component of time derivative (velocity) of the
2248  /// generalised position, dx(k,i)/dt at local node n.
2249  /// `Type': k; Coordinate direction: i.
2250  double dnodal_position_gen_dt(const unsigned &n,
2251  const unsigned &k,
2252  const unsigned &i) const
2253  {return node_pt(n)->dposition_gen_dt(k,i);}
2254 
2255  /// \short i-th component of j-th time derivative of the
2256  /// generalised position, dx(k,i)/dt at local node n.
2257  /// `Type': k; Coordinate direction: i.
2258  double dnodal_position_gen_dt(const unsigned &j,
2259  const unsigned &n,
2260  const unsigned &k,
2261  const unsigned &i) const
2262  {return node_pt(n)->dposition_gen_dt(j,k,i);}
2263 
2264 
2265  /// \short Compute derivatives of elemental residual vector with respect
2266  /// to nodal coordinates. Default implementation by FD can be overwritten
2267  /// for specific elements.
2268  /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
2270  dresidual_dnodal_coordinates);
2271 
2272  /// \short This is an empty function that establishes a uniform
2273  /// interface for all (derived) elements that involve time-derivatives.
2274  /// Such elements are/should be implemented in ALE form to allow
2275  /// mesh motions. The additional expense associated with the
2276  /// computation of the mesh velocities is, of course, superfluous
2277  /// if the elements are used in problems in which the mesh is
2278  /// stationary. This function should therefore be overloaded
2279  /// in all derived elements that are formulated in ALE form
2280  /// to suppress the computation of the mesh velocities.
2281  /// The user disables the ALE functionality at his/her own risk!
2282  /// If the mesh does move after all, then the results will be wrong.
2283  /// Here we simply issue a warning message stating that the empty
2284  /// function has been called.
2285  virtual void disable_ALE()
2286  {
2287  std::ostringstream warn_message;
2288  warn_message
2289  << "Warning: You have just called the default (empty) function \n\n"
2290  << " FiniteElement::disable_ALE() \n\n"
2291  << "This suggests that you've either tried to call it for an element\n"
2292  << "that \n"
2293  << "(1) does not involve time-derivatives (e.g. a Poisson element) \n"
2294  << "(2) an element for which the time-derivatives aren't implemented \n"
2295  << " in ALE form \n"
2296  << "(3) an element for which the ALE form of the equations can't be \n"
2297  << " be disabled (yet).\n";
2298  OomphLibWarning(warn_message.str(),
2299  "Problem::disable_ALE()",
2300  OOMPH_EXCEPTION_LOCATION);
2301  }
2302 
2303 
2304  /// \short (Re-)enable ALE, i.e. take possible mesh motion into account
2305  /// when evaluating the time-derivative. This function is empty
2306  /// and simply establishes a common interface for all derived
2307  /// elements that are formulated in ALE form.
2308  virtual void enable_ALE()
2309  {
2310  std::ostringstream warn_message;
2311  warn_message
2312  << "Warning: You have just called the default (empty) function \n\n"
2313  << " FiniteElement::enable_ALE() \n\n"
2314  << "This suggests that you've either tried to call it for an element\n"
2315  << "that \n"
2316  << "(1) does not involve time-derivatives (e.g. a Poisson element) \n"
2317  << "(2) an element for which the time-derivatives aren't implemented \n"
2318  << " in ALE form \n"
2319  << "(3) an element for which the ALE form of the equations can't be \n"
2320  << " be disabled (yet)\n"
2321  << "(4) an element for which this function has not been (properly) \n "
2322  << " implemented. This is likely to be a bug!\n ";
2323  OomphLibWarning(warn_message.str(),
2324  "Problem::enable_ALE()",
2325  OOMPH_EXCEPTION_LOCATION);
2326  }
2327 
2328  /// \short Number of values that must be stored at local node n by
2329  /// the element. The default is 0, until over-ridden by a particular element.
2330  /// For example, a Poisson equation requires only one value to be stored
2331  /// at each node; 2D Navier--Stokes equations require two values (velocity
2332  /// components) to be stored at each Node
2333  /// (provided that the pressure interpolation is discontinuous).
2334  virtual unsigned required_nvalue(const unsigned &n) const
2335  {return Default_Initial_Nvalue;}
2336 
2337  /// \short Return the number of coordinate types that the element requires
2338  /// to interpolate the geometry between
2339  /// the nodes. For Lagrange elements it is 1.
2340  unsigned nnodal_position_type() const {return Nnodal_position_type;}
2341 
2342  /// \short Return boolean to indicate if any of the element's nodes
2343  /// are geometrically hanging.
2344  bool has_hanging_nodes() const
2345  {
2346  unsigned nnod=nnode();
2347  for(unsigned j=0;j<nnod;j++)
2348  {
2349  if(node_pt(j)->is_hanging())
2350  {
2351  return true;
2352  }
2353  }
2354  return false;
2355  }
2356 
2357  /// Return the required Eulerian dimension of the nodes in this element
2358  unsigned nodal_dimension() const {return Nodal_dimension;}
2359 
2360  /// Return the number of vertex nodes in this element. Broken virtual
2361  /// function in "pure" finite elements.
2362  virtual unsigned nvertex_node() const
2363  {
2364  std::string error_msg = "Not implemented for FiniteElement.";
2365  throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
2366  OOMPH_EXCEPTION_LOCATION);
2367  }
2368 
2369  /// \short Pointer to the j-th vertex node in the element. Broken virtual
2370  /// function in "pure" finite elements.
2371  virtual Node* vertex_node_pt(const unsigned& j) const
2372  {
2373  std::string error_msg = "Not implemented for FiniteElement.";
2374  throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
2375  OOMPH_EXCEPTION_LOCATION);
2376  }
2377 
2378  /// \short Construct the local node n and return a pointer to the newly
2379  /// created node object.
2380  virtual Node* construct_node(const unsigned &n)
2381  {
2382  //Create a node and assign it to the local node pointer
2383  //The dimension and number of values are taken from internal element data
2385  required_nvalue(n));
2386  //Now return a pointer to the node, so that the mesh can use it
2387  return node_pt(n);
2388  }
2389 
2390  /// \short Construct the local node n, including
2391  /// storage for history values required by timestepper, and return a pointer
2392  /// to the newly created node object.
2393  virtual Node* construct_node(const unsigned &n,
2394  TimeStepper* const &time_stepper_pt)
2395  {
2396  //Create a node and assign it to the local node pointer.
2397  //The dimension and number of values are taken from internal element data
2398  node_pt(n) = new Node(time_stepper_pt, Nodal_dimension,
2400  required_nvalue(n));
2401  //Now return a pointer to the node, so that the mesh can find it
2402  return node_pt(n);
2403  }
2404 
2405  /// \short Construct the local node n as a boundary node;
2406  /// that is a node that MAY be placed on a mesh boundary
2407  /// and return a pointer to the newly created node object.
2408  virtual Node* construct_boundary_node(const unsigned &n)
2409  {
2410  //Create a node and assign it to the local node pointer
2411  //The dimension and number of values are taken from internal element data
2414  required_nvalue(n));
2415  //Now return a pointer to the node, so that the mesh can use it
2416  return node_pt(n);
2417  }
2418 
2419  /// \short Construct the local node n, including
2420  /// storage for history values required by timestepper,
2421  /// as a boundary node; that is a node that MAY be placed on a mesh
2422  /// boundary and return a pointer to the newly created node object.
2423  virtual Node* construct_boundary_node(const unsigned &n,
2424  TimeStepper* const &time_stepper_pt)
2425  {
2426  //Create a node and assign it to the local node pointer.
2427  //The dimension and number of values are taken from internal element data
2430  required_nvalue(n));
2431  //Now return a pointer to the node, so that the mesh can find it
2432  return node_pt(n);
2433  }
2434 
2435  /// \short Return the number of the node *node_pt if this node
2436  /// is in the element, else return -1;
2437  int get_node_number(Node* const &node_pt) const;
2438 
2439 
2440  /// \short If there is a node at this local coordinate, return the pointer to
2441  /// the node
2442  virtual Node* get_node_at_local_coordinate(const Vector<double> &s) const;
2443 
2444  /// \short Return the i-th value stored at local node n
2445  /// but do NOT take hanging nodes into account
2446  double raw_nodal_value(const unsigned &n, const unsigned &i) const
2447  {return node_pt(n)->raw_value(i);}
2448 
2449  /// \short Return the i-th value stored at local node n, at time level t
2450  /// (t=0: present; t>0 previous timesteps), but do NOT take hanging nodes
2451  /// into account
2452  double raw_nodal_value(const unsigned &t, const unsigned &n,
2453  const unsigned &i)
2454  const {return node_pt(n)->raw_value(t,i);}
2455 
2456  /// \short Return the i-th value stored at local node n.
2457  /// Produces suitably interpolated values for hanging nodes.
2458  double nodal_value(const unsigned &n, const unsigned &i) const
2459  {return node_pt(n)->value(i);}
2460 
2461  /// \short Return the i-th value stored at local node n, at time level t
2462  /// (t=0: present; t>0 previous timesteps). Produces suitably interpolated
2463  /// values for hanging nodes.
2464  double nodal_value(const unsigned &t, const unsigned &n, const unsigned &i)
2465  const {return node_pt(n)->value(t,i);}
2466 
2467  /// \short Return the spatial dimension of the element, i.e.
2468  /// the number of local coordinates required to parametrise its
2469  /// geometry.
2470  unsigned dim() const {return Elemental_dimension;}
2471 
2472  /// Return the geometry type of the element (either Q or T usually).
2474  {
2475  std::string err = "Broken virtual function.";
2476  throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
2477  OOMPH_EXCEPTION_LOCATION);
2478  }
2479 
2480  /// Return FE interpolated coordinate x[i] at local coordinate s
2481  virtual double interpolated_x(const Vector<double> &s,
2482  const unsigned &i) const;
2483 
2484  ///\short Return FE interpolated coordinate x[i] at local coordinate s
2485  ///at previous timestep t (t=0: present; t>0: previous timestep)
2486  virtual double interpolated_x(const unsigned& t, const Vector<double> &s,
2487  const unsigned &i) const;
2488 
2489  /// Return FE interpolated position x[] at local coordinate s as Vector
2490  virtual void interpolated_x(const Vector<double> &s,Vector<double>& x) const;
2491 
2492  /// \short Return FE interpolated position x[] at local coordinate s
2493  /// at previous timestep t as Vector (t=0: present; t>0: previous timestep)
2494  virtual void interpolated_x(const unsigned& t, const Vector<double> &s,
2495  Vector<double>& x) const;
2496 
2497  /// \short Return t-th time-derivative of the
2498  /// i-th FE-interpolated Eulerian coordinate at
2499  /// local coordinate s.
2500  virtual double interpolated_dxdt(const Vector<double> &s,
2501  const unsigned &i,
2502  const unsigned &t);
2503 
2504  /// \short Compte t-th time-derivative of the
2505  /// FE-interpolated Eulerian coordinate vector at
2506  /// local coordinate s.
2507  virtual void interpolated_dxdt(const Vector<double> &s,
2508  const unsigned &t,
2509  Vector<double>& dxdt);
2510 
2511  /// \short A standard FiniteElement is fixed, so there are no geometric
2512  /// data when viewed in its GeomObject incarnation
2513  inline unsigned ngeom_data() const {return 0;}
2514 
2515  /// \short A standard FiniteElement is fixed, so there are no geometric
2516  /// data when viewed in its GeomObject incarnation
2517  inline Data* geom_data_pt(const unsigned &j) {return 0;}
2518 
2519  /// \short Return the parametrised position of the FiniteElement in
2520  /// its incarnation as a GeomObject, r(zeta).
2521  /// The position is given by the Eulerian coordinate and the intrinsic
2522  /// coordinate (zeta) is the local coordinate of the element (s).
2523  void position(const Vector<double>& zeta, Vector<double>& r)
2524  const {this->interpolated_x(zeta,r);}
2525 
2526  /// \short Return the parametrised position of the FiniteElement
2527  /// in its GeomObject incarnation:
2528  /// r(zeta). The position is given by the Eulerian coordinate and the
2529  /// intrinsic coordinate (zeta) is the local coordinate of the element (s)
2530  /// This version of the function returns the position as a function
2531  /// of time t=0: current time; t>0: previous
2532  /// timestep. Works for t=0 but needs to be overloaded
2533  /// if genuine time-dependence is required.
2534  void position(const unsigned& t, const Vector<double>& zeta,
2535  Vector<double>& r) const
2536  {this->interpolated_x(t,zeta,r);}
2537 
2538  /// \short Return the t-th time derivative of the
2539  /// parametrised position of the FiniteElement
2540  /// in its GeomObject incarnation:
2541  /// \f$ \frac{d^{t} dr(zeta)}{d t^{t}} \f$.
2542  /// Call the t-th time derivative of the FE-interpolated Eulerian coordinate
2543  void dposition_dt(const Vector<double> &zeta, const unsigned &t,
2544  Vector<double> &drdt)
2545  {this->interpolated_dxdt(zeta,t,drdt);}
2546 
2547  /// \short Specify the values of the "global" intrinsic coordinate, zeta,
2548  /// of a compound geometric object (a mesh of elements) when
2549  /// the element is viewied as a sub-geometric object.
2550  /// The default assumption is that the element will be
2551  /// treated as a sub-geometric object in a bulk Mesh of other elements
2552  /// (geometric objects). The "global" coordinate of the compound geometric
2553  /// object is simply the Eulerian coordinate, x.
2554  /// The second default assumption is that the coordinate zeta will
2555  /// be stored at the nodes and interpolated using the shape functions
2556  /// of the element. This function returns the value of zeta stored at
2557  /// local node n, where k is the type of coordinate and i is the
2558  /// coordinate direction. The function is virtual so that it can
2559  /// be overloaded by different types of element: FaceElements
2560  /// and SolidFiniteElements
2561  virtual double zeta_nodal(const unsigned &n, const unsigned &k,
2562  const unsigned &i) const
2563  {
2564  //By default return the value for nodal_position_gen
2565  return nodal_position_gen(n,k,i);
2566  }
2567 
2568 
2569  /// \short Calculate the interpolated value of zeta, the intrinsic coordinate
2570  /// of the element when viewed as a compound geometric object within a Mesh
2571  /// as a function of the local coordinate of the element, s. The default
2572  /// assumption is the zeta is interpolated using the shape functions of
2573  /// the element with the values given by zeta_nodal().
2574  /// A MacroElement representation of the intrinsic coordinate parametrised
2575  /// by the local coordinate s is used if available.
2576  /// Choosing the MacroElement representation of zeta (Eulerian x by default)
2577  /// allows a correspondence to be established between elements on different
2578  /// Meshes covering the same curvilinear domain in cases where one element
2579  /// is much coarser than the other.
2580  void interpolated_zeta(const Vector<double> &s,
2581  Vector<double> &zeta) const;
2582 
2583  /// \short For a given value of zeta, the "global" intrinsic coordinate of
2584  /// a mesh of FiniteElements represented as a compound geometric object,
2585  /// find the local coordinate in this element that corresponds to the
2586  /// requested value of zeta.
2587  /// If zeta cannot be located in this element, geom_object_pt is set
2588  /// to NULL. If zeta is located in this element, we return its "this"
2589  /// pointer.
2590  /// By default don't use any value passed in to the local coordinate s
2591  /// as the initial guess in the Newton method
2592  void locate_zeta(const Vector<double> &zeta,
2593  GeomObject*& geom_object_pt, Vector<double> &s,
2594  const bool& use_coordinate_as_initial_guess=false);
2595 
2596 
2597  /// \short Update the positions of all nodes in the element using
2598  /// each node update function. The default implementation may
2599  /// be overloaded so that more efficient versions can be written
2600  virtual void node_update();
2601 
2602  /// \short The purpose of this function is to identify all possible
2603  /// Data that can affect the fields interpolated by the FiniteElement.
2604  /// The information will typically be used in interaction problems in
2605  /// which the FiniteElement provides a forcing term for an
2606  /// ElementWithExternalElement. The Data must be provided as
2607  /// \c paired_load data containing
2608  /// - the pointer to a Data object
2609  /// and
2610  /// - the index of the value in that Data object
2611  /// .
2612  /// The generic implementation (should be overloaded in more specific
2613  /// applications) is to include all nodal and internal Data stored in
2614  /// the FiniteElement. The geometric data, which includes the positions
2615  /// of SolidNodes, is treated separately by the function
2616  /// \c identify_geometric_data()
2618  std::set<std::pair<Data*,unsigned> > &paired_field_data);
2619 
2620 
2621  /// \short The purpose of this
2622  /// function is to identify all \c Data objects that affect the elements'
2623  /// geometry. This function is implemented as an empty virtual
2624  /// function since it can only be implemented in conjunction
2625  /// with a node-update strategy. A specific implementation
2626  /// is provided in the ElementWithMovingNodes class.
2627  virtual void identify_geometric_data(std::set<Data*> &geometric_data_pt){}
2628 
2629 
2630  /// Min value of local coordinate
2631  virtual double s_min() const
2632  {
2633  throw OomphLibError(
2634  "s_min() isn't implemented for this element\n",
2635  OOMPH_CURRENT_FUNCTION,
2636  OOMPH_EXCEPTION_LOCATION);
2637  // Dummy return
2638  return 0.0;
2639  }
2640 
2641  /// Max. value of local coordinate
2642  virtual double s_max() const
2643  {
2644  throw OomphLibError(
2645  "s_max() isn't implemented for this element\n",
2646  OOMPH_CURRENT_FUNCTION,
2647  OOMPH_EXCEPTION_LOCATION);
2648  // Dummy return
2649  return 0.0;
2650  }
2651 
2652  /// Calculate the size of the element (length, area, volume,...)
2653  /// in Eulerian computational
2654  /// coordinates. Use suitably overloaded compute_physical_size()
2655  /// function to compute the actual size (taking into account
2656  /// factors such as 2pi or radii the integrand) -- such function
2657  /// can only be implemented on an equation-by-equation basis.
2658  double size() const;
2659 
2660 
2661  /// \short Broken virtual
2662  /// function to compute the actual size (taking into account
2663  /// factors such as 2pi or radii the integrand) -- such function
2664  /// can only be implemented on an equation-by-equation basis.
2665  virtual double compute_physical_size() const
2666  {
2667  throw OomphLibError(
2668  "compute_physical_size() isn't implemented for this element\n",
2669  OOMPH_CURRENT_FUNCTION,
2670  OOMPH_EXCEPTION_LOCATION);
2671  // Dummy return
2672  return 0.0;
2673  }
2674 
2675  /// \short Virtual function to write the double precision numbers that
2676  /// appear in a single line of output into the data vector. Empty virtual,
2677  /// can be overloaded for specific elements; used e.g. by LineVisualiser.
2678  virtual void point_output_data(const Vector<double> &s, Vector<double>& data)
2679  {}
2680 
2681  /// \short Output solution (as defined by point_output_data())
2682  /// at local cordinates s
2683  void point_output(std::ostream &outfile, const Vector<double> &s)
2684  {
2685  // Get point data
2686  Vector<double> data;
2687  this->point_output_data(s,data);
2688 
2689  // Output
2690  unsigned n=data.size();
2691  for (unsigned i=0;i<n;i++)
2692  {
2693  outfile << data[i] << " ";
2694  }
2695  }
2696 
2697  /// \short Return the number of actual plot points for paraview
2698  /// plot with parameter nplot. Broken virtual; can be overloaded
2699  /// in specific elements.
2700  virtual unsigned nplot_points_paraview(const unsigned& nplot) const
2701  {
2702  throw OomphLibError(
2703  "This function hasn't been implemented for this element",
2704  OOMPH_CURRENT_FUNCTION,
2705  OOMPH_EXCEPTION_LOCATION);
2706 
2707  // Dummy unsigned
2708  return 0;
2709  }
2710 
2711  /// \short Return the number of local sub-elements for paraview plot with
2712  /// parameter nplot. Broken virtual; can be overloaded
2713  /// in specific elements.
2714  virtual unsigned nsub_elements_paraview(const unsigned& nplot) const
2715  {
2716  throw OomphLibError(
2717  "This function hasn't been implemented for this element",
2718  OOMPH_CURRENT_FUNCTION,
2719  OOMPH_EXCEPTION_LOCATION);
2720 
2721  // Dummy unsigned
2722  return 0;
2723  }
2724 
2725  /// \short Paraview output -- this outputs the coordinates at the plot
2726  /// points (for parameter nplot) to specified output file.
2727  void output_paraview(std::ofstream& file_out, const unsigned& nplot) const
2728  {
2729  //Decide the dimensions of the nodes
2730  unsigned nnod=nnode();
2731  if (nnod==0) return;
2732  unsigned n=node_pt(0)->ndim();
2733 
2734  // Vector for local coordinates
2735  Vector<double> s(n,0.0);
2736 
2737  //Vector for cartesian coordinates
2738  Vector<double> x(n,0.0);
2739 
2740  //Determine the total number of plotpoints
2741  unsigned plot=nplot_points_paraview(nplot);
2742 
2743  // Loop over the local points
2744  for(unsigned j=0;j<plot;j++)
2745  {
2746  // Determine where in the local element the point is
2747  this->get_s_plot(j,nplot,s);
2748 
2749  // Update the cartesian coordinates vector
2750  this->interpolated_x(s,x);
2751 
2752  // Print the global coordinates. Note: no whitespace after last
2753  // coordinate or paraview is very unhappy.
2754  for(unsigned i=0;i<n-1;i++)
2755  {
2756  file_out << x[i] << " ";
2757  }
2758  file_out << x[n-1];
2759 
2760  // Since unstructured grid always needs 3 components for each
2761  // point, output 0's by default
2762  switch(n)
2763  {
2764  case 1:
2765  file_out << " 0" << " 0"<< std::endl;
2766  break;
2767 
2768  case 2:
2769  file_out << " 0"<< std::endl;
2770  break;
2771 
2772  case 3:
2773  file_out << std::endl;
2774  break;
2775 
2776  // Paraview can't handle more than 3 dimensions, output error
2777  default:
2778  throw OomphLibError(
2779  "Printing PlotPoint to .vtu failed; it has >3 dimensions.",
2780  OOMPH_CURRENT_FUNCTION,
2781  OOMPH_EXCEPTION_LOCATION);
2782  }
2783  }
2784  }
2785 
2786  /// \short Fill in the offset information for paraview plot. Broken virtual.
2787  /// Needs to be implemented for each new geometric element type; see
2788  /// http://www.vtk.org/VTK/img/file-formats.pdf
2789  virtual void write_paraview_output_offset_information(std::ofstream& file_out,
2790  const unsigned& nplot,
2791  unsigned& counter) const
2792  {
2793  throw OomphLibError(
2794  "This function hasn't been implemented for this element",
2795  OOMPH_CURRENT_FUNCTION,
2796  OOMPH_EXCEPTION_LOCATION);
2797  }
2798 
2799  /// \short Return the paraview element type. Broken virtual.
2800  /// Needs to be implemented for each new geometric element type; see
2801  /// http://www.vtk.org/VTK/img/file-formats.pdf
2802  virtual void write_paraview_type(std::ofstream& file_out,
2803  const unsigned& nplot) const
2804  {
2805  throw OomphLibError(
2806  "This function hasn't been implemented for this element",
2807  OOMPH_CURRENT_FUNCTION,
2808  OOMPH_EXCEPTION_LOCATION);
2809  }
2810 
2811  /// \short Return the offsets for the paraview sub-elements. Broken
2812  /// virtual. Needs to be implemented for each new geometric element type; see
2813  /// http://www.vtk.org/VTK/img/file-formats.pdf
2814  virtual void write_paraview_offsets(std::ofstream& file_out,
2815  const unsigned& nplot,
2816  unsigned& offset_sum) const
2817  {
2818  throw OomphLibError(
2819  "This function hasn't been implemented for this element",
2820  OOMPH_CURRENT_FUNCTION,
2821  OOMPH_EXCEPTION_LOCATION);
2822  }
2823 
2824  /// \short Number of scalars/fields output by this element. Broken
2825  /// virtual. Needs to be implemented for each new specific element type.
2826  virtual unsigned nscalar_paraview() const
2827  {
2828  throw OomphLibError(
2829  "This function hasn't been implemented for this element",
2830  OOMPH_CURRENT_FUNCTION,
2831  OOMPH_EXCEPTION_LOCATION);
2832 
2833  // Dummy unsigned
2834  return 0;
2835  }
2836 
2837  /// \short Write values of the i-th scalar field at the plot points. Broken
2838  /// virtual. Needs to be implemented for each new specific element type.
2839  virtual void scalar_value_paraview(std::ofstream& file_out,
2840  const unsigned& i,
2841  const unsigned& nplot) const
2842  {
2843  throw OomphLibError(
2844  "This function hasn't been implemented for this element",
2845  OOMPH_CURRENT_FUNCTION,
2846  OOMPH_EXCEPTION_LOCATION);
2847  }
2848 
2849  /// \short Write values of the i-th scalar field at the plot points. Broken
2850  /// virtual. Needs to be implemented for each new specific element type.
2851  virtual void scalar_value_fct_paraview(std::ofstream& file_out,
2852  const unsigned& i,
2853  const unsigned& nplot,
2854  FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
2855  {
2856  throw OomphLibError(
2857  "This function hasn't been implemented for this element",
2858  OOMPH_CURRENT_FUNCTION,
2859  OOMPH_EXCEPTION_LOCATION);
2860  }
2861 
2862  /// \short Name of the i-th scalar field. Default implementation
2863  /// returns V1 for the first one, V2 for the second etc. Can (should!) be
2864  /// overloaded with more meaningful names in specific elements.
2865  virtual std::string scalar_name_paraview(const unsigned& i) const
2866  {
2867  return "V"+StringConversion::to_string(i);
2868  }
2869 
2870  /// \short Output the element data --- typically the values at the
2871  /// nodes in a format suitable for post-processing.
2872  virtual void output(std::ostream &outfile)
2873  {
2874  outfile << "Output function undefined" << std::endl;
2875  }
2876 
2877  /// \short Output the element data --- pass (some measure of)
2878  /// the number of plot points per element
2879  virtual void output(std::ostream &outfile, const unsigned &n_plot)
2880  {
2881  outfile << "Output function undefined" << std::endl;
2882  }
2883 
2884  /// \short Output the element data at time step t. This is const because
2885  /// it is newly added and so can be done easily. Really all the
2886  /// output(...) functions should be const!
2887  virtual void output(const unsigned& t,
2888  std::ostream &outfile,
2889  const unsigned &n_plot) const
2890  {
2891  outfile << "Output function undefined" << std::endl;
2892  }
2893 
2894  /// \short Output the element data --- typically the values at the
2895  /// nodes in a format suitable for post-processing.
2896  /// (C style output)
2897  virtual void output(FILE* file_pt)
2898  {
2899  fprintf(file_pt,"C_style output function undefined\n");
2900  }
2901 
2902  /// \short Output the element data --- pass (some measure of)
2903  /// the number of plot points per element
2904  /// (C style output)
2905  virtual void output(FILE* file_pt, const unsigned &n_plot)
2906  {
2907  fprintf(file_pt,"C_style output function undefined\n");
2908  }
2909 
2910  /// Output an exact solution over the element.
2911  virtual void output_fct(std::ostream &outfile, const unsigned &n_plot,
2913  {
2914  outfile << "Output function undefined for exact solution \n";
2915  }
2916 
2917 
2918  /// Output a time-dependent exact solution over the element.
2919  virtual void output_fct(std::ostream &outfile, const unsigned &n_plot,
2920  const double& time,
2922  {
2923  outfile << "Output function undefined for exact solution \n";
2924  }
2925 
2926  /// Output a time-dependent exact solution over the element.
2927  virtual void output_fct(std::ostream &outfile, const unsigned &n_plot,
2928  const double& time,
2929  const SolutionFunctorBase& exact_soln) const
2930  {
2931  outfile << "Output function undefined for exact solution \n";
2932  }
2933 
2934 
2935 
2936  /// \short Get cector of local coordinates of plot point i (when plotting
2937  /// nplot points in each "coordinate direction").
2938  virtual void get_s_plot(const unsigned& i, const unsigned& nplot,
2939  Vector<double>& s) const //senare
2940  {
2941  throw OomphLibError(
2942  "get_s_plot(...) is not implemented for this element\n",
2943  OOMPH_CURRENT_FUNCTION,
2944  OOMPH_EXCEPTION_LOCATION);
2945  };
2946 
2947  /// \short Return string for tecplot zone header (when plotting
2948  /// nplot points in each "coordinate direction")
2949  virtual std::string tecplot_zone_string(const unsigned& nplot) const
2950  {
2951  throw OomphLibError(
2952  "tecplot_zone_string(...) is not implemented for this element\n",
2953  OOMPH_CURRENT_FUNCTION,
2954  OOMPH_EXCEPTION_LOCATION);
2955  return "dummy return";
2956  }
2957 
2958  /// \short Add tecplot zone "footer" to output stream (when plotting
2959  /// nplot points in each "coordinate direction").
2960  /// Empty by default -- can be used, e.g., to add FE connectivity
2961  /// lists to elements that need it.
2962  virtual void write_tecplot_zone_footer(std::ostream& outfile,
2963  const unsigned& nplot)
2964  const {};
2965 
2966  /// \short Add tecplot zone "footer" to C-style output. (when plotting
2967  /// nplot points in each "coordinate direction").
2968  /// Empty by default -- can be used, e.g., to add FE connectivity
2969  /// lists to elements that need it.
2970  virtual void write_tecplot_zone_footer(FILE* file_pt,
2971  const unsigned& nplot)
2972  const {};
2973 
2974  /// \short Return total number of plot points (when plotting
2975  /// nplot points in each "coordinate direction")
2976  virtual unsigned nplot_points(const unsigned& nplot) const
2977  {
2978  throw OomphLibError(
2979  "nplot_points(...) is not implemented for this element",
2980  OOMPH_CURRENT_FUNCTION,
2981  OOMPH_EXCEPTION_LOCATION);
2982  // Dummy return
2983  return 0;
2984  }
2985 
2986 
2987 
2988  /// \short Plot the error when compared
2989  /// against a given exact solution \f$ {\bf f}({\bf x}) \f$.
2990  /// Also calculates the norm of the
2991  /// error and that of the exact solution.
2992  virtual void compute_error(std::ostream &outfile,
2994  double& error, double& norm)
2995  {
2996  std::string error_message = "compute_error undefined for this element \n";
2997  outfile << error_message;
2998 
2999  throw OomphLibError(error_message,
3000  OOMPH_CURRENT_FUNCTION,
3001  OOMPH_EXCEPTION_LOCATION);
3002 
3003  }
3004 
3005  /// \short Plot the error when compared
3006  /// against a given time-dependent exact solution \f$ {\bf f}(t,{\bf x}) \f$.
3007  /// Also calculates the norm of the error and that of the exact solution.
3008  virtual void compute_error(std::ostream &outfile,
3010  const double& time, double& error, double& norm)
3011  {
3012  std::string error_message =
3013  "time-dependent compute_error undefined for this element \n";
3014  outfile << error_message;
3015 
3016  throw OomphLibError(error_message,
3017  OOMPH_CURRENT_FUNCTION,
3018  OOMPH_EXCEPTION_LOCATION);
3019  }
3020 
3021  /// \short Plot the error when compared
3022  /// against a given exact solution \f$ {\bf f}({\bf x}) \f$.
3023  /// Also calculates the norm of the
3024  /// error and that of the exact solution.
3025  /// Version with vectors of norms and errors so that different variables' norms
3026  /// and errors can be returned individually
3027  virtual void compute_error(std::ostream &outfile,
3029  exact_soln_pt,
3030  Vector<double>& error, Vector<double>& norm)
3031  {
3032  std::string error_message = "compute_error undefined for this element \n";
3033  outfile << error_message;
3034 
3035  throw OomphLibError(error_message,
3036  "FiniteElement::compute_error()",
3037  OOMPH_EXCEPTION_LOCATION);
3038 
3039  }
3040 
3041  /// \short Plot the error when compared
3042  /// against a given time-dependent exact solution \f$ {\bf f}(t,{\bf x}) \f$.
3043  /// Also calculates the norm of the error and that of the exact solution.
3044  /// Version with vectors of norms and errors so that different variables' norms
3045  /// and errors can be returned individually
3046  virtual void compute_error(std::ostream &outfile,
3048  exact_soln_pt,
3049  const double& time,
3050  Vector<double>& error,
3051  Vector<double>& norm)
3052  {
3053  std::string error_message =
3054  "time-dependent compute_error undefined for this element \n";
3055  outfile << error_message;
3056 
3057  throw OomphLibError(error_message,
3058  "FiniteElement::compute_error()",
3059  OOMPH_EXCEPTION_LOCATION);
3060 
3061  }
3062 
3063  /// \short Plot the error when compared
3064  /// against a given exact solution \f$ {\bf f}({\bf x}) \f$.
3065  /// Also calculates the maximum absolute error
3066  virtual void compute_abs_error(std::ostream &outfile,
3068  exact_soln_pt,
3069  double& error)
3070  {
3071  std::string error_message =
3072  "compute_abs_error undefined for this element \n";
3073  outfile << error_message;
3074 
3075  throw OomphLibError(error_message,
3076  OOMPH_CURRENT_FUNCTION,
3077  OOMPH_EXCEPTION_LOCATION);
3078  }
3079 
3080  /// \short Evaluate integral of a Vector-valued function
3081  /// \f$ {\bf f}({\bf x}) \f$ over the element.
3083  Vector<double>& integral);
3084 
3085 
3086  /// \short Evaluate integral of a Vector-valued, time-dependent function
3087  /// \f$ {\bf f}(t,{\bf x}) \f$ over the element.
3089  const double& time,
3090  Vector<double>& integral);
3091 
3092  /// \short Function for building a lower dimensional FaceElement
3093  /// on the specified face of the FiniteElement.
3094  /// The arguments are the index of the face, an integer whose value
3095  /// depends on the particular element type, and a pointer to the
3096  /// FaceElement.
3097  virtual void build_face_element(const int &face_index,
3098  FaceElement* face_element_pt);
3099 
3100 
3101  /// \short Self-test: Check inversion of element & do self-test for
3102  /// GeneralisedElement. Return 0 if OK.
3103  virtual unsigned self_test();
3104 
3105  /// Get the number of the ith node on face face_index (in the bulk node
3106  /// vector).
3107  virtual unsigned get_bulk_node_number(const int& face_index,
3108  const unsigned& i) const
3109  {
3110  std::string err = "Not implemented for this element.";
3111  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
3112  OOMPH_CURRENT_FUNCTION);
3113  }
3114 
3115  /// Get the sign of the outer unit normal on the face given by face_index.
3116  virtual int face_outer_unit_normal_sign(const int& face_index) const
3117  {
3118  std::string err = "Not implemented for this element.";
3119  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
3120  OOMPH_CURRENT_FUNCTION);
3121  }
3122 
3123  virtual unsigned nnode_on_face() const
3124  {
3125  std::string err = "Not implemented for this element.";
3126  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
3127  OOMPH_CURRENT_FUNCTION);
3128  }
3129 
3130  /// Range check for face node numbers
3131  void face_node_number_error_check(const unsigned& i) const
3132  {
3133 #ifdef RANGE_CHECKING
3134  if(i > nnode_on_face())
3135  {
3136  std::string err = "Face node index i out of range on face.";
3137  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
3138  OOMPH_CURRENT_FUNCTION);
3139  }
3140 #endif
3141  }
3142 
3143  /// Get a pointer to the function mapping face coordinates to bulk coordinates
3145  (const int& face_index) const
3146  {
3147  std::string err = "Not implemented for this element.";
3148  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
3149  OOMPH_CURRENT_FUNCTION);
3150  }
3151 
3152  /// Get a pointer to the derivative of the mapping from face to bulk
3153  /// coordinates.
3155  (const int& face_index) const
3156  {
3157  std::string err = "Not implemented for this element.";
3158  throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
3159  OOMPH_CURRENT_FUNCTION);
3160  }
3161 
3162 };
3163 
3164 
3165 //////////////////////////////////////////////////////////////////////////
3166 //////////////////////////////////////////////////////////////////////////
3167 //////////////////////////////////////////////////////////////////////////
3168 
3169 
3170 
3171 //=======================================================================
3172 /// Point element has just a single node and a single shape function
3173 /// which is identically equal to one.
3174 //=======================================================================
3175 class PointElement : public virtual FiniteElement
3176 {
3177 
3178  public:
3179 
3180  /// Constructor
3182  {
3183  /// Allocate storage for the pointers to the nodes,
3184  /// There is only one nodes
3185  this->set_n_node(1);
3186 
3187  // Set the integration scheme
3189  }
3190 
3191  /// Broken copy constructor
3193  {
3194  BrokenCopy::broken_copy("PointElement");
3195  }
3196 
3197  /// Broken assignment operator
3198  /*void operator=(const PointElement&)
3199  {
3200  BrokenCopy::broken_assign("PointElement");
3201  }*/
3202 
3203  /// Calculate the geometric shape functions at local coordinate s
3204  void shape(const Vector<double> &s, Shape &psi) const;
3205 
3206  /// Get local coordinates of node j in the element; vector sets its own size
3207  void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
3208  {
3209  s.resize(0);
3210  }
3211 
3212  /// \short Return spatial dimension of element (=number of local coordinates
3213  /// needed to parametrise the element)
3214  //unsigned dim() const {return 0;}
3215 
3216  private:
3217 
3218  /// \short Default integration scheme
3220 
3221 
3222 };
3223 
3224 
3225 
3226 //////////////////////////////////////////////////////////////////////////
3227 //////////////////////////////////////////////////////////////////////////
3228 //////////////////////////////////////////////////////////////////////////
3229 
3230 
3231 class GeomObject;
3232 
3233 //=======================================================================
3234 /// \short A class to specify the initial conditions for a solid body.
3235 /// Solid bodies are often discretised with
3236 /// Hermite-type elements, for which the assignment
3237 /// of the generalised nodal values is nontrivial since they represent
3238 /// derivatives w.r.t. to the local coordinates. A SolidInitialCondition
3239 /// object specifies initial position (i.e. shape), velocity and acceleration
3240 /// of the structure with a geometric object.
3241 /// An integer specifies which time-derivative derivative is currently
3242 /// assigned. See example codes for a demonstration of its use.
3243 //=======================================================================
3245 {
3246 
3247  public:
3248 
3249  /// Constructor: Pass geometric object; initialise time deriv to 0
3251  Geom_object_pt(geom_object_pt), IC_time_deriv(0)
3252  {};
3253 
3254 
3255  /// Broken copy constructor
3257  {
3258  BrokenCopy::broken_copy("SolidInitialCondition");
3259  }
3260 
3261  /// Broken assignment operator
3263  {
3264  BrokenCopy::broken_assign("SolidInitialCondition");
3265  }
3266 
3267 
3268  /// (Reference to) pointer to geom object that specifies the initial condition
3270  {
3271  return Geom_object_pt;
3272  }
3273 
3274  /// Which time derivative are we currently assigning?
3275  unsigned& ic_time_deriv()
3276  {
3277  return IC_time_deriv;
3278  }
3279 
3280 
3281  private:
3282 
3283  /// \short Pointer to the GeomObject that specifies the initial
3284  /// condition (shape, veloc and accel)
3286 
3287  /// \short Which time derivative (0,1,2) are we currently assigning
3288  unsigned IC_time_deriv;
3289 
3290 };
3291 
3292 
3293 /////////////////////////////////////////////////////////////////////////
3294 /////////////////////////////////////////////////////////////////////////
3295 /////////////////////////////////////////////////////////////////////////
3296 
3297 
3298 //============================================================================
3299 /// \short SolidFiniteElement class.
3300 ///
3301 /// Solid elements are elements whose nodal positions are unknowns
3302 /// in the problem -- their nodes are SolidNodes. In such elements,
3303 /// the nodes not only have a variable (Eulerian) but also a fixed
3304 /// (Lagrangian) position. The positional variables have their own
3305 /// local equation numbering scheme which is set up with
3306 /// \code assign_solid_local_eqn_numbers () \endcode
3307 /// The derivatives of the
3308 /// `solid equations' (i.e. the equations that determine the
3309 /// nodal positions) with respect to the nodal positions, required
3310 /// in the Jacobian matrix, are determined by finite differencing.
3311 ///
3312 /// In the present form, the SolidFiniteElement represents
3313 /// a fully functional base class for `proper' solid mechanics elements,
3314 /// but it can also be combined
3315 /// (via inheritance) with elements that solve additional equations.
3316 /// This is particularly useful in cases where the solid equations
3317 /// are merely used to update the nodal positions in a moving mesh
3318 /// problem, although this can prove costly in practice.
3319 //============================================================================
3320 class SolidFiniteElement : public virtual FiniteElement
3321 {
3322 
3323 public:
3324 
3325  /// \short Set the lagrangian dimension of the element --- the number
3326  /// of lagrangian coordinates stored at the nodes in the element.
3329 
3330  /// \short Return whether there is internal solid data (e.g. discontinuous
3331  /// solid pressure). At present, this is used to report an error in
3332  /// setting initial conditions for ElasticProblems which can't
3333  /// handle such cases. The default is false.
3334  virtual bool has_internal_solid_data() {return false;}
3335 
3336  /// \short Pointer to function that computes the "multiplier" for the
3337  /// inertia terms in the consistent determination of the initial
3338  /// conditions for Newmark timestepping.
3339  typedef double (*MultiplierFctPt)(const Vector<double>& xi);
3340 
3341  ///Constructor: Set defaults
3347  {}
3348 
3349  ///Destructor to clean up any allocated memory
3350  virtual ~SolidFiniteElement();
3351 
3352  /// Broken copy constructor
3354  {
3355  BrokenCopy::broken_copy("SolidFiniteElement");
3356  }
3357 
3358 
3359  /// Broken assignment operator
3360  /*void operator=(const SolidFiniteElement&)
3361  {
3362  BrokenCopy::broken_assign("SolidFiniteElement");
3363  }*/
3364 
3365  ///\short The number of geometric data affecting a SolidFiniteElemnet is
3366  ///the same as the number of nodes (one variable position data per node)
3367  inline unsigned ngeom_data() const {return nnode();}
3368 
3369  /// \short Return pointer to the j-th Data item that the object's
3370  /// shape depends on. (Redirects to the nodes' positional Data).
3371  inline Data* geom_data_pt(const unsigned& j)
3372  {return static_cast<SolidNode*>(node_pt(j))->variable_position_pt();}
3373 
3374  /// \short Specify Data that affects the geometry of the element
3375  /// by adding the position Data to the set that's passed in.
3376  /// (This functionality is required in FSI problems; set is used to
3377  /// avoid double counting).
3378  void identify_geometric_data(std::set<Data*> &geometric_data_pt)
3379  {
3380  //Loop over the node update data and add to the set
3381  const unsigned n_node=this->nnode();
3382  for(unsigned n=0;n<n_node;n++)
3383  {
3384  geometric_data_pt.insert(dynamic_cast<SolidNode*>(this->node_pt(n))
3385  ->variable_position_pt());
3386  }
3387  }
3388 
3389 
3390  /// \short In a SolidFiniteElement, the "global" intrinsic coordinate
3391  /// of the element when viewed as part of a compound geometric
3392  /// object (a Mesh) is, by default, the Lagrangian coordinate
3393  /// Note the assumption here is that we are always using isoparameteric
3394  /// elements in which the lagrangian coordinate is interpolated by the
3395  /// same shape functions as the eulerian coordinate.
3396  inline double zeta_nodal(const unsigned &n, const unsigned &k,
3397  const unsigned &i) const
3398  {return lagrangian_position_gen(n,k,i);}
3399 
3400  /// \short Eulerian and Lagrangian coordinates as function of the
3401  /// local coordinates: The Eulerian position is returned in
3402  /// FE-interpolated form (\c x_fe) and then in the form obtained
3403  /// from the "current" MacroElement representation (if it exists -- if not,
3404  /// \c x is the same as \c x_fe). This allows the Domain/MacroElement-
3405  /// based representation to be used to apply displacement boundary
3406  /// conditions exactly. Ditto for the Lagrangian coordinates returned
3407  /// in xi_fe and xi.
3408  /// (Broken virtual -- overload in specific geometric element class
3409  /// if you want to use this functionality.)
3410  virtual void get_x_and_xi(const Vector<double>& s,
3411  Vector<double>& x_fe,
3412  Vector<double>& x,
3413  Vector<double>& xi_fe,
3414  Vector<double>& xi) const
3415  {
3416  throw OomphLibError(
3417  "get_x_and_xi(...) is not implemented for this element\n",
3418  OOMPH_CURRENT_FUNCTION,
3419  OOMPH_EXCEPTION_LOCATION);
3420  }
3421 
3422  /// \short Set pointer to MacroElement -- overloads generic version
3423  /// and uses the MacroElement
3424  /// also as the default for the "undeformed" configuration.
3425  /// This assignment must be overwritten with
3426  /// set_undeformed_macro_elem_pt(...) if the deformation of
3427  /// the solid body is driven by a deformation of the
3428  /// "current" Domain/MacroElement representation of it's boundary.
3429  /// Can be overloaded in derived classes to perform additional
3430  /// tasks
3432  {
3435  }
3436 
3437  /// \short Set pointers to "current" and "undeformed" MacroElements.
3438  /// Can be overloaded in derived classes to perform additional
3439  /// tasks
3442  {
3445  }
3446 
3447  /// \short Set pointer to "undeformed" macro element.
3448  /// Can be overloaded in derived classes to perform additional
3449  /// tasks.
3453 
3454 
3455  /// Access function to pointer to "undeformed" macro element
3457  {return Undeformed_macro_elem_pt;}
3458 
3459  /// \short Calculate shape functions and derivatives w.r.t. Lagrangian
3460  /// coordinates at local coordinate s. Returns the Jacobian of the mapping
3461  /// from Lagrangian to local coordinates.
3462  double dshape_lagrangian(const Vector<double> &s, Shape &psi,
3463  DShape &dpsidxi) const;
3464 
3465  /// \short Return the geometric shape functions and also first
3466  /// derivatives w.r.t. Lagrangian coordinates at ipt-th integration point.
3467  virtual double dshape_lagrangian_at_knot(const unsigned &ipt,
3468  Shape &psi, DShape &dpsidxi) const;
3469 
3470  /// \short Compute the geometric shape functions and also first
3471  /// and second derivatives w.r.t. Lagrangian coordinates at
3472  /// local coordinate s;
3473  /// Returns Jacobian of mapping from Lagrangian to local coordinates.
3474  /// Numbering:
3475  /// \b 1D:
3476  /// d2pidxi(i,0) = \f$ d^2 \psi_j / d \xi^2 \f$
3477  /// \b 2D:
3478  /// d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial \xi_0^2 \f$
3479  /// d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial \xi_1^2 \f$
3480  /// d2psidxi(i,2) = \f$ \partial^2 \psi_j / \partial \xi_0 \partial \xi_1 \f$
3481  /// \b 3D:
3482  /// d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial \xi_0^2 \f$
3483  /// d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial \xi_1^2 \f$
3484  /// d2psidxi(i,2) = \f$ \partial^2 \psi_j / \partial \xi_2^2 \f$
3485  /// d2psidxi(i,3) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_1 \f$
3486  /// d2psidxi(i,4) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_2 \f$
3487  /// d2psidxi(i,5) = \f$ \partial^2 \psi_j/\partial \xi_1 \partial \xi_2 \f$
3488  double d2shape_lagrangian(const Vector<double> &s, Shape &psi,
3489  DShape &dpsidxi, DShape &d2psidxi) const;
3490 
3491  /// \short Return the geometric shape functions and also first
3492  /// and second derivatives w.r.t. Lagrangian coordinates at
3493  /// the ipt-th integration point.
3494  /// Returns Jacobian of mapping from Lagrangian to local coordinates.
3495  /// Numbering:
3496  /// \b 1D:
3497  /// d2pidxi(i,0) = \f$ d^2 \psi_j / d^2 \xi^2 \f$
3498  /// \b 2D:
3499  /// d2psidxi(i,0) = \f$ \partial^2 \psi_j/\partial^2 \xi_0^2 \f$
3500  /// d2psidxi(i,1) = \f$ \partial^2 \psi_j/\partial^2 \xi_1^2 \f$
3501  /// d2psidxi(i,2) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_1 \f$
3502  /// \b 3D:
3503  /// d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial^2 \xi_0^2 \f$
3504  /// d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial^2 \xi_1^2 \f$
3505  /// d2psidxi(i,2) = \f$ \partial^2 \psi_j / \partial^2 \xi_2^2 \f$
3506  /// d2psidxi(i,3) = \f$ \partial^2 \psi_j / \partial \xi_0 \partial \xi_1 \f$
3507  /// d2psidxi(i,4) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_2 \f$
3508  /// d2psidxi(i,5) = \f$ \partial^2 \psi_j/\partial \xi_1 \partial \xi_2 \f$
3509  virtual double d2shape_lagrangian_at_knot(const unsigned &ipt,
3510  Shape &psi,
3511  DShape &dpsidxi,
3512  DShape &d2psidxi) const;
3513 
3514  /// \short Return the number of Lagrangian coordinates that the
3515  /// element requires at all nodes.
3516  /// This is by default the elemental dimension. If we ever need any
3517  /// other case, it can be implemented.
3518  unsigned lagrangian_dimension() const {return Lagrangian_dimension;}
3519 
3520  /// \short Return the number of types of (generalised)
3521  /// nodal Lagrangian coordinates required to
3522  /// interpolate the Lagrangian coordinates in the element. (E.g.
3523  /// 1 for Lagrange-type elements; 2 for Hermite beam elements;
3524  /// 4 for Hermite shell elements). Default value is 1. Needs to
3525  /// be overloaded for any other element.
3527 
3528  /// Construct the local node n and return a pointer to it.
3529  Node* construct_node(const unsigned &n)
3530  {
3531  //Construct a solid node and assign it to the local node pointer vector.
3532  //The dimension and number of values are taken from internal element data
3533  //The number of solid pressure dofs are also taken from internal data
3534  //The number of timesteps to be stored comes from the problem!
3537  nodal_dimension(),
3539  required_nvalue(n));
3540  //Now return a pointer to the node, so that the mesh can find it
3541  return node_pt(n);
3542  }
3543 
3544  ///\short Construct the local node n and return
3545  /// a pointer to it. Additionally, create storage for `history'
3546  /// values as required by timestepper
3547  Node* construct_node(const unsigned &n,
3548  TimeStepper* const &time_stepper_pt)
3549  {
3550  //Construct a solid node and assign it to the local node pointer vector
3551  //The dimension and number of values are taken from internal element data
3552  //The number of solid pressure dofs are also taken from internal data
3553  node_pt(n) = new SolidNode(time_stepper_pt,
3556  nodal_dimension(),
3558  required_nvalue(n));
3559  //Now return a pointer to the node, so that the mesh can find it
3560  return node_pt(n);
3561  }
3562 
3563  /// \short Construct the local node n and return a pointer to it.
3564  /// in the case when it is a boundary node; that is it MAY be
3565  /// located on a Mesh boundary
3566  Node* construct_boundary_node(const unsigned &n)
3567  {
3568  //Construct a solid node and assign it to the local node pointer vector.
3569  //The dimension and number of values are taken from internal element data
3570  //The number of solid pressure dofs are also taken from internal data
3571  //The number of timesteps to be stored comes from the problem!
3572  node_pt(n) =
3575  nodal_dimension(),
3577  required_nvalue(n));
3578  //Now return a pointer to the node, so that the mesh can find it
3579  return node_pt(n);
3580  }
3581 
3582  ///\short Construct the local node n and return
3583  /// a pointer to it, in the case when the node MAY be located
3584  /// on a boundary. Additionally, create storage for `history'
3585  /// values as required by timestepper
3586  Node* construct_boundary_node(const unsigned &n,
3587  TimeStepper* const &time_stepper_pt)
3588  {
3589  //Construct a solid node and assign it to the local node pointer vector
3590  //The dimension and number of values are taken from internal element data
3591  //The number of solid pressure dofs are also taken from internal data
3592  node_pt(n) =
3596  nodal_dimension(),
3598  required_nvalue(n));
3599  //Now return a pointer to the node, so that the mesh can find it
3600  return node_pt(n);
3601  }
3602 
3603  /// \short Overload assign_all_generic_local_equation numbers to
3604  /// include the data associated with solid dofs.
3605  /// It remains virtual so that it can be overloaded
3606  /// by RefineableSolidElements. If the boolean argument is true
3607  /// then the degrees of freedom are stored in Dof_pt
3609  const bool &store_local_dof_pt)
3610  {
3611  //Call the standard finite element equation numbering
3612  //(internal, external and nodal data).
3614  //Assign the numbering for the solid dofs
3615  assign_solid_local_eqn_numbers(store_local_dof_pt);
3616  }
3617 
3618  /// \short Function to describe the local dofs of the element. The ostream
3619  /// specifies the output stream to which the description
3620  /// is written; the string stores the currently
3621  /// assembled output that is ultimately written to the
3622  /// output stream by Data::describe_dofs(...); it is typically
3623  /// built up incrementally as we descend through the
3624  /// call hierarchy of this function when called from
3625  /// Problem::describe_dofs(...)
3626  void describe_local_dofs(std::ostream& out,
3627  const std::string& current_string) const;
3628 
3629  /// \short Return i-th Lagrangian coordinate at local node n without using
3630  /// the hanging representation
3631  double raw_lagrangian_position(const unsigned &n, const unsigned &i) const
3632  {return static_cast<SolidNode*>(node_pt(n))->xi(i);}
3633 
3634  /// \short Return Generalised Lagrangian coordinate at local node n.
3635  /// `Direction' i, `Type' k. Does not use the hanging node representation
3636  double raw_lagrangian_position_gen(const unsigned &n, const unsigned &k,
3637  const unsigned &i) const
3638  {return static_cast<SolidNode*>(node_pt(n))->xi_gen(k,i);}
3639 
3640  /// Return i-th Lagrangian coordinate at local node n
3641  double lagrangian_position(const unsigned &n, const unsigned &i) const
3642  {return static_cast<SolidNode*>(node_pt(n))->lagrangian_position(i);}
3643 
3644  /// \short Return Generalised Lagrangian coordinate at local node n.
3645  /// `Direction' i, `Type' k.
3646  double lagrangian_position_gen(const unsigned &n, const unsigned &k,
3647  const unsigned &i) const
3648  {return static_cast<SolidNode*>(node_pt(n))->lagrangian_position_gen(k,i);}
3649 
3650  /// \short Return i-th FE-interpolated Lagrangian coordinate xi[i] at
3651  /// local coordinate s.
3652  virtual double interpolated_xi(const Vector<double> &s,
3653  const unsigned &i) const;
3654 
3655  /// \short Compute FE interpolated Lagrangian coordinate vector xi[] at
3656  /// local coordinate s as Vector
3657  virtual void interpolated_xi(const Vector<double> &s,
3658  Vector<double>& xi) const;
3659 
3660  /// \short Compute derivatives of FE-interpolated Lagrangian coordinates xi
3661  /// with respect to local coordinates: dxids[i][j]=dxi_i/ds_j.
3662  virtual void interpolated_dxids(const Vector<double> &s,
3663  DenseMatrix<double> &dxids) const;
3664 
3665  /// \short Return the Jacobian of mapping from local to Lagrangian
3666  /// coordinates at local position s. NOT YET IMPLEMENTED
3667  virtual void J_lagrangian(const Vector<double> &s) const
3668  {
3669  // Must be implemented and overloaded in FaceElements
3670  throw OomphLibError("Function not implemented yet",
3671  OOMPH_CURRENT_FUNCTION,
3672  OOMPH_EXCEPTION_LOCATION);
3673  }
3674 
3675  /// \short Return the Jacobian of the mapping from local to Lagrangian
3676  /// coordinates at the ipt-th integration point. NOT YET IMPLEMENTED
3677  virtual double J_lagrangian_at_knot(const unsigned &ipt) const
3678  {
3679  // Must be implemented and overloaded in FaceElements
3680  throw OomphLibError("Function not implemented yet",
3681  OOMPH_CURRENT_FUNCTION,
3682  OOMPH_EXCEPTION_LOCATION);
3683  }
3684 
3685  /// \short Pointer to object that describes the initial condition.
3687  {
3688  return Solid_ic_pt;
3689  }
3690 
3691  /// \short Set to alter the problem being solved when
3692  /// assigning the initial conditions for time-dependent problems:
3693  /// solve for the history value that
3694  /// corresponds to the acceleration in the Newmark scheme by demanding
3695  /// that the PDE is satisifed at the initial time. In this case
3696  /// the Jacobian is replaced by the mass matrix.
3699 
3700  /// \short Set to reset the problem being solved to be the standard problem
3703 
3704  /// \short Access function: Pointer to multiplicator function for assignment
3705  /// of consistent assignement of initial conditions for Newmark scheme
3707  {
3708  return Multiplier_fct_pt;
3709  }
3710 
3711 
3712  /// \short Access function: Pointer to multiplicator function for assignment
3713  /// of consistent assignement of initial conditions for Newmark scheme
3714  /// (const version)
3716  {
3717  return Multiplier_fct_pt;
3718  }
3719 
3720 
3721  /// \short Compute the residuals for the setup of an initial condition.
3722  /// The global equations are:
3723  /// \f[
3724  /// 0 = \int \left( \sum_{j=1}^N \sum_{k=1}^K X_{ijk} \psi_{jk}(\xi_n)
3725  /// - \frac{\partial^D R^{(IC)}_i(\xi_n)}{\partial t^D}
3726  /// \right) \psi_{lm}(\xi_n) \ dv
3727  /// \mbox{ \ \ \ \ for \ \ \ $l=1,...,N, \ \ m=1,...,K$}
3728  /// \f]
3729  /// where \f$ N \f$ is the number of nodes in the mesh and \f$ K \f$
3730  /// the number of generalised nodal coordinates. The initial shape
3731  /// of the solid body, \f$ {\bf R}^{(IC)},\f$ and its time-derivatives
3732  /// are specified via the \c GeomObject that is stored in the
3733  /// \c SolidFiniteElement::SolidInitialCondition object. The latter also
3734  /// stores the order of the time-derivative \f$ D \f$ to be assigned.
3736  {
3737  residuals.initialise(0.0);
3738  fill_in_residuals_for_solid_ic(residuals);
3739  }
3740 
3741  /// \short Fill in the residuals for the setup of an initial condition.
3742  /// The global equations are:
3743  /// \f[
3744  /// 0 = \int \left( \sum_{j=1}^N \sum_{k=1}^K X_{ijk} \psi_{jk}(\xi_n)
3745  /// - \frac{\partial^D R^{(IC)}_i(\xi_n)}{\partial t^D}
3746  /// \right) \psi_{lm}(\xi_n) \ dv
3747  /// \mbox{ \ \ \ \ for \ \ \ $l=1,...,N, \ \ m=1,...,K$}
3748  /// \f]
3749  /// where \f$ N \f$ is the number of nodes in the mesh and \f$ K \f$
3750  /// the number of generalised nodal coordinates. The initial shape
3751  /// of the solid body, \f$ {\bf R}^{(IC)},\f$ and its time-derivatives
3752  /// are specified via the \c GeomObject that is stored in the
3753  /// \c SolidFiniteElement::SolidInitialCondition object. The latter also
3754  /// stores the order of the time-derivative \f$ D \f$ to be assigned.
3756  {
3757  //Call the generic residuals function with flag set to 0
3758  //using a dummy matrix argument
3760  residuals,GeneralisedElement::Dummy_matrix,0);
3761  }
3762 
3763  /// \short Fill in the residuals and Jacobian for the setup of an
3764  /// initial condition. The global equations are:
3765  /// \f[
3766  /// 0 = \int \left( \sum_{j=1}^N \sum_{k=1}^K X_{ijk} \psi_{jk}(\xi_n)
3767  /// - \frac{\partial^D R^{(IC)}_i(\xi_n)}{\partial t^D}
3768  /// \right) \psi_{lm}(\xi_n) \ dv
3769  /// \mbox{ \ \ \ \ for \ \ \ $l=1,...,N, \ \ m=1,...,K$}
3770  /// \f]
3771  /// where \f$ N \f$ is the number of nodes in the mesh and \f$ K \f$
3772  /// the number of generalised nodal coordinates. The initial shape
3773  /// of the solid body, \f$ {\bf R}^{(IC)},\f$ and its time-derivatives
3774  /// are specified via the \c GeomObject that is stored in the
3775  /// \c SolidFiniteElement::SolidInitialCondition object. The latter also
3776  /// stores the order of the time-derivative \f$ D \f$ to be assigned.
3778  DenseMatrix<double> &jacobian)
3779  {
3780  //Call the generic routine with the flag set to 1
3781  fill_in_generic_jacobian_for_solid_ic(residuals,jacobian,1);
3782  }
3783 
3784 
3785  /// \short Fill in the contributions of the Jacobian matrix
3786  /// for the consistent assignment of the initial "accelerations" in
3787  /// Newmark scheme. In this case the Jacobian is the mass matrix.
3789 
3790 
3791 
3792  /// \short Calculate the L2 norm of the displacement u=R-r to overload the
3793  /// compute_norm function in the GeneralisedElement base class
3794  void compute_norm(double& el_norm);
3795 
3796  protected:
3797 
3798  /// \short Helper function to fill in the residuals and (if flag==1) the
3799  /// Jacobian for the setup of an initial condition. The global equations are:
3800  /// \f[
3801  /// 0 = \int \left( \sum_{j=1}^N \sum_{k=1}^K X_{ijk} \psi_{jk}(\xi_n)
3802  /// - \frac{\partial^D R^{(IC)}_i(\xi_n)}{\partial t^D}
3803  /// \right) \psi_{lm}(\xi_n) \ dv
3804  /// \mbox{ \ \ \ \ for \ \ \ $l=1,...,N, \ \ m=1,...,K$}
3805  /// \f]
3806  /// where \f$ N \f$ is the number of nodes in the mesh and \f$ K \f$
3807  /// the number of generalised nodal coordinates. The initial shape
3808  /// of the solid body, \f$ {\bf R}^{(IC)},\f$ and its time-derivatives
3809  /// are specified via the \c GeomObject that is stored in the
3810  /// \c SolidFiniteElement::SolidInitialCondition object. The latter also
3811  /// stores the order of the time-derivative \f$ D \f$ to be assigned.
3813  DenseMatrix<double> &jacobian,
3814  const unsigned& flag);
3815 
3816  /// \short Set the number of types required to interpolate the
3817  /// Lagrangian coordinates
3818  void set_nnodal_lagrangian_type(const unsigned &nlagrangian_type)
3819  {Nnodal_lagrangian_type = nlagrangian_type;}
3820 
3821  /// Pointer to the element's "undeformed" macro element (NULL by default)
3823 
3824  /// \short Calculate the mapping from local to lagrangian coordinates,
3825  /// given the derivatives of the shape functions w.r.t. local coorindates.
3826  /// Return the determinant of the jacobian, the jacobian and inverse jacobian
3828  const DShape &dpsids, DenseMatrix<double> &jacobian,
3829  DenseMatrix<double> &inverse_jacobian) const
3830  {
3831  //Assemble the jacobian
3832  assemble_local_to_lagrangian_jacobian(dpsids,jacobian);
3833  //Invert the jacobian
3834  return invert_jacobian_mapping(jacobian,inverse_jacobian);
3835  }
3836 
3837  /// \short Calculate the mapping from local to lagrangian coordinates,
3838  /// given the derivatives of the shape functions w.r.t. local coordinates,
3839  /// Return only the determinant of the jacobian and the inverse of the
3840  /// mapping (ds/dx)
3841  double local_to_lagrangian_mapping(const DShape &dpsids,
3842  DenseMatrix<double> &inverse_jacobian) const
3843  {
3844  //Find the dimension of the element
3845  unsigned el_dim = dim();
3846  //Assign memory for the jacobian
3847  DenseMatrix<double> jacobian(el_dim);
3848  //Calculate the jacobian and inverse
3849  return local_to_lagrangian_mapping(dpsids,jacobian,inverse_jacobian);
3850  }
3851 
3852  /// \short Calculate the mapping from local to Lagrangian coordinates given
3853  /// the derivatives of the shape functions w.r.t the local coorindates.
3854  /// assuming that the coordinates are aligned in the direction of the local
3855  /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
3856  /// This function returns the determinant of the jacobian, the jacobian
3857  /// and the inverse jacobian.
3858  virtual double local_to_lagrangian_mapping_diagonal(
3859  const DShape &dpsids,DenseMatrix<double> &jacobian,
3860  DenseMatrix<double> &inverse_jacobian) const;
3861 
3862 
3863  /// \short Assigns local equation numbers for the generic solid
3864  /// local equation numbering schemes.
3865  /// If the boolean flag is true the the degrees of freedom are stored in
3866  /// Dof_pt
3867  virtual void assign_solid_local_eqn_numbers(const bool &store_local_dof);
3868 
3869  /// \short Classifies dofs locally for solid specific aspects.
3870  void describe_solid_local_dofs(std::ostream& out,
3871  const std::string& current_string) const;
3872 
3873  /// \short Pointer to object that specifies the initial condition
3875 
3876  public:
3877 
3878  /// \short Access function that returns the local equation number that
3879  /// corresponds to the j-th coordinate of the k-th position-type at the
3880  /// n-th local node.
3881  inline int position_local_eqn(const unsigned &n, const unsigned &k,
3882  const unsigned &j) const
3883  {
3884 #ifdef RANGE_CHECKING
3885  std::ostringstream error_message;
3886  bool error=false;
3887  if(n >= nnode())
3888  {
3889  error = true;
3890  error_message << "Range Error: Nodal number " << n
3891  << " is not in the range (0,"
3892  << nnode() << ")";
3893  }
3894 
3895  if(k >= nnodal_position_type())
3896  {
3897  error = true;
3898  error_message << "Range Error: Position type " << k
3899  << " is not in the range (0,"
3900  << nnodal_position_type() << ")";
3901  }
3902 
3903  if(j >= nodal_dimension())
3904  {
3905  error=true;
3906  error_message << "Range Error: Nodal coordinate " << j
3907  << " is not in the range (0,"
3908  << nodal_dimension() << ")";
3909  }
3910 
3911  if(error)
3912  {
3913  //Throw the error
3914  throw OomphLibError(error_message.str(),
3915  OOMPH_CURRENT_FUNCTION,
3916  OOMPH_EXCEPTION_LOCATION);
3917  }
3918 #endif
3919 
3920  //Return the value
3922  + j];
3923  }
3924 
3925  protected:
3926 
3927 
3928  /// \short Overload the fill_in_contribution_to_jacobian() function to
3929  /// use finite
3930  /// differences to calculate the solid residuals by default
3932  DenseMatrix<double> &jacobian)
3933  {
3934  //Add the contribution to the residuals
3936 
3937  //Solve for the consistent acceleration in Newmark scheme?
3939  {
3941  return;
3942  }
3943 
3944  //Allocate storage for the full residuals (residuals of entire element)
3945  unsigned n_dof = ndof();
3946  Vector<double> full_residuals(n_dof);
3947  //Get the residuals for the entire element
3948  get_residuals(full_residuals);
3949  //Get the solid entries in the jacobian using finite differences
3950  fill_in_jacobian_from_solid_position_by_fd(full_residuals,jacobian);
3951  //There could be internal data
3952  //(finite-difference the lot by default)
3953  fill_in_jacobian_from_internal_by_fd(full_residuals,jacobian,true);
3954  //There could also be external data
3955  //(finite-difference the lot by default)
3956  fill_in_jacobian_from_external_by_fd(full_residuals,jacobian,true);
3957  //There could also be nodal data
3958  fill_in_jacobian_from_nodal_by_fd(full_residuals,jacobian);
3959  }
3960 
3961 
3962 
3963  /// \short Use finite differences to calculate the Jacobian entries
3964  /// corresponding to the solid positions. This version assumes
3965  /// that the residuals vector has already been computed
3966 virtual void
3968  DenseMatrix<double> &jacobian);
3969 
3970 
3971  /// \short Use finite differences to calculate the Jacobian entries
3972  /// corresponding to the solid positions.
3974  {
3975  //Allocate storage for a residuals vector and initialise to zero
3976  unsigned n_dof = ndof();
3977  Vector<double> residuals(n_dof,0.0);
3978  //Get the residuals for the entire element
3979  get_residuals(residuals);
3980  //Call the jacobian calculation
3981  fill_in_jacobian_from_solid_position_by_fd(residuals,jacobian);
3982  }
3983 
3984  /// \short Function that is called before the finite differencing of
3985  /// any solid position data. This may be overloaded to update any slaved
3986  /// data before finite differencing takes place.
3987  virtual inline void update_before_solid_position_fd() { }
3988 
3989  /// \short Function that is call after the finite differencing of
3990  /// the solid position data. This may be overloaded to reset any slaved
3991  /// variables that may have changed during the finite differencing.
3992  virtual inline void reset_after_solid_position_fd() { }
3993 
3994  /// \short Function called within the finite difference loop for
3995  /// the solid position dat after a change in any values in the n-th
3996  /// node.
3997  virtual inline void update_in_solid_position_fd(const unsigned &i) { }
3998 
3999  /// \short Function called within the finite difference loop for
4000  /// solid position data after the values in the i-th node
4001  /// are reset. The default behaviour is to call the update function.
4002  virtual inline void reset_in_solid_position_fd(const unsigned &i)
4004 
4005 
4006 private:
4007 
4008  /// \short Assemble the jacobian matrix for the mapping from local
4009  /// to lagrangian coordinates, given the derivatives of the shape function
4011  const DShape &dpsids, DenseMatrix<double> &jacobian) const;
4012 
4013 
4014  /// \short Assemble the the "jacobian" matrix of second derivatives, given
4015  /// the second derivatives of the shape functions w.r.t. local coordinates
4017  const DShape &d2psids, DenseMatrix<double> &jacobian2) const;
4018 
4019  /// \short Pointer to function that computes the "multiplier" for the
4020  /// inertia terms in the consistent determination of the initial
4021  /// conditions for Newmark timestepping.
4023 
4024  /// \short Array to hold the local equation number information for the
4025  /// solid equations, whatever they may be.
4026  // ALH: (This is here so that the numbering can be done automatically)
4028 
4029 /// \short The Lagrangian dimension of the nodes stored in the element,
4030 //// i.e. the number of Lagrangian coordinates
4032 
4033  /// \short The number of coordinate types requried to intepolate the
4034  /// Lagrangian coordinates in the element. For Lagrange elements
4035  /// it is 1 (the default). It must be over-ridden by using
4036  /// the set_nlagrangian_type() function in the constructors of elements
4037  /// that use generalised coordinate, e.g. for 1D Hermite elements
4038  /// Nnodal_position_types =2.
4040 
4041 protected:
4042 
4043 /// \short Flag to indicate which system of equations to solve
4044 /// when assigning initial conditions for time-dependent problems.
4045 /// If true, solve for the history value that
4046 /// corresponds to the acceleration in the Newmark scheme by demanding
4047 /// that the PDE is satisifed at the initial time. In this case
4048 /// the Jacobian is replaced by the mass matrix.
4050 
4051 private:
4052 
4053  /// \short Access to the "multiplier" for the
4054  /// inertia terms in the consistent determination of the initial
4055  /// conditions for Newmark timestepping.
4056  inline double multiplier(const Vector<double>& xi)
4057  {
4058  //If no function has been set, return 1
4059  if(Multiplier_fct_pt==0)
4060  {
4061  return 1.0;
4062  }
4063  else
4064  {
4065  // Evaluate function pointer
4066  return (*Multiplier_fct_pt)(xi);
4067  }
4068  }
4069 
4070 
4071 };
4072 
4073 
4074 
4075 
4076 //========================================================================
4077 /// FaceElements are elements that coincide with the faces of
4078 /// higher-dimensional "bulk" elements. They are used on
4079 /// boundaries where additional non-trivial boundary conditions need to be
4080 /// applied. Examples include free surfaces, and applied traction conditions.
4081 /// In many cases, FaceElements need to evaluate to quantities
4082 /// in the associated bulk elements. For instance, the evaluation
4083 /// of a shear stresses on 2D FaceElement requires the evaluation
4084 /// of velocity derivatives in the associated 3D volume element etc.
4085 /// Therefore we store a pointer to the associated bulk element,
4086 /// and information about the relation between the local
4087 /// coordinates in the face and bulk elements.
4088 //========================================================================
4089 class FaceElement: public virtual FiniteElement
4090 {
4091  /// \short Typedef for the function that translates the face coordinate
4092  /// to the coordinate in the bulk element
4093  typedef void (*CoordinateMappingFctPt)(const Vector<double> &s,
4094  Vector<double> &s_bulk);
4095 
4096  /// \short Typedef for the function that returns the partial derivative
4097  /// of the local coordinates in the bulk element
4098  /// with respect to the coordinates along the face.
4099  /// In addition this function returns an index of one of the
4100  /// bulk local coordinates that varies away from the edge
4102  const Vector<double> &s, DenseMatrix<double> &ds_bulk_dsface,
4103  unsigned &interior_direction);
4104 
4105  private:
4106 
4107  /// \short Pointer to a function that translates the face coordinate
4108  /// to the coordinate in the bulk element
4110 
4111  /// \short Pointer to a function that returns the derivatives of the local
4112  /// "bulk" coordinates with respect to the local face coordinates.
4114 
4115  /// \short Vector holding integers to translate additional position types
4116  /// to those of the "bulk" element; e.g. Bulk_position_type(1) is
4117  /// the position type of the "bulk" element that corresponds to face position
4118  /// type 1. This is required in QHermiteElements, where the slope in the
4119  /// direction of the 1D element is face position type 1, but may be
4120  /// position type 1 or 2 in the "bulk" element, depending upon which face,
4121  /// we are located.
4123 
4124  /// \short Sign of outer unit normal (relative to cross-products of tangent
4125  /// vectors in the corresponding "bulk" element.
4127 
4128  /// \short Index of the face
4130 
4131  /// \short Ignores the warning when the tangent vectors from
4132  /// continuous_tangent_and_outer_unit_normal may not be continuous
4133  /// as a result of the unit normal being aligned with the z axis.
4134  /// This can be avoided by supplying a general direction vector for
4135  /// the tangent vector via set_tangent_direction(...).
4137 
4138  protected:
4139 
4140  /// The boundary number in the bulk mesh to which this element is attached
4142 
4143 #ifdef PARANOID
4144 
4145  /// \short Has the Boundary_number_in_bulk_mesh been set? Only included if
4146  /// compiled with PARANOID switched on.
4148 
4149 #endif
4150 
4151  /// \short Pointer to the associated higher-dimensional "bulk" element
4153 
4154  /// \short List of indices of the local node numbers in the "bulk" element
4155  /// that correspond to the local node numbers in the FaceElement
4157 
4158  /// \short A vector that will hold the number of data values at the nodes that
4159  /// are associated with the "bulk" element. i.e. not including any additional
4160  /// degrees of freedom that might be required for extra equations
4161  /// that are being solved by
4162  /// the FaceElement.
4163  // NOTE: This breaks if the nodes have already been resized
4164  // before calling the face element, i.e. two separate sets of equations on the
4165  // face!
4167 
4168  /// \short A general direction pointer for the tangent vectors.
4169  /// This is used in the function continuous_tangent_and_outer_unit_normal()
4170  /// for creating continuous tangent vectors in spatial dimensions.
4171  /// The general direction is projected on to the surface. This technique is
4172  /// not required in two spatial dimensions.
4174 
4175  /// \short Helper function adding additional values for the unknowns
4176  /// associated with the FaceElement. This function also sets the map
4177  /// containing the position of the first entry of this face element's
4178  /// additional values.The inputs are the number of additional values
4179  /// and the face element's ID. Note the same number of additonal values
4180  /// are allocated at ALL nodes.
4181  void add_additional_values(const Vector<unsigned> &nadditional_values,
4182  const unsigned &id)
4183  {
4184  // How many nodes?
4185  const unsigned n_node=nnode();
4186 
4187  // loop over the nodes
4188  for (unsigned n=0;n<n_node;n++)
4189  {
4190  //Assign the required number of additional nodes to the node
4191  dynamic_cast<BoundaryNodeBase*>(
4192  this->node_pt(n))->assign_additional_values_with_face_id(
4193  nadditional_values[n],id);
4194  }
4195  }
4196 
4197 
4198  public:
4199 
4200  /// Constructor: Initialise all appropriate member data
4204  {
4205  // Check whether things have been set
4206 #ifdef PARANOID
4208 #endif
4209 
4210  // Bulk_position_type[0] is always 0 (the position)
4211  Bulk_position_type.push_back(0);
4212  }
4213 
4214  /// Empty virtual destructor
4215  virtual ~FaceElement() {}
4216 
4217 
4218  /// Broken copy constructor
4220  {
4221  BrokenCopy::broken_copy("FaceElement");
4222  }
4223 
4224  /// Broken assignment operator
4225  /*void operator=(const FaceElement&)
4226  {
4227  BrokenCopy::broken_assign("FaceElement");
4228  }*/
4229 
4230  /// Access function for the boundary number in bulk mesh
4231  inline const unsigned& boundary_number_in_bulk_mesh() const
4233 
4234 
4235  /// Set function for the boundary number in bulk mesh
4236  inline void set_boundary_number_in_bulk_mesh(const unsigned& b)
4237  {
4239 #ifdef PARANOID
4241 #endif
4242  }
4243 
4244  /// \short In a FaceElement, the "global" intrinsic coordinate
4245  /// of the element along the boundary, when viewed as part of
4246  /// a compound geometric object is specified using the
4247  /// boundary coordinate defined by the mesh.
4248  /// Note: Boundary coordinates will have been set up when
4249  /// creating the underlying mesh, and their values will have
4250  /// been stored at the nodes.
4251  double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i)
4252  const
4253  {
4254  //Vector in which to hold the intrinsic coordinate
4255  Vector<double> zeta(this->dim());
4256 
4257  //Get the k-th generalised boundary coordinate at node n
4260 
4261  //Return the individual coordinate
4262  return zeta[i];
4263  }
4264 
4265  /// \short Return the Jacobian of mapping from local to global
4266  /// coordinates at local position s.
4267  /// Overloaded from FiniteElement.
4268  double J_eulerian(const Vector<double> &s) const;
4269 
4270  /// \short Return the Jacobian of the mapping from local to global
4271  /// coordinates at the ipt-th integration point
4272  /// Overloaded from FiniteElement.
4273  double J_eulerian_at_knot(const unsigned &ipt) const;
4274 
4275  /// \short Return FE interpolated coordinate x[i] at local coordinate s.
4276  /// Overloaded to get information from bulk.
4277  double interpolated_x(const Vector<double> &s, const unsigned &i) const
4278  {
4279  // Local coordinates in bulk element
4280  Vector<double> s_bulk(dim()+1);
4281  s_bulk=local_coordinate_in_bulk(s);
4282 
4283  // Return Eulerian coordinate as computed by bulk
4284  return bulk_element_pt()->interpolated_x(s_bulk,i);
4285  }
4286 
4287  ///\short Return FE interpolated coordinate x[i] at local coordinate s
4288  ///at previous timestep t (t=0: present; t>0: previous timestep).
4289  /// Overloaded to get information from bulk.
4290  double interpolated_x(const unsigned& t, const Vector<double> &s,
4291  const unsigned &i) const
4292  {
4293  // Local coordinates in bulk element
4294  Vector<double> s_bulk(dim()+1);
4295  s_bulk=local_coordinate_in_bulk(s);
4296 
4297  // Return Eulerian coordinate as computed by bulk
4298  return bulk_element_pt()->interpolated_x(t,s_bulk,i);
4299  }
4300 
4301  /// \short Return FE interpolated position x[] at local coordinate s as Vector
4302  /// Overloaded to get information from bulk.
4304  {
4305  // Local coordinates in bulk element
4306  Vector<double> s_bulk(dim()+1);
4307  s_bulk=local_coordinate_in_bulk(s);
4308 
4309  // Get Eulerian position vector
4310  bulk_element_pt()->interpolated_x(s_bulk,x);
4311  }
4312 
4313  /// \short Return FE interpolated position x[] at local coordinate s
4314  /// at previous timestep t as Vector (t=0: present; t>0: previous timestep).
4315  /// Overloaded to get information from bulk.
4316  void interpolated_x(const unsigned& t, const Vector<double> &s,
4317  Vector<double>& x) const
4318  {
4319  // Local coordinates in bulk element
4320  Vector<double> s_bulk(dim()+1);
4321  s_bulk=local_coordinate_in_bulk(s);
4322 
4323  // Get Eulerian position vector
4324  bulk_element_pt()->interpolated_x(t,s_bulk,x);
4325  }
4326 
4327  /// \short Return t-th time-derivative of the
4328  /// i-th FE-interpolated Eulerian coordinate at
4329  /// local coordinate s. Overloaded to get information from bulk.
4331  const unsigned &i,
4332  const unsigned &t)
4333  {
4334  // Local coordinates in bulk element
4335  Vector<double> s_bulk(dim()+1);
4336  s_bulk=local_coordinate_in_bulk(s);
4337 
4338  // Return Eulerian coordinate as computed by bulk
4339  return bulk_element_pt()->interpolated_dxdt(s_bulk,i,t);
4340  }
4341 
4342  /// \short Compte t-th time-derivative of the
4343  /// FE-interpolated Eulerian coordinate vector at
4344  /// local coordinate s. Overloaded to get information from bulk.
4346  const unsigned &t,
4347  Vector<double>& dxdt)
4348  {
4349  // Local coordinates in bulk element
4350  Vector<double> s_bulk(dim()+1);
4351  s_bulk=local_coordinate_in_bulk(s);
4352 
4353  // Get Eulerian position vector
4354  bulk_element_pt()->interpolated_dxdt(s_bulk,t,dxdt);
4355  }
4356 
4357  /// \short Sign of outer unit normal (relative to cross-products of tangent
4358  /// vectors in the corresponding "bulk" element.
4359  int& normal_sign(){return Normal_sign;}
4360 
4361  /// \short Return sign of outer unit normal (relative to cross-products
4362  /// of tangent vectors in the corresponding "bulk" element. (const version)
4363  int normal_sign() const {return Normal_sign;}
4364 
4365  /// \short Index of the face (a number that uniquely identifies the face
4366  /// in the element)
4367  int& face_index(){return Face_index;}
4368 
4369  /// \short Index of the face (a number that uniquely identifies the face
4370  /// in the element) (const version)
4371  int face_index() const {return Face_index;}
4372 
4373  /// \short Public access function for the tangent direction pointer.
4375  {
4376  return Tangent_direction_pt;
4377  }
4378 
4379  /// \short Set the tangent direction vector.
4381  {
4382 #ifdef PARANOID
4383  // Check that tangent_direction_pt is not null.
4384  if(tangent_direction_pt == 0)
4385  {
4386  std::ostringstream error_message;
4387  error_message
4388  << "The pointer tangent_direction_pt is null.\n";
4389  throw OomphLibError(error_message.str(),
4390  OOMPH_CURRENT_FUNCTION,
4391  OOMPH_EXCEPTION_LOCATION);
4392  }
4393 
4394  // Check that the vector is the correct size.
4395  // The size of the tangent vector.
4396  unsigned tang_dir_size = tangent_direction_pt->size();
4397  unsigned spatial_dimension = this->nodal_dimension();
4398  if(tang_dir_size != spatial_dimension)
4399  {
4400  std::ostringstream error_message;
4401  error_message
4402  << "The tangent direction vector has size " << tang_dir_size << "\n"
4403  << "but this element has spatial dimension " << spatial_dimension
4404  << ".\n";
4405  throw OomphLibError(error_message.str(),
4406  OOMPH_CURRENT_FUNCTION,
4407  OOMPH_EXCEPTION_LOCATION);
4408  }
4409 
4410  if(tang_dir_size == 2)
4411  {
4412  std::ostringstream warning_message;
4413  warning_message
4414  << "The spatial dimension is " << spatial_dimension << ".\n"
4415  << "I do not need a tangent direction vector to create \n"
4416  << "continuous tangent vectors in two spatial dimensions.";
4417  OomphLibWarning(warning_message.str(),
4418  OOMPH_CURRENT_FUNCTION,
4419  OOMPH_EXCEPTION_LOCATION);
4420 
4421  }
4422 #endif
4423 
4424  // Set the direction vector for the tangent.
4426  }
4427 
4428  /// \short Turn on warning for when there may be discontinuous tangent vectors
4429  /// from continuous_tangent_and_outer_unit_normal(...)
4431  {
4433  }
4434 
4435  /// \short Turn off warning for when there may be discontinuous tangent
4436  /// vectors from continuous_tangent_and_outer_unit_normal(...)
4438  {
4440  }
4441 
4442  /// \short Compute the tangent vector(s) and the outer unit normal
4443  /// vector at the specified local coordinate.
4444  /// In two spatial dimensions, a "tangent direction" is not required.
4445  /// In three spatial dimensions, a tangent direction is required
4446  /// (set via set_tangent_direction(...)), and we project the tanent direction
4447  /// on to the surface. The second tangent vector is taken to be the cross
4448  /// product of the projection and the unit normal.
4450  (const Vector<double> &s, Vector<Vector<double> > &tang_vec,
4451  Vector<double> &unit_normal) const;
4452 
4453  /// \short Compute the tangent vector(s) and the outer unit normal
4454  /// vector at the ipt-th integration point. This is a wrapper around
4455  /// continuous_tangent_and_outer_unit_normal(...) with the integration points
4456  /// converted into local coordinates.
4458  (const unsigned &ipt, Vector<Vector<double> > &tang_vec,
4459  Vector<double> &unit_normal) const;
4460 
4461  /// \short Compute outer unit normal at the specified local coordinate
4462  void outer_unit_normal(const Vector<double> &s,
4463  Vector<double> &unit_normal) const;
4464 
4465  /// \short Compute outer unit normal at ipt-th integration point
4466  void outer_unit_normal(const unsigned &ipt,
4467  Vector<double> &unit_normal) const;
4468 
4469  /// Pointer to higher-dimensional "bulk" element
4471 
4472 
4473  /// Pointer to higher-dimensional "bulk" element (const version)
4475 
4476  //Clang specific pragma's
4477 #ifdef __clang__
4478 #pragma clang diagnostic push
4479 #pragma clang diagnostic ignored "-Woverloaded-virtual"
4480 #endif
4481 
4482  /// \short Return the pointer to the function that maps the face
4483  /// coordinate to the bulk coordinate
4486 
4487  /// \short Return the pointer to the function that maps the face
4488  /// coordinate to the bulk coordinate (const version)
4491 
4492 
4493  /// \short Return the pointer to the function that returns the derivatives
4494  /// of the bulk coordinates wrt the face coordinates
4497 
4498  /// \short Return the pointer to the function that returns the derivatives
4499  /// of the bulk coordinates wrt the face coordinates (const version)
4502 
4503 #ifdef __clang__
4504 #pragma clang diagnostic pop
4505 #endif
4506 
4507  /// \short Return vector of local coordinates in bulk element,
4508  /// given the local coordinates in this FaceElement
4510 
4511  /// \short Calculate the vector of local coordinate in the bulk element
4512  /// given the local coordinates in this FaceElement
4514  Vector<double> &s_bulk) const;
4515 
4516  /// \short Calculate the derivatives of the local coordinates in the
4517  /// bulk element with respect to the local coordinates in this FaceElement.
4518  /// In addition return the index of a bulk local coordinate that varies away
4519  /// from the face.
4520  void get_ds_bulk_ds_face(const Vector<double> &s,
4521  DenseMatrix<double> &dsbulk_dsface,
4522  unsigned &interior_direction) const;
4523 
4524  /// \short Return the position type in the "bulk" element that corresponds
4525  /// to position type i on the FaceElement.
4526  unsigned& bulk_position_type(const unsigned &i)
4527  {return Bulk_position_type[i];}
4528 
4529  /// \short Return the position type in the "bulk" element that corresponds
4530  /// to the position type i on the FaceElement. Const version
4531  const unsigned &bulk_position_type(const unsigned &i) const
4532  {return Bulk_position_type[i];}
4533 
4534  /// \short Resize the storage for the bulk node numbers
4535  void bulk_node_number_resize(const unsigned &i) {Bulk_node_number.resize(i);}
4536 
4537  /// \short Return the bulk node number that corresponds to the n-th
4538  /// local node number
4539  unsigned &bulk_node_number(const unsigned &n) {return Bulk_node_number[n];}
4540 
4541  /// \short Return the bulk node number that corresponds to the n-th
4542  /// local node number (const version)
4543  const unsigned &bulk_node_number(const unsigned &n) const
4544  {return Bulk_node_number[n];}
4545 
4546  /// \short Resize the storage for bulk_position_type to i entries.
4547  void bulk_position_type_resize(const unsigned &i)
4548  {Bulk_position_type.resize(i);}
4549 
4550  /// \short Return the number of values originally stored at local node n
4551  /// (before the FaceElement added additional values to it (if it did))
4552  unsigned& nbulk_value(const unsigned &n) {return Nbulk_value[n];}
4553 
4554  /// \short Return the number of values originally stored at local node n
4555  /// (before the FaceElement added additional values to it (if it did))
4556  /// (const version)
4557  unsigned nbulk_value(const unsigned &n) const
4558  {return Nbulk_value[n];}
4559 
4560  /// \short Resize the storage for the number of values originally
4561  /// stored at the local nodes to i entries.
4562  void nbulk_value_resize(const unsigned &i) {Nbulk_value.resize(i);}
4563 
4564  /// \short Provide additional storage for a specified number of
4565  /// values at the nodes of the FaceElement. (This is needed, for
4566  /// instance, in free-surface elements, if the non-penetration
4567  /// condition is imposed by Lagrange multipliers whose values
4568  /// are only stored at the surface nodes but not in the interior of
4569  /// the bulk element). \c nadditional_data_values[n] specifies the
4570  /// number of additional values required at node \c n of the
4571  /// FaceElement. \b Note: Since this function is executed separately
4572  /// for each FaceElement, nodes that are common to multiple elements
4573  /// might be resized repeatedly. To avoid this, we only allow
4574  /// a single resize operation by comparing the number of values stored
4575  /// at each node to the number of values the node had when it
4576  /// was simply a member of the associated "bulk"
4577  /// element. <b>There are cases where this will break!</b> -- e.g. if
4578  /// a node is common to two FaceElements which require
4579  /// additional storage for distinct quantities. Such cases
4580  /// need to be handled by "hand-crafted" face elements.
4581  void resize_nodes(Vector<unsigned> &nadditional_data_values)
4582  {
4583  //Locally cache the number of node
4584  unsigned n_node = nnode();
4585  //Resize the storage for values at the nodes:
4586  for(unsigned l=0;l<n_node;l++)
4587  {
4588  //Find number of values stored at the node
4589  unsigned Initial_Nvalue = node_pt(l)->nvalue();
4590  //Read out the number of additional values
4591  unsigned Nadditional = nadditional_data_values[l];
4592  //If the node has not already been resized, resize it
4593  if((Initial_Nvalue == Nbulk_value[l]) && (Nadditional > 0))
4594  {
4595  //Resize the node according to the number of additional values
4596  node_pt(l)->resize(Nbulk_value[l]+Nadditional);
4597  }
4598  } //End of loop over nodes
4599  }
4600 
4601  /// Output boundary coordinate zeta
4602  void output_zeta(std::ostream &outfile, const unsigned& nplot);
4603 
4604 };
4605 
4606 
4607 
4608 //========================================================================
4609 /// SolidFaceElements combine FaceElements and SolidFiniteElements
4610 /// and overload various functions so they work properly in the
4611 /// FaceElement context
4612 //========================================================================
4613 class SolidFaceElement: public virtual FaceElement,
4614  public virtual SolidFiniteElement
4615 {
4616 
4617  public:
4618 
4619  /// \short The "global" intrinsic coordinate of the element when
4620  /// viewed as part of a geometric object should be given by
4621  /// the FaceElement representation, by default
4622  /// This final over-ride is required because both SolidFiniteElements
4623  /// and FaceElements overload zeta_nodal
4624  double zeta_nodal(const unsigned &n, const unsigned &k,
4625  const unsigned &i) const
4626  {return FaceElement::zeta_nodal(n,k,i);}
4627 
4628 
4629 
4630  /// \short Return i-th FE-interpolated Lagrangian coordinate xi[i] at
4631  /// local coordinate s. Overloaded from SolidFiniteElement. Note that
4632  /// the Lagrangian coordinates are those defined in the bulk!
4633  /// For instance, in a 1D FaceElement that is aligned with
4634  /// the Lagrangian coordinate line xi_0=const, only xi_1 will vary
4635  /// in the FaceElement. This may confuse you if you (wrongly!) believe that
4636  /// in a 1D SolidElement there should only a single Lagrangian
4637  /// coordinate, namely xi_0!
4639  const unsigned &i) const
4640  {
4641  // Local coordinates in bulk element
4642  Vector<double> s_bulk(dim()+1);
4643  s_bulk=local_coordinate_in_bulk(s);
4644 
4645  // Return Lagrangian coordinate as computed by bulk
4646  return dynamic_cast<SolidFiniteElement*>(bulk_element_pt())->
4647  interpolated_xi(s_bulk,i);
4648  }
4649 
4650 
4651  /// \short Compute FE interpolated Lagrangian coordinate vector xi[] at
4652  /// local coordinate s as Vector. Overloaded from SolidFiniteElement. Note
4653  /// that the Lagrangian coordinates are those defined in the bulk!
4654  /// For instance, in a 1D FaceElement that is aligned with
4655  /// the Lagrangian coordinate line xi_0=const, only xi_1 will vary
4656  /// in the FaceElement. This may confuse you if you (wrongly!) believe that
4657  /// in a 1D SolidElement there should only a single Lagrangian
4658  /// coordinate, namely xi_0!
4660  Vector<double>& xi) const
4661  {
4662  // Local coordinates in bulk element
4663  Vector<double> s_bulk(dim()+1);
4664  s_bulk=local_coordinate_in_bulk(s);
4665 
4666  // Get Lagrangian position vector
4667  dynamic_cast<SolidFiniteElement*>(bulk_element_pt())->
4668  interpolated_xi(s_bulk,xi);
4669  }
4670 
4671 
4672 };
4673 
4674 
4675 
4676 ///////////////////////////////////////////////////////////////////////
4677 ///////////////////////////////////////////////////////////////////////
4678 ///////////////////////////////////////////////////////////////////////
4679 
4680 
4681 
4682 //=======================================================================
4683 /// Solid point element
4684 //=======================================================================
4685  class SolidPointElement : public virtual SolidFiniteElement,
4686  public virtual PointElement
4687 {
4688 };
4689 
4690 
4691 ///////////////////////////////////////////////////////////////////////
4692 ///////////////////////////////////////////////////////////////////////
4693 ///////////////////////////////////////////////////////////////////////
4694 
4695 
4696 
4697 
4698 
4699 //=======================================================================
4700 /// FaceGeometry class definition: This policy class is used to allow
4701 /// construction of face elements that solve arbitrary equations without
4702 /// having to tamper with the corresponding "bulk" elements.
4703 /// The geometrical information for the face element must be specified
4704 /// by each "bulk" element using an explicit specialisation of this class.
4705 //=======================================================================
4706 template<class ELEMENT>
4708 {
4709 
4710 };
4711 
4712 
4713 ///////////////////////////////////////////////////////////////////////
4714 ///////////////////////////////////////////////////////////////////////
4715 ///////////////////////////////////////////////////////////////////////
4716 
4717 
4718 
4719 //======================================================================
4720 /// Dummy FaceElement for use with purely geometric operations
4721 /// such as mesh generation.
4722 //======================================================================
4723 template <class ELEMENT>
4724 class DummyFaceElement : public virtual FaceGeometry<ELEMENT>,
4725 public virtual FaceElement
4726 {
4727 
4728 public:
4729 
4730  /// \short Constructor, which takes a "bulk" element and the
4731  /// face index
4732  DummyFaceElement(FiniteElement* const &element_pt,
4733  const int &face_index) :
4734  FaceGeometry<ELEMENT>(), FaceElement()
4735  {
4736  //Attach the geometrical information to the element. N.B. This function
4737  //also assigns nbulk_value from the required_nvalue of the bulk element
4738  element_pt->build_face_element(face_index,this);
4739  }
4740 
4741 
4742  /// \short Constructor
4744  FaceGeometry<ELEMENT>(), FaceElement()
4745  {
4746  }
4747 
4748 
4749  /// \short The "global" intrinsic coordinate of the element when
4750  /// viewed as part of a geometric object should be given by
4751  /// the FaceElement representation, by default
4752  double zeta_nodal(const unsigned &n, const unsigned &k,
4753  const unsigned &i) const
4754  {return FaceElement::zeta_nodal(n,k,i);}
4755 
4756  /// Output nodal coordinates
4757  void output(std::ostream &outfile)
4758  {
4759  outfile << "ZONE" << std::endl;
4760  unsigned nnod=nnode();
4761  for (unsigned j=0;j<nnod;j++)
4762  {
4763  Node* nod_pt=node_pt(j);
4764  unsigned dim=nod_pt->ndim();
4765  for (unsigned i=0;i<dim;i++)
4766  {
4767  outfile << nod_pt->x(i) << " ";
4768  }
4769  outfile << std::endl;
4770  }
4771  }
4772 
4773  /// Output at n_plot points
4774  void output(std::ostream &outfile, const unsigned &n_plot)
4775  {FiniteElement::output(outfile,n_plot);}
4776 
4777  /// C-style output
4778  void output(FILE* file_pt)
4779  {FiniteElement::output(file_pt);}
4780 
4781  /// C_style output at n_plot points
4782  void output(FILE* file_pt, const unsigned &n_plot)
4783  {FiniteElement::output(file_pt,n_plot);}
4784 
4785 };
4786 
4787 ///////////////////////////////////////////////////////////////////////
4788 ///////////////////////////////////////////////////////////////////////
4789 ///////////////////////////////////////////////////////////////////////
4790 
4791 
4792 //=========================================================================
4793 /// Base class for elements that can specify a drag and torque
4794 /// (about the origin) -- typically used for immersed particle
4795 /// computations.
4796 //=========================================================================
4798 {
4799 
4800 public:
4801 
4802  /// Empty constructor
4804 
4805  /// Empty virtual destructor
4807 
4808  /// \short Function that specifies the drag force and the torque about
4809  /// the origin
4810  virtual void get_drag_and_torque(Vector<double>& drag_force,
4811  Vector<double>& drag_torque)=0;
4812 
4813 };
4814 
4815 
4816 ///////////////////////////////////////////////////////////////////////
4817 ///////////////////////////////////////////////////////////////////////
4818 ///////////////////////////////////////////////////////////////////////
4819 
4820 
4821 //======================================================================
4822 /// Basic-ified FaceElement, without any of the functionality of
4823 /// of actual FaceElements -- it's just a surface element of the
4824 /// same geometric type as the FaceGeometry associated with
4825 /// bulk element specified by the template parameter. The element
4826 /// can be used to represent boundaries without actually being
4827 /// attached to a bulk element. Used mainly during unstructured
4828 /// mesh generation.
4829 //======================================================================
4830 template <class ELEMENT>
4831 class FreeStandingFaceElement : public virtual FaceGeometry<ELEMENT>
4832 {
4833 
4834 public:
4835 
4836  /// \short Constructor
4839  {
4840 
4841  //Check whether things have been set
4842 #ifdef PARANOID
4844 #endif
4845 
4846  }
4847 
4848  /// Access function for the boundary number in bulk mesh
4849  inline const unsigned& boundary_number_in_bulk_mesh() const
4851 
4852 
4853  /// Set function for the boundary number in bulk mesh
4854  inline void set_boundary_number_in_bulk_mesh(const unsigned& b)
4855  {
4857 #ifdef PARANOID
4859 #endif
4860  }
4861 
4862  /// \short In a FaceElement, the "global" intrinsic coordinate
4863  /// of the element along the boundary, when viewed as part of
4864  /// a compound geometric object is specified using the
4865  /// boundary coordinate defined by the mesh.
4866  /// Note: Boundary coordinates will have been set up when
4867  /// creating the underlying mesh, and their values will have
4868  /// been stored at the nodes.
4869  double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i)
4870  const
4871  {
4872  //Vector in which to hold the intrinsic coordinate
4873  Vector<double> zeta(this->dim());
4874 
4875  //Get the k-th generalised boundary coordinate at node n
4876  this->node_pt(n)->get_coordinates_on_boundary(
4878 
4879  //Return the individual coordinate
4880  return zeta[i];
4881  }
4882 
4883  protected:
4884 
4885  /// The boundary number in the bulk mesh to which this element is attached
4887 
4888 #ifdef PARANOID
4889 
4890  /// \short Has the Boundary_number_in_bulk_mesh been set? Only included if
4891  /// compiled with PARANOID switched on.
4893 
4894 #endif
4895 
4896 };
4897 
4898 
4899 
4900 ///////////////////////////////////////////////////////////////////////
4901 ///////////////////////////////////////////////////////////////////////
4902 ///////////////////////////////////////////////////////////////////////
4903 
4904 
4905 //======================================================================
4906 /// Pure virtual base class for elements that can be used with
4907 /// PressureBasedSolidLSCPreconditioner.
4908 //======================================================================
4910 {
4911 
4912  public:
4913 
4914  /// Empty constructor
4916 
4917  /// Virtual destructor
4919 
4920  /// Broken copy constructor
4923  {
4925  "SolidElementWithDiagonalMassMatrix");
4926  }
4927 
4928  /// Broken assignment operator
4929  void operator=(const
4931  {
4933  "SolidElementWithDiagonalMassMatrix");
4934  }
4935 
4936  /// \short Get the diagonal of whatever represents the mass matrix
4937  /// in the specific preconditionable element. For Navier-Stokes
4938  /// elements this is the velocity mass matrix; for incompressible solids it's
4939  /// the mass matrix; ...
4940  virtual void get_mass_matrix_diagonal(Vector<double> &mass_diag)=0;
4941 
4942 };
4943 
4944 
4945 
4946 
4947 ///////////////////////////////////////////////////////////////////////
4948 ///////////////////////////////////////////////////////////////////////
4949 ///////////////////////////////////////////////////////////////////////
4950 
4951 
4952 
4953 //======================================================================
4954 /// Pure virtual base class for elements that can be used with
4955 /// Navier-Stokes Schur complement preconditioner and provide the diagonal
4956 /// of their velocity and pressure mass matrices -- needs to be defined here
4957 /// (in generic) because this applies to a variety of Navier-Stokes
4958 /// elements (cartesian, cylindrical polar, ...)
4959 /// that can be preconditioned effectively by the Navier Stokes (!)
4960 /// preconditioners in the (cartesian) Navier-Stokes directory.
4961 //======================================================================
4963 {
4964 
4965  public:
4966 
4967  /// Empty constructor
4969 
4970  /// Virtual destructor
4972 
4973  /// Broken copy constructor
4976  {
4978  "NavierStokesElementWithDiagonalMassMatrices");
4979  }
4980 
4981  /// Broken assignment operator
4982  void operator=(const
4984  {
4986  "NavierStokesElementWithDiagonalMassMatrices");
4987  }
4988 
4989  /// \short Compute the diagonal of the velocity/pressure mass matrices.
4990  /// If which one=0, both are computed, otherwise only the pressure
4991  /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
4992  /// LSC version of the preconditioner only needs that one)
4994  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
4995  const unsigned& which_one=0)=0;
4996 
4997 };
4998 
4999 
5000 
5001 
5002 ///////////////////////////////////////////////////////////////////////
5003 ///////////////////////////////////////////////////////////////////////
5004 ///////////////////////////////////////////////////////////////////////
5005 
5006 } // end of oomph-lib namespace
5007 
5008 #endif
5009 
5010 
5011 
5012 
5013 
5014 
5015 
5016 
unsigned Boundary_number_in_bulk_mesh
The boundary number in the bulk mesh to which this element is attached.
Definition: elements.h:4886
bool external_data_fd(const unsigned &i) const
Return the status of the boolean flag indicating whether the external data is included in the finite ...
Definition: elements.h:759
virtual void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative. This function is empty and simply establishes a common interface for all derived elements that are formulated in ALE form.
Definition: elements.h:2308
A Generalised Element class.
Definition: elements.h:76
PointElement(const PointElement &)
Broken copy constructor.
Definition: elements.h:3192
virtual void get_residuals(Vector< double > &residuals)
Calculate the vector of residuals of the equations in the element. By default initialise the vector t...
Definition: elements.h:975
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4869
virtual void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
Plot the error when compared against a given exact solution . Also calculates the norm of the error a...
Definition: elements.h:3027
virtual void assemble_eulerian_base_vectors(const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
Assemble the covariant Eulerian base vectors, assuming that the derivatives of the shape functions wi...
Definition: elements.cc:1987
virtual void reset_after_internal_fd()
Function that is call after the finite differencing of the internal data. This may be overloaded to r...
Definition: elements.h:460
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data...
Definition: elements.h:268
void operator=(const GeneralisedElement &)
Broken assignment operator.
Definition: elements.h:617
DummyFaceElement()
Constructor.
Definition: elements.h:4743
virtual void point_output_data(const Vector< double > &s, Vector< double > &data)
Virtual function to write the double precision numbers that appear in a single line of output into th...
Definition: elements.h:2678
virtual unsigned nnode_on_face() const
Definition: elements.h:3123
double interpolated_x(const unsigned &t, const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s at previous timestep t (t=0: present; t>...
Definition: elements.h:4290
virtual void compute_abs_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error)
Plot the error when compared against a given exact solution . Also calculates the maximum absolute er...
Definition: elements.h:3066
void clear_global_eqn_numbers()
Clear the storage for the global equation numbers and pointers to dofs (if stored) ...
Definition: elements.h:213
virtual double compute_physical_size() const
Broken virtual function to compute the actual size (taking into account factors such as 2pi or radii ...
Definition: elements.h:2665
void read_internal_data_values_from_vector(const Vector< double > &vector_of_values, unsigned &index)
Read all internal data and time history values from the vector starting from index. On return the index will be set to the value at the end of the data that has been read in.
Definition: elements.cc:612
unsigned lagrangian_dimension() const
Return the number of Lagrangian coordinates that the element requires at all nodes. This is by default the elemental dimension. If we ever need any other case, it can be implemented.
Definition: elements.h:3518
virtual void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
Definition: elements.cc:545
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
virtual void reset_in_nodal_fd(const unsigned &i)
Function called within the finite difference loop for nodal data after the i-th nodal values is reset...
Definition: elements.h:1684
virtual double s_max() const
Max. value of local coordinate.
Definition: elements.h:2642
virtual void get_mass_matrix_diagonal(Vector< double > &mass_diag)=0
Get the diagonal of whatever represents the mass matrix in the specific preconditionable element...
virtual void assign_additional_local_eqn_numbers()
Setup any additional look-up schemes for local equation numbers. Examples of use include using local ...
Definition: elements.h:264
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
void set_must_be_kept_as_halo()
Insist that this element be kept as a halo element during a distribute?
Definition: elements.h:1151
virtual void get_x_and_xi(const Vector< double > &s, Vector< double > &x_fe, Vector< double > &x, Vector< double > &xi_fe, Vector< double > &xi) const
Eulerian and Lagrangian coordinates as function of the local coordinates: The Eulerian position is re...
Definition: elements.h:3410
virtual void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the elemental contribution to the residuals vector. Note that this function will NOT initialise t...
Definition: elements.h:360
virtual double dshape_lagrangian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidxi) const
Return the geometric shape functions and also first derivatives w.r.t. Lagrangian coordinates at ipt-...
Definition: elements.cc:6411
virtual void get_residuals_for_solid_ic(Vector< double > &residuals)
Compute the residuals for the setup of an initial condition. The global equations are: where is the...
Definition: elements.h:3735
void initialise(const T &val)
Initialize all values in the matrix to val.
Definition: matrices.h:523
void dposition_dt(const Vector< double > &zeta, const unsigned &t, Vector< double > &drdt)
Return the t-th time derivative of the parametrised position of the FiniteElement in its GeomObject i...
Definition: elements.h:2543
std::vector< bool > Data_fd
Storage for boolean flags of size Ninternal_data + Nexternal_data that correspond to the data used in...
Definition: elements.h:127
void output_zeta(std::ostream &outfile, const unsigned &nplot)
Output boundary coordinate zeta.
Definition: elements.cc:4991
virtual void assign_all_generic_local_eqn_numbers(const bool &store_local_dof_pt)
Overloaded version of the calculation of the local equation numbers. If the boolean argument is true ...
Definition: elements.h:2087
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
Definition: elements.h:4624
virtual void reset_in_internal_fd(const unsigned &i)
Function called within the finite difference loop for internal data after the values in the i-th exte...
Definition: elements.h:470
Node * construct_boundary_node(const unsigned &n, TimeStepper *const &time_stepper_pt)
Construct the local node n and return a pointer to it, in the case when the node MAY be located on a ...
Definition: elements.h:3586
virtual void assign_all_generic_local_eqn_numbers(const bool &store_local_dof_pt)
Overload assign_all_generic_local_equation numbers to include the data associated with solid dofs...
Definition: elements.h:3608
void add_internal_value_pt_to_map(std::map< unsigned, double * > &map_of_value_pt)
Add pointers to the internal data values to map indexed by the global equation number.
Definition: elements.cc:583
virtual void disable_ALE()
This is an empty function that establishes a uniform interface for all (derived) elements that involv...
Definition: elements.h:2285
double nodal_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Return the value of the k-th type of the i-th positional variable at the local node n...
Definition: elements.h:2237
Data *const & external_data_pt(const unsigned &i) const
Return a pointer to i-th external data object (const version)
Definition: elements.h:681
void include_internal_data_fd(const unsigned &i)
Set the boolean flag to include the internal datum in the finite difference loop when computing the j...
Definition: elements.h:193
unsigned Nnodal_lagrangian_type
The number of coordinate types requried to intepolate the Lagrangian coordinates in the element...
Definition: elements.h:4039
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation...
Definition: nodes.cc:2314
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object's shape depends on. (Redirects to the nodes' pos...
Definition: elements.h:3371
void get_x(const Vector< double > &s, Vector< double > &x) const
Global coordinates as function of local coordinates. Either via FE representation or via macro-elemen...
Definition: elements.h:1833
void resize_nodes(Vector< unsigned > &nadditional_data_values)
Provide additional storage for a specified number of values at the nodes of the FaceElement. (This is needed, for instance, in free-surface elements, if the non-penetration condition is imposed by Lagrange multipliers whose values are only stored at the surface nodes but not in the interior of the bulk element). nadditional_data_values[n] specifies the number of additional values required at node n of the FaceElement. Note: Since this function is executed separately for each FaceElement, nodes that are common to multiple elements might be resized repeatedly. To avoid this, we only allow a single resize operation by comparing the number of values stored at each node to the number of values the node had when it was simply a member of the associated "bulk" element. There are cases where this will break! – e.g. if a node is common to two FaceElements which require additional storage for distinct quantities. Such cases need to be handled by "hand-crafted" face elements.
Definition: elements.h:4581
static bool Ignore_discontinuous_tangent_warning
Ignores the warning when the tangent vectors from continuous_tangent_and_outer_unit_normal may not be...
Definition: elements.h:4136
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5121
double raw_lagrangian_position(const unsigned &n, const unsigned &i) const
Return i-th Lagrangian coordinate at local node n without using the hanging representation.
Definition: elements.h:3631
Node *const & node_pt(const unsigned &n) const
Return a pointer to the local node n (const version)
Definition: elements.h:2115
double dnodal_position_dt(const unsigned &n, const unsigned &j, const unsigned &i) const
Return the i-th component of j-th derivative of nodal position: d^jx/dt^j at node n...
Definition: elements.h:2231
Vector< unsigned > Bulk_position_type
Vector holding integers to translate additional position types to those of the "bulk" element; e...
Definition: elements.h:4122
A class to specify the initial conditions for a solid body. Solid bodies are often discretised with H...
Definition: elements.h:3244
virtual void describe_nodal_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element[s]. The ostream specifies the output stream to whi...
Definition: elements.cc:1699
virtual void get_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable".
Definition: elements.h:984
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
virtual void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)=0
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed...
void operator=(const SolidElementWithDiagonalMassMatrix &)
Broken assignment operator.
Definition: elements.h:4929
unsigned long * Eqn_number
Storage for the global equation numbers represented by the degrees of freedom in the element...
Definition: elements.h:82
NavierStokesElementWithDiagonalMassMatrices(const NavierStokesElementWithDiagonalMassMatrices &)
Broken copy constructor.
Definition: elements.h:4974
unsigned & ic_time_deriv()
Which time derivative are we currently assigning?
Definition: elements.h:3275
Data * geom_data_pt(const unsigned &j)
A standard FiniteElement is fixed, so there are no geometric data when viewed in its GeomObject incar...
Definition: elements.h:2517
virtual void get_djacobian_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Calculate the derivatives of the elemental Jacobian matrix and residuals with respect to a parameter...
Definition: elements.h:1040
void fill_in_jacobian_from_internal_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
Calculate the contributions to the jacobian from the internal degrees of freedom using finite differe...
Definition: elements.cc:1064
unsigned ngeom_data() const
Broken assignment operator.
Definition: elements.h:3367
void continuous_tangent_and_outer_unit_normal(const Vector< double > &s, Vector< Vector< double > > &tang_vec, Vector< double > &unit_normal) const
Compute the tangent vector(s) and the outer unit normal vector at the specified local coordinate...
Definition: elements.cc:5211
virtual void write_paraview_type(std::ofstream &file_out, const unsigned &nplot) const
Return the paraview element type. Broken virtual. Needs to be implemented for each new geometric elem...
Definition: elements.h:2802
void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
Definition: elements.cc:6197
SolidInitialCondition *& solid_ic_pt()
Pointer to object that describes the initial condition.
Definition: elements.h:3686
void output_paraview(std::ofstream &file_out, const unsigned &nplot) const
Paraview output – this outputs the coordinates at the plot points (for parameter nplot) to specified ...
Definition: elements.h:2727
unsigned Max_newton_iterations
Maximum number of newton iterations.
Definition: elements.cc:1627
GeomObject *& geom_object_pt()
(Reference to) pointer to geom object that specifies the initial condition
Definition: elements.h:3269
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
cstr elem_len * i
Definition: cfortran.h:607
virtual std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.
Definition: elements.h:2865
virtual double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point...
Definition: elements.cc:4049
double dshape_lagrangian(const Vector< double > &s, Shape &psi, DShape &dpsidxi) const
Calculate shape functions and derivatives w.r.t. Lagrangian coordinates at local coordinate s...
Definition: elements.cc:6386
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
double raw_nodal_position(const unsigned &t, const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n, at time level t (t=0: present; t>0: previous time level)...
Definition: elements.h:2164
virtual unsigned nsub_elements_paraview(const unsigned &nplot) const
Return the number of local sub-elements for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2714
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3148
virtual void update_in_solid_position_fd(const unsigned &i)
Function called within the finite difference loop for the solid position dat after a change in any va...
Definition: elements.h:3997
int local_eqn_number(const unsigned long &ieqn_global) const
Return the local equation number corresponding to the ieqn_global-th global equation number...
Definition: elements.h:731
virtual void get_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
Calculate the residuals and the elemental "mass" matrix, the matrix that multiplies the time derivati...
Definition: elements.h:997
void get_ds_bulk_ds_face(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction) const
Calculate the derivatives of the local coordinates in the bulk element with respect to the local coor...
Definition: elements.cc:6089
virtual void fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
Add the elemental contribution to the derivative of the jacobian matrix, mass matrix and the residual...
Definition: elements.cc:1416
virtual void fill_in_contribution_to_hessian_vector_products(Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Fill in contribution to the product of the Hessian (derivative of Jacobian with respect to all variab...
Definition: elements.cc:1469
virtual void fill_in_contribution_to_inner_products(Vector< std::pair< unsigned, unsigned > > const &history_index, Vector< double > &inner_product)
Fill in the contribution to the inner products between given pairs of history values.
Definition: elements.cc:1511
static bool Accept_negative_jacobian
Boolean that if set to true allows a negative jacobian in the transform between global and local coor...
Definition: elements.h:1736
double raw_dnodal_position_dt(const unsigned &n, const unsigned &j, const unsigned &i) const
Return the i-th component of j-th derivative of nodal position: d^jx/dt^j at node n...
Definition: elements.h:2175
void output(std::ostream &outfile)
Output nodal coordinates.
Definition: elements.h:4757
virtual double J_lagrangian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to Lagrangian coordinates at the ipt-th integration poi...
Definition: elements.h:3677
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
Definition: elements.cc:4247
unsigned Boundary_number_in_bulk_mesh
The boundary number in the bulk mesh to which this element is attached.
Definition: elements.h:4141
void get_x(const unsigned &t, const Vector< double > &s, Vector< double > &x)
Global coordinates as function of local coordinates at previous time "level" t (t=0: present; t>0: pr...
Definition: elements.h:1846
unsigned Elemental_dimension
The spatial dimension of the element, i.e. the number of local coordinates used to parametrize it...
Definition: elements.h:1290
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
void(* BulkCoordinateDerivativesFctPt)(const Vector< double > &s, DenseMatrix< double > &ds_bulk_dsface, unsigned &interior_direction)
Typedef for the function that returns the partial derivative of the local coordinates in the bulk ele...
Definition: elements.h:1250
NavierStokesElementWithDiagonalMassMatrices()
Empty constructor.
Definition: elements.h:4968
virtual void update_in_external_fd(const unsigned &i)
Function called within the finite difference loop for external data after a change in any values in t...
Definition: elements.h:487
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
virtual ~NavierStokesElementWithDiagonalMassMatrices()
Virtual destructor.
Definition: elements.h:4971
virtual void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Plot the error when compared against a given time-dependent exact solution . Also calculates the norm...
Definition: elements.h:3008
BulkCoordinateDerivativesFctPt & bulk_coordinate_derivatives_fct_pt()
Return the pointer to the function that returns the derivatives of the bulk coordinates wrt the face ...
Definition: elements.h:4495
virtual void fill_in_contribution_to_djacobian_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Add the elemental contribution to the derivatives of the elemental Jacobian matrix and residuals with...
Definition: elements.cc:1363
A general Finite Element class.
Definition: elements.h:1271
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition: elements.h:2226
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2340
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:3981
GeneralisedElement() GeneralisedElement(const GeneralisedElement &)
Constructor: Initialise all pointers and all values to zero.
Definition: elements.h:611
void set_internal_data_time_stepper(const unsi