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-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// 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
47namespace 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
108
109 /// Number of external 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),
601 Ndof(0),
604#ifdef OOMPH_HAS_MPI
605 ,
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
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
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 {
1242 T
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
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 {
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.
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
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$.
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>&,
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
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
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}
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 {
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.
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,
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,
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
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,
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.
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,
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,
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
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),
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 {
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 {
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!
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
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!
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,
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
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);
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.
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
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
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)
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
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.
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 //=======================================================================
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