stored_shape_function_elements.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Non-inline member functions for generic elements
27 
29 #include "shape.h"
30 #include "integral.h"
31 
32 
33 namespace oomph
34 {
35  /// ////////////////////////////////////////////////////////////////////////
36  /// ////////////////////////////////////////////////////////////////////////
37  // Functions for finite elements
38  /// ////////////////////////////////////////////////////////////////////////
39  /// ////////////////////////////////////////////////////////////////////////
40 
41 
42  //=========================================================================
43  /// Delete all the objects stored in the vectors:
44  /// Shape_stored, DShape_local_stored, and D2Shape_local_stored
45  //========================================================================
47  {
51  }
52 
53 
54  //=========================================================================
55  /// Delete stored shape functions
56  //========================================================================
58  {
59  // Can this element delete the stored data?
61  {
62  // Find the number of entries in the shape function storage vector
63  // N.B. This *should* be the same for all three vectors
64  unsigned n_intpt = Shape_stored_pt->size();
65  // Loop over the entries of each vector, in reverse and delete them
66  for (unsigned ipt = n_intpt; ipt > 0; ipt--)
67  {
68  delete (*Shape_stored_pt)[ipt - 1];
69  (*Shape_stored_pt)[ipt - 1] = 0;
70  }
71  // Delete the vector itself
72  delete Shape_stored_pt;
73  }
74  // Reset the pointer to zero {Must do this even if copied}
75  Shape_stored_pt = 0;
76  }
77 
78 
79  //=========================================================================
80  /// Delete stored derivatives of shape functions w.r.t. to local
81  /// coordinates
82  //========================================================================
84  {
85  // Can this element delete the stored data?
87  {
88  unsigned n_intpt = DShape_local_stored_pt->size();
89  for (unsigned ipt = n_intpt; ipt > 0; ipt--)
90  {
91  delete (*DShape_local_stored_pt)[ipt - 1];
92  (*DShape_local_stored_pt)[ipt - 1] = 0;
93  }
94  // Delete the vector itself
96  }
97  // Reset the pointer to zero, must do this, even if copied
99  }
100 
101 
102  //=========================================================================
103  /// Delete stored second derivatives of shape functions w.r.t. to local
104  /// coordinates
105  //========================================================================
107  {
108  // Can this element delete the stored data?
110  {
111  unsigned n_intpt = D2Shape_local_stored_pt->size();
112  for (unsigned ipt = n_intpt; ipt > 0; ipt--)
113  {
114  delete (*D2Shape_local_stored_pt)[ipt - 1];
115  (*D2Shape_local_stored_pt)[ipt - 1] = 0;
116  }
117  // Delete the vector itself
119  }
120  // Reset the pointer to zero, must do it, even if copied
122  }
123 
124 
125  //=========================================================================
126  /// Delete all stored quantities related to derivatives of shape
127  /// fcts w.r.t. to global Eulerian coordinates
128  //========================================================================
130  {
134  }
135 
136  //=========================================================================
137  /// Delete stored derivatives w.r.t. Eulerian coordinates
138  //========================================================================
140  {
141  // If the storage has been allocated and we can delete it
143  {
144  // Find the number of entries in the first vector
145  unsigned n_intpt = DShape_eulerian_stored_pt->size();
146  // Loop over the entries of the vectors, in reverse and delete them
147  for (unsigned ipt = n_intpt; ipt > 0; ipt--)
148  {
149  delete (*DShape_eulerian_stored_pt)[ipt - 1];
150  (*DShape_eulerian_stored_pt)[ipt - 1] = 0;
151  }
152  // Delete the vector itself
154  }
155  // Reset the pointer to zero, even if copied
157  }
158 
159 
160  //=========================================================================
161  /// Delete stored 2nd derivatives w.r.t. Eulerian coordinates
162  //========================================================================
164  {
165  // If the storage has been allocated
167  {
168  // Find the number of entries in the second vector
169  unsigned n_intpt = D2Shape_eulerian_stored_pt->size();
170  // Loop over the entries in reverse and delete them
171  for (unsigned ipt = n_intpt; ipt > 0; ipt--)
172  {
173  delete (*D2Shape_eulerian_stored_pt)[ipt - 1];
174  (*D2Shape_eulerian_stored_pt)[ipt - 1] = 0;
175  }
176  // Delete the vector itself
178  }
179  // Reset the pointer to zero, even if copied
181  }
182 
183 
184  //=========================================================================
185  /// Delete stored Jacobian of mapping between local and global Eulerian
186  /// coordinates
187  //========================================================================
189  {
190  // If the element originally allocated the storage, delete it
192  {
193  // Delete the stored Jacobians
195  }
196  // Reset the pointer to zero, even if copied
198  }
199 
200 
201  //======================================================================
202  /// The destructor cleans up the memory allocated
203  /// for shape function storage.
204  //=======================================================================
206  {
207  // Merely need to call the private functions to clean up storages
210  }
211 
212  //=======================================================================
213  /// Set the spatial integration scheme and also calculate the values of the
214  /// shape functions and their derivatives w.r.t. the local coordinates,
215  /// placing the values into storage so that they may be re-used,
216  /// without recalculation
217  //=======================================================================
219  Integral* const& integral_pt)
220  {
221  // Assign the integration scheme
223 
224  // If we are storing the shape functions and first and second derivatives
225  if (D2Shape_local_stored_pt != 0)
226  {
228  }
229  // If we are storing the shape functions and first derivatives
230  else if (DShape_local_stored_pt != 0)
231  {
233  }
234  // If we are storing the shape functions
235  else if (Shape_stored_pt != 0)
236  {
238  }
239 
240  // If we are storing Eulerian first and second derivatives, recompute them
242  {
244  }
245  // If we are storing Eulerian first derivatives, recompute them
246  else if (DShape_eulerian_stored_pt != 0)
247  {
249  }
250  // If we are storing the Jacobian of the mapping from local to Eulerian
251  // coordinates, recompute it
252  else if (Jacobian_eulerian_stored_pt != 0)
253  {
255  }
256  }
257 
258  //========================================================================
259  /// Calculate the shape functions at the integration points and store in
260  /// internal storage of the element
261  //========================================================================
263  {
264  // Find the number of nodes in the element
265  unsigned n_node = nnode();
266 #ifdef PARANOID
267  if (n_node == 0)
268  {
269  std::string error_message =
270  "FiniteElement::Node_pt must be resized to a value greater than\n";
271  error_message += "zero before calling pre_compute_shape_at_knots()";
272 
273  throw OomphLibError(
274  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
275  }
276 #endif
277  // Find number of interpolated position dofs
278  unsigned n_position_type = nnodal_position_type();
279 
280  // In case we have an exisiting integration scheme, wipe the stored data
282  // Element is now in charge of deleting its own stored data again
284 
285  // Allocate internal storage for the shape functions
287 
288  // Storage for the shape functions and their local derivatives
289  Shape psi(n_node, n_position_type);
290 
291  // Loop over the integration points
292  unsigned n_intpt = integral_pt()->nweight();
293  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
294  {
295  // Get the shape functions
297 
298  // Set up local storage for the shape functions and derivatives
299  Shape* psi_pt = new Shape(n_node, n_position_type);
300 
301  // Copy the values of the shape functions and their local derivatives
302  // into the element storage
303  for (unsigned n = 0; n < n_node; n++)
304  {
305  for (unsigned k = 0; k < n_position_type; k++)
306  {
307  (*psi_pt)(n, k) = psi(n, k);
308  }
309  }
310 
311  // Add the pointers to the shape functions to the internal storage
312  Shape_stored_pt->push_back(psi_pt);
313  } // End of loop over integration points
314  }
315 
316  //========================================================================
317  /// Calculate the shape functions and thir first derivatives
318  /// w.r.t. local coordinates at the integration points and store in
319  /// internal storage of the element
320  //========================================================================
322  {
323  // Find the number of nodes in the element
324  unsigned n_node = nnode();
325 #ifdef PARANOID
326  if (n_node == 0)
327  {
328  std::string error_message =
329  "FiniteElement::Node_pt must be resized to a value greater than\n";
330  error_message +=
331  "zero before calling pre_compute_dshape_local_at_knots()";
332 
333  throw OomphLibError(
334  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
335  }
336 #endif
337  // Find number of interpolated position dofs
338  unsigned n_position_type = nnodal_position_type();
339  // Find spatial dimension of the element
340  unsigned Dim = dim();
341 
342  // In case we have an exisiting integration scheme, wipe the stored data
345  // Element is now in charge of deleting its own stored data again
347 
348  // Allocate internal storage for the shape functions
351 
352  // Storage for the shape functions and their local derivatives
353  Shape psi(n_node, n_position_type);
354  DShape dpsids(n_node, n_position_type, Dim);
355 
356  // Loop over the integration points
357  unsigned n_intpt = integral_pt()->nweight();
358  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
359  {
360  // Get the shape functions and local derivatives at the integration point
361  FiniteElement::dshape_local_at_knot(ipt, psi, dpsids);
362 
363  // Set up local storage for the shape functions and derivatives
364  Shape* psi_pt = new Shape(n_node, n_position_type);
365  DShape* dpsids_pt = new DShape(n_node, n_position_type, Dim);
366 
367  // Copy the values of the shape functions and their local derivatives
368  // into the element storage
369  for (unsigned n = 0; n < n_node; n++)
370  {
371  for (unsigned k = 0; k < n_position_type; k++)
372  {
373  (*psi_pt)(n, k) = psi(n, k);
374 
375  for (unsigned i = 0; i < Dim; i++)
376  {
377  (*dpsids_pt)(n, k, i) = dpsids(n, k, i);
378  }
379  }
380  }
381 
382  // Add the pointers to the shape functions and derivatives to the internal
383  // storage
384  Shape_stored_pt->push_back(psi_pt);
385  DShape_local_stored_pt->push_back(dpsids_pt);
386  } // End of loop over integration points
387  }
388 
389  //========================================================================
390  /// Calculate the shape functions and thir first and second derivatives
391  /// w.r.t. local coordinates at the integration points and store in
392  /// internal storage of the element
393  //========================================================================
395  {
396  // Find the number of nodes in the element
397  unsigned n_node = nnode();
398 #ifdef PARANOID
399  if (n_node == 0)
400  {
401  std::string error_message =
402  "FiniteElement::Node_pt must be resized to a value greater than\n";
403  error_message +=
404  "zero before calling pre_compute_d2shape_local_at_knots()";
405 
406  throw OomphLibError(
407  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
408  }
409 #endif
410  // Find number of interpolated position dofs
411  unsigned n_position_type = nnodal_position_type();
412  // Find spatial dimension of the element
413  unsigned Dim = dim();
414 
415  // Find the number of second derivatives required
416  // N.B. We are assuming that the mixed derivatives are symmetric here
417  unsigned n_deriv = 0;
418  switch (Dim)
419  {
420  case 1:
421  n_deriv = 1;
422  break;
423  case 2:
424  n_deriv = 3;
425  break;
426  case 3:
427  n_deriv = 6;
428  break;
429  default:
430  oomph_info << "Really more than 3 dimensions?" << std::endl;
431  break;
432  }
433 
434  // In case we have an exisiting integration scheme, wipe the stored data
438  // Element is now in charge of deleting its own stored data again
440 
441  // Allocate internal storage for the shape functions
445 
446  // Storage for the shape functions and their local derivatives
447  Shape psi(n_node, n_position_type);
448  DShape dpsids(n_node, n_position_type, Dim);
449  DShape d2psids(n_node, n_position_type, n_deriv);
450 
451  // Loop over the integration points
452  unsigned n_intpt = integral_pt()->nweight();
453  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
454  {
455  // Get the shape functions and local derivatives at the integration point
456  FiniteElement::d2shape_local_at_knot(ipt, psi, dpsids, d2psids);
457 
458  // Set up local storage for the shape functions and derivatives
459  Shape* psi_pt = new Shape(n_node, n_position_type);
460  DShape* dpsids_pt = new DShape(n_node, n_position_type, Dim);
461  DShape* d2psids_pt = new DShape(n_node, n_position_type, n_deriv);
462 
463  // Copy the values of the shape functions and their local derivatives
464  // into the element storage
465  for (unsigned n = 0; n < n_node; n++)
466  {
467  for (unsigned k = 0; k < n_position_type; k++)
468  {
469  (*psi_pt)(n, k) = psi(n, k);
470 
471  for (unsigned i = 0; i < Dim; i++)
472  {
473  (*dpsids_pt)(n, k, i) = dpsids(n, k, i);
474  }
475 
476  for (unsigned i = 0; i < n_deriv; i++)
477  {
478  (*d2psids_pt)(n, k, i) = d2psids(n, k, i);
479  }
480  }
481  }
482 
483  // Add the pointers to the shape functions and derivatives to the internal
484  // storage
485  Shape_stored_pt->push_back(psi_pt);
486  DShape_local_stored_pt->push_back(dpsids_pt);
487  D2Shape_local_stored_pt->push_back(d2psids_pt);
488  } // End of loop over integration points
489  }
490 
491  //=======================================================================
492  /// Calculate the value of the Jacobian of the mapping from local to
493  /// global coordinates at the integration points and store internally
494  //=======================================================================
496  {
497  // Delete previously existing storage
499  // Now we're in change of deletion again
501 
502  // Allocate storage for the stored Jacobian values
504 
505  // Loop over the integration points
506  unsigned n_intpt = integral_pt()->nweight();
507  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
508  {
509  // Add the value of the Jacobian to the internally stored vector
510  Jacobian_eulerian_stored_pt->push_back(
512  }
513  }
514 
515  //========================================================================
516  /// Calculate the values of the derivatives of the shape functions at the
517  /// integration points and store in the internal storage of the element
518  //========================================================================
520  {
521  // Pre-compute the basic shape functions
523 
524  // Find the number of nodes
525  unsigned n_node = nnode();
526  // Get the number of position types and the dimension from the element
527  unsigned n_position_type = this->nnodal_position_type();
528  unsigned n_dim = this->nodal_dimension();
529 
530  // Delete the exisiting stored objects
533  // Now we're in change of deletion again
535 
536  // Allocate storage for the stored shape function derivatives
539 
540  // Assign local variables for the shape function and derivatives
541  Shape psi(n_node, n_position_type);
542  DShape dpsidx(n_node, n_position_type, n_dim);
543 
544  // Loop over the integration points
545  unsigned n_intpt = integral_pt()->nweight();
546  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
547  {
548  // Get the values of the shape function and derivatives at the
549  // integration point and add to the value of the Jacobian to the
550  // internally stored vector
551  Jacobian_eulerian_stored_pt->push_back(
552  FiniteElement::dshape_eulerian_at_knot(ipt, psi, dpsidx));
553 
554  // Set up local storage for the shape function derivatives
555  DShape* dpsidx_pt = new DShape(n_node, n_position_type, n_dim);
556 
557  // Now copy the values over
558  for (unsigned l = 0; l < n_node; l++)
559  {
560  for (unsigned k = 0; k < n_position_type; k++)
561  {
562  // First derivatives
563  for (unsigned i = 0; i < n_dim; i++)
564  {
565  (*dpsidx_pt)(l, k, i) = dpsidx(l, k, i);
566  }
567  }
568  }
569 
570  // Add the pointer to the vector of stored DShape objects
571  DShape_eulerian_stored_pt->push_back(dpsidx_pt);
572  } // End of loop over integration points
573  }
574 
575  //========================================================================
576  /// Calculate the values of the first and second derivatives of the shape
577  /// functions at the integration points and store in the internal storage
578  /// of the element
579  //========================================================================
581  {
582  // Pre-compute the basic shape functions
584 
585  // Find the number of nodes
586  unsigned n_node = nnode();
587  // Get the number of position types and the dimension from element
588  unsigned n_position_type = this->nnodal_position_type();
589  unsigned n_dim = this->nodal_dimension();
590 
591  // Find the number of second derivatives required
592  // N.B. We are assuming that the mixed derivatives are symmetric here
593  unsigned n_deriv = 0;
594  switch (n_dim)
595  {
596  case 1:
597  n_deriv = 1;
598  break;
599  case 2:
600  n_deriv = 3;
601  break;
602  case 3:
603  n_deriv = 6;
604  break;
605  default:
606  oomph_info << "Really more than 3 dimensions?" << std::endl;
607  break;
608  }
609 
610  // Delete the existing objects, if there are any
614  // Now we're in change of deletion again
616 
617  // Allocate storage for the stored shape function derivatives
621 
622  // Assign local variables for the shape function and derivatives
623  Shape psi(n_node, n_position_type);
624  DShape dpsidx(n_node, n_position_type, n_dim);
625  DShape d2psidx(n_node, n_position_type, n_deriv);
626 
627  // Loop over the integration points
628  unsigned n_intpt = integral_pt()->nweight();
629  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
630  {
631  // Get the values of the shape function and derivatives at the
632  // integration point and assign the value of the Jacobian to the
633  // internally stored vector
634  Jacobian_eulerian_stored_pt->push_back(
635  FiniteElement::d2shape_eulerian_at_knot(ipt, psi, dpsidx, d2psidx));
636 
637  // Set up local storage for the shape function derivatives
638  DShape* dpsidx_pt = new DShape(n_node, n_position_type, n_dim);
639  DShape* d2psidx_pt = new DShape(n_node, n_position_type, n_deriv);
640 
641  // Now copy the values over
642  for (unsigned l = 0; l < n_node; l++)
643  {
644  for (unsigned k = 0; k < n_position_type; k++)
645  {
646  // First derivatives
647  for (unsigned i = 0; i < n_dim; i++)
648  {
649  (*dpsidx_pt)(l, k, i) = dpsidx(l, k, i);
650  }
651 
652  // Second derivatives
653  for (unsigned i = 0; i < n_deriv; i++)
654  {
655  (*d2psidx_pt)(l, k, i) = d2psidx(l, k, i);
656  }
657  }
658  }
659 
660  // Add the pointers to the shape function derivatives to the internal
661  // storage
662  DShape_eulerian_stored_pt->push_back(dpsidx_pt);
663  D2Shape_eulerian_stored_pt->push_back(d2psidx_pt);
664  } // End of loop over the shape functions
665  }
666 
667  //=========================================================================
668  /// Return the shape function stored at the ipt-th integration
669  /// point.
670  //=========================================================================
671  void StorableShapeElementBase::shape_at_knot(const unsigned& ipt,
672  Shape& psi) const
673  {
674  // If we are not storing the shape functions, calculate the values
675  if (Shape_stored_pt == 0)
676  {
678  }
679  else
680  {
681  // Read out the stored shape functions
682  // Note this will copy the values by pointer (fast)
683  psi = (*Shape_stored_pt)[ipt];
684  }
685  }
686 
687 
688  //=========================================================================
689  /// Return the shape function and its derivatives w.r.t. the local
690  /// coordinates at the ipt-th integration point.
691  //=========================================================================
693  Shape& psi,
694  DShape& dpsids) const
695  {
696  // If we are not storing the first derivatives, calculate them
697  if (DShape_local_stored_pt == 0)
698  {
699  FiniteElement::dshape_local_at_knot(ipt, psi, dpsids);
700  }
701  else
702  {
703  // Read out the stored shape functions
704  // Set the internal pointers in psi and dpsids
705  psi = (*Shape_stored_pt)[ipt];
706  dpsids = (*DShape_local_stored_pt)[ipt];
707  }
708  }
709 
710  //=========================================================================
711  /// Return the shape function and its first and second derivatives
712  /// w.r.t. the local coordinates at the ipt-th integration point.
713  //=========================================================================
715  Shape& psi,
716  DShape& dpsids,
717  DShape& d2psids) const
718  {
719  // If we are not storing the second derivatives, calculate them on the fly
720  if (D2Shape_local_stored_pt == 0)
721  {
722  FiniteElement::d2shape_local_at_knot(ipt, psi, dpsids, d2psids);
723  }
724  else
725  {
726  // Read out the stored shape functions
727  // Set the internal pointers in psi, dpsids, and d2psids
728  psi = (*Shape_stored_pt)[ipt];
729  dpsids = (*DShape_local_stored_pt)[ipt];
730  d2psids = (*D2Shape_local_stored_pt)[ipt];
731  }
732  }
733 
734 
735  //==========================================================================
736  /// Compute the geometric shape functions, and
737  /// derivatives w.r.t eulerian coordinates at the ipt-th integration point.
738  /// If the values have already been computed, return the stored values.
739  //==========================================================================
741  Shape& psi,
742  DShape& dpsidx) const
743  {
744  // If we are not storing the values, return the calculated values
745  if (DShape_eulerian_stored_pt == 0)
746  {
747  return FiniteElement::dshape_eulerian_at_knot(ipt, psi, dpsidx);
748  }
749  else
750  {
751  // Set internal pointers in the shape functions.
752  psi = (*Shape_stored_pt)[ipt];
753  dpsidx = (*DShape_eulerian_stored_pt)[ipt];
754 
755  // Return the stored value of the jacobian
756  return ((*Jacobian_eulerian_stored_pt)[ipt]);
757  }
758  }
759 
760  //==========================================================================
761  /// Return the geometric shape functions, first and second
762  /// derivatives w.r.t eulerian coordinates at the ipt-th integration point.
763  /// If the values have already been computed, return the stored values.
764  //==========================================================================
766  const unsigned& ipt, Shape& psi, DShape& dpsidx, DShape& d2psidx) const
767  {
768  // If we are not storing the values return the calculated values
770  {
771  return FiniteElement::d2shape_eulerian_at_knot(ipt, psi, dpsidx, d2psidx);
772  }
773  else
774  {
775  // Set internal pointers in the shape functions
776  psi = (*Shape_stored_pt)[ipt];
777  dpsidx = (*DShape_eulerian_stored_pt)[ipt];
778  d2psidx = (*D2Shape_eulerian_stored_pt)[ipt];
779 
780  // Return the stored value of the jacobian
781  return ((*Jacobian_eulerian_stored_pt)[ipt]);
782  }
783  }
784 
785  //==========================================================================
786  /// Return the Jacobian of the mapping between the local and global
787  /// coordinates. If the value has been precomputed return that
788  //==========================================================================
789  double StorableShapeElementBase::J_eulerian_at_knot(const unsigned& ipt) const
790  {
791  // If we are not storing the values, return the calculated values
793  {
795  }
796  else
797  {
798  // Return the stored value of the jacobian
799  return ((*Jacobian_eulerian_stored_pt)[ipt]);
800  }
801  }
802 
803  //========================================================================
804  /// Set the shape functions referenced in the internal vectors to be
805  /// those stored in the StorableShapeElementBase pointed to by element_pt.
806  /// Using this function will allow a saving in the storage required
807  /// for integration schemes in the (most common)
808  /// case when a large number of elements have the same integration scheme
809  //=======================================================================
811  StorableShapeElementBase* const& element_pt)
812  {
813 #ifdef PARANOID
814  // Check that we aren't nulling out
815  if (element_pt->shape_stored_pt() == 0)
816  {
817  std::string error_message =
818  "Element does not have stored shape functions\n";
819 
820  throw OomphLibError(
821  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
822  }
823 #endif
824 
825  // Only do this if the referenced Shape objects are not already the same
826  // Assume that if the stored shape functions are the same, the rest will be
827  if (Shape_stored_pt != element_pt->shape_stored_pt())
828  {
829  // Delete the existing data
831  // Now this element can no longer delete the data pointed at by
832  // the internal vectors
834 
835  // Assign the pointers
836  Shape_stored_pt = element_pt->shape_stored_pt();
839  }
840  }
841 
842  //========================================================================
843  /// Set the stored derivatives of shape functions w.r.t Eulerian coordinates
844  /// to be
845  /// those stored in the StorableShapeElementBase pointed to by element_pt.
846  /// Using this function will allow a saving in the storage required
847  /// for integration schemes in the (most common)
848  /// case when a large number of elements have the same integration scheme
849  //=======================================================================
851  StorableShapeElementBase* const& element_pt)
852  {
854 
855  // Only do this if the referenced Shape objects are not already the same
856  // Assume that if the stored shape functions are the same, the rest will be
858  {
859  // Delete the existing data
861  // Now this element can no longer delete the data pointed at by
862  // the internal vectors
864 
865  // Assign the pointers
869  }
870  }
871 
872 
873  /// ////////////////////////////////////////////////////////////////////////
874  /// ////////////////////////////////////////////////////////////////////////
875  // Functions for solid elements with stored shape functions
876  /// ////////////////////////////////////////////////////////////////////////
877  /// ////////////////////////////////////////////////////////////////////////
878 
879  //=========================================================================
880  /// Delete all the objects and storage associated with stored derivatives
881  /// of shape functions with respect to Lagrangian coordinates
882  //=========================================================================
884  {
888  }
889 
890 
891  //=========================================================================
892  /// Delete all the objects stored in the vectors:
893  /// DShape_lagrangian_stored
894  //========================================================================
896  {
897  // If storage has been allocated
899  {
900  // Find the number of entries in the shape function storage vector
901  unsigned n_intpt = DShape_lagrangian_stored_pt->size();
902  // Loop over the entries of the vectors, in reverse and delete them
903  for (unsigned ipt = n_intpt; ipt > 0; ipt--)
904  {
905  delete (*DShape_lagrangian_stored_pt)[ipt - 1];
906  (*DShape_lagrangian_stored_pt)[ipt - 1] = 0;
907  }
908  // Delete the actual vector
910  }
911  // Reset the pointer to zero
913  }
914 
915 
916  //=========================================================================
917  /// Delete all the objects stored in the vectors:
918  /// D2Shape_lagrangian_stored
919  //========================================================================
921  {
922  // If storage has been allocated
924  {
925  // Now find the number of entries in the second vector
926  unsigned n_intpt = D2Shape_lagrangian_stored_pt->size();
927  // Loop over the entries in reverse and delete them
928  for (unsigned ipt = n_intpt; ipt > 0; ipt--)
929  {
930  delete (*D2Shape_lagrangian_stored_pt)[ipt - 1];
931  (*D2Shape_lagrangian_stored_pt)[ipt - 1] = 0;
932  }
933  // Delete the actual vector
935  }
936  // Reset the pointer to zero
938  }
939 
940  //=======================================================================
941  /// Delete the stored Jacobian of the mapping between the Lagrangian and
942  /// local coordinates.
943  //=======================================================================
945  {
946  // If we allocated the storage, delete it
948  {
949  // Delete the stored Jacobian
951  }
952  // Reset the pointer to zero
954  }
955 
956  //========================================================================
957  /// Calculate the values of the derivatives of the shape functions at the
958  /// integration points and store in the internal storage of the element
959  //========================================================================
961  {
962  // Pre-compute the basic shape functions
964 
965  // Find the number of nodes
966  unsigned n_node = nnode();
967  // Get the number of position types and the dimension from first node
968  // N.B. Assume that it is the same for all nodes
969  unsigned n_lagrangian_type =
970  static_cast<SolidNode*>(node_pt(0))->nlagrangian_type();
971  unsigned n_lagrangian = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
972 
973  // Delete the exisiting stored objects
976  // Now we're in charge of deletion again
978 
979  // Allocate storage for the stored shape function derivatives
982 
983  // Assign local variables for the shape function and derivatives
984  Shape psi(n_node, n_lagrangian_type);
985  DShape dpsidxi(n_node, n_lagrangian_type, n_lagrangian);
986 
987  // Loop over the integration points
988  unsigned n_intpt = integral_pt()->nweight();
989  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
990  {
991  // Get the values of the shape function and derivatives at the
992  // integration point and add to the value of the Jacobian to the
993  // internally stored vector
996 
997  // Set up local storage for the shape function derivatives
998  DShape* dpsidxi_pt = new DShape(n_node, n_lagrangian_type, n_lagrangian);
999 
1000  // Now copy the values over
1001  for (unsigned l = 0; l < n_node; l++)
1002  {
1003  for (unsigned k = 0; k < n_lagrangian_type; k++)
1004  {
1005  // First derivatives
1006  for (unsigned i = 0; i < n_lagrangian; i++)
1007  {
1008  (*dpsidxi_pt)(l, k, i) = dpsidxi(l, k, i);
1009  }
1010  }
1011  }
1012 
1013  // Add the pointer to the vector of stored DShape objects
1014  DShape_lagrangian_stored_pt->push_back(dpsidxi_pt);
1015  } // End of loop over integration points
1016  }
1017 
1018  //========================================================================
1019  /// Calculate the values of the first and second derivatives of the shape
1020  /// functions at the integration points and store in the internal storage
1021  /// of the element
1022  //========================================================================
1024  {
1025  // Pre-compute the basic shape functions
1027 
1028  // Find the number of nodes
1029  unsigned n_node = nnode();
1030  // Get the number of position types and the dimension from first node
1031  // N.B. Assume that it is the same for all nodes
1032  unsigned n_lagrangian_type =
1033  static_cast<SolidNode*>(node_pt(0))->nlagrangian_type();
1034  unsigned n_lagrangian = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
1035 
1036  // Find the number of second derivatives required
1037  // N.B. We are assuming that the mixed derivatives are symmetric here
1038  unsigned n_deriv = 0;
1039  switch (n_lagrangian)
1040  {
1041  case 1:
1042  n_deriv = 1;
1043  break;
1044  case 2:
1045  n_deriv = 3;
1046  break;
1047  case 3:
1048  n_deriv = 6;
1049  break;
1050  default:
1051  oomph_info << "Really more than 3 dimensions?" << std::endl;
1052  break;
1053  }
1054 
1055  // Delete the existing objects, if there are any
1059  // We are in charge of deleting again
1061 
1062  // Allocate storage for the stored shape function derivatives
1066 
1067  // Assign local variables for the shape function and derivatives
1068  Shape psi(n_node, n_lagrangian_type);
1069  DShape dpsidxi(n_node, n_lagrangian_type, n_lagrangian);
1070  DShape d2psidxi(n_node, n_lagrangian_type, n_deriv);
1071 
1072  // Loop over the integration points
1073  unsigned n_intpt = integral_pt()->nweight();
1074  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1075  {
1076  // Get the values of the shape function and derivatives at the
1077  // integration point and assign the value of the Jacobian to the
1078  // internally stored vector
1079  Jacobian_lagrangian_stored_pt->push_back(
1081  ipt, psi, dpsidxi, d2psidxi));
1082 
1083  // Set up local storage for the shape function derivatives
1084  DShape* dpsidxi_pt = new DShape(n_node, n_lagrangian_type, n_lagrangian);
1085  DShape* d2psidxi_pt = new DShape(n_node, n_lagrangian_type, n_deriv);
1086 
1087  // Now copy the values over
1088  for (unsigned l = 0; l < n_node; l++)
1089  {
1090  for (unsigned k = 0; k < n_lagrangian_type; k++)
1091  {
1092  // First derivatives
1093  for (unsigned i = 0; i < n_lagrangian; i++)
1094  {
1095  (*dpsidxi_pt)(l, k, i) = dpsidxi(l, k, i);
1096  }
1097 
1098  // Second derivatives
1099  for (unsigned i = 0; i < n_deriv; i++)
1100  {
1101  (*d2psidxi_pt)(l, k, i) = d2psidxi(l, k, i);
1102  }
1103  }
1104  }
1105 
1106  // Add the pointers to the shape function derivatives to the internal
1107  // storage
1108  DShape_lagrangian_stored_pt->push_back(dpsidxi_pt);
1109  D2Shape_lagrangian_stored_pt->push_back(d2psidxi_pt);
1110  } // End of loop over the shape functions
1111  }
1112 
1113 
1114  //==========================================================================
1115  /// Compute the geometric shape functions, and
1116  /// derivatives w.r.t Lagrangian coordinates at the ipt-th integration point.
1117  /// If the values have already been computed, return the stored values.
1118  //==========================================================================
1120  const unsigned& ipt, Shape& psi, DShape& dpsidxi) const
1121  {
1122  // If we are not storing the values, return the calculated values
1123  if (DShape_lagrangian_stored_pt == 0)
1124  {
1125  return SolidFiniteElement::dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
1126  }
1127  else
1128  {
1129  // Set the internal pointers in the shape functions
1130  psi = shape_stored_pt(ipt);
1131  dpsidxi = (*DShape_lagrangian_stored_pt)[ipt];
1132 
1133  // Return the stored value of the jacobian
1134  return ((*Jacobian_lagrangian_stored_pt)[ipt]);
1135  }
1136  }
1137 
1138  //==========================================================================
1139  /// Compute the geometric shape functions, first and second
1140  /// derivatives w.r.t Lagrangian coordinates at the ipt-th integration point.
1141  /// If the values have already been computed, return the stored values.
1142  //==========================================================================
1144  const unsigned& ipt, Shape& psi, DShape& dpsidxi, DShape& d2psidxi) const
1145  {
1146  // If we are not storing the values return the calculated values
1148  {
1150  ipt, psi, dpsidxi, d2psidxi);
1151  }
1152  else
1153  {
1154  // Set the internal values of the pointers in the Shape objects
1155  psi = shape_stored_pt(ipt);
1156  dpsidxi = (*DShape_lagrangian_stored_pt)[ipt];
1157  d2psidxi = (*D2Shape_lagrangian_stored_pt)[ipt];
1158 
1159  // Return the stored value of the jacobian
1160  return ((*Jacobian_lagrangian_stored_pt)[ipt]);
1161  }
1162  }
1163 
1164 
1165  //========================================================================
1166  /// Set the stored derivatives of shape functions w.r.t Lagrangian coordinates
1167  /// to be
1168  /// those stored in the StorableShapeElementBase pointed to by element_pt.
1169  /// Using this function will allow a saving in the storage required
1170  /// for integration schemes in the (most common)
1171  /// case when a large number of elements have the same integration scheme
1172  //=======================================================================
1174  StorableShapeSolidElementBase* const& element_pt)
1175  {
1177 
1178  // Only do this if the referenced Shape objects are not already the same
1179  // Assume that if the stored shape functions are the same, the rest will be
1181  element_pt->dshape_lagrangian_stored_pt())
1182  {
1183  // Delete the existing data
1185  // Now this element can no longer delete the data pointed at by
1186  // the internal vectors
1188 
1189  // Assign the pointers
1193  element_pt->jacobian_lagrangian_stored_pt();
1194  }
1195  }
1196 
1197 
1198 } // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
virtual void dshape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids) const
Return the geometric shape function and its derivative w.r.t. the local coordinates at the ipt-th int...
Definition: elements.cc:3269
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2179
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3355
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2467
virtual void d2shape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
Return the geometric shape function and its first and second derivatives w.r.t. the local coordinates...
Definition: elements.cc:3304
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2615
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2214
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1967
virtual double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point.
Definition: elements.cc:4198
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2488
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3240
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3250
virtual double d2shape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
Return the geometric shape functions and also first and second derivatives w.r.t. global coordinates ...
Definition: elements.cc:3532
unsigned nlagrangian() const
Access function to # of Lagrangian coordinates.
Definition: geom_objects.h:171
Generic class for numerical integration schemes:
Definition: integral.h:49
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
An OomphLibError object which should be thrown when an run-time error is encountered....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
virtual double d2shape_lagrangian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidxi, DShape &d2psidxi) const
Return the geometric shape functions and also first and second derivatives w.r.t. Lagrangian coordina...
Definition: elements.cc:6862
virtual double dshape_lagrangian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidxi) const
Return the geometric shape functions and also first derivatives w.r.t. Lagrangian coordinates at ipt-...
Definition: elements.cc:6767
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
Base class for elements that allow storage of precomputed shape functions and their derivatives w....
Vector< DShape * > *& dshape_eulerian_stored_pt()
Return a pointer to the vector of pointers to the stored first derivatives of the shape functions w....
void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme – overloaded from the finite element base class since a change in ...
void delete_all_dshape_eulerian_stored()
Delete all storage related to deriatives of shape fcts w.r.t. to global Eulerian coords.
Vector< DShape * > * DShape_local_stored_pt
Pointer to storage for the pointers to the derivatives of the nodal shape functions w....
Vector< DShape * > * D2Shape_local_stored_pt
Pointer to storage for the pointers to the second derivatives of the nodal shape functions w....
virtual ~StorableShapeElementBase()
The destructor cleans up the static memory allocated for shape function storage. Internal and externa...
Vector< DShape * > * D2Shape_eulerian_stored_pt
Pointer to storage for the second derivatives of the shape functions w.r.t. global coordinates at int...
void delete_dshape_local_stored()
Delete stored derivatives of shape functions w.r.t. to local coordinates.
void delete_shape_local_stored()
Delete stored shape functions.
Vector< DShape * > *& dshape_local_stored_pt()
Return a pointer to the vector of pointers to the stored first derivatives of the shape functions w....
void delete_d2shape_local_stored()
Delete stored 2nd derivatives of shape functions w.r.t. to local coordinates.
Vector< DShape * > *& d2shape_eulerian_stored_pt()
Return a pointer to the vector of pointers to the stored second derivatives of the shape functions w....
double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Vector< DShape * > *& d2shape_local_stored_pt()
Return a pointer to the vector of pointers to the stored second derivatives of the shape functions w....
Vector< Shape * > * Shape_stored_pt
Pointer to storage for the pointers to the nodal shape functions at the integration points (knots)
void pre_compute_d2shape_eulerian_at_knots()
Calculate the first and second derivatives of the shape functions w.r.t global coordinates at the int...
Vector< double > * Jacobian_eulerian_stored_pt
Pointer to storage for the Jacobian of the element w.r.t global coordinates.
void delete_d2shape_eulerian_stored()
Delete stored second derivatives of shape functions w.r.t. to global Eulerian coordinates.
void dshape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids) const
Return the geometric shape function and its derivative w.r.t. the local coordinates at the ipt-th int...
void pre_compute_dshape_eulerian_at_knots()
Calculate the first derivatives of the shape functions w.r.t the global coordinates at the integratio...
void set_shape_local_stored_from_element(StorableShapeElementBase *const &element_pt)
Set the shape functions pointed to internally to be those pointed to by the FiniteElement element_pt ...
void pre_compute_J_eulerian_at_knots()
Calculate the Jacobian of the mapping from local to global coordinates at the integration points and ...
double d2shape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
Return the geometric shape functions and also first and second derivatives w.r.t. global coordinates ...
void set_dshape_eulerian_stored_from_element(StorableShapeElementBase *const &element_pt)
Set the derivatives of stored shape functions with respect to the global coordinates to be the same a...
void delete_J_eulerian_stored()
Delete stored Jacobian of mapping between local and global (Eulerian) coordinates.
void pre_compute_dshape_local_at_knots()
Calculate the shape functions and first derivatives w.r.t. local coordinatess at the integration poin...
void delete_dshape_eulerian_stored()
Delete stored deriatives of shape fcts w.r.t. to global Eulerian coords.
bool Can_delete_dshape_eulerian_stored
Boolean to determine whether the element can delete the stored derivatives of shape functions w....
Vector< double > *& jacobian_eulerian_stored_pt()
Return a pointer to the vector of Jacobians of the mapping between the local and global (eulerian) co...
bool Can_delete_shape_local_stored
Boolean to determine whether the element can delete the stored local shape functions.
void pre_compute_shape_at_knots()
Calculate the shape functions at the integration points and store the results internally.
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point.
Vector< DShape * > * DShape_eulerian_stored_pt
Pointer to storage for the derivatives of the shape functions w.r.t. global coordinates at integratio...
void delete_all_shape_local_stored()
Delete all the objects stored in the [...]_local_stored_pt vectors and delete the vectors themselves.
void pre_compute_d2shape_local_at_knots()
Calculate the second derivatives of the shape functions w.r.t. local coordinates at the integration p...
void d2shape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
Return the geometric shape function and its first and second derivatives w.r.t. the local coordinates...
Vector< Shape * > *& shape_stored_pt()
Return a pointer to the vector of pointers to the stored shape functions.
Base class for solid elements that allow storage of precomputed shape functions and their derivatives...
Vector< double > * Jacobian_lagrangian_stored_pt
Pointer to storage for the Jacobian of the mapping between the local and the global Lagrangian coordi...
Vector< DShape * > *& dshape_lagrangian_stored_pt()
Return a pointer to the vector of pointers to the stored first derivatives of the shape functions w....
void pre_compute_d2shape_lagrangian_at_knots()
Calculate the first and second derivatives of the shape functions w.r.t Lagrangian coordinates at the...
void pre_compute_dshape_lagrangian_at_knots()
Calculate the first derivatives of the shape functions w.r.t Lagrangian coordinates at the integratio...
double d2shape_lagrangian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidxi, DShape &d2psidxi) const
Return the geometric shape functions and also first and second derivatives w.r.t. Lagrangian coordina...
Vector< DShape * > * D2Shape_lagrangian_stored_pt
Pointer to storage for the pointers to the second derivatives of the shape functions w....
double dshape_lagrangian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidxi) const
Return the geometric shape functions and also first derivatives w.r.t. Lagrangian coordinates at ipt-...
void delete_J_lagrangian_stored()
Delete stored Jaocbian of mapping between local and Lagrangian coordinates.
Vector< DShape * > *& d2shape_lagrangian_stored_pt()
Return a pointer to the vector of pointers to the stored second derivatives of the shape functions w....
void delete_d2shape_lagrangian_stored()
Delete stored second derivatives of shape functions w.r.t. Lagrangian coordinates.
bool Can_delete_dshape_lagrangian_stored
Boolean to determine whether the element can delete the stored shape function derivatives w....
void delete_all_dshape_lagrangian_stored()
Delete all the objects stored in the [...]_lagrangian_stored_pt vectors and delete the vectors themse...
void delete_dshape_lagrangian_stored()
Delete all the objects stored in the Lagrangian_stored vectors.
Vector< DShape * > * DShape_lagrangian_stored_pt
Pointer to storage for the pointers to the derivatives of the shape functions w.r....
Vector< double > *& jacobian_lagrangian_stored_pt()
Return a pointer to the vector of Jacobians of the mapping between the local and global (eulerian) co...
void set_dshape_lagrangian_stored_from_element(StorableShapeSolidElementBase *const &element_pt)
Set the derivatives of stored shape functions with respect to the lagrangian coordinates to be the sa...
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:60
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...