shape.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //This header file includes generic shape function classes
31 
32 //Include guards to prevent multiple inclusions of the file
33 #ifndef OOMPH_SHAPE_HEADER
34 #define OOMPH_SHAPE_HEADER
35 
36 
37 // Config header generated by autoconfig
38 #ifdef HAVE_CONFIG_H
39  #include <oomph-lib-config.h>
40 #endif
41 
42 #ifdef OOMPH_HAS_MPI
43 #include "mpi.h"
44 #endif
45 
46 //oomph-lib includes
47 #include "Vector.h"
48 #include "matrices.h"
49 #include "orthpoly.h"
50 
51 namespace oomph
52 {
53 
54 
55 //========================================================================
56 /// A Class for shape functions. In simple cases, the shape functions
57 /// have only one index that can be thought of as corresponding to the
58 /// nodal points. In general, however, when quantities and
59 /// their gradients are interpolated separately, the shape function have
60 /// two indices: one corresponding to the nodal points, and the other
61 /// to the "type" of quantity being interpolated: function, derivative, &c
62 /// The second index can also represent the vector coordinate for
63 /// vector-valued (Nedelec) shape functions.
64 ///
65 /// The implementation of Shape functions is designed to permit fast
66 /// copying of entire sets of values by resetting the internal pointer
67 /// to the data, Psi;
68 /// functionality that is required, for example,
69 /// when setting the test functions
70 /// in Galerkin elements and when reading pre-computed values of the shape
71 /// functions.
72 /// In general, we cannot know at construction time whether the pointer to
73 /// the values will be reset or not and, therefore,
74 /// whether the storage for values should be allocated by the object.
75 /// We choose to allocate storage on construction and store an
76 /// additional pointer Allocated_data that \b always addresses the storage
77 /// allocated by the object. If the Psi pointer is reset then this storage
78 /// will be "wasted", but only for the lifetime of the object. The cost for
79 /// non-copied Shape functions is one additional pointer.
80 //=========================================================================
81 class Shape
82 {
83 protected:
84 
85  /// \short Pointer that addresses the storage that will be used to read and
86  /// set the shape functions. The shape functions are packed into
87  /// a flat array of doubles.
88  double *Psi;
89 
90  /// \short Pointer that addresses the storage allocated by the object on
91  /// construction. This will be the same as Psi if the object is not
92  /// copied.
94 
95  ///Size of the first index of the shape function
96  unsigned Index1;
97 
98  ///Size of the second index of the shape function
99  unsigned Index2;
100 
101  ///Private function that checks whether the index is in range
102  void range_check(const unsigned &i, const unsigned &j) const
103  {
104  //If an index is out of range, throw an error
105  if((i >= Index1) || (j >= Index2))
106  {
107  std::ostringstream error_stream;
108  error_stream << "Range Error: ";
109  if(i >= Index1)
110  {
111  error_stream << i << " is not in the range (0," << Index1-1 << ")"
112  << std::endl;
113  }
114  if(j >= Index2)
115  {
116  error_stream << j << " is not in the range (0," << Index2-1 << ")"
117  << std::endl;
118  }
119  throw OomphLibError(error_stream.str(),
120  OOMPH_CURRENT_FUNCTION,
121  OOMPH_EXCEPTION_LOCATION);
122  }
123  }
124 
125 public:
126 
127  /// Constructor for a single-index set of shape functions.
128  Shape(const unsigned &N) : Index1(N), Index2(1)
129  {Allocated_storage = new double[N]; Psi = Allocated_storage;}
130 
131  /// Constructor for a two-index set of shape functions.
132  Shape(const unsigned &N, const unsigned &M) : Index1(N), Index2(M)
133  {Allocated_storage = new double[N*M]; Psi = Allocated_storage;}
134 
135  /// Broken copy constructor
137 
138  /// Default constructor - just assigns a null pointers and zero index
139  /// sizes.
140  Shape() : Psi(0), Allocated_storage(0), Index1(0), Index2(0) {}
141 
142  /// The assignment operator does a shallow copy
143  /// (resets the pointer to the data)
144  void operator=(const Shape &shape)
145  {
146 #ifdef PARANOID
147  //Check the dimensions
148  if((shape.Index1 != Index1) ||
149  (shape.Index2 != Index2))
150  {
151  std::ostringstream error_stream;
152  error_stream << "Cannot assign Shape object:" << std::endl
153  << "Indices do not match "
154  << "LHS: " << Index1 << " " << Index2
155  << ", RHS: " << shape.Index1 << " " << shape.Index2
156  << std::endl;
157  throw OomphLibError(error_stream.str(),
158  OOMPH_CURRENT_FUNCTION,
159  OOMPH_EXCEPTION_LOCATION);
160  }
161 #endif
162  Psi = shape.Psi;
163  }
164 
165  /// The assignment operator does a shallow copy
166  /// (resets the pointer to the data)
167  void operator=(Shape* const &shape_pt)
168  {
169 #ifdef PARANOID
170  //Check the dimensions
171  if((shape_pt->Index1 != Index1) ||
172  (shape_pt->Index2 != Index2))
173  {
174  std::ostringstream error_stream;
175  error_stream << "Cannot assign Shape object:" << std::endl
176  << "Indices do not match "
177  << "LHS: " << Index1 << " " << Index2
178  << ", RHS: " << shape_pt->Index1 << " "
179  << shape_pt->Index2
180  << std::endl;
181  throw OomphLibError(error_stream.str(),
182  OOMPH_CURRENT_FUNCTION,
183  OOMPH_EXCEPTION_LOCATION);
184  }
185 #endif
186  Psi = shape_pt->Psi;
187  }
188 
189  /// Destructor, clear up the memory allocated by the object
191 
192  /// Change the size of the storage
193  void resize(const unsigned& N, const unsigned& M=1)
194  {
195  // Clear old storage
196  delete[] Allocated_storage; Allocated_storage = 0;
197  Psi = 0;
198 
199  // Allocate new storage
200  Index1 = N;
201  Index2 = M;
202  Allocated_storage = new double[N*M];
204  }
205 
206  /// Overload the bracket operator to provide access to values.
207  inline double & operator[](const unsigned &i)
208  {
209 #ifdef RANGE_CHECKING
210  range_check(i,0);
211 #endif
212  return Psi[i*Index2];
213  }
214 
215  /// Overload the bracket operator (const version)
216  inline const double & operator[](const unsigned &i) const
217  {
218 #ifdef RANGE_CHECKING
219  range_check(i,0);
220 #endif
221  return Psi[i*Index2];
222  }
223 
224  /// Overload the round bracket operator to provide access to values.
225  inline double &operator()(const unsigned &i)
226  {
227 #ifdef RANGE_CHECKING
228  range_check(i,0);
229 #endif
230  return Psi[i*Index2];
231  }
232 
233  /// Overload the round bracket operator (const version)
234  inline const double &operator()(const unsigned &i) const
235  {
236 #ifdef RANGE_CHECKING
237  range_check(i,0);
238 #endif
239  return Psi[i*Index2];
240  }
241 
242  /// Overload the round bracket operator, allowing for two indices
243  inline double &operator()(const unsigned &i, const unsigned &j)
244  {
245 #ifdef RANGE_CHECKING
246  range_check(i,j);
247 #endif
248  return Psi[i*Index2 + j];
249  }
250 
251  ///\short Overload the round bracket operator, allowing for two indices
252  /// (const version)
253  inline const double &operator()(const unsigned &i, const unsigned &j)
254  const{
255 #ifdef RANGE_CHECKING
256  range_check(i,j);
257 #endif
258  return Psi[i*Index2 + j];
259  }
260 
261  /// Return the range of index 1 of the shape function object
262  inline unsigned nindex1() const {return Index1;}
263 
264  /// Return the range of index 2 of the shape function object
265  inline unsigned nindex2() const {return Index2;}
266 
267 };
268 
269 //================================================================
270 /// A Class for the derivatives of shape functions
271 /// The class design is essentially the same as Shape, but there is
272 /// on additional index that is used to indicate the coordinate direction in
273 /// which the derivative is taken.
274 //================================================================
275 class DShape
276 {
277  private:
278 
279  /// \short Pointer that addresses the storage that will be used to read and
280  /// set the shape-function derivatives. The values are packed into
281  /// a flat array of doubles.
282  double *DPsi;
283 
284  /// \short Pointer that addresses the storage allocated by the object on
285  /// construction. This will be the same as DPsi if the object is not
286  /// copied.
288 
289  ///Size of the first index of the shape function
290  unsigned Index1;
291 
292  ///Size of the second index of the shape function
293  unsigned Index2;
294 
295  ///Size of the third index of the shape function
296  unsigned Index3;
297 
298  /// Private function that checks whether the indices are in range
299  void range_check(const unsigned &i, const unsigned &j,
300  const unsigned &k) const
301  {
302  //Check the first index
303  if((i >= Index1) || (j >= Index2) || (k >= Index3))
304  {
305  std::ostringstream error_stream;
306  error_stream << "Range Error: ";
307  if(i >= Index1)
308  {
309  error_stream << i
310  << " is not in the range (0," << Index1-1 << ")"
311  << std::endl;
312  }
313  if(j >= Index2)
314  {
315  error_stream << j
316  << " is not in the range (0," << Index2-1 << ")"
317  << std::endl;
318  }
319  if(k >= Index3)
320  {
321  error_stream << k
322  << " is not in the range (0," << Index3-1 << ")"
323  << std::endl;
324  }
325  throw OomphLibError(error_stream.str(),
326  OOMPH_CURRENT_FUNCTION,
327  OOMPH_EXCEPTION_LOCATION);
328  }
329  }
330 
331 
332  public:
333 
334  /// Constructor with two parameters: a single-index shape function
335  DShape(const unsigned &N, const unsigned &P) : Index1(N), Index2(1),
336  Index3(P)
337  {Allocated_storage = new double[N*P]; DPsi = Allocated_storage;}
338 
339  /// Constructor with three paramters: a two-index shape function
340  DShape(const unsigned &N, const unsigned &M, const unsigned &P) :
341  Index1(N), Index2(M), Index3(P)
342  {Allocated_storage = new double[N*M*P]; DPsi = Allocated_storage;}
343 
344  /// Default constructor - just assigns a null pointers and zero index
345  /// sizes.
346  DShape() : DPsi(0), Allocated_storage(0), Index1(0), Index2(0), Index3(0) {}
347 
348  /// Broken copy constructor
350 
351  /// The assignment operator does a shallow copy
352  /// (resets the pointer to the data)
353  void operator=(const DShape &dshape)
354  {
355 #ifdef PARANOID
356  //Check the dimensions
357  if((dshape.Index1 != Index1) ||
358  (dshape.Index2 != Index2) ||
359  (dshape.Index3 != Index3))
360  {
361  std::ostringstream error_stream;
362  error_stream << "Cannot assign DShape object:" << std::endl
363  << "Indices do not match "
364  << "LHS: " << Index1 << " " << Index2 << " "
365  << Index3
366  << ", RHS: " << dshape.Index1 << " "
367  << dshape.Index2 << " " << dshape.Index3
368  << std::endl;
369  throw OomphLibError(error_stream.str(),
370  OOMPH_CURRENT_FUNCTION,
371  OOMPH_EXCEPTION_LOCATION);
372  }
373 #endif
374  DPsi = dshape.DPsi;
375  }
376 
377  /// The assignment operator does a shallow copy
378  /// (resets the pointer to the data)
379  void operator=(DShape* const &dshape_pt)
380  {
381 #ifdef PARANOID
382  //Check the dimensions
383  if((dshape_pt->Index1 != Index1) ||
384  (dshape_pt->Index2 != Index2) ||
385  (dshape_pt->Index3 != Index3))
386  {
387  std::ostringstream error_stream;
388  error_stream << "Cannot assign Shape object:" << std::endl
389  << "Indices do not match "
390  << "LHS: " << Index1 << " " << Index2
391  << " " << Index3
392  << ", RHS: " << dshape_pt->Index1 << " "
393  << dshape_pt->Index2 << " "
394  << dshape_pt->Index3
395  << std::endl;
396  throw OomphLibError(error_stream.str(),
397  OOMPH_CURRENT_FUNCTION,
398  OOMPH_EXCEPTION_LOCATION);
399  }
400 #endif
401  DPsi = dshape_pt->DPsi;
402  }
403 
404 
405 
406  /// Destructor, clean up the memory allocated by this object
408 
409  /// Change the size of the storage. Note that (for some strange reason)
410  /// index2 is the "optional" index, to conform with the existing
411  /// constructor.
412  void resize(const unsigned& N, const unsigned& P, const unsigned& M=1)
413  {
414  // Clear old storage
415  delete[] Allocated_storage; Allocated_storage = 0;
416  DPsi = 0;
417 
418  // Allocate new storage
419  Index1 = N;
420  Index2 = M;
421  Index3 = P;
422  Allocated_storage = new double[N*M*P];
424  }
425 
426  /// Overload the round bracket operator for access to the data
427  inline double &operator()(const unsigned &i, const unsigned &k)
428  {
429 #ifdef RANGE_CHECKING
430  range_check(i,0,k);
431 #endif
432  return DPsi[i*Index2*Index3 + k];
433  }
434 
435  /// Overload the round bracket operator (const version)
436  inline const double &operator()(const unsigned &i, const unsigned &k) const
437  {
438 #ifdef RANGE_CHECKING
439  range_check(i,0,k);
440 #endif
441  return DPsi[i*Index2*Index3 + k];
442  }
443 
444  /// Overload the round bracket operator, with 3 indices
445  inline double &operator()(const unsigned &i, const unsigned &j,
446  const unsigned &k)
447  {
448 #ifdef RANGE_CHECKING
449  range_check(i,j,k);
450 #endif
451  return DPsi[(i*Index2 + j)*Index3 + k];
452  }
453 
454  /// Overload the round bracket operator (const version)
455  inline const double &operator()(const unsigned &i, const unsigned &j,
456  const unsigned &k) const
457  {
458 #ifdef RANGE_CHECKING
459  range_check(i,j,k);
460 #endif
461  return DPsi[(i*Index2 + j)*Index3 + k];
462  }
463 
464  /// \short Direct access to internal storage of data in flat-packed C-style
465  /// column-major format. WARNING: Only for experienced users. Only
466  /// use this if raw speed is of the essence, as in the solid mechanics
467  /// problems.
468  inline double& raw_direct_access(const unsigned long &i)
469  {return DPsi[i];}
470 
471  /// \short Direct access to internal storage of data in flat-packed C-style
472  /// column-major format. WARNING: Only for experienced users. Only
473  /// use this if raw speed is of the essence, as in the solid mechanics
474  /// problems.
475  inline const double& raw_direct_access(const unsigned long &i) const
476  {return DPsi[i];}
477 
478  /// \short Caculate the offset in flat-packed C-style, column-major format,
479  /// required for a given i,j. WARNING: Only for experienced users. Only
480  /// use this if raw speed is of the essence, as in the solid mechanics
481  /// problems.
482  unsigned offset(const unsigned long &i,
483  const unsigned long &j) const
484  {return (i*Index2 + j)*Index3 + 0;}
485 
486 
487 
488  /// Return the range of index 1 of the derivatives of the shape functions
489  inline unsigned long nindex1() const {return Index1;}
490 
491  /// Return the range of index 2 of the derivatives of the shape functions
492  inline unsigned long nindex2() const {return Index2;}
493 
494  /// Return the range of index 3 of the derivatives of the shape functions
495  inline unsigned long nindex3() const {return Index3;}
496 
497 };
498 
499 //======================================================================
500 /// A shape function with a deep copy constructor. This allows for use with stl
501 /// operations (e.g. manipulating vectors of shape functions). A seperate class
502 /// is needed because the basic shape function uses a shallow copy.
503 //======================================================================
504 class ShapeWithDeepCopy : public Shape
505 {
506 public:
507 
508  /// Constructor for a single-index set of shape functions.
509  ShapeWithDeepCopy(const unsigned &N) : Shape(N) {}
510 
511  /// Constructor for a two-index set of shape functions.
512  ShapeWithDeepCopy(const unsigned &N, const unsigned &M) : Shape(N,M) {}
513 
514  /// Default constructor
516 
517  /// Deep copy constructor
519  Shape(old_shape.Index1, old_shape.Index2)
520  {
521  for(unsigned i=0; i<Index1*Index2; i++)
522  {
523  Psi[i] = old_shape.Psi[i];
524  }
525  }
526 
527  /// Broken assignment operator
528  void operator=(const ShapeWithDeepCopy &old_shape)
529  {BrokenCopy::broken_assign("ShapeWithDeepCopy");}
530 
531  /// Destructor, clear up the memory allocated by the object
533 
534 };
535 
536 ////////////////////////////////////////////////////////////////////
537 //
538 // One dimensional shape functions and derivatives.
539 // empty -- simply establishes the template parameters.
540 //
541 ////////////////////////////////////////////////////////////////////
542 
543 namespace OneDimLagrange
544 {
545  /// \short Definition for 1D Lagrange shape functions. The
546  /// value of all the shape functions at the local coordinate s
547  /// are returned in the array Psi.
548  template<unsigned NNODE_1D>
549  void shape(const double &s, double *Psi)
550  {
551  std::ostringstream error_stream;
552  error_stream << "One dimensional Lagrange shape functions "
553  << "have not been defined "
554  << "for " << NNODE_1D << " nodes." << std::endl;
555  throw OomphLibError(error_stream.str(),
556  OOMPH_CURRENT_FUNCTION,
557  OOMPH_EXCEPTION_LOCATION);
558  }
559 
560  /// \short Definition for derivatives of 1D Lagrange shape functions. The
561  /// value of all the shape function derivatives at the local coordinate s
562  /// are returned in the array DPsi.
563  template<unsigned NNODE_1D>
564  void dshape(const double &s, double *DPsi)
565  {
566  std::ostringstream error_stream;
567  error_stream << "One dimensional Lagrange shape function derivatives "
568  << "have not been defined "
569  << "for " << NNODE_1D << " nodes." << std::endl;
570  throw OomphLibError(error_stream.str(),
571  OOMPH_CURRENT_FUNCTION,
572  OOMPH_EXCEPTION_LOCATION);
573  }
574 
575  /// \short Definition for second derivatives of
576  /// 1D Lagrange shape functions. The
577  /// value of all the shape function derivatives at the local coordinate s
578  /// are returned in the array DPsi.
579  template<unsigned NNODE_1D>
580  void d2shape(const double &s, double *DPsi)
581  {
582  std::ostringstream error_stream;
583  error_stream << "One dimensional Lagrange shape function "
584  << "second derivatives "
585  << "have not been defined "
586  << "for " << NNODE_1D << " nodes." << std::endl;
587  throw OomphLibError(error_stream.str(),
588  OOMPH_CURRENT_FUNCTION,
589  OOMPH_EXCEPTION_LOCATION);
590  }
591 
592  /// \short 1D shape functions specialised to linear order (2 Nodes)
593  // Note that the numbering is such that shape[0] is at s = -1.0.
594  // and shape[1] is at s = 1.0
595  template<>
596  inline void shape<2>(const double &s,double *Psi)
597  {
598  Psi[0] = 0.5*(1.0-s);
599  Psi[1] = 0.5*(1.0+s);
600  }
601 
602  /// Derivatives of 1D shape functions specialised to linear order (2 Nodes)
603  template<>
604  inline void dshape<2>(const double &s, double *DPsi)
605  {
606  DPsi[0] = -0.5;
607  DPsi[1] = 0.5;
608  }
609 
610  /// \short Second Derivatives of 1D shape functions,
611  /// specialised to linear order (2 Nodes)
612  template<>
613  inline void d2shape<2>(const double &s, double *DPsi)
614  {
615  DPsi[0] = 0.0;
616  DPsi[1] = 0.0;
617  }
618 
619  /// \short 1D shape functions specialised to quadratic order (3 Nodes)
620  // Note that the numbering is such that shape[0] is at s = -1.0,
621  // shape[1] is at s = 0.0 and shape[2] is at s = 1.0.
622  template<>
623  inline void shape<3>(const double &s, double *Psi)
624  {
625  Psi[0] = 0.5*s*(s-1.0);
626  Psi[1] = 1.0 - s*s;
627  Psi[2] = 0.5*s*(s+1.0);
628  }
629 
630  /// Derivatives of 1D shape functions specialised to quadratic order (3 Nodes)
631  template<>
632  inline void dshape<3>(const double &s, double *DPsi)
633  {
634  DPsi[0] = s - 0.5;
635  DPsi[1] = -2.0*s;
636  DPsi[2] = s + 0.5;
637  }
638 
639 
640  /// Second Derivatives of 1D shape functions, specialised to quadratic order
641  /// (3 Nodes)
642  template<>
643  inline void d2shape<3>(const double &s, double *DPsi)
644  {
645  DPsi[0] = 1.0;
646  DPsi[1] = -2.0;
647  DPsi[2] = 1.0;
648  }
649 
650  /// 1D shape functions specialised to cubic order (4 Nodes)
651  template<>
652  inline void shape<4>(const double &s, double *Psi)
653  {
654  //Output from Maple
655  double t1 = s*s;
656  double t2 = t1*s;
657  double t3 = 0.5625*t2;
658  double t4 = 0.5625*t1;
659  double t5 = 0.625E-1*s;
660  double t7 = 0.16875E1*t2;
661  double t8 = 0.16875E1*s;
662  Psi[0] = -t3+t4+t5-0.625E-1;
663  Psi[1] = t7-t4-t8+0.5625;
664  Psi[2] = -t7-t4+t8+0.5625;
665  Psi[3] = t3+t4-t5-0.625E-1;
666  }
667 
668 
669  /// Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
670  template<>
671  inline void dshape<4>(const double &s, double *DPsi)
672  {
673  //Output from Maple
674  double t1 = s*s;
675  double t2 = 0.16875E1*t1;
676  double t3 = 0.1125E1*s;
677  double t5 = 0.50625E1*t1;
678  DPsi[0] = -t2+t3+0.625E-1;
679  DPsi[1] = t5-t3-0.16875E1;
680  DPsi[2] = -t5-t3+0.16875E1;
681  DPsi[3] = t2+t3-0.625E-1;
682  }
683 
684  /// Second Derivatives of 1D shape functions specialised to cubic
685  /// order (4 Nodes)
686  template<>
687  inline void d2shape<4>(const double &s, double *DPsi)
688  {
689  //Output from Maple (modified by ALH, CHECK IT)
690  double t1 = 2.0*s;
691  double t2 = 0.16875E1*t1;
692  double t5 = 0.50625E1*t1;
693  DPsi[0] = -t2+0.1125E1;
694  DPsi[1] = t5-0.1125E1;
695  DPsi[2] = -t5-0.1125E1;
696  DPsi[3] = t2+0.1125E1;
697  }
698 
699 };
700 
701 //===============================================================
702 /// One Dimensional Hermite shape functions
703 //===============================================================
704 namespace OneDimHermite
705 {
706  //Convention for polynomial numbering scheme
707  //Type 0 is position, 1 is slope
708  //Node 0 is at s=0 and 1 is s=1
709 
710  ///Constructor sets the values of the shape functions at the position s.
711  inline void shape(const double &s, double Psi[2][2])
712  {
713  //Node 0
714  Psi[0][0] = 0.25*(s*s*s - 3.0*s + 2.0);
715  Psi[0][1] = 0.25*(s*s*s - s*s - s + 1.0);
716  //Node 1
717  Psi[1][0] = 0.25*(2.0 + 3.0*s - s*s*s);
718  Psi[1][1] = 0.25*(s*s*s + s*s - s - 1.0);
719  }
720 
721 
722 /// Derivatives of 1D Hermite shape functions
723  inline void dshape(const double &s, double DPsi[2][2])
724  {
725  //Node 0
726  DPsi[0][0] = 0.75*(s*s - 1.0);
727  DPsi[0][1] = 0.25*(3.0*s*s - 2.0*s - 1.0);
728  //Node 1
729  DPsi[1][0] = 0.75*(1.0 - s*s);
730  DPsi[1][1] = 0.25*(3.0*s*s + 2.0*s - 1.0);
731  }
732 
733 /// Second derivatives of the Hermite shape functions
734  inline void d2shape(const double &s, double DPsi[2][2])
735  {
736  //Node 0
737  DPsi[0][0] = 1.5*s;
738  DPsi[0][1] = 0.5*(3.0*s - 1.0);
739  //Node 1
740  DPsi[1][0] = -1.5*s;
741  DPsi[1][1] = 0.5*(3.0*s + 1.0);
742 
743  }
744 
745 };
746 
747 //=====================================================================
748 /// Class that returns the shape functions associated with legendre
749 //=====================================================================
750 template<unsigned NNODE_1D>
752 {
753  static bool Nodes_calculated;
754 
755  public:
757 
758  /// Static function used to populate the stored positions
759  static inline void calculate_nodal_positions()
760  {
761  if(!Nodes_calculated)
762  {
763  Orthpoly::gll_nodes(NNODE_1D,z);
764  Nodes_calculated=true;
765  }
766  }
767 
768  static inline double nodal_position(const unsigned &n)
769  {return z[n];}
770 
771  /// Constructor
772  OneDimensionalLegendreShape(const double &s) : Shape(NNODE_1D)
773  {
774  using namespace Orthpoly;
775 
776  unsigned p = NNODE_1D-1;
777  //Now populate the shape function
778  for(unsigned i=0;i<NNODE_1D;i++)
779  {
780  //If we're at one of the nodes, the value must be 1.0
781  if(std::fabs(s-z[i]) < Orthpoly::eps) {(*this)[i] = 1.0;}
782  //Otherwise use the lagrangian interpolant
783  else
784  {
785  (*this)[i] = (1.0 - s*s)*dlegendre(p,s)/
786  (p*(p+1)*legendre(p,z[i])*(z[i] - s));
787  }
788  }
789  }
790 };
791 
792 template<unsigned NNODE_1D>
793 Vector<double> OneDimensionalLegendreShape<NNODE_1D>::z;
794 
795 template<unsigned NNODE_1D>
796 bool OneDimensionalLegendreShape<NNODE_1D>::Nodes_calculated=false;
797 
798 
799 template <unsigned NNODE_1D>
801 {
802  public:
803  // Constructor
804  OneDimensionalLegendreDShape(const double &s) : Shape(NNODE_1D)
805  {
806  unsigned p = NNODE_1D - 1;
808 
809 
810  bool root=false;
811 
812  for(unsigned i=0;i<NNODE_1D;i++)
813  {
814  unsigned rootnum=0;
815  for(unsigned j=0;j<NNODE_1D;j++)
816  { // Loop over roots to check if
817  if(std::fabs(s-z[j])<10*Orthpoly::eps)
818  { // s happens to be a root.
819  root=true;
820  break;
821  }
822  rootnum+=1;
823  }
824  if(root==true){
825  if(i==rootnum && i==0)
826  {
827  (*this)[i]=-(1.0+p)*p/4.0;
828  }
829  else if(i==rootnum && i==p)
830  {
831  (*this)[i]=(1.0+p)*p/4.0;
832  }
833  else if(i==rootnum)
834  {
835  (*this)[i]=0.0;
836  }
837  else
838  {
839  (*this)[i]=Orthpoly::legendre(p,z[rootnum])
840  /Orthpoly::legendre(p,z[i])/(z[rootnum]-z[i]);
841  }
842  }
843  else
844  {
845  (*this)[i]=((1+s*(s-2*z[i]))/(s-z[i])* Orthpoly::dlegendre(p,s)
846  -(1-s*s)* Orthpoly::ddlegendre(p,s))
847  /p/(p+1.0)/Orthpoly::legendre(p,z[i])/(s-z[i]);
848  }
849  root = false;
850  }
851 
852 
853  }
854 
855 };
856 
857 
858 
859 //=====================================================================
860 /// Non-templated class that returns modal hierachical shape functions
861 /// based on Legendre polynomials
862 //=====================================================================
864 {
865  public:
866  /// Constructor
867  OneDimensionalModalShape(const unsigned p_order, const double &s)
868  : Shape(p_order)
869  {
870  //Populate the shape functions
871  (*this)[0] = 0.5*(1.0 - s);
872  (*this)[1] = 0.5*(1.0 + s);
873  for(unsigned i=2;i<p_order;i++)
874  {
875  (*this)[i] = (0.5*(1.0 - s))*(0.5*(1.0 + s))
876  *Orthpoly::legendre(i-2,s);
877  }
878  }
879 };
880 
882 {
883  public:
884  // Constructor
885  OneDimensionalModalDShape(const unsigned p_order, const double &s)
886  : Shape(p_order)
887  {
888  //Populate the shape functions
889  (*this)[0] = -0.5;
890  (*this)[1] = 0.5;
891  for(unsigned i=2;i<p_order;i++)
892  {
893  (*this)[i] = (0.5*(1.0 - s))*(0.5*(1.0 + s))
894  *Orthpoly::dlegendre(i-2,s)
895  - 0.5*s*Orthpoly::legendre(i-2,s);
896  }
897  }
898 
899 };
900 
901 }
902 
903 #endif
unsigned Index2
Size of the second index of the shape function.
Definition: shape.h:99
double & operator()(const unsigned &i, const unsigned &k)
Overload the round bracket operator for access to the data.
Definition: shape.h:427
void shape(const double &s, double Psi[2][2])
Constructor sets the values of the shape functions at the position s.
Definition: shape.h:711
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void d2shape< 4 >(const double &s, double *DPsi)
Definition: shape.h:687
double dlegendre(const unsigned &p, const double &x)
Calculates first derivative of Legendre polynomial of degree p at x using three term recursive formul...
Definition: orthpoly.h:118
OneDimensionalModalDShape(const unsigned p_order, const double &s)
Definition: shape.h:885
unsigned long nindex3() const
Return the range of index 3 of the derivatives of the shape functions.
Definition: shape.h:495
void dshape< 2 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:604
DShape(const unsigned &N, const unsigned &M, const unsigned &P)
Constructor with three paramters: a two-index shape function.
Definition: shape.h:340
Shape(const Shape &shape)
Broken copy constructor.
Definition: shape.h:136
cstr elem_len * i
Definition: cfortran.h:607
double & raw_direct_access(const unsigned long &i)
Direct access to internal storage of data in flat-packed C-style column-major format. WARNING: Only for experienced users. Only use this if raw speed is of the essence, as in the solid mechanics problems.
Definition: shape.h:468
double * Allocated_storage
Pointer that addresses the storage allocated by the object on construction. This will be the same as ...
Definition: shape.h:93
void operator=(const ShapeWithDeepCopy &old_shape)
Broken assignment operator.
Definition: shape.h:528
unsigned offset(const unsigned long &i, const unsigned long &j) const
Caculate the offset in flat-packed C-style, column-major format, required for a given i...
Definition: shape.h:482
void operator=(const DShape &dshape)
Definition: shape.h:353
void shape< 3 >(const double &s, double *Psi)
1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:623
double & operator[](const unsigned &i)
Overload the bracket operator to provide access to values.
Definition: shape.h:207
void d2shape(const double &s, double DPsi[2][2])
Second derivatives of the Hermite shape functions.
Definition: shape.h:734
void d2shape< 2 >(const double &s, double *DPsi)
Second Derivatives of 1D shape functions, specialised to linear order (2 Nodes)
Definition: shape.h:613
static Vector< double > z
Definition: shape.h:756
unsigned long nindex2() const
Return the range of index 2 of the derivatives of the shape functions.
Definition: shape.h:492
ShapeWithDeepCopy()
Default constructor.
Definition: shape.h:515
double legendre(const unsigned &p, const double &x)
Calculates Legendre polynomial of degree p at x using the three term recurrence relation ...
Definition: orthpoly.h:62
const double & raw_direct_access(const unsigned long &i) const
Direct access to internal storage of data in flat-packed C-style column-major format. WARNING: Only for experienced users. Only use this if raw speed is of the essence, as in the solid mechanics problems.
Definition: shape.h:475
static double nodal_position(const unsigned &n)
Definition: shape.h:768
void shape< 4 >(const double &s, double *Psi)
1D shape functions specialised to cubic order (4 Nodes)
Definition: shape.h:652
const double & operator()(const unsigned &i, const unsigned &j, const unsigned &k) const
Overload the round bracket operator (const version)
Definition: shape.h:455
void operator=(Shape *const &shape_pt)
Definition: shape.h:167
double * Allocated_storage
Pointer that addresses the storage allocated by the object on construction. This will be the same as ...
Definition: shape.h:287
ShapeWithDeepCopy(const ShapeWithDeepCopy &old_shape)
Deep copy constructor.
Definition: shape.h:518
~ShapeWithDeepCopy()
Destructor, clear up the memory allocated by the object.
Definition: shape.h:532
void range_check(const unsigned &i, const unsigned &j, const unsigned &k) const
Private function that checks whether the indices are in range.
Definition: shape.h:299
double * Psi
Pointer that addresses the storage that will be used to read and set the shape functions. The shape functions are packed into a flat array of doubles.
Definition: shape.h:88
void gll_nodes(const unsigned &Nnode, Vector< double > &x)
Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1.
Definition: orthpoly.cc:38
double ddlegendre(const unsigned &p, const double &x)
Calculates second derivative of Legendre polynomial of degree p at x using three term recursive formu...
Definition: orthpoly.h:137
void range_check(const unsigned &i, const unsigned &j) const
Private function that checks whether the index is in range.
Definition: shape.h:102
OneDimensionalModalShape(const unsigned p_order, const double &s)
Constructor.
Definition: shape.h:867
ShapeWithDeepCopy(const unsigned &N, const unsigned &M)
Constructor for a two-index set of shape functions.
Definition: shape.h:512
DShape(const DShape &dshape)
Broken copy constructor.
Definition: shape.h:349
void resize(const unsigned &N, const unsigned &M=1)
Change the size of the storage.
Definition: shape.h:193
const double & operator()(const unsigned &i) const
Overload the round bracket operator (const version)
Definition: shape.h:234
static char t char * s
Definition: cfortran.h:572
void dshape< 3 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:632
void d2shape< 3 >(const double &s, double *DPsi)
Definition: shape.h:643
Shape(const unsigned &N)
Constructor for a single-index set of shape functions.
Definition: shape.h:128
DShape(const unsigned &N, const unsigned &P)
Constructor with two parameters: a single-index shape function.
Definition: shape.h:335
OneDimensionalLegendreDShape(const double &s)
Definition: shape.h:804
unsigned Index2
Size of the second index of the shape function.
Definition: shape.h:293
Shape(const unsigned &N, const unsigned &M)
Constructor for a two-index set of shape functions.
Definition: shape.h:132
void shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
Definition: shape.h:549
ShapeWithDeepCopy(const unsigned &N)
Constructor for a single-index set of shape functions.
Definition: shape.h:509
unsigned nindex1() const
Return the range of index 1 of the shape function object.
Definition: shape.h:262
void dshape(const double &s, double DPsi[2][2])
Derivatives of 1D Hermite shape functions.
Definition: shape.h:723
unsigned Index3
Size of the third index of the shape function.
Definition: shape.h:296
double * DPsi
Pointer that addresses the storage that will be used to read and set the shape-function derivatives...
Definition: shape.h:282
const double & operator()(const unsigned &i, const unsigned &j) const
Overload the round bracket operator, allowing for two indices (const version)
Definition: shape.h:253
~Shape()
Destructor, clear up the memory allocated by the object.
Definition: shape.h:190
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
void dshape(const double &s, double *DPsi)
Definition for derivatives of 1D Lagrange shape functions. The value of all the shape function deriva...
Definition: shape.h:564
OneDimensionalLegendreShape(const double &s)
Constructor.
Definition: shape.h:772
void dshape< 4 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
Definition: shape.h:671
const double & operator()(const unsigned &i, const unsigned &k) const
Overload the round bracket operator (const version)
Definition: shape.h:436
static void calculate_nodal_positions()
Static function used to populate the stored positions.
Definition: shape.h:759
double & operator()(const unsigned &i, const unsigned &j)
Overload the round bracket operator, allowing for two indices.
Definition: shape.h:243
const double & operator[](const unsigned &i) const
Overload the bracket operator (const version)
Definition: shape.h:216
void d2shape(const double &s, double *DPsi)
Definition for second derivatives of 1D Lagrange shape functions. The value of all the shape function...
Definition: shape.h:580
~DShape()
Destructor, clean up the memory allocated by this object.
Definition: shape.h:407
void resize(const unsigned &N, const unsigned &P, const unsigned &M=1)
Definition: shape.h:412
unsigned long nindex1() const
Return the range of index 1 of the derivatives of the shape functions.
Definition: shape.h:489
double & operator()(const unsigned &i)
Overload the round bracket operator to provide access to values.
Definition: shape.h:225
const double eps
Definition: orthpoly.h:57
void operator=(DShape *const &dshape_pt)
Definition: shape.h:379
unsigned Index1
Size of the first index of the shape function.
Definition: shape.h:96
unsigned Index1
Size of the first index of the shape function.
Definition: shape.h:290
unsigned nindex2() const
Return the range of index 2 of the shape function object.
Definition: shape.h:265
double & operator()(const unsigned &i, const unsigned &j, const unsigned &k)
Overload the round bracket operator, with 3 indices.
Definition: shape.h:445
void operator=(const Shape &shape)
Definition: shape.h:144
Class that returns the shape functions associated with legendre.
Definition: shape.h:751
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:596