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