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// 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// This header file includes generic shape function classes
27
28// Include guards to prevent multiple inclusions of the file
29#ifndef OOMPH_SHAPE_HEADER
30#define OOMPH_SHAPE_HEADER
31
32
33// Config header generated by autoconfig
34#ifdef HAVE_CONFIG_H
35#include <oomph-lib-config.h>
36#endif
37
38#ifdef OOMPH_HAS_MPI
39#include "mpi.h"
40#endif
41
42// oomph-lib includes
43#include "Vector.h"
44#include "matrices.h"
45#include "orthpoly.h"
46
47namespace oomph
48{
49 //========================================================================
50 /// A Class for shape functions. In simple cases, the shape functions
51 /// have only one index that can be thought of as corresponding to the
52 /// nodal points. In general, however, when quantities and
53 /// their gradients are interpolated separately, the shape function have
54 /// two indices: one corresponding to the nodal points, and the other
55 /// to the "type" of quantity being interpolated: function, derivative, &c
56 /// The second index can also represent the vector coordinate for
57 /// vector-valued (Nedelec) shape functions.
58 ///
59 /// The implementation of Shape functions is designed to permit fast
60 /// copying of entire sets of values by resetting the internal pointer
61 /// to the data, Psi;
62 /// functionality that is required, for example,
63 /// when setting the test functions
64 /// in Galerkin elements and when reading pre-computed values of the shape
65 /// functions.
66 /// In general, we cannot know at construction time whether the pointer to
67 /// the values will be reset or not and, therefore,
68 /// whether the storage for values should be allocated by the object.
69 /// We choose to allocate storage on construction and store an
70 /// additional pointer Allocated_data that \b always addresses the storage
71 /// allocated by the object. If the Psi pointer is reset then this storage
72 /// will be "wasted", but only for the lifetime of the object. The cost for
73 /// non-copied Shape functions is one additional pointer.
74 //=========================================================================
75 class Shape
76 {
77 protected:
78 /// Pointer that addresses the storage that will be used to read and
79 /// set the shape functions. The shape functions are packed into
80 /// a flat array of doubles.
81 double* Psi;
82
83 /// Pointer that addresses the storage allocated by the object on
84 /// construction. This will be the same as Psi if the object is not
85 /// copied.
87
88 /// Size of the first index of the shape function
89 unsigned Index1;
90
91 /// Size of the second index of the shape function
92 unsigned Index2;
93
94 /// Private function that checks whether the index is in range
95 void range_check(const unsigned& i, const unsigned& j) const
96 {
97 // If an index is out of range, throw an error
98 if ((i >= Index1) || (j >= Index2))
99 {
100 std::ostringstream error_stream;
101 error_stream << "Range Error: ";
102 if (i >= Index1)
103 {
104 error_stream << i << " is not in the range (0," << Index1 - 1 << ")"
105 << std::endl;
106 }
107 if (j >= Index2)
108 {
109 error_stream << j << " is not in the range (0," << Index2 - 1 << ")"
110 << std::endl;
111 }
112 throw OomphLibError(
113 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
114 }
115 }
116
117 public:
118 /// Constructor for a single-index set of shape functions.
119 Shape(const unsigned& N) : Index1(N), Index2(1)
120 {
121 Allocated_storage = new double[N];
123 }
124
125 /// Constructor for a two-index set of shape functions.
126 Shape(const unsigned& N, const unsigned& M) : Index1(N), Index2(M)
127 {
128 Allocated_storage = new double[N * M];
130 }
131
132 /// Broken copy constructor
133 Shape(const Shape& shape) = delete;
134
135 /// Default constructor - just assigns a null pointers and zero index
136 /// sizes.
138
139 /// The assignment operator does a shallow copy
140 /// (resets the pointer to the data)
141 void operator=(const Shape& shape)
142 {
143#ifdef PARANOID
144 // Check the dimensions
145 if ((shape.Index1 != Index1) || (shape.Index2 != Index2))
146 {
147 std::ostringstream error_stream;
148 error_stream << "Cannot assign Shape object:" << std::endl
149 << "Indices do not match "
150 << "LHS: " << Index1 << " " << Index2
151 << ", RHS: " << shape.Index1 << " " << shape.Index2
152 << std::endl;
153 throw OomphLibError(
154 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
155 }
156#endif
157 Psi = shape.Psi;
158 }
159
160 /// The assignment operator does a shallow copy
161 /// (resets the pointer to the data)
162 void operator=(Shape* const& shape_pt)
163 {
164#ifdef PARANOID
165 // Check the dimensions
166 if ((shape_pt->Index1 != Index1) || (shape_pt->Index2 != Index2))
167 {
168 std::ostringstream error_stream;
169 error_stream << "Cannot assign Shape object:" << std::endl
170 << "Indices do not match "
171 << "LHS: " << Index1 << " " << Index2
172 << ", RHS: " << shape_pt->Index1 << " " << shape_pt->Index2
173 << std::endl;
174 throw OomphLibError(
175 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
176 }
177#endif
178 Psi = shape_pt->Psi;
179 }
180
181 /// Destructor, clear up the memory allocated by the object
183 {
184 delete[] Allocated_storage;
186 }
187
188 /// Change the size of the storage
189 void resize(const unsigned& N, const unsigned& M = 1)
190 {
191 // Clear old storage
192 delete[] Allocated_storage;
194 Psi = 0;
195
196 // Allocate new storage
197 Index1 = N;
198 Index2 = M;
199 Allocated_storage = new double[N * M];
201 }
202
203 /// Overload the bracket operator to provide access to values.
204 inline double& operator[](const unsigned& i)
205 {
206#ifdef RANGE_CHECKING
207 range_check(i, 0);
208#endif
209 return Psi[i * Index2];
210 }
211
212 /// Overload the bracket operator (const version)
213 inline const double& operator[](const unsigned& i) const
214 {
215#ifdef RANGE_CHECKING
216 range_check(i, 0);
217#endif
218 return Psi[i * Index2];
219 }
220
221 /// Overload the round bracket operator to provide access to values.
222 inline double& operator()(const unsigned& i)
223 {
224#ifdef RANGE_CHECKING
225 range_check(i, 0);
226#endif
227 return Psi[i * Index2];
228 }
229
230 /// Overload the round bracket operator (const version)
231 inline const double& operator()(const unsigned& i) const
232 {
233#ifdef RANGE_CHECKING
234 range_check(i, 0);
235#endif
236 return Psi[i * Index2];
237 }
238
239 /// Overload the round bracket operator, allowing for two indices
240 inline double& operator()(const unsigned& i, const unsigned& j)
241 {
242#ifdef RANGE_CHECKING
243 range_check(i, j);
244#endif
245 return Psi[i * Index2 + j];
246 }
247
248 /// Overload the round bracket operator, allowing for two indices
249 /// (const version)
250 inline const double& operator()(const unsigned& i, const unsigned& j) const
251 {
252#ifdef RANGE_CHECKING
253 range_check(i, j);
254#endif
255 return Psi[i * Index2 + j];
256 }
257
258 /// Return the range of index 1 of the shape function object
259 inline unsigned nindex1() const
260 {
261 return Index1;
262 }
263
264 /// Return the range of index 2 of the shape function object
265 inline unsigned nindex2() const
266 {
267 return Index2;
268 }
269 };
270
271 //================================================================
272 /// A Class for the derivatives of shape functions
273 /// The class design is essentially the same as Shape, but there is
274 /// on additional index that is used to indicate the coordinate direction in
275 /// which the derivative is taken.
276 //================================================================
277 class DShape
278 {
279 private:
280 /// Pointer that addresses the storage that will be used to read and
281 /// set the shape-function derivatives. The values are packed into
282 /// a flat array of doubles.
283 double* DPsi;
284
285 /// Pointer that addresses the storage allocated by the object on
286 /// construction. This will be the same as DPsi if the object is not
287 /// copied.
289
290 /// Size of the first index of the shape function
291 unsigned Index1;
292
293 /// Size of the second index of the shape function
294 unsigned Index2;
295
296 /// Size of the third index of the shape function
297 unsigned Index3;
298
299 /// Private function that checks whether the indices are in range
300 void range_check(const unsigned& i,
301 const unsigned& j,
302 const unsigned& k) const
303 {
304 // Check the first index
305 if ((i >= Index1) || (j >= Index2) || (k >= Index3))
306 {
307 std::ostringstream error_stream;
308 error_stream << "Range Error: ";
309 if (i >= Index1)
310 {
311 error_stream << i << " is not in the range (0," << Index1 - 1 << ")"
312 << std::endl;
313 }
314 if (j >= Index2)
315 {
316 error_stream << j << " is not in the range (0," << Index2 - 1 << ")"
317 << std::endl;
318 }
319 if (k >= Index3)
320 {
321 error_stream << k << " is not in the range (0," << Index3 - 1 << ")"
322 << std::endl;
323 }
324 throw OomphLibError(
325 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
326 }
327 }
328
329
330 public:
331 /// Constructor with two parameters: a single-index shape function
332 DShape(const unsigned& N, const unsigned& P)
333 : Index1(N), Index2(1), Index3(P)
334 {
335 Allocated_storage = new double[N * P];
337 }
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 {
343 Allocated_storage = new double[N * M * P];
345 }
346
347 /// Default constructor - just assigns a null pointers and zero index
348 /// sizes.
350
351 /// Broken copy constructor
352 DShape(const DShape& dshape) = delete;
353
354 /// The assignment operator does a shallow copy
355 /// (resets the pointer to the data)
357 {
358#ifdef PARANOID
359 // Check the dimensions
360 if ((dshape.Index1 != Index1) || (dshape.Index2 != Index2) ||
361 (dshape.Index3 != Index3))
362 {
363 std::ostringstream error_stream;
364 error_stream << "Cannot assign DShape object:" << std::endl
365 << "Indices do not match "
366 << "LHS: " << Index1 << " " << Index2 << " " << Index3
367 << ", RHS: " << dshape.Index1 << " " << dshape.Index2
368 << " " << dshape.Index3 << std::endl;
369 throw OomphLibError(
370 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
371 }
372#endif
373 DPsi = dshape.DPsi;
374 }
375
376 /// The assignment operator does a shallow copy
377 /// (resets the pointer to the data)
378 void operator=(DShape* const& dshape_pt)
379 {
380#ifdef PARANOID
381 // Check the dimensions
382 if ((dshape_pt->Index1 != Index1) || (dshape_pt->Index2 != Index2) ||
383 (dshape_pt->Index3 != Index3))
384 {
385 std::ostringstream error_stream;
386 error_stream << "Cannot assign Shape object:" << std::endl
387 << "Indices do not match "
388 << "LHS: " << Index1 << " " << Index2 << " " << Index3
389 << ", RHS: " << dshape_pt->Index1 << " "
390 << dshape_pt->Index2 << " " << dshape_pt->Index3
391 << std::endl;
392 throw OomphLibError(
393 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
394 }
395#endif
396 DPsi = dshape_pt->DPsi;
397 }
398
399
400 /// Destructor, clean up the memory allocated by this object
402 {
403 delete[] Allocated_storage;
405 }
406
407 /// Change the size of the storage. Note that (for some strange reason)
408 /// index2 is the "optional" index, to conform with the existing
409 /// constructor.
410 void resize(const unsigned& N, const unsigned& P, const unsigned& M = 1)
411 {
412 // Clear old storage
413 delete[] Allocated_storage;
415 DPsi = 0;
416
417 // Allocate new storage
418 Index1 = N;
419 Index2 = M;
420 Index3 = P;
421 Allocated_storage = new double[N * M * P];
423 }
424
425 /// Overload the round bracket operator for access to the data
426 inline double& operator()(const unsigned& i, const unsigned& k)
427 {
428#ifdef RANGE_CHECKING
429 range_check(i, 0, k);
430#endif
431 return DPsi[i * Index2 * Index3 + k];
432 }
433
434 /// Overload the round bracket operator (const version)
435 inline const double& operator()(const unsigned& i, const unsigned& k) const
436 {
437#ifdef RANGE_CHECKING
438 range_check(i, 0, k);
439#endif
440 return DPsi[i * Index2 * Index3 + k];
441 }
442
443 /// Overload the round bracket operator, with 3 indices
444 inline double& operator()(const unsigned& i,
445 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,
456 const unsigned& j,
457 const unsigned& k) const
458 {
459#ifdef RANGE_CHECKING
460 range_check(i, j, k);
461#endif
462 return DPsi[(i * Index2 + j) * Index3 + k];
463 }
464
465 /// Direct access to internal storage of data in flat-packed C-style
466 /// column-major format. WARNING: Only for experienced users. Only
467 /// use this if raw speed is of the essence, as in the solid mechanics
468 /// problems.
469 inline double& raw_direct_access(const unsigned long& i)
470 {
471 return DPsi[i];
472 }
473
474 /// Direct access to internal storage of data in flat-packed C-style
475 /// column-major format. WARNING: Only for experienced users. Only
476 /// use this if raw speed is of the essence, as in the solid mechanics
477 /// problems.
478 inline const double& raw_direct_access(const unsigned long& i) const
479 {
480 return DPsi[i];
481 }
482
483 /// Caculate the offset in flat-packed C-style, column-major format,
484 /// required for a given i,j. WARNING: Only for experienced users. Only
485 /// use this if raw speed is of the essence, as in the solid mechanics
486 /// problems.
487 unsigned offset(const unsigned long& i, const unsigned long& j) const
488 {
489 return (i * Index2 + j) * Index3 + 0;
490 }
491
492
493 /// Return the range of index 1 of the derivatives of the shape functions
494 inline unsigned long nindex1() const
495 {
496 return Index1;
497 }
498
499 /// Return the range of index 2 of the derivatives of the shape functions
500 inline unsigned long nindex2() const
501 {
502 return Index2;
503 }
504
505 /// Return the range of index 3 of the derivatives of the shape functions
506 inline unsigned long nindex3() const
507 {
508 return Index3;
509 }
510 };
511
512 //======================================================================
513 /// A shape function with a deep copy constructor. This allows for use with
514 /// stl
515 /// operations (e.g. manipulating vectors of shape functions). A seperate
516 /// class is needed because the basic shape function uses a shallow copy.
517 //======================================================================
519 {
520 public:
521 /// Constructor for a single-index set of shape functions.
522 ShapeWithDeepCopy(const unsigned& N) : Shape(N) {}
523
524 /// Constructor for a two-index set of shape functions.
525 ShapeWithDeepCopy(const unsigned& N, const unsigned& M) : Shape(N, M) {}
526
527 /// Default constructor
529
530 /// Deep copy constructor
532 : Shape(old_shape.Index1, old_shape.Index2)
533 {
534 for (unsigned i = 0; i < Index1 * Index2; i++)
535 {
536 Psi[i] = old_shape.Psi[i];
537 }
538 }
539
540 /// Broken assignment operator
541 void operator=(const ShapeWithDeepCopy& old_shape) = delete;
542
543 /// Destructor, clear up the memory allocated by the object
545 {
546 delete[] Allocated_storage;
548 }
549 };
550
551 /// /////////////////////////////////////////////////////////////////
552 //
553 // One dimensional shape functions and derivatives.
554 // empty -- simply establishes the template parameters.
555 //
556 /// /////////////////////////////////////////////////////////////////
557
558 namespace OneDimLagrange
559 {
560 /// Definition for 1D Lagrange shape functions. The
561 /// value of all the shape functions at the local coordinate s
562 /// are returned in the array Psi.
563 template<unsigned NNODE_1D>
564 void shape(const double& s, double* Psi)
565 {
566 std::ostringstream error_stream;
567 error_stream << "One dimensional Lagrange shape functions "
568 << "have not been defined "
569 << "for " << NNODE_1D << " nodes." << std::endl;
570 throw OomphLibError(
571 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
572 }
573
574 /// Definition for derivatives of 1D Lagrange shape functions. The
575 /// value of all the shape function derivatives at the local coordinate s
576 /// are returned in the array DPsi.
577 template<unsigned NNODE_1D>
578 void dshape(const double& s, double* DPsi)
579 {
580 std::ostringstream error_stream;
581 error_stream << "One dimensional Lagrange shape function derivatives "
582 << "have not been defined "
583 << "for " << NNODE_1D << " nodes." << std::endl;
584 throw OomphLibError(
585 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
586 }
587
588 /// Definition for second derivatives of
589 /// 1D Lagrange shape functions. The
590 /// value of all the shape function derivatives at the local coordinate s
591 /// are returned in the array DPsi.
592 template<unsigned NNODE_1D>
593 void d2shape(const double& s, double* DPsi)
594 {
595 std::ostringstream error_stream;
596 error_stream << "One dimensional Lagrange shape function "
597 << "second derivatives "
598 << "have not been defined "
599 << "for " << NNODE_1D << " nodes." << std::endl;
600 throw OomphLibError(
601 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
602 }
603
604 /// 1D shape functions specialised to linear order (2 Nodes)
605 // Note that the numbering is such that shape[0] is at s = -1.0.
606 // and shape[1] is at s = 1.0
607 template<>
608 inline void shape<2>(const double& s, double* Psi)
609 {
610 Psi[0] = 0.5 * (1.0 - s);
611 Psi[1] = 0.5 * (1.0 + s);
612 }
613
614 /// Derivatives of 1D shape functions specialised to linear order (2 Nodes)
615 template<>
616 inline void dshape<2>(const double& s, double* DPsi)
617 {
618 DPsi[0] = -0.5;
619 DPsi[1] = 0.5;
620 }
621
622 /// Second Derivatives of 1D shape functions,
623 /// specialised to linear order (2 Nodes)
624 template<>
625 inline void d2shape<2>(const double& s, double* DPsi)
626 {
627 DPsi[0] = 0.0;
628 DPsi[1] = 0.0;
629 }
630
631 /// 1D shape functions specialised to quadratic order (3 Nodes)
632 // Note that the numbering is such that shape[0] is at s = -1.0,
633 // shape[1] is at s = 0.0 and shape[2] is at s = 1.0.
634 template<>
635 inline void shape<3>(const double& s, double* Psi)
636 {
637 Psi[0] = 0.5 * s * (s - 1.0);
638 Psi[1] = 1.0 - s * s;
639 Psi[2] = 0.5 * s * (s + 1.0);
640 }
641
642 /// Derivatives of 1D shape functions specialised to quadratic order (3
643 /// Nodes)
644 template<>
645 inline void dshape<3>(const double& s, double* DPsi)
646 {
647 DPsi[0] = s - 0.5;
648 DPsi[1] = -2.0 * s;
649 DPsi[2] = s + 0.5;
650 }
651
652
653 /// Second Derivatives of 1D shape functions, specialised to quadratic order
654 /// (3 Nodes)
655 template<>
656 inline void d2shape<3>(const double& s, double* DPsi)
657 {
658 DPsi[0] = 1.0;
659 DPsi[1] = -2.0;
660 DPsi[2] = 1.0;
661 }
662
663 /// 1D shape functions specialised to cubic order (4 Nodes)
664 template<>
665 inline void shape<4>(const double& s, double* Psi)
666 {
667 // Output from Maple
668 double t1 = s * s;
669 double t2 = t1 * s;
670 double t3 = 0.5625 * t2;
671 double t4 = 0.5625 * t1;
672 double t5 = 0.625E-1 * s;
673 double t7 = 0.16875E1 * t2;
674 double t8 = 0.16875E1 * s;
675 Psi[0] = -t3 + t4 + t5 - 0.625E-1;
676 Psi[1] = t7 - t4 - t8 + 0.5625;
677 Psi[2] = -t7 - t4 + t8 + 0.5625;
678 Psi[3] = t3 + t4 - t5 - 0.625E-1;
679 }
680
681
682 /// Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
683 template<>
684 inline void dshape<4>(const double& s, double* DPsi)
685 {
686 // Output from Maple
687 double t1 = s * s;
688 double t2 = 0.16875E1 * t1;
689 double t3 = 0.1125E1 * s;
690 double t5 = 0.50625E1 * t1;
691 DPsi[0] = -t2 + t3 + 0.625E-1;
692 DPsi[1] = t5 - t3 - 0.16875E1;
693 DPsi[2] = -t5 - t3 + 0.16875E1;
694 DPsi[3] = t2 + t3 - 0.625E-1;
695 }
696
697 /// Second Derivatives of 1D shape functions specialised to cubic
698 /// order (4 Nodes)
699 template<>
700 inline void d2shape<4>(const double& s, double* DPsi)
701 {
702 // Output from Maple (modified by ALH, CHECK IT)
703 double t1 = 2.0 * s;
704 double t2 = 0.16875E1 * t1;
705 double t5 = 0.50625E1 * t1;
706 DPsi[0] = -t2 + 0.1125E1;
707 DPsi[1] = t5 - 0.1125E1;
708 DPsi[2] = -t5 - 0.1125E1;
709 DPsi[3] = t2 + 0.1125E1;
710 }
711
712 }; // namespace OneDimLagrange
713
714 //======================================================================
715 /// One dimensional shape functions and derivatives. Empty -- simply
716 /// establishes the template parameters.
717 //======================================================================
718 namespace OneDimDiscontinuousGalerkin
719 {
720 /// Definition for 1D Lagrange shape functions. The
721 /// value of all the shape functions at the local coordinate s
722 /// are returned in the array Psi.
723 template<unsigned NNODE_1D>
724 void shape(const double& s, double* Psi)
725 {
726 // Create an output stream
727 std::ostringstream error_stream;
728
729 // Create the error message
730 error_stream << "One dimensional Lagrange shape functions "
731 << "have not been defined "
732 << "for " << NNODE_1D << " nodes." << std::endl;
733
734 // Throw the error message
735 throw OomphLibError(
736 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
737 }
738
739 /// Definition for derivatives of 1D Lagrange shape functions. The
740 /// value of all the shape function derivatives at the local coordinate s
741 /// are returned in the array DPsi.
742 template<unsigned NNODE_1D>
743 void dshape(const double& s, double* DPsi)
744 {
745 // Create an output stream
746 std::ostringstream error_stream;
747
748 // Create the error message
749 error_stream << "One dimensional Lagrange shape function derivatives "
750 << "have not been defined "
751 << "for " << NNODE_1D << " nodes." << std::endl;
752
753 // Throw the error message
754 throw OomphLibError(
755 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
756 }
757
758 /// Definition for second derivatives of
759 /// 1D Lagrange shape functions. The
760 /// value of all the shape function derivatives at the local coordinate s
761 /// are returned in the array DPsi.
762 template<unsigned NNODE_1D>
763 void d2shape(const double& s, double* DPsi)
764 {
765 // Create an output stream
766 std::ostringstream error_stream;
767
768 // Create the error message
769 error_stream << "One dimensional Lagrange shape function "
770 << "second derivatives "
771 << "have not been defined "
772 << "for " << NNODE_1D << " nodes." << std::endl;
773
774 // Throw the error message
775 throw OomphLibError(
776 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
777 }
778
779 /// 1D shape functions specialised to linear order (2 Nodes)
780 // Note that the numbering is such that shape[0] is at s = -1.0.
781 // and shape[1] is at s = 1.0
782 template<>
783 inline void shape<2>(const double& s, double* Psi)
784 {
785 Psi[0] = 0.0;
786 Psi[1] = 1.0;
787 }
788
789 /// Derivatives of 1D shape functions specialised to linear order (2 Nodes)
790 template<>
791 inline void dshape<2>(const double& s, double* DPsi)
792 {
793 DPsi[0] = 0.0;
794 DPsi[1] = 0.0;
795 }
796
797 /// Second Derivatives of 1D shape functions,
798 /// specialised to linear order (2 Nodes)
799 template<>
800 inline void d2shape<2>(const double& s, double* DPsi)
801 {
802 DPsi[0] = 0.0;
803 DPsi[1] = 0.0;
804 }
805
806 /// 1D shape functions specialised to quadratic order (3 Nodes)
807 // Note that the numbering is such that shape[0] is at s = -1.0,
808 // shape[1] is at s = 0.0 and shape[2] is at s = 1.0.
809 template<>
810 inline void shape<3>(const double& s, double* Psi)
811 {
812 Psi[0] = 0.0;
813 Psi[1] = 0.5 * (1.0 - s);
814 Psi[2] = 0.5 * (1.0 + s);
815 }
816
817 /// Derivatives of 1D shape functions specialised to quadratic order (3
818 /// Nodes)
819 template<>
820 inline void dshape<3>(const double& s, double* DPsi)
821 {
822 DPsi[0] = 0.0;
823 DPsi[1] = -0.5;
824 DPsi[2] = 0.5;
825 }
826
827 /// Second Derivatives of 1D shape functions, specialised to quadratic order
828 /// (3 Nodes)
829 template<>
830 inline void d2shape<3>(const double& s, double* DPsi)
831 {
832 DPsi[0] = 0.0;
833 DPsi[1] = 0.0;
834 DPsi[2] = 0.0;
835 }
836
837 /// 1D shape functions specialised to cubic order (4 Nodes)
838 template<>
839 inline void shape<4>(const double& s, double* Psi)
840 {
841 Psi[0] = 0.0;
842 Psi[1] = 0.5 * s * (s - 1.0);
843 Psi[2] = 1.0 - s * s;
844 Psi[3] = 0.5 * s * (s + 1.0);
845 }
846
847
848 /// Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
849 template<>
850 inline void dshape<4>(const double& s, double* DPsi)
851 {
852 DPsi[0] = 0.0;
853 DPsi[1] = s - 0.5;
854 DPsi[2] = -2.0 * s;
855 DPsi[3] = s + 0.5;
856 }
857
858 /// Second derivatives of 1D shape functions specialised to cubic
859 /// order (4 nodes)
860 template<>
861 inline void d2shape<4>(const double& s, double* DPsi)
862 {
863 DPsi[0] = 0.0;
864 DPsi[1] = 1.0;
865 DPsi[2] = -2.0;
866 DPsi[3] = 1.0;
867 }
868 }; // namespace OneDimDiscontinuousGalerkin
869
870
871 //======================================================================
872 /// One dimensional shape functions and derivatives. Empty -- simply
873 /// establishes the template parameters.
874 //======================================================================
875 namespace OneDimDiscontinuousGalerkinMixedOrderBasis
876 {
877 /// Definition for 1D Lagrange shape functions. The
878 /// value of all the shape functions at the local coordinate s
879 /// are returned in the array Psi.
880 template<unsigned NNODE_1D>
881 void shape(const double& s, double* Psi)
882 {
883 // Create an output stream
884 std::ostringstream error_stream;
885
886 // Create the error message
887 error_stream << "One dimensional Lagrange shape functions "
888 << "have not been defined "
889 << "for " << NNODE_1D << " nodes." << std::endl;
890
891 // Throw the error message
892 throw OomphLibError(
893 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
894 }
895
896 /// Definition for derivatives of 1D Lagrange shape functions. The
897 /// value of all the shape function derivatives at the local coordinate s
898 /// are returned in the array DPsi.
899 template<unsigned NNODE_1D>
900 void dshape(const double& s, double* DPsi)
901 {
902 // Create an output stream
903 std::ostringstream error_stream;
904
905 // Create the error message
906 error_stream << "One dimensional Lagrange shape function derivatives "
907 << "have not been defined "
908 << "for " << NNODE_1D << " nodes." << std::endl;
909
910 // Throw the error message
911 throw OomphLibError(
912 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
913 }
914
915 /// Definition for second derivatives of
916 /// 1D Lagrange shape functions. The
917 /// value of all the shape function derivatives at the local coordinate s
918 /// are returned in the array DPsi.
919 template<unsigned NNODE_1D>
920 void d2shape(const double& s, double* DPsi)
921 {
922 // Create an output stream
923 std::ostringstream error_stream;
924
925 // Create the error message
926 error_stream << "One dimensional Lagrange shape function "
927 << "second derivatives "
928 << "have not been defined "
929 << "for " << NNODE_1D << " nodes." << std::endl;
930
931 // Throw the error message
932 throw OomphLibError(
933 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
934 }
935
936 /// 1D shape functions specialised to linear order (2 Nodes)
937 // Note that the numbering is such that shape[0] is at s = -1.0.
938 // and shape[1] is at s = 1.0
939 template<>
940 inline void shape<2>(const double& s, double* Psi)
941 {
942 Psi[0] = 0.5 * (1.0 - s);
943 Psi[1] = 0.5 * (1.0 + s);
944 }
945
946 /// Derivatives of 1D shape functions specialised to linear order (2 Nodes)
947 template<>
948 inline void dshape<2>(const double& s, double* DPsi)
949 {
950 DPsi[0] = -0.5;
951 DPsi[1] = 0.5;
952 }
953
954 /// Second Derivatives of 1D shape functions,
955 /// specialised to linear order (2 Nodes)
956 template<>
957 inline void d2shape<2>(const double& s, double* DPsi)
958 {
959 DPsi[0] = 0.0;
960 DPsi[1] = 0.0;
961 }
962
963 /// 1D shape functions specialised to quadratic order (3 Nodes)
964 // Note that the numbering is such that shape[0] is at s = -1.0,
965 // shape[1] is at s = 0.0 and shape[2] is at s = 1.0.
966 template<>
967 inline void shape<3>(const double& s, double* Psi)
968 {
969 Psi[0] = 0.5 * (1.0 - s);
970 Psi[1] = 0.0;
971 Psi[2] = 0.5 * (1.0 + s);
972 }
973
974 /// Derivatives of 1D shape functions specialised to quadratic order (3
975 /// Nodes)
976 template<>
977 inline void dshape<3>(const double& s, double* DPsi)
978 {
979 DPsi[0] = -0.5;
980 DPsi[1] = 0.0;
981 DPsi[2] = 0.5;
982 }
983
984 /// Second Derivatives of 1D shape functions, specialised to quadratic order
985 /// (3 Nodes)
986 template<>
987 inline void d2shape<3>(const double& s, double* DPsi)
988 {
989 DPsi[0] = 0.0;
990 DPsi[1] = 0.0;
991 DPsi[2] = 0.0;
992 }
993
994 /// 1D shape functions specialised to cubic order (4 Nodes)
995 template<>
996 inline void shape<4>(const double& s, double* Psi)
997 {
998 Psi[0] = 0.5 * (1.0 - s);
999 Psi[1] = 0.0;
1000 Psi[2] = 0.0;
1001 Psi[3] = 0.5 * (1.0 + s);
1002 }
1003
1004
1005 /// Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
1006 template<>
1007 inline void dshape<4>(const double& s, double* DPsi)
1008 {
1009 DPsi[0] = -0.5;
1010 DPsi[1] = 0.0;
1011 DPsi[2] = 0.0;
1012 DPsi[3] = 0.5;
1013 }
1014
1015 /// Second derivatives of 1D shape functions specialised to cubic
1016 /// order (4 nodes)
1017 template<>
1018 inline void d2shape<4>(const double& s, double* DPsi)
1019 {
1020 DPsi[0] = 0.0;
1021 DPsi[1] = 0.0;
1022 DPsi[2] = 0.0;
1023 DPsi[3] = 0.0;
1024 }
1025 }; // namespace OneDimDiscontinuousGalerkinMixedOrderBasis
1026
1027
1028 //======================================================================
1029 /// One dimensional shape functions and derivatives. Empty -- simply
1030 /// establishes the template parameters.
1031 //======================================================================
1032 namespace OneDimDiscontinuousGalerkinMixedOrderTest
1033 {
1034 /// Definition for 1D Lagrange shape functions. The
1035 /// value of all the shape functions at the local coordinate s
1036 /// are returned in the array Psi.
1037 template<unsigned NNODE_1D>
1038 void shape(const double& s, double* Psi)
1039 {
1040 // Create an output stream
1041 std::ostringstream error_stream;
1042
1043 // Create the error message
1044 error_stream << "One dimensional Lagrange shape functions "
1045 << "have not been defined "
1046 << "for " << NNODE_1D << " nodes." << std::endl;
1047
1048 // Throw the error message
1049 throw OomphLibError(
1050 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1051 }
1052
1053 /// Definition for derivatives of 1D Lagrange shape functions. The
1054 /// value of all the shape function derivatives at the local coordinate s
1055 /// are returned in the array DPsi.
1056 template<unsigned NNODE_1D>
1057 void dshape(const double& s, double* DPsi)
1058 {
1059 // Create an output stream
1060 std::ostringstream error_stream;
1061
1062 // Create the error message
1063 error_stream << "One dimensional Lagrange shape function derivatives "
1064 << "have not been defined "
1065 << "for " << NNODE_1D << " nodes." << std::endl;
1066
1067 // Throw the error message
1068 throw OomphLibError(
1069 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1070 }
1071
1072 /// Definition for second derivatives of
1073 /// 1D Lagrange shape functions. The
1074 /// value of all the shape function derivatives at the local coordinate s
1075 /// are returned in the array DPsi.
1076 template<unsigned NNODE_1D>
1077 void d2shape(const double& s, double* DPsi)
1078 {
1079 // Create an output stream
1080 std::ostringstream error_stream;
1081
1082 // Create the error message
1083 error_stream << "One dimensional Lagrange shape function "
1084 << "second derivatives "
1085 << "have not been defined "
1086 << "for " << NNODE_1D << " nodes." << std::endl;
1087
1088 // Throw the error message
1089 throw OomphLibError(
1090 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1091 }
1092
1093 /// 1D shape functions specialised to linear order (2 Nodes)
1094 // Note that the numbering is such that shape[0] is at s = -1.0.
1095 // and shape[1] is at s = 1.0
1096 template<>
1097 inline void shape<2>(const double& s, double* Psi)
1098 {
1099 Psi[0] = 0.0;
1100 Psi[1] = 1.0;
1101 }
1102
1103 /// Derivatives of 1D shape functions specialised to linear order (2 Nodes)
1104 template<>
1105 inline void dshape<2>(const double& s, double* DPsi)
1106 {
1107 DPsi[0] = 0.0;
1108 DPsi[1] = 0.0;
1109 }
1110
1111 /// Second Derivatives of 1D shape functions,
1112 /// specialised to linear order (2 Nodes)
1113 template<>
1114 inline void d2shape<2>(const double& s, double* DPsi)
1115 {
1116 DPsi[0] = 0.0;
1117 DPsi[1] = 0.0;
1118 }
1119
1120 /// 1D shape functions specialised to quadratic order (3 Nodes)
1121 // Note that the numbering is such that shape[0] is at s = -1.0,
1122 // shape[1] is at s = 0.0 and shape[2] is at s = 1.0.
1123 template<>
1124 inline void shape<3>(const double& s, double* Psi)
1125 {
1126 Psi[0] = 0.0;
1127 Psi[1] = 0.0;
1128 Psi[2] = 1.0;
1129 }
1130
1131 /// Derivatives of 1D shape functions specialised to quadratic order (3
1132 /// Nodes)
1133 template<>
1134 inline void dshape<3>(const double& s, double* DPsi)
1135 {
1136 DPsi[0] = 0.0;
1137 DPsi[1] = 0.0;
1138 DPsi[2] = 0.0;
1139 }
1140
1141 /// Second Derivatives of 1D shape functions, specialised to quadratic order
1142 /// (3 Nodes)
1143 template<>
1144 inline void d2shape<3>(const double& s, double* DPsi)
1145 {
1146 DPsi[0] = 0.0;
1147 DPsi[1] = 0.0;
1148 DPsi[2] = 0.0;
1149 }
1150
1151 /// 1D shape functions specialised to cubic order (4 Nodes)
1152 template<>
1153 inline void shape<4>(const double& s, double* Psi)
1154 {
1155 Psi[0] = 0.0;
1156 Psi[1] = 0.0;
1157 Psi[2] = 0.0;
1158 Psi[3] = 1.0;
1159 }
1160
1161
1162 /// Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
1163 template<>
1164 inline void dshape<4>(const double& s, double* DPsi)
1165 {
1166 DPsi[0] = 0.0;
1167 DPsi[1] = 0.0;
1168 DPsi[2] = 0.0;
1169 DPsi[3] = 0.0;
1170 }
1171
1172 /// Second derivatives of 1D shape functions specialised to cubic
1173 /// order (4 nodes)
1174 template<>
1175 inline void d2shape<4>(const double& s, double* DPsi)
1176 {
1177 DPsi[0] = 0.0;
1178 DPsi[1] = 0.0;
1179 DPsi[2] = 0.0;
1180 DPsi[3] = 0.0;
1181 }
1182 }; // namespace OneDimDiscontinuousGalerkinMixedOrderTest
1183
1184 //===============================================================
1185 /// One Dimensional Hermite shape functions
1186 //===============================================================
1187 namespace OneDimHermite
1188 {
1189 // Convention for polynomial numbering scheme
1190 // Type 0 is position, 1 is slope
1191 // Node 0 is at s=0 and 1 is s=1
1192
1193 /// Constructor sets the values of the shape functions at the position s.
1194 inline void shape(const double& s, double Psi[2][2])
1195 {
1196 // Node 0
1197 Psi[0][0] = 0.25 * (s * s * s - 3.0 * s + 2.0);
1198 Psi[0][1] = 0.25 * (s * s * s - s * s - s + 1.0);
1199 // Node 1
1200 Psi[1][0] = 0.25 * (2.0 + 3.0 * s - s * s * s);
1201 Psi[1][1] = 0.25 * (s * s * s + s * s - s - 1.0);
1202 }
1203
1204
1205 /// Derivatives of 1D Hermite shape functions
1206 inline void dshape(const double& s, double DPsi[2][2])
1207 {
1208 // Node 0
1209 DPsi[0][0] = 0.75 * (s * s - 1.0);
1210 DPsi[0][1] = 0.25 * (3.0 * s * s - 2.0 * s - 1.0);
1211 // Node 1
1212 DPsi[1][0] = 0.75 * (1.0 - s * s);
1213 DPsi[1][1] = 0.25 * (3.0 * s * s + 2.0 * s - 1.0);
1214 }
1215
1216 /// Second derivatives of the Hermite shape functions
1217 inline void d2shape(const double& s, double DPsi[2][2])
1218 {
1219 // Node 0
1220 DPsi[0][0] = 1.5 * s;
1221 DPsi[0][1] = 0.5 * (3.0 * s - 1.0);
1222 // Node 1
1223 DPsi[1][0] = -1.5 * s;
1224 DPsi[1][1] = 0.5 * (3.0 * s + 1.0);
1225 }
1226
1227 }; // namespace OneDimHermite
1228
1229 //=====================================================================
1230 /// Class that returns the shape functions associated with legendre
1231 //=====================================================================
1232 template<unsigned NNODE_1D>
1234 {
1235 static bool Nodes_calculated;
1236
1237 public:
1239
1240 /// Static function used to populate the stored positions
1241 static inline void calculate_nodal_positions()
1242 {
1243 if (!Nodes_calculated)
1244 {
1245 Orthpoly::gll_nodes(NNODE_1D, z);
1246 Nodes_calculated = true;
1247 }
1248 }
1249
1250 static inline double nodal_position(const unsigned& n)
1251 {
1252 return z[n];
1253 }
1254
1255 /// Constructor
1256 OneDimensionalLegendreShape(const double& s) : Shape(NNODE_1D)
1257 {
1258 using namespace Orthpoly;
1259
1260 unsigned p = NNODE_1D - 1;
1261 // Now populate the shape function
1262 for (unsigned i = 0; i < NNODE_1D; i++)
1263 {
1264 // If we're at one of the nodes, the value must be 1.0
1265 if (std::fabs(s - z[i]) < Orthpoly::eps)
1266 {
1267 (*this)[i] = 1.0;
1268 }
1269 // Otherwise use the lagrangian interpolant
1270 else
1271 {
1272 (*this)[i] = (1.0 - s * s) * dlegendre(p, s) /
1273 (p * (p + 1) * legendre(p, z[i]) * (z[i] - s));
1274 }
1275 }
1276 }
1277 };
1278
1279 template<unsigned NNODE_1D>
1281
1282 template<unsigned NNODE_1D>
1284
1285
1286 template<unsigned NNODE_1D>
1288 {
1289 public:
1290 // Constructor
1291 OneDimensionalLegendreDShape(const double& s) : Shape(NNODE_1D)
1292 {
1293 unsigned p = NNODE_1D - 1;
1295
1296
1297 bool root = false;
1298
1299 for (unsigned i = 0; i < NNODE_1D; i++)
1300 {
1301 unsigned rootnum = 0;
1302 for (unsigned j = 0; j < NNODE_1D; j++)
1303 { // Loop over roots to check if
1304 if (std::fabs(s - z[j]) < 10 * Orthpoly::eps)
1305 { // s happens to be a root.
1306 root = true;
1307 break;
1308 }
1309 rootnum += 1;
1310 }
1311 if (root == true)
1312 {
1313 if (i == rootnum && i == 0)
1314 {
1315 (*this)[i] = -(1.0 + p) * p / 4.0;
1316 }
1317 else if (i == rootnum && i == p)
1318 {
1319 (*this)[i] = (1.0 + p) * p / 4.0;
1320 }
1321 else if (i == rootnum)
1322 {
1323 (*this)[i] = 0.0;
1324 }
1325 else
1326 {
1327 (*this)[i] = Orthpoly::legendre(p, z[rootnum]) /
1328 Orthpoly::legendre(p, z[i]) / (z[rootnum] - z[i]);
1329 }
1330 }
1331 else
1332 {
1333 (*this)[i] =
1334 ((1 + s * (s - 2 * z[i])) / (s - z[i]) * Orthpoly::dlegendre(p, s) -
1335 (1 - s * s) * Orthpoly::ddlegendre(p, s)) /
1336 p / (p + 1.0) / Orthpoly::legendre(p, z[i]) / (s - z[i]);
1337 }
1338 root = false;
1339 }
1340 }
1341 };
1342
1343
1344 //=====================================================================
1345 /// Non-templated class that returns modal hierachical shape functions
1346 /// based on Legendre polynomials
1347 //=====================================================================
1349 {
1350 public:
1351 /// Constructor
1352 OneDimensionalModalShape(const unsigned p_order, const double& s)
1353 : Shape(p_order)
1354 {
1355 // Populate the shape functions
1356 (*this)[0] = 0.5 * (1.0 - s);
1357 (*this)[1] = 0.5 * (1.0 + s);
1358 for (unsigned i = 2; i < p_order; i++)
1359 {
1360 (*this)[i] =
1361 (0.5 * (1.0 - s)) * (0.5 * (1.0 + s)) * Orthpoly::legendre(i - 2, s);
1362 }
1363 }
1364 };
1365
1367 {
1368 public:
1369 // Constructor
1370 OneDimensionalModalDShape(const unsigned p_order, const double& s)
1371 : Shape(p_order)
1372 {
1373 // Populate the shape functions
1374 (*this)[0] = -0.5;
1375 (*this)[1] = 0.5;
1376 for (unsigned i = 2; i < p_order; i++)
1377 {
1378 (*this)[i] = (0.5 * (1.0 - s)) * (0.5 * (1.0 + s)) *
1379 Orthpoly::dlegendre(i - 2, s) -
1380 0.5 * s * Orthpoly::legendre(i - 2, s);
1381 }
1382 }
1383 };
1384
1385} // namespace oomph
1386
1387#endif
static char t char * s
Definition: cfortran.h:568
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
double & operator()(const unsigned &i, const unsigned &j, const unsigned &k)
Overload the round bracket operator, with 3 indices.
Definition: shape.h:444
const double & operator()(const unsigned &i, const unsigned &j, const unsigned &k) const
Overload the round bracket operator (const version)
Definition: shape.h:455
unsigned long nindex3() const
Return the range of index 3 of the derivatives of the shape functions.
Definition: shape.h:506
unsigned long nindex2() const
Return the range of index 2 of the derivatives of the shape functions.
Definition: shape.h:500
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:487
double & raw_direct_access(const unsigned long &i)
Direct access to internal storage of data in flat-packed C-style column-major format....
Definition: shape.h:469
unsigned Index2
Size of the second index of the shape function.
Definition: shape.h:294
const double & operator()(const unsigned &i, const unsigned &k) const
Overload the round bracket operator (const version)
Definition: shape.h:435
~DShape()
Destructor, clean up the memory allocated by this object.
Definition: shape.h:401
void operator=(DShape *const &dshape_pt)
The assignment operator does a shallow copy (resets the pointer to the data)
Definition: shape.h:378
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....
Definition: shape.h:478
void operator=(const DShape &dshape)
The assignment operator does a shallow copy (resets the pointer to the data)
Definition: shape.h:356
unsigned long nindex1() const
Return the range of index 1 of the derivatives of the shape functions.
Definition: shape.h:494
DShape()
Default constructor - just assigns a null pointers and zero index sizes.
Definition: shape.h:349
double * Allocated_storage
Pointer that addresses the storage allocated by the object on construction. This will be the same as ...
Definition: shape.h:288
void resize(const unsigned &N, const unsigned &P, const unsigned &M=1)
Change the size of the storage. Note that (for some strange reason) index2 is the "optional" index,...
Definition: shape.h:410
unsigned Index1
Size of the first index of the shape function.
Definition: shape.h:291
unsigned Index3
Size of the third index of the shape function.
Definition: shape.h:297
DShape(const unsigned &N, const unsigned &M, const unsigned &P)
Constructor with three paramters: a two-index shape function.
Definition: shape.h:340
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:300
DShape(const unsigned &N, const unsigned &P)
Constructor with two parameters: a single-index shape function.
Definition: shape.h:332
DShape(const DShape &dshape)=delete
Broken copy constructor.
double * DPsi
Pointer that addresses the storage that will be used to read and set the shape-function derivatives....
Definition: shape.h:283
double & operator()(const unsigned &i, const unsigned &k)
Overload the round bracket operator for access to the data.
Definition: shape.h:426
OneDimensionalLegendreDShape(const double &s)
Definition: shape.h:1291
Class that returns the shape functions associated with legendre.
Definition: shape.h:1234
OneDimensionalLegendreShape(const double &s)
Constructor.
Definition: shape.h:1256
static double nodal_position(const unsigned &n)
Definition: shape.h:1250
static Vector< double > z
Definition: shape.h:1238
static void calculate_nodal_positions()
Static function used to populate the stored positions.
Definition: shape.h:1241
OneDimensionalModalDShape(const unsigned p_order, const double &s)
Definition: shape.h:1370
Non-templated class that returns modal hierachical shape functions based on Legendre polynomials.
Definition: shape.h:1349
OneDimensionalModalShape(const unsigned p_order, const double &s)
Constructor.
Definition: shape.h:1352
An OomphLibError object which should be thrown when an run-time error is encountered....
A shape function with a deep copy constructor. This allows for use with stl operations (e....
Definition: shape.h:519
ShapeWithDeepCopy()
Default constructor.
Definition: shape.h:528
ShapeWithDeepCopy(const unsigned &N)
Constructor for a single-index set of shape functions.
Definition: shape.h:522
~ShapeWithDeepCopy()
Destructor, clear up the memory allocated by the object.
Definition: shape.h:544
ShapeWithDeepCopy(const unsigned &N, const unsigned &M)
Constructor for a two-index set of shape functions.
Definition: shape.h:525
void operator=(const ShapeWithDeepCopy &old_shape)=delete
Broken assignment operator.
ShapeWithDeepCopy(const ShapeWithDeepCopy &old_shape)
Deep copy constructor.
Definition: shape.h:531
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
const double & operator[](const unsigned &i) const
Overload the bracket operator (const version)
Definition: shape.h:213
const double & operator()(const unsigned &i, const unsigned &j) const
Overload the round bracket operator, allowing for two indices (const version)
Definition: shape.h:250
double * Allocated_storage
Pointer that addresses the storage allocated by the object on construction. This will be the same as ...
Definition: shape.h:86
unsigned nindex1() const
Return the range of index 1 of the shape function object.
Definition: shape.h:259
const double & operator()(const unsigned &i) const
Overload the round bracket operator (const version)
Definition: shape.h:231
double & operator()(const unsigned &i)
Overload the round bracket operator to provide access to values.
Definition: shape.h:222
unsigned Index1
Size of the first index of the shape function.
Definition: shape.h:89
double & operator[](const unsigned &i)
Overload the bracket operator to provide access to values.
Definition: shape.h:204
void resize(const unsigned &N, const unsigned &M=1)
Change the size of the storage.
Definition: shape.h:189
double & operator()(const unsigned &i, const unsigned &j)
Overload the round bracket operator, allowing for two indices.
Definition: shape.h:240
void operator=(const Shape &shape)
The assignment operator does a shallow copy (resets the pointer to the data)
Definition: shape.h:141
unsigned Index2
Size of the second index of the shape function.
Definition: shape.h:92
Shape(const unsigned &N, const unsigned &M)
Constructor for a two-index set of shape functions.
Definition: shape.h:126
double * Psi
Pointer that addresses the storage that will be used to read and set the shape functions....
Definition: shape.h:81
void operator=(Shape *const &shape_pt)
The assignment operator does a shallow copy (resets the pointer to the data)
Definition: shape.h:162
unsigned nindex2() const
Return the range of index 2 of the shape function object.
Definition: shape.h:265
Shape(const Shape &shape)=delete
Broken copy constructor.
Shape(const unsigned &N)
Constructor for a single-index set of shape functions.
Definition: shape.h:119
void range_check(const unsigned &i, const unsigned &j) const
Private function that checks whether the index is in range.
Definition: shape.h:95
Shape()
Default constructor - just assigns a null pointers and zero index sizes.
Definition: shape.h:137
~Shape()
Destructor, clear up the memory allocated by the object.
Definition: shape.h:182
void dshape< 4 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
Definition: shape.h:1007
void d2shape< 2 >(const double &s, double *DPsi)
Second Derivatives of 1D shape functions, specialised to linear order (2 Nodes)
Definition: shape.h:957
void dshape< 3 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:977
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:940
void d2shape< 3 >(const double &s, double *DPsi)
Second Derivatives of 1D shape functions, specialised to quadratic order (3 Nodes)
Definition: shape.h:987
void d2shape< 4 >(const double &s, double *DPsi)
Second derivatives of 1D shape functions specialised to cubic order (4 nodes)
Definition: shape.h:1018
void shape< 3 >(const double &s, double *Psi)
1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:967
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:900
void shape< 4 >(const double &s, double *Psi)
1D shape functions specialised to cubic order (4 Nodes)
Definition: shape.h:996
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:920
void dshape< 2 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:948
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:881
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:1077
void d2shape< 3 >(const double &s, double *DPsi)
Second Derivatives of 1D shape functions, specialised to quadratic order (3 Nodes)
Definition: shape.h:1144
void dshape< 3 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:1134
void shape< 3 >(const double &s, double *Psi)
1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:1124
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:1057
void dshape< 2 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:1105
void d2shape< 4 >(const double &s, double *DPsi)
Second derivatives of 1D shape functions specialised to cubic order (4 nodes)
Definition: shape.h:1175
void dshape< 4 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
Definition: shape.h:1164
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:1038
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:1097
void d2shape< 2 >(const double &s, double *DPsi)
Second Derivatives of 1D shape functions, specialised to linear order (2 Nodes)
Definition: shape.h:1114
void shape< 4 >(const double &s, double *Psi)
1D shape functions specialised to cubic order (4 Nodes)
Definition: shape.h:1153
void dshape< 3 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:820
void dshape< 4 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
Definition: shape.h:850
void shape< 3 >(const double &s, double *Psi)
1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:810
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:743
void d2shape< 3 >(const double &s, double *DPsi)
Second Derivatives of 1D shape functions, specialised to quadratic order (3 Nodes)
Definition: shape.h:830
void dshape< 2 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:791
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:783
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:763
void shape< 4 >(const double &s, double *Psi)
1D shape functions specialised to cubic order (4 Nodes)
Definition: shape.h:839
void d2shape< 4 >(const double &s, double *DPsi)
Second derivatives of 1D shape functions specialised to cubic order (4 nodes)
Definition: shape.h:861
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:724
void d2shape< 2 >(const double &s, double *DPsi)
Second Derivatives of 1D shape functions, specialised to linear order (2 Nodes)
Definition: shape.h:800
void dshape(const double &s, double DPsi[2][2])
Derivatives of 1D Hermite shape functions.
Definition: shape.h:1206
void d2shape(const double &s, double DPsi[2][2])
Second derivatives of the Hermite shape functions.
Definition: shape.h:1217
void shape(const double &s, double Psi[2][2])
Constructor sets the values of the shape functions at the position s.
Definition: shape.h:1194
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:564
void d2shape< 2 >(const double &s, double *DPsi)
Second Derivatives of 1D shape functions, specialised to linear order (2 Nodes)
Definition: shape.h:625
void dshape< 3 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:645
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:593
void dshape< 2 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:616
void d2shape< 3 >(const double &s, double *DPsi)
Second Derivatives of 1D shape functions, specialised to quadratic order (3 Nodes)
Definition: shape.h:656
void d2shape< 4 >(const double &s, double *DPsi)
Second Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
Definition: shape.h:700
void dshape< 4 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
Definition: shape.h:684
void shape< 3 >(const double &s, double *Psi)
1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:635
void shape< 4 >(const double &s, double *Psi)
1D shape functions specialised to cubic order (4 Nodes)
Definition: shape.h:665
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:608
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:578
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:121
const double eps
Definition: orthpoly.h:52
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:144
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:57
void gll_nodes(const unsigned &Nnode, Vector< double > &x)
Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1.
Definition: orthpoly.cc:33
//////////////////////////////////////////////////////////////////// ////////////////////////////////...