Telements.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-2025 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header functions for classes that define Telements
27#ifndef OOMPH_TELEMENT_HEADER
28#define OOMPH_TELEMENT_HEADER
29
30// Config header
31#ifdef HAVE_CONFIG_H
32#include <oomph-lib-config.h>
33#endif
34
35#ifdef OOMPH_HAS_MPI
36#include "mpi.h"
37#endif
38
39// oomph-lib headers
40#include "Vector.h"
41#include "shape.h"
42#include "integral.h"
43#include "elements.h"
44
45namespace oomph
46{
47 ///////////////////////////////////////////////////////////////////////////
48 ///////////////////////////////////////////////////////////////////////////
49 ///////////////////////////////////////////////////////////////////////////
50
51
52 //===================================================
53 /// Triangular Face class
54 //===================================================
55 class TFace
56 {
57 public:
58 /// Constructor: Pass in the three vertex nodes
59 TFace(Node* node1_pt, Node* node2_pt, Node* node3_pt)
60 {
61 if ((node1_pt == node2_pt) || (node2_pt == node3_pt) ||
62 (node1_pt == node3_pt))
63 {
64#ifdef PARANOID
65 std::ostringstream error_stream;
66 error_stream << "TFace cannot have two identical vertex nodes\n";
67 throw OomphLibError(
68 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
69#endif
70 }
71
72 // Sort lexicographically based on pointer address of nodes
73 std::set<Node*> nodes;
74 nodes.insert(node1_pt);
75 nodes.insert(node2_pt);
76 nodes.insert(node3_pt);
77 std::set<Node*>::iterator it = nodes.begin();
78 Node1_pt = (*it);
79 it++;
80 Node2_pt = (*it);
81 it++;
82 Node3_pt = (*it);
83 it++;
84 }
85
86
87 /// Access to the first vertex node
88 Node* node1_pt() const
89 {
90 return Node1_pt;
91 }
92
93 /// Access to the second vertex node
94 Node* node2_pt() const
95 {
96 return Node2_pt;
97 }
98
99 /// Access to the third vertex node
100 Node* node3_pt() const
101 {
102 return Node3_pt;
103 }
104
105 /// Comparison operator
106 bool operator==(const TFace& other) const
107 {
108 if ((Node1_pt == other.node1_pt()) && (Node2_pt == other.node2_pt()) &&
109 (Node3_pt == other.node3_pt()))
110 {
111 return true;
112 }
113 else
114 {
115 return false;
116 }
117 }
118
119
120 /// Less-than operator
121 bool operator<(const TFace& other) const
122 {
123 if (Node1_pt < other.node1_pt())
124 {
125 return true;
126 }
127 else if (Node1_pt == other.node1_pt())
128 {
129 if (Node2_pt < other.node2_pt())
130 {
131 return true;
132 }
133 else if (Node2_pt == other.node2_pt())
134 {
135 if (Node3_pt < other.node3_pt())
136 {
137 return true;
138 }
139 else
140 {
141 return false;
142 }
143 }
144 else
145 {
146 return false;
147 }
148 }
149 else
150 {
151 return false;
152 }
153 }
154
155
156 /// Test whether the face lies on a boundary. Relatively simple
157 /// test, based on all vertices lying on (some) boundary.
158 bool is_on_boundary() const
159 {
160 return (Node1_pt->is_on_boundary() && Node2_pt->is_on_boundary() &&
161 Node3_pt->is_on_boundary());
162 }
163
164
165 /// Test whether the face is a boundary face, i.e. does it
166 /// connnect three boundary nodes?
167 bool is_boundary_face() const
168 {
169 return ((dynamic_cast<BoundaryNodeBase*>(Node1_pt) != 0) &&
170 (dynamic_cast<BoundaryNodeBase*>(Node2_pt) != 0) &&
171 (dynamic_cast<BoundaryNodeBase*>(Node3_pt) != 0));
172 }
173
174 /// Access to pointer to set of mesh boundaries that this
175 /// face occupies; NULL if the node is not on any boundary.
176 /// Construct via set intersection of the boundary sets for the
177 /// associated vertex nodes
178 void get_boundaries_pt(std::set<unsigned>*& boundaries_pt)
179 {
180 std::set<unsigned> set1;
181 std::set<unsigned>* set1_pt = &set1;
182 Node1_pt->get_boundaries_pt(set1_pt);
183 std::set<unsigned> set2;
184 std::set<unsigned>* set2_pt = &set2;
185 Node2_pt->get_boundaries_pt(set2_pt);
186 std::set<unsigned> set3;
187 std::set<unsigned>* set3_pt = &set3;
188 Node3_pt->get_boundaries_pt(set3_pt);
189 std::set<unsigned> aux_set;
190 set_intersection((*set1_pt).begin(),
191 (*set1_pt).end(),
192 (*set2_pt).begin(),
193 (*set2_pt).end(),
194 inserter(aux_set, aux_set.begin()));
195 set_intersection(aux_set.begin(),
196 aux_set.end(),
197 (*set3_pt).begin(),
198 (*set3_pt).end(),
199 inserter((*boundaries_pt), (*boundaries_pt).begin()));
200 }
201
202
203 private:
204 /// First vertex node
205 Node* Node1_pt;
206
207 /// Second vertex node
208 Node* Node2_pt;
209
210 /// Third vertex node
211 Node* Node3_pt;
212 };
213
214
215 ///////////////////////////////////////////////////////////////////////
216 ///////////////////////////////////////////////////////////////////////
217 ///////////////////////////////////////////////////////////////////////
218
219
220 //========================================================================
221 /// A class for those member functions that must be fully specialised
222 /// for the Telements. The fact that member functions of partially
223 /// specialised classes cannot necessarily be fully specialised
224 /// means that we must either fully specialise every class, or use this
225 /// base class to fully specialize only those functions that are required.
226 //========================================================================
227 template<unsigned DIM, unsigned NNODE_1D>
228 class TElementShape
229 {
230 };
231
232 /////////////////////////////////////////////////////////////////////////
233 /// TElementShape inline functions:
234 /////////////////////////////////////////////////////////////////////////
235 template<>
236 class TElementShape<1, 2>
237 {
238 public:
239 //=======================================================================
240 /// Return local coordinates of node j
241 //=======================================================================
242 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
243 {
244 s.resize(1);
245 switch (j)
246 {
247 case 0:
248 s[0] = 0.0;
249 break;
250
251 case 1:
252 s[0] = 1.0;
253 break;
254
255 default:
256 std::ostringstream error_message;
257 error_message
258 << "Element only has two nodes; called with node number " << j
259 << std::endl;
260
261 throw OomphLibError(error_message.str(),
262 OOMPH_CURRENT_FUNCTION,
263 OOMPH_EXCEPTION_LOCATION);
264 }
265 }
266
267
268 //=======================================================================
269 /// Shape function for specific TElement<1,2>
270 //=======================================================================
271 void shape(const Vector<double>& s, Shape& psi) const
272 {
273 psi[0] = 1.0 - s[0];
274 psi[1] = s[0];
275 }
276
277
278 //=======================================================================
279 /// Derivatives of shape functions for specific TElement<2,2>
280 //=======================================================================
281 void dshape_local(const Vector<double>& s, Shape& psi, DShape& dpsids) const
282 {
283 this->shape(s, psi);
284
285 // Derivatives
286 dpsids(0, 0) = -1.0;
287 dpsids(1, 0) = 1.0;
288 }
289
290
291 //=======================================================================
292 /// Second derivatives of shape functions for specific TElement<1,2>:
293 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
294 //=======================================================================
295 void d2shape_local(const Vector<double>& s,
296 Shape& psi,
297 DShape& dpsids,
298 DShape& d2psids) const
299 {
300 this->dshape_local(s, psi, dpsids);
301
302 d2psids(0, 0) = 0.0;
303 d2psids(1, 0) = 0.0;
304 }
305 };
306
307
308 template<>
309 class TElementShape<1, 3>
310 {
311 public:
312 //=======================================================================
313 /// Return local coordinates of node j
314 //=======================================================================
315 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
316 {
317 s.resize(1);
318 switch (j)
319 {
320 case 0:
321 s[0] = 0.0;
322 break;
323
324 case 1:
325 s[0] = 0.5;
326 break;
327
328 case 2:
329 s[0] = 1.0;
330 break;
331
332 default:
333 std::ostringstream error_message;
334 error_message
335 << "Element only has three nodes; called with node number " << j
336 << std::endl;
337
338 throw OomphLibError(error_message.str(),
339 OOMPH_CURRENT_FUNCTION,
340 OOMPH_EXCEPTION_LOCATION);
341 }
342 }
343
344
345 //=======================================================================
346 /// Shape function for specific TElement<1,3>
347 //=======================================================================
348 void shape(const Vector<double>& s, Shape& psi) const
349 {
350 psi[0] = 2.0 * (s[0] - 1.0) * (s[0] - 0.5);
351 psi[1] = 4.0 * (1.0 - s[0]) * s[0];
352 psi[2] = 2.0 * (s[0] - 0.5) * s[0];
353 }
354
355
356 //=======================================================================
357 /// Derivatives of shape functions for specific TElement<1,3>
358 //=======================================================================
359 void dshape_local(const Vector<double>& s, Shape& psi, DShape& dpsids) const
360 {
361 // ALH: Don't know why object qualifier is needed
362 this->shape(s, psi);
363
364 dpsids(0, 0) = 4.0 * s[0] - 3.0;
365 dpsids(1, 0) = 4.0 - 8.0 * s[0];
366 dpsids(2, 0) = 4.0 * s[0] - 1.0;
367 }
368
369
370 //=======================================================================
371 /// Second derivatives of shape functions for specific TElement<1,3>:
372 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
373 //=======================================================================
374 void d2shape_local(const Vector<double>& s,
375 Shape& psi,
376 DShape& dpsids,
377 DShape& d2psids) const
378 {
379 // ALH: Don't know why object qualifier is needed
380 this->dshape_local(s, psi, dpsids);
381
382 d2psids(0, 0) = 4.0;
383 d2psids(1, 0) = -8.0;
384 d2psids(2, 0) = 4.0;
385 }
386 };
387
388 template<>
389 class TElementShape<1, 4>
390 {
391 public:
392 //=======================================================================
393 /// Return local coordinates of node j
394 //=======================================================================
395 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
396 {
397 s.resize(1);
398 switch (j)
399 {
400 case 0:
401 s[0] = 0.0;
402 break;
403
404 case 1:
405 s[0] = (1.0 / 3.0);
406 break;
407
408 case 2:
409 s[0] = (2.0 / 3.0);
410 break;
411
412 case 3:
413 s[0] = 1.0;
414 break;
415
416 default:
417 std::ostringstream error_message;
418 error_message
419 << "Element only has four nodes; called with node number " << j
420 << std::endl;
421
422 throw OomphLibError(error_message.str(),
423 OOMPH_CURRENT_FUNCTION,
424 OOMPH_EXCEPTION_LOCATION);
425 }
426 }
427
428
429 //=======================================================================
430 /// Shape function for specific TElement<1,4>
431 //=======================================================================
432 void shape(const Vector<double>& s, Shape& psi) const
433 {
434 psi[0] = 0.5 * (1.0 - s[0]) * (3.0 * s[0] - 2.0) * (3.0 * s[0] - 1.0);
435 psi[1] = -4.5 * s[0] * (1.0 - s[0]) * (3.0 * s[0] - 2.0);
436 psi[2] = 4.5 * s[0] * (1.0 - s[0]) * (3.0 * s[0] - 1.0);
437 psi[3] = 0.5 * s[0] * (3.0 * s[0] - 2.0) * (3.0 * s[0] - 1.0);
438 }
439
440 //=======================================================================
441 /// Derivatives of shape functions for specific TElement<1,4>
442 //=======================================================================
443 void dshape_local(const Vector<double>& s, Shape& psi, DShape& dpsids) const
444 {
445 // ALH: Don't know why object qualifier is needed
446 this->shape(s, psi);
447
448 dpsids(0, 0) = -13.5 * s[0] * s[0] + 18.0 * s[0] - 5.5;
449 dpsids(1, 0) = 40.5 * s[0] * s[0] - 45.0 * s[0] + 9.0;
450 dpsids(2, 0) = -40.5 * s[0] * s[0] + 36.0 * s[0] - 4.5;
451 dpsids(3, 0) = 13.5 * s[0] * s[0] - 9.0 * s[0] + 1.0;
452 }
453
454 //=======================================================================
455 /// Second derivatives of shape functions for specific TElement<2,4>:
456 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
457 //=======================================================================
458 void d2shape_local(const Vector<double>& s,
459 Shape& psi,
460 DShape& dpsids,
461 DShape& d2psids) const
462 {
463 throw OomphLibError(
464 "Not checked yet", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
465
466 // ALH: Don't know why object qualifier is needed
467 this->dshape_local(s, psi, dpsids);
468
469 d2psids(0, 0) = 0.0;
470 d2psids(1, 0) = 0.0;
471 d2psids(2, 0) = 0.0;
472 d2psids(3, 0) = 0.0;
473 }
474 };
475
476
477 template<>
478 class TElementShape<2, 2>
479 {
480 public:
481 //=======================================================================
482 /// Return local coordinates of node j
483 //=======================================================================
484 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
485 {
486 s.resize(2);
487
488 switch (j)
489 {
490 case 0:
491 s[0] = 1.0;
492 s[1] = 0.0;
493 break;
494
495 case 1:
496 s[0] = 0.0;
497 s[1] = 1.0;
498 break;
499
500 case 2:
501 s[0] = 0.0;
502 s[1] = 0.0;
503 break;
504
505 default:
506 std::ostringstream error_message;
507 error_message
508 << "Element only has three nodes; called with node number " << j
509 << std::endl;
510
511 throw OomphLibError(error_message.str(),
512 OOMPH_CURRENT_FUNCTION,
513 OOMPH_EXCEPTION_LOCATION);
514 }
515 }
516
517
518 //=======================================================================
519 /// Shape function for specific TElement<2,2>
520 //=======================================================================
521 void shape(const Vector<double>& s, Shape& psi) const
522 {
523 psi[0] = s[0];
524 psi[1] = s[1];
525 psi[2] = 1.0 - s[0] - s[1];
526 }
527
528
529 //=======================================================================
530 /// Derivatives of shape functions for specific TElement<2,2>
531 //=======================================================================
532 void dshape_local(const Vector<double>& s, Shape& psi, DShape& dpsids) const
533 {
534 this->shape(s, psi);
535
536 // Derivatives
537 dpsids(0, 0) = 1.0;
538 dpsids(0, 1) = 0.0;
539 dpsids(1, 0) = 0.0;
540 dpsids(1, 1) = 1.0;
541 dpsids(2, 0) = -1.0;
542 dpsids(2, 1) = -1.0;
543 }
544
545
546 //=======================================================================
547 /// Second derivatives of shape functions for specific TElement<2,2>:
548 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
549 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
550 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
551 //=======================================================================
552 void d2shape_local(const Vector<double>& s,
553 Shape& psi,
554 DShape& dpsids,
555 DShape& d2psids) const
556 {
557 this->dshape_local(s, psi, dpsids);
558
559 for (unsigned i = 0; i < 3; i++)
560 {
561 d2psids(i, 0) = 0.0;
562 d2psids(i, 1) = 0.0;
563 d2psids(i, 2) = 0.0;
564 }
565 }
566 };
567
568 template<>
569 class TElementShape<2, 3>
570 {
571 public:
572 //=======================================================================
573 /// Return local coordinates of node j
574 //=======================================================================
575 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
576 {
577 s.resize(2);
578
579 switch (j)
580 {
581 case 0:
582 s[0] = 1.0;
583 s[1] = 0.0;
584 break;
585
586 case 1:
587 s[0] = 0.0;
588 s[1] = 1.0;
589 break;
590
591 case 2:
592 s[0] = 0.0;
593 s[1] = 0.0;
594 break;
595
596 case 3:
597 s[0] = 0.5;
598 s[1] = 0.5;
599 break;
600
601 case 4:
602 s[0] = 0.0;
603 s[1] = 0.5;
604 break;
605
606 case 5:
607 s[0] = 0.5;
608 s[1] = 0.0;
609 break;
610
611 default:
612 std::ostringstream error_message;
613 error_message
614 << "Element only has six nodes; called with node number " << j
615 << std::endl;
616
617 throw OomphLibError(error_message.str(),
618 OOMPH_CURRENT_FUNCTION,
619 OOMPH_EXCEPTION_LOCATION);
620 }
621 }
622
623
624 //=======================================================================
625 /// Shape function for specific TElement<2,3>
626 //=======================================================================
627 void shape(const Vector<double>& s, Shape& psi) const
628 {
629 // Reconstruct the third area coordinate
630 double s_2 = 1.0 - s[0] - s[1];
631
632 // note that s[2] needs replacing by s_2 everywhere since only
633 // two independent variables s[0],s[1] and s_2 is expressed in terms of
634 // those later.
635 psi[0] = 2.0 * s[0] * (s[0] - 0.5);
636 psi[1] = 2.0 * s[1] * (s[1] - 0.5);
637 psi[2] = 2.0 * s_2 * (s_2 - 0.5);
638 psi[3] = 4.0 * s[0] * s[1];
639 psi[4] = 4.0 * s[1] * s_2;
640 psi[5] = 4.0 * s_2 * s[0];
641 }
642
643
644 //=======================================================================
645 /// Derivatives of shape functions for specific TElement<2,3>
646 //=======================================================================
647 void dshape_local(const Vector<double>& s, Shape& psi, DShape& dpsids) const
648 {
649 // ALH: Don't know why object qualifier is needed
650 this->shape(s, psi);
651
652 dpsids(0, 0) = 4.0 * s[0] - 1.0;
653 dpsids(0, 1) = 0.0;
654 dpsids(1, 0) = 0.0;
655 dpsids(1, 1) = 4.0 * s[1] - 1.0;
656 dpsids(2, 0) = 2.0 * (2.0 * s[0] - 1.5 + 2.0 * s[1]);
657 dpsids(2, 1) = 2.0 * (2.0 * s[0] - 1.5 + 2.0 * s[1]);
658 dpsids(3, 0) = 4.0 * s[1];
659 dpsids(3, 1) = 4.0 * s[0];
660 dpsids(4, 0) = -4.0 * s[1];
661 dpsids(4, 1) = 4.0 * (1.0 - s[0] - 2.0 * s[1]);
662 dpsids(5, 0) = 4.0 * (1.0 - 2.0 * s[0] - s[1]);
663 dpsids(5, 1) = -4.0 * s[0];
664 }
665
666
667 //=======================================================================
668 /// Second derivatives of shape functions for specific TElement<2,3>:
669 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
670 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
671 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
672 //=======================================================================
673 void d2shape_local(const Vector<double>& s,
674 Shape& psi,
675 DShape& dpsids,
676 DShape& d2psids) const
677 {
678 // ALH: Don't know why object qualifier is needed
679 this->dshape_local(s, psi, dpsids);
680
681 d2psids(0, 0) = 4.0;
682 d2psids(0, 1) = 0.0;
683 d2psids(0, 2) = 0.0;
684
685 d2psids(1, 0) = 0.0;
686 d2psids(1, 1) = 4.0;
687 d2psids(1, 2) = 0.0;
688
689 d2psids(2, 0) = 4.0;
690 d2psids(2, 1) = 4.0;
691 d2psids(2, 2) = 4.0;
692
693 d2psids(3, 0) = 0.0;
694 d2psids(3, 1) = 0.0;
695 d2psids(3, 2) = 4.0;
696
697 d2psids(4, 0) = 0.0;
698 d2psids(4, 1) = -8.0;
699 d2psids(4, 2) = -4.0;
700
701 d2psids(5, 0) = -8.0;
702 d2psids(5, 1) = 0.0;
703 d2psids(5, 2) = -4.0;
704 }
705 };
706
707 template<>
708 class TElementShape<2, 4>
709 {
710 public:
711 //=======================================================================
712 /// Return local coordinates of node j
713 //=======================================================================
714 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
715 {
716 s.resize(2);
717
718 switch (j)
719 {
720 case 0:
721 s[0] = 1.0;
722 s[1] = 0.0;
723 break;
724
725 case 1:
726 s[0] = 0.0;
727 s[1] = 1.0;
728 break;
729
730 case 2:
731 s[0] = 0.0;
732 s[1] = 0.0;
733 break;
734
735 case 3:
736 s[0] = 2.0 / 3.0;
737 s[1] = 1.0 / 3.0;
738 break;
739
740 case 4:
741 s[0] = 1.0 / 3.0;
742 s[1] = 2.0 / 3.0;
743 break;
744
745 case 5:
746 s[0] = 0.0;
747 s[1] = 2.0 / 3.0;
748 break;
749
750 case 6:
751 s[0] = 0.0;
752 s[1] = 1.0 / 3.0;
753 break;
754
755 case 8:
756 s[0] = 2.0 / 3.0;
757 s[1] = 0.0;
758 break;
759
760 case 7:
761 s[0] = 1.0 / 3.0;
762 s[1] = 0.0;
763 break;
764
765 case 9:
766 s[0] = 1.0 / 3.0;
767 s[1] = 1.0 / 3.0;
768 break;
769
770 default:
771 std::ostringstream error_message;
772 error_message
773 << "Element only has ten nodes; called with node number " << j
774 << std::endl;
775
776 throw OomphLibError(error_message.str(),
777 OOMPH_CURRENT_FUNCTION,
778 OOMPH_EXCEPTION_LOCATION);
779 }
780 }
781
782
783 //=======================================================================
784 /// Shape function for specific TElement<2,4>
785 //=======================================================================
786 void shape(const Vector<double>& s, Shape& psi) const
787 {
788 psi[0] = 0.5 * s[0] * (3.0 * s[0] - 2.0) * (3.0 * s[0] - 1.0);
789 psi[1] = 0.5 * s[1] * (3.0 * s[1] - 2.0) * (3.0 * s[1] - 1.0);
790 psi[2] = 0.5 * (1.0 - s[0] - s[1]) * (1.0 - 3.0 * s[0] - 3.0 * s[1]) *
791 (2.0 - 3.0 * s[0] - 3.0 * s[1]);
792 psi[3] = 4.5 * s[0] * s[1] * (3.0 * s[0] - 1.0);
793 psi[4] = 4.5 * s[0] * s[1] * (3.0 * s[1] - 1.0);
794 psi[5] = 4.5 * s[1] * (1.0 - s[0] - s[1]) * (3.0 * s[1] - 1.0);
795 psi[6] =
796 4.5 * s[1] * (1.0 - s[0] - s[1]) * (3.0 * (1.0 - s[0] - s[1]) - 1.0);
797 psi[7] = 4.5 * s[0] * (1.0 - s[0] - s[1]) * (2.0 - 3 * s[0] - 3 * s[1]);
798 psi[8] = 4.5 * s[0] * (1.0 - s[0] - s[1]) * (3.0 * s[0] - 1.0);
799 psi[9] = 27.0 * s[0] * s[1] * (1.0 - s[0] - s[1]);
800 }
801
802 //=======================================================================
803 /// Derivatives of shape functions for specific TElement<2,4>
804 //=======================================================================
805 void dshape_local(const Vector<double>& s, Shape& psi, DShape& dpsids) const
806 {
807 // ALH: Don't know why object qualifier is needed
808 this->shape(s, psi);
809
810 dpsids(0, 0) = 13.5 * s[0] * s[0] - 9.0 * s[0] + 1.0;
811 dpsids(0, 1) = 0.0;
812 dpsids(1, 0) = 0.0;
813 dpsids(1, 1) = 13.5 * s[1] * s[1] - 9.0 * s[1] + 1.0;
814 dpsids(2, 0) = 0.5 * (36.0 * s[0] + 36.0 * s[1] - 27.0 * s[0] * s[0] -
815 27.0 * s[1] * s[1] - 54.0 * s[0] * s[1] - 11.0);
816 dpsids(2, 1) = 0.5 * (36.0 * s[0] + 36.0 * s[1] - 27.0 * s[0] * s[0] -
817 27.0 * s[1] * s[1] - 54.0 * s[0] * s[1] - 11.0);
818 dpsids(3, 0) = 27.0 * s[0] * s[1] - 4.5 * s[1];
819 dpsids(3, 1) = 4.5 * s[0] * (3.0 * s[0] - 1.0);
820 dpsids(4, 0) = 4.5 * s[1] * (3.0 * s[1] - 1.0);
821 dpsids(4, 1) = 27.0 * s[0] * s[1] - 4.5 * s[0];
822 dpsids(5, 0) = 4.5 * (s[1] - 3.0 * s[1] * s[1]);
823 dpsids(5, 1) =
824 4.5 * (s[0] - 6.0 * s[0] * s[1] - 9.0 * s[1] * s[1] + 8 * s[1] - 1.0);
825 dpsids(6, 0) = 4.5 * (6.0 * s[0] * s[1] - 5.0 * s[1] + 6.0 * s[1] * s[1]);
826 dpsids(6, 1) =
827 4.5 * (2.0 - 5.0 * s[0] + 3.0 * s[0] * s[0] + 12.0 * s[0] * s[1] -
828 10.0 * s[1] + 9.0 * s[1] * s[1]);
829 dpsids(7, 0) =
830 4.5 * (2.0 - 10.0 * s[0] + 9.0 * s[0] * s[0] + 12.0 * s[0] * s[1] -
831 5.0 * s[1] + 3.0 * s[1] * s[1]);
832 dpsids(7, 1) = 4.5 * (6.0 * s[0] * s[0] - 5.0 * s[0] + 6.0 * s[0] * s[1]);
833 dpsids(8, 0) =
834 4.5 * (s[1] - 6.0 * s[0] * s[1] - 9.0 * s[0] * s[0] + 8 * s[0] - 1.0);
835 dpsids(8, 1) = 4.5 * (s[0] - 3.0 * s[0] * s[0]);
836 dpsids(9, 0) = 27.0 * s[1] - 54.0 * s[0] * s[1] - 27.0 * s[1] * s[1];
837 dpsids(9, 1) = 27.0 * s[0] - 54.0 * s[0] * s[1] - 27.0 * s[0] * s[0];
838 }
839
840 //=======================================================================
841 /// Second derivatives of shape functions for specific TElement<2,4>:
842 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
843 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
844 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
845 //=======================================================================
846 void d2shape_local(const Vector<double>& s,
847 Shape& psi,
848 DShape& dpsids,
849 DShape& d2psids) const
850 {
851 throw OomphLibError(
852 "Not checked yet", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
853
854 // ALH: Don't know why object qualifier is needed
855 this->dshape_local(s, psi, dpsids);
856
857 d2psids(0, 0) = 9.0 * (3.0 * s[0] - 1.0);
858 d2psids(0, 1) = 0.0;
859 d2psids(0, 2) = 0.0;
860 d2psids(1, 0) = 0.0;
861 d2psids(1, 1) = 0.0;
862 d2psids(1, 2) = 9.0 * (3.0 * s[1] - 1.0);
863 d2psids(2, 0) = 9.0 * (2.0 - 3.0 * s[0] - 3.0 * s[1]);
864 d2psids(2, 1) = 9.0 * (2.0 - 3.0 * s[0] - 3.0 * s[1]);
865 d2psids(2, 2) = 9.0 * (2.0 - 3.0 * s[0] - 3.0 * s[1]);
866 d2psids(3, 0) = 27.0 * s[1];
867 d2psids(3, 1) = 0.0;
868 d2psids(3, 2) = 27.0 * s[0] - 4.5;
869 d2psids(4, 0) = 0.0;
870 d2psids(4, 1) = 27.0 * s[0];
871 d2psids(4, 2) = 27.0 * s[1] - 4.5;
872 d2psids(5, 0) = 0.0;
873 d2psids(5, 1) = 9.0 * (4.0 - 3.0 * s[0] - 9.0 * s[1]);
874 d2psids(5, 2) = 4.5 * (1.0 - 6.0 * s[1]);
875 d2psids(6, 0) = 27.0 * s[1];
876 d2psids(6, 1) = 9.0 * (6.0 * s[0] + 9.0 * s[1] - 5.0);
877 d2psids(6, 2) = 4.5 * (6.0 * s[0] + 12.0 * s[1] - 5.0);
878 d2psids(8, 0) = 9.0 * (4.0 - 9.0 * s[0] - 3.0 * s[1]);
879 d2psids(8, 1) = 0.0;
880 d2psids(8, 2) = 4.5 * (1.0 - 6.0 * s[0]);
881 d2psids(7, 0) = 9.0 * (9.0 * s[0] + 6.0 * s[1] - 5.0);
882 d2psids(7, 1) = 27.0 * s[0];
883 d2psids(7, 2) = 4.5 * (12.0 * s[0] + 6.0 * s[1] - 5.0);
884 d2psids(9, 0) = -54.0 * s[1];
885 d2psids(9, 1) = -54.0 * s[0];
886 d2psids(9, 2) = 27.0 - 54.0 * s[0] - 54.0 * s[1];
887 }
888 };
889
890
891 //========================================================================
892 /// A class for those member functions that must be fully specialised
893 /// for Telements that are enriched by bubbble functions.
894 /// The fact that member functions of partially
895 /// specialised classes cannot necessarily be fully specialised
896 /// means that we must either fully specialise every class, or use this
897 /// base class to fully specialize only those functions that are required.
898 //========================================================================
899 template<unsigned DIM, unsigned NNODE_1D>
900 class TBubbleEnrichedElementShape
901 {
902 };
903
904
905 ///////////////////////////////////////////////////////////////////////
906 /// Specific Enriched TElementShape inline functions
907 //////////////////////////////////////////////////////////////////////
908
909 //===============================================================
910 /// Standard quadratic shape functions enriched by the addition
911 /// of a cubic bubble, which consists of adding a single node
912 /// at the centroid
913 //=============================================================
914 template<>
915 class TBubbleEnrichedElementShape<2, 3>
916 {
917 public:
918 //=====================================================================
919 /// Return the number of nodes required for enrichement
920 //====================================================================
921 unsigned n_enriched_nodes()
922 {
923 return 1;
924 }
925
926 //=======================================================================
927 /// Return local coordinates of node j
928 //=======================================================================
929 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
930 {
931 s.resize(2);
932
933 switch (j)
934 {
935 case 0:
936 s[0] = 1.0;
937 s[1] = 0.0;
938 break;
939
940 case 1:
941 s[0] = 0.0;
942 s[1] = 1.0;
943 break;
944
945 case 2:
946 s[0] = 0.0;
947 s[1] = 0.0;
948 break;
949
950 case 3:
951 s[0] = 0.5;
952 s[1] = 0.5;
953 break;
954
955 case 4:
956 s[0] = 0.0;
957 s[1] = 0.5;
958 break;
959
960 case 5:
961 s[0] = 0.5;
962 s[1] = 0.0;
963 break;
964
965 // Add the centroid as the enriched node
966 case 6:
967 s[0] = 1.0 / 3.0;
968 s[1] = 1.0 / 3.0;
969 break;
970
971 default:
972 std::ostringstream error_message;
973 error_message
974 << "Element only has seven nodes; called with node number " << j
975 << std::endl;
976
977 throw OomphLibError(error_message.str(),
978 OOMPH_CURRENT_FUNCTION,
979 OOMPH_EXCEPTION_LOCATION);
980 }
981 }
982
983
984 //=======================================================================
985 /// Shape function for specific TBubbleEnrichedElement<2,3>
986 //=======================================================================
987 void shape(const Vector<double>& s, Shape& psi) const
988 {
989 // Reconstruct the third area coordinate
990 const double s_2 = 1.0 - s[0] - s[1];
991
992 // Calculate the enrichment function
993 const double cubic_bubble = s[0] * s[1] * s_2;
994 // The appropriate amount of the cubic bubble function is
995 // added/subtracted to each original quadratic shape function to ensure
996 // that it is zero at the centroid (1/3,1/3).
997
998 // note that s[2] needs replacing by s_2 everywhere since only
999 // two independent variables s[0],s[1] and s_2 is expressed in terms of
1000 // those later.
1001 psi[0] = 2.0 * s[0] * (s[0] - 0.5) + 3.0 * cubic_bubble;
1002 psi[1] = 2.0 * s[1] * (s[1] - 0.5) + 3.0 * cubic_bubble;
1003 psi[2] = 2.0 * s_2 * (s_2 - 0.5) + 3.0 * cubic_bubble;
1004 psi[3] = 4.0 * s[0] * s[1] - 12.0 * cubic_bubble;
1005 psi[4] = 4.0 * s[1] * s_2 - 12.0 * cubic_bubble;
1006 psi[5] = 4.0 * s_2 * s[0] - 12.0 * cubic_bubble;
1007 // The bubble function scaled to have magnitude one at (1/3,1/3)
1008 psi[6] = 27.0 * cubic_bubble;
1009 }
1010
1011
1012 //=======================================================================
1013 /// Derivatives of shape functions for specific TBubbleElement<2,3>
1014 //=======================================================================
1015 void dshape_local(const Vector<double>& s, Shape& psi, DShape& dpsids) const
1016 {
1017 // ALH: Don't know why object qualifier is needed
1018 this->shape(s, psi);
1019
1020 // Calculate derivatives of bubble functions
1021 const double d_bubble_ds0 = s[1] * (1.0 - s[1] - 2.0 * s[0]);
1022 const double d_bubble_ds1 = s[0] * (1.0 - s[0] - 2.0 * s[1]);
1023
1024 // Add the appropriate derivatives to the shape functions
1025 dpsids(0, 0) = 4.0 * s[0] - 1.0 + 3.0 * d_bubble_ds0;
1026 dpsids(0, 1) = 0.0 + 3.0 * d_bubble_ds1;
1027 dpsids(1, 0) = 0.0 + 3.0 * d_bubble_ds0;
1028 dpsids(1, 1) = 4.0 * s[1] - 1.0 + 3.0 * d_bubble_ds1;
1029 dpsids(2, 0) = 2.0 * (2.0 * s[0] - 1.5 + 2.0 * s[1]) + 3.0 * d_bubble_ds0;
1030 dpsids(2, 1) = 2.0 * (2.0 * s[0] - 1.5 + 2.0 * s[1]) + 3.0 * d_bubble_ds1;
1031 dpsids(3, 0) = 4.0 * s[1] - 12.0 * d_bubble_ds0;
1032 dpsids(3, 1) = 4.0 * s[0] - 12.0 * d_bubble_ds1;
1033 dpsids(4, 0) = -4.0 * s[1] - 12.0 * d_bubble_ds0;
1034 dpsids(4, 1) = 4.0 * (1.0 - s[0] - 2.0 * s[1]) - 12.0 * d_bubble_ds1;
1035 dpsids(5, 0) = 4.0 * (1.0 - 2.0 * s[0] - s[1]) - 12.0 * d_bubble_ds0;
1036 dpsids(5, 1) = -4.0 * s[0] - 12.0 * d_bubble_ds1;
1037 dpsids(6, 0) = 27.0 * d_bubble_ds0;
1038 dpsids(6, 1) = 27.0 * d_bubble_ds1;
1039 }
1040
1041
1042 //=======================================================================
1043 /// Second derivatives of shape functions for specific
1044 /// TBubbleElement<2,3>:
1045 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
1046 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
1047 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
1048 //=======================================================================
1049 void d2shape_local(const Vector<double>& s,
1050 Shape& psi,
1051 DShape& dpsids,
1052 DShape& d2psids) const
1053 {
1054 // ALH: Don't know why object qualifier is needed
1055 this->dshape_local(s, psi, dpsids);
1056
1057 // Calculate derivatives of bubble functions
1058 const double d2_bubble_ds0 = -2.0 * s[1];
1059 const double d2_bubble_ds1 = -2.0 * s[0];
1060 const double d2_bubble_ds2 = 1.0 - 2.0 * s[0] - 2.0 * s[1];
1061
1062 d2psids(0, 0) = 4.0 + 3.0 * d2_bubble_ds0;
1063 d2psids(0, 1) = 0.0 + 3.0 * d2_bubble_ds1;
1064 d2psids(0, 2) = 0.0 + 3.0 * d2_bubble_ds2;
1065
1066 d2psids(1, 0) = 0.0 + 3.0 * d2_bubble_ds0;
1067 d2psids(1, 1) = 4.0 + 3.0 * d2_bubble_ds1;
1068 d2psids(1, 2) = 0.0 + 3.0 * d2_bubble_ds2;
1069
1070 d2psids(2, 0) = 4.0 + 3.0 * d2_bubble_ds0;
1071 d2psids(2, 1) = 4.0 + 3.0 * d2_bubble_ds1;
1072 d2psids(2, 2) = 4.0 + 3.0 * d2_bubble_ds2;
1073
1074 d2psids(3, 0) = 0.0 - 12.0 * d2_bubble_ds0;
1075 d2psids(3, 1) = 0.0 - 12.0 * d2_bubble_ds1;
1076 d2psids(3, 2) = 4.0 - 12.0 * d2_bubble_ds2;
1077
1078 d2psids(4, 0) = 0.0 - 12.0 * d2_bubble_ds0;
1079 d2psids(4, 1) = -8.0 - 12.0 * d2_bubble_ds1;
1080 d2psids(4, 2) = -4.0 - 12.0 * d2_bubble_ds2;
1081
1082 d2psids(5, 0) = -8.0 - 12.0 * d2_bubble_ds0;
1083 d2psids(5, 1) = 0.0 - 12.0 * d2_bubble_ds1;
1084 d2psids(5, 2) = -4.0 - 12.0 * d2_bubble_ds2;
1085
1086 d2psids(6, 0) = 27.0 * d2_bubble_ds0;
1087 d2psids(6, 1) = 27.0 * d2_bubble_ds1;
1088 d2psids(6, 2) = 27.0 * d2_bubble_ds2;
1089 }
1090 };
1091
1092 //////////////////////////////////////////////////////////////////////
1093 //////////////////////////////////////////////////////////////////////
1094 //////////////////////////////////////////////////////////////////////
1095
1096 //========================================================================
1097 /// Empty base class for Telements (created so that
1098 /// we can use dynamic_cast<>() to figure out if a an element
1099 /// is a Telement (from a purely geometric point of view).
1100 //========================================================================
1101 class TElementGeometricBase : public virtual FiniteElement
1102 {
1103 public:
1104 /// Empty default constructor
1105 TElementGeometricBase() {}
1106
1107 /// Broken copy constructor
1108 TElementGeometricBase(const TElementGeometricBase&) = delete;
1109
1110 /// Broken assignment operator
1111 // Commented out broken assignment operator because this can lead to a
1112 // conflict warning when used in the virtual inheritence hierarchy.
1113 // Essentially the compiler doesn't realise that two separate
1114 // implementations of the broken function are the same and so, quite
1115 // rightly, it shouts.
1116 /*void operator=(const TElementGeometricBase&) = delete;*/
1117 };
1118
1119
1120 //////////////////////////////////////////////////////////////////////
1121 //////////////////////////////////////////////////////////////////////
1122 //////////////////////////////////////////////////////////////////////
1123
1124 //========================================================================
1125 /// Empty base class for Telements (created so that
1126 /// we can use dynamic_cast<>() to figure out if a an element
1127 /// is a Telement).
1128 //========================================================================
1129 class TElementBase : public virtual TElementGeometricBase
1130 {
1131 public:
1132 /// Empty default constructor
1133 TElementBase() {}
1134
1135 /// Broken copy constructor
1136 TElementBase(const TElementBase&) = delete;
1137
1138 /// Broken assignment operator
1139 /*void operator=(const TElementBase&) = delete;*/
1140
1141 /// It's a T element!
1146
1147 /// Check whether the local coordinates are valid or not
1148 bool local_coord_is_valid(const Vector<double>& s)
1149 {
1150 // Check coordinates
1151 unsigned ncoord = dim();
1152 double sum = 0.0;
1153 for (unsigned i = 0; i < ncoord; i++)
1154 {
1155 // Each local coordinate must be positive
1156 if (s[i] < 0.0)
1157 {
1158 return false;
1159 }
1160 sum += s[i];
1161 }
1162
1163 // Sum must be less than 1
1164 if (sum <= 1.0)
1165 {
1166 return true;
1167 }
1168
1169 // We're outside...
1170 return false;
1171 }
1172
1173 /// Adjust local coordinates so that they're located inside
1174 /// the element
1175 void move_local_coord_back_into_element(Vector<double>& s) const
1176 {
1177 // Check coordinates
1178 unsigned ncoord = dim();
1179 double sum = 0.0;
1180 for (unsigned i = 0; i < ncoord; i++)
1181 {
1182 // Each coordinate must be positive individually
1183 if (s[i] < 0.0) s[i] = 0.0;
1184 sum += s[i];
1185 }
1186
1187 // Sum must be less than 1
1188 double excess = sum - 1.0;
1189 if (excess > 0.0)
1190 {
1191 // Subtract excess equally from all coordinates
1192 double sub = excess / double(ncoord);
1193 for (unsigned i = 0; i < ncoord; i++)
1194 {
1195 s[i] -= sub;
1196 }
1197 }
1198 }
1199 };
1200
1201 //=======================================================================
1202 /// General TElement class
1203 ///
1204 /// Empty, just establishes the template parameters
1205 //=======================================================================
1206 template<unsigned DIM, unsigned NNODE_1D>
1207 class TElement
1208 {
1209 };
1210
1211
1212 //=======================================================================
1213 /// General TElement class specialised to one spatial dimensions
1214 /// Ordering of nodes is 0 at local coordinate s[0] = 0, 1 at local
1215 /// coordinate s[0] = 1 and then filling in the intermediate values
1216 /// from s[0]=0 to 1.
1217 //=======================================================================
1218 template<unsigned NNODE_1D>
1219 class TElement<1, NNODE_1D> : public virtual TElementBase,
1220 public TElementShape<1, NNODE_1D>
1221 {
1222 private:
1223 /// Default integration rule: Gaussian integration of same 'order' as
1224 /// the element
1225 // This is sort of optimal, because it means that the integration is exact
1226 // for the shape functions. Can overwrite this in specific element
1227 // defintion.
1228 static TGauss<1, NNODE_1D> Default_integration_scheme;
1229
1230 public:
1231 /// Constructor
1232 TElement()
1233 {
1234 // Number of nodes
1235 switch (NNODE_1D)
1236 {
1237 case 2:
1238 case 3:
1239 case 4:
1240 break;
1241
1242 default:
1243 std::string error_message =
1244 "One-dimensional TElements are currently only implemented for\n";
1245 error_message += "three and six nodes, i.e. NNODE_1D=2 , 3 , 4\n";
1246
1247 throw OomphLibError(
1248 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1249 }
1250
1251 // Set the number of nodes
1252 unsigned n_node = NNODE_1D;
1253 this->set_n_node(n_node);
1254
1255 // Set the elemental and nodal dimension
1256 set_dimension(1);
1257
1258 // Assign default (full) spatial integration scheme
1259 set_integration_scheme(&Default_integration_scheme);
1260 }
1261
1262
1263 /// Broken copy constructor
1264 TElement(const TElement&) = delete;
1265
1266 /// Broken assignment operator
1267 /*void operator=(const TElement&) = delete;*/
1268
1269
1270 /// Destructor
1271 ~TElement() {}
1272
1273 /// Number of nodes along each element edge
1274 unsigned nnode_1d() const
1275 {
1276 return NNODE_1D;
1277 }
1278
1279
1280 /// Number of vertex nodes in the element: One more
1281 /// than spatial dimension
1282 unsigned nvertex_node() const
1283 {
1284 return 2;
1285 }
1286
1287 /// Pointer to the j-th vertex node in the element
1288 Node* vertex_node_pt(const unsigned& j) const
1289 {
1290 switch (j)
1291 {
1292 case 0:
1293
1294 return node_pt(0);
1295 break;
1296
1297 case 1:
1298
1299 return node_pt(NNODE_1D - 1);
1300 break;
1301
1302 default:
1303
1304 std::ostringstream error_message;
1305 error_message
1306 << "Element only has two vertex nodes; called with node number "
1307 << j << std::endl;
1308 throw OomphLibError(error_message.str(),
1309 OOMPH_CURRENT_FUNCTION,
1310 OOMPH_EXCEPTION_LOCATION);
1311 }
1312 }
1313
1314 /// Calculate the geometric shape functions at local coordinate s
1315 inline void shape(const Vector<double>& s, Shape& psi) const
1316 {
1317 TElementShape<1, NNODE_1D>::shape(s, psi);
1318 }
1319
1320 /// Compute the geometric shape functions and
1321 /// derivatives w.r.t. local coordinates at local coordinate s
1322 inline void dshape_local(const Vector<double>& s,
1323 Shape& psi,
1324 DShape& dpsids) const
1325 {
1326 TElementShape<1, NNODE_1D>::dshape_local(s, psi, dpsids);
1327 }
1328
1329 /// Compute the geometric shape functions, derivatives and
1330 /// second derivatives w.r.t local coordinates at local coordinate s
1331 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
1332 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
1333 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
1334 inline void d2shape_local(const Vector<double>& s,
1335 Shape& psi,
1336 DShape& dpsids,
1337 DShape& d2psids) const
1338 {
1339 TElementShape<1, NNODE_1D>::d2shape_local(s, psi, dpsids, d2psids);
1340 }
1341
1342 /// Overload the template-free interface for the calculation of
1343 /// inverse jacobian matrix. This is a one dimensional element, so use
1344 /// the 1D version.
1345 double invert_jacobian_mapping(const DenseMatrix<double>& jacobian,
1346 DenseMatrix<double>& inverse_jacobian) const
1347 {
1348 return FiniteElement::invert_jacobian<1>(jacobian, inverse_jacobian);
1349 }
1350
1351 /// Min. value of local coordinate
1352 double s_min() const
1353 {
1354 return 0.0;
1355 }
1356
1357 /// Max. value of local coordinate
1358 double s_max() const
1359 {
1360 return 1.0;
1361 }
1362
1363 /// Return local coordinates of node j
1364 inline void local_coordinate_of_node(const unsigned& j,
1365 Vector<double>& s) const
1366 {
1367 TElementShape<1, NNODE_1D>::local_coordinate_of_node(j, s);
1368 }
1369
1370 /// Return the number of actual plot points for paraview
1371 /// plot with parameter nplot.
1372 unsigned nplot_points_paraview(const unsigned& nplot) const
1373 {
1374 return nplot;
1375 }
1376
1377 /// Return the number of local sub-elements for paraview plot with
1378 /// parameter nplot.
1379 unsigned nsub_elements_paraview(const unsigned& nplot) const
1380 {
1381 return (nplot - 1);
1382 }
1383
1384 /// Fill in the offset information for paraview plot.
1385 /// Needs to be implemented for each new geometric element type; see
1386 /// http://www.vtk.org/VTK/img/file-formats.pdf
1387 void write_paraview_output_offset_information(std::ofstream& file_out,
1388 const unsigned& nplot,
1389 unsigned& counter) const
1390 {
1391 // Number of local elements we want to plot over
1392 unsigned plot = nsub_elements_paraview(nplot);
1393
1394 // loops over the i-th local element in parent element
1395 for (unsigned i = 0; i < plot; i++)
1396 {
1397 file_out << i + counter << " " << i + 1 + counter << std::endl;
1398 }
1399 counter += nplot_points_paraview(nplot);
1400 }
1401
1402 /// Return the paraview element type.
1403 /// Needs to be implemented for each new geometric element type; see
1404 /// http://www.vtk.org/VTK/img/file-formats.pdf
1405 void write_paraview_type(std::ofstream& file_out,
1406 const unsigned& nplot) const
1407 {
1408 unsigned local_loop = nsub_elements_paraview(nplot);
1409 for (unsigned i = 0; i < local_loop; i++)
1410 {
1411 file_out << "3" << std::endl;
1412 }
1413 }
1414
1415 /// Return the offsets for the paraview sub-elements. Needs
1416 /// to be implemented for each new geometric element type; see
1417 /// http://www.vtk.org/VTK/img/file-formats.pdf
1418 void write_paraview_offsets(std::ofstream& file_out,
1419 const unsigned& nplot,
1420 unsigned& offset_sum) const
1421 {
1422 // Loop over all local elements and add its offset to the overall
1423 // offset_sum
1424 unsigned local_loop = nsub_elements_paraview(nplot);
1425 for (unsigned i = 0; i < local_loop; i++)
1426 {
1427 offset_sum += 2;
1428 file_out << offset_sum << std::endl;
1429 }
1430 }
1431
1432 /// Output
1433 void output(std::ostream& output);
1434
1435 /// Output at specified number of plot points
1436 void output(std::ostream& outfile, const unsigned& nplot);
1437
1438 /// C-style output
1439 void output(FILE* file_pt);
1440
1441 /// C_style output at n_plot points
1442 void output(FILE* file_pt, const unsigned& n_plot);
1443
1444 /// Get vector of local coordinates of plot point i (when plotting
1445 /// nplot points in each "coordinate direction").
1446 void get_s_plot(
1447 const unsigned& i,
1448 const unsigned& nplot,
1449 Vector<double>& s,
1450 const bool& use_equally_spaced_interior_sample_points = false) const
1451 {
1452 if (nplot > 1)
1453 {
1454 s[0] = double(i) / double(nplot - 1);
1455
1456 if (use_equally_spaced_interior_sample_points)
1457 {
1458 double range = 1.0;
1459 double dx_new = range / double(nplot);
1460 double range_new = double(nplot - 1) * dx_new;
1461 s[0] = 0.5 * dx_new + range_new * s[0] / range;
1462 }
1463 }
1464 else
1465 {
1466 s[0] = 0.5;
1467 }
1468 }
1469
1470 /// Return string for tecplot zone header (when plotting
1471 /// nplot points in each "coordinate direction)
1472 std::string tecplot_zone_string(const unsigned& nplot) const
1473 {
1474 std::ostringstream header;
1475 header << "ZONE I=" << nplot << "\n";
1476 return header.str();
1477 }
1478
1479 /// Return total number of plot points (when plotting
1480 /// nplot points in each "coordinate direction)
1481 unsigned nplot_points(const unsigned& nplot) const
1482 {
1483 return nplot;
1484 }
1485
1486 /// Build the lower-dimensional FaceElement (an element of type
1487 /// PointElement). The face index takes two values
1488 /// corresponding to the two possible faces:
1489 /// -1 (Left) s[0] = -1.0
1490 /// +1 (Right) s[0] = 1.0
1491 void build_face_element(const int& face_index,
1492 FaceElement* face_element_pt);
1493 };
1494
1495
1496 //=======================================================================
1497 /// General TElement class specialised to two spatial dimensions
1498 /// Ordering of nodes as in Zienkiwizc sketches: vertex nodes
1499 /// 0 - 1 - 2 anticlockwise. Midside nodes filled in progressing
1500 /// along the consecutive edges. Central node(s) come(s) last.
1501 //=======================================================================
1502 template<unsigned NNODE_1D>
1503 class TElement<2, NNODE_1D> : public virtual TElementBase,
1504 public TElementShape<2, NNODE_1D>
1505 {
1506 private:
1507 /// Nodal translation scheme for use when generating face elements
1508 static const unsigned Node_on_face[3][NNODE_1D];
1509
1510 /// Default integration rule: Gaussian integration of same 'order' as
1511 /// the element
1512 // This is sort of optimal, because it means that the integration is exact
1513 // for the shape functions. Can overwrite this in specific element
1514 // defintion.
1515 static TGauss<2, NNODE_1D> Default_integration_scheme;
1516
1517
1518 public:
1519 /// Constructor
1520 TElement()
1521 {
1522 // Number of nodes
1523 switch (NNODE_1D)
1524 {
1525 case 2:
1526 case 3:
1527 case 4:
1528 break;
1529
1530 default:
1531 std::string error_message =
1532 "Triangles are currently only implemented for\n";
1533 error_message += "three and six nodes, i.e. NNODE_1D=2 , 3 , 4\n";
1534
1535 throw OomphLibError(
1536 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1537 }
1538
1539 // Set the number of nodes
1540 unsigned n_node = (NNODE_1D * (NNODE_1D + 1)) / 2;
1541 this->set_n_node(n_node);
1542
1543 // Set the elemental and nodal dimension
1544 set_dimension(2);
1545
1546 // Assign default (full) spatial integration scheme
1547 set_integration_scheme(&Default_integration_scheme);
1548 }
1549
1550 /// Alternative constructor
1551 TElement(const bool& allow_high_order)
1552 {
1553 // Check if we are overriding the restriction on NNODE_1D
1554 if (!allow_high_order)
1555 {
1556 // Call the default constructor
1557 TElement<2, NNODE_1D>();
1558 }
1559 else
1560 {
1561 // Set the number of nodes
1562 unsigned n_node = (NNODE_1D * (NNODE_1D + 1)) / 2;
1563 this->set_n_node(n_node);
1564
1565 // Set the elemental and nodal dimension
1566 set_dimension(2);
1567
1568 // Assign default (full) spatial integration scheme
1569 set_integration_scheme(&Default_integration_scheme);
1570 }
1571 }
1572
1573
1574 /// Broken copy constructor
1575 TElement(const TElement&) = delete;
1576
1577 /// Broken assignment operator
1578 /*void operator=(const TElement&) = delete;*/
1579
1580
1581 /// Destructor
1582 ~TElement() {}
1583
1584 /// Number of nodes along each element edge
1585 unsigned nnode_1d() const
1586 {
1587 return NNODE_1D;
1588 }
1589
1590 /// Number of vertex nodes in the element: One more
1591 /// than spatial dimension
1592 unsigned nvertex_node() const
1593 {
1594 return 3;
1595 }
1596
1597 /// Public access function for Node_on_face.
1598 unsigned get_bulk_node_number(const int& face_index,
1599 const unsigned& i) const
1600 {
1601 return Node_on_face[face_index][i];
1602 }
1603
1604 /// Pointer to the j-th vertex node in the element
1605 Node* vertex_node_pt(const unsigned& j) const
1606 {
1607 // Vertex nodes come first:
1608#ifdef PARANOID
1609 if (j > 2)
1610 {
1611 std::ostringstream error_message;
1612 error_message
1613 << "Element only has three vertex nodes; called with node number "
1614 << j << std::endl;
1615 throw OomphLibError(error_message.str(),
1616 OOMPH_CURRENT_FUNCTION,
1617 OOMPH_EXCEPTION_LOCATION);
1618 }
1619#endif
1620 return node_pt(j);
1621 }
1622
1623 /// Calculate the geometric shape functions at local coordinate s
1624 inline void shape(const Vector<double>& s, Shape& psi) const
1625 {
1626 TElementShape<2, NNODE_1D>::shape(s, psi);
1627 }
1628
1629 /// Compute the geometric shape functions and
1630 /// derivatives w.r.t. local coordinates at local coordinate s
1631 inline void dshape_local(const Vector<double>& s,
1632 Shape& psi,
1633 DShape& dpsids) const
1634 {
1635 TElementShape<2, NNODE_1D>::dshape_local(s, psi, dpsids);
1636 }
1637
1638 /// Compute the geometric shape functions, derivatives and
1639 /// second derivatives w.r.t local coordinates at local coordinate s
1640 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
1641 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
1642 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
1643 inline void d2shape_local(const Vector<double>& s,
1644 Shape& psi,
1645 DShape& dpsids,
1646 DShape& d2psids) const
1647 {
1648 TElementShape<2, NNODE_1D>::d2shape_local(s, psi, dpsids, d2psids);
1649 }
1650
1651 /// Overload the template-free interface for the calculation of
1652 /// inverse jacobian matrix. This is a two dimensional element, so use
1653 /// the 2D version.
1654 double invert_jacobian_mapping(const DenseMatrix<double>& jacobian,
1655 DenseMatrix<double>& inverse_jacobian) const
1656 {
1657 return FiniteElement::invert_jacobian<2>(jacobian, inverse_jacobian);
1658 }
1659
1660 /// Min. value of local coordinate
1661 double s_min() const
1662 {
1663 return 0.0;
1664 }
1665
1666 /// Max. value of local coordinate
1667 double s_max() const
1668 {
1669 return 1.0;
1670 }
1671
1672 /// Return local coordinates of node j
1673 inline void local_coordinate_of_node(const unsigned& j,
1674 Vector<double>& s) const
1675 {
1676 TElementShape<2, NNODE_1D>::local_coordinate_of_node(j, s);
1677 }
1678
1679 /// Return the number of actual plot points for paraview
1680 /// plot with parameter nplot.
1681 unsigned nplot_points_paraview(const unsigned& nplot) const
1682 {
1683 unsigned node_sum = 0;
1684 for (unsigned i = 1; i <= nplot; i++)
1685 {
1686 node_sum += i;
1687 }
1688 return node_sum;
1689 }
1690
1691 /// Return the number of local sub-elements for paraview plot with
1692 /// parameter nplot.
1693 unsigned nsub_elements_paraview(const unsigned& nplot) const
1694 {
1695 unsigned local_sum = 0;
1696 for (unsigned i = 1; i < nplot; i++)
1697 {
1698 local_sum += 2 * (nplot - i - 1) + 1;
1699 }
1700 return local_sum;
1701 }
1702
1703 /// Fill in the offset information for paraview plot.
1704 /// Needs to be implemented for each new geometric element type; see
1705 /// http://www.vtk.org/VTK/img/file-formats.pdf
1706 void write_paraview_output_offset_information(std::ofstream& file_out,
1707 const unsigned& nplot,
1708 unsigned& counter) const
1709 {
1710 // Outputs list of connectivity of Paraview elements,
1711 // whilst remembering the overall ordering
1712
1713 unsigned node_count = 0;
1714 for (unsigned i = 0; i < nplot - 1; i++)
1715 {
1716 for (unsigned j = 0; j < nplot - i - 1; j++)
1717 {
1718 file_out << j + node_count + counter << " "
1719 << j + node_count + 1 + counter << " "
1720 << j + nplot + node_count - i + counter << std::endl;
1721
1722 if (j < nplot - i - 2)
1723 {
1724 file_out << j + node_count + 1 + counter << " "
1725 << j + nplot + node_count - i + 1 + counter << " "
1726 << j + nplot + node_count - i + counter << std::endl;
1727 }
1728 }
1729 node_count += (nplot - i);
1730 }
1731 counter += nplot_points_paraview(nplot);
1732 }
1733
1734 /// Return the paraview element type.
1735 /// Needs to be implemented for each new geometric element type; see
1736 /// http://www.vtk.org/VTK/img/file-formats.pdf
1737 void write_paraview_type(std::ofstream& file_out,
1738 const unsigned& nplot) const
1739 {
1740 unsigned local_loop = nsub_elements_paraview(nplot);
1741
1742 // Loop over all the local elements and print its paraview type
1743 for (unsigned i = 0; i < local_loop; i++)
1744 {
1745 file_out << "5" << std::endl;
1746 }
1747 }
1748
1749 /// Return the offsets for the paraview sub-elements. Needs
1750 /// to be implemented for each new geometric element type; see
1751 /// http://www.vtk.org/VTK/img/file-formats.pdf
1752 void write_paraview_offsets(std::ofstream& file_out,
1753 const unsigned& nplot,
1754 unsigned& offset_sum) const
1755 {
1756 unsigned local_loop = nsub_elements_paraview(nplot);
1757
1758 // Loop over all local elements and add its offset to the overall
1759 // offset_sum
1760 for (unsigned i = 0; i < local_loop; i++)
1761 {
1762 offset_sum += 3;
1763 file_out << offset_sum << std::endl;
1764 }
1765 }
1766
1767 /// Output
1768 void output(std::ostream& output);
1769
1770 /// Output at specified number of plot points
1771 void output(std::ostream& outfile, const unsigned& nplot);
1772
1773 /// C-style output
1774 void output(FILE* file_pt);
1775
1776 /// C_style output at n_plot points
1777 void output(FILE* file_pt, const unsigned& n_plot);
1778
1779 /// Get vector of local coordinates of plot point i (when plotting
1780 /// nplot points in each "coordinate direction").
1781 void get_s_plot(
1782 const unsigned& iplot,
1783 const unsigned& nplot,
1784 Vector<double>& s,
1785 const bool& use_equally_spaced_interior_sample_points = false) const
1786 {
1787 if (nplot > 1)
1788 {
1789 unsigned np = 0, i, j;
1790 for (i = 0; i < nplot; ++i)
1791 {
1792 for (j = 0; j < nplot - i; ++j)
1793 {
1794 if (np == iplot)
1795 {
1796 s[0] = double(j) / double(nplot - 1);
1797 s[1] = double(i) / double(nplot - 1);
1798 if (use_equally_spaced_interior_sample_points)
1799 {
1800 double range = 1.0;
1801 double dx_new = range / (double(nplot) + 0.5);
1802 double range_new = double(nplot - 1) * dx_new;
1803 s[0] = 0.5 * dx_new + range_new * s[0] / range;
1804 s[1] = 0.5 * dx_new + range_new * s[1] / range;
1805 }
1806 return;
1807 }
1808 ++np;
1809 }
1810 }
1811 }
1812 else
1813 {
1814 s[0] = 1.0 / 3.0;
1815 s[1] = 1.0 / 3.0;
1816 }
1817 }
1818
1819 /// Return string for tecplot zone header (when plotting
1820 /// nplot points in each "coordinate direction)
1821 std::string tecplot_zone_string(const unsigned& nplot) const
1822 {
1823 std::ostringstream header;
1824 unsigned nel = 0;
1825 for (unsigned i = 1; i < nplot; i++)
1826 {
1827 nel += 2 * i - 1;
1828 }
1829 header << "ZONE N=" << nplot_points(nplot) << ", E=" << nel
1830 << ", F=FEPOINT, ET=TRIANGLE\n";
1831 return header.str();
1832 }
1833
1834 /// Add tecplot zone "footer" to output stream (when plotting
1835 /// nplot points in each "coordinate direction).
1836 /// Empty by default -- can be used, e.g., to add FE connectivity
1837 /// lists to elements that need it.
1838 void write_tecplot_zone_footer(std::ostream& outfile,
1839 const unsigned& nplot) const
1840 {
1841 // Output node lists for sub elements for Tecplot (node index
1842 // must start at 1)
1843 unsigned nod_count = 1;
1844 for (unsigned i = 0; i < nplot; i++)
1845 {
1846 for (unsigned j = 0; j < nplot - i; j++)
1847 {
1848 if (j < nplot - i - 1)
1849 {
1850 outfile << nod_count << " " << nod_count + 1 << " "
1851 << nod_count + nplot - i << std::endl;
1852 if (j < nplot - i - 2)
1853 {
1854 outfile << nod_count + 1 << " " << nod_count + nplot - i + 1
1855 << " " << nod_count + nplot - i << std::endl;
1856 }
1857 }
1858 ++nod_count;
1859 }
1860 }
1861 }
1862
1863 /// Add tecplot zone "footer" to C-style output. (when plotting
1864 /// nplot points in each "coordinate direction).
1865 /// Empty by default -- can be used, e.g., to add FE connectivity
1866 /// lists to elements that need it.
1867 void write_tecplot_zone_footer(FILE* file_pt, const unsigned& nplot) const
1868 {
1869 // Output node lists for sub elements for Tecplot (node index
1870 // must start at 1)
1871 unsigned nod_count = 1;
1872 for (unsigned i = 0; i < nplot; i++)
1873 {
1874 for (unsigned j = 0; j < nplot - i; j++)
1875 {
1876 if (j < nplot - i - 1)
1877 {
1878 fprintf(file_pt,
1879 "%i %i %i \n",
1880 nod_count,
1881 nod_count + 1,
1882 nod_count + nplot - i);
1883 if (j < nplot - i - 2)
1884 {
1885 fprintf(file_pt,
1886 "%i %i %i \n",
1887 nod_count + 1,
1888 nod_count + nplot - i + 1,
1889 nod_count + nplot - i);
1890 }
1891 }
1892 ++nod_count;
1893 }
1894 }
1895 }
1896
1897 /// Return total number of plot points (when plotting
1898 /// nplot points in each "coordinate direction)
1899 unsigned nplot_points(const unsigned& nplot) const
1900 {
1901 unsigned np = 0;
1902 for (unsigned i = 1; i <= nplot; i++)
1903 {
1904 np += i;
1905 }
1906 return np;
1907 }
1908
1909
1910 /// Build the lower-dimensional FaceElement (an element of type
1911 /// TElement<1,NNODE_1D>). The face index takes three possible values:
1912 /// 0 (Left) s[0] = 0.0
1913 /// 1 (Bottom) s[1] = 0.0
1914 /// 2 (Sloping face) s[0] = 1.0 - s[1]
1915 void build_face_element(const int& face_index,
1916 FaceElement* face_element_pt);
1917 };
1918
1919
1920 ///////////////////////////////////////////////////////////////////////
1921 ///////////////////////////////////////////////////////////////////////
1922 ///////////////////////////////////////////////////////////////////////
1923
1924
1925 //=======================================================================
1926 /// Return local coordinates of node j
1927 //=======================================================================
1928 template<>
1929 class TElementShape<3, 2>
1930 {
1931 public:
1932 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
1933 {
1934 s.resize(3);
1935
1936 switch (j)
1937 {
1938 case 0:
1939 s[0] = 1.0;
1940 s[1] = 0.0;
1941 s[2] = 0.0;
1942 break;
1943
1944 case 1:
1945 s[0] = 0.0;
1946 s[1] = 1.0;
1947 s[2] = 0.0;
1948 break;
1949
1950 case 2:
1951 s[0] = 0.0;
1952 s[1] = 0.0;
1953 s[2] = 1.0;
1954 break;
1955
1956 case 3:
1957 s[0] = 0.0;
1958 s[1] = 0.0;
1959 s[2] = 0.0;
1960 break;
1961
1962 default:
1963 std::ostringstream error_message;
1964 error_message
1965 << "Element only has four nodes; called with node number " << j
1966 << std::endl;
1967
1968 throw OomphLibError(error_message.str(),
1969 OOMPH_CURRENT_FUNCTION,
1970 OOMPH_EXCEPTION_LOCATION);
1971 }
1972 }
1973
1974
1975 //=======================================================================
1976 /// Shape function for specific TElement<3,2>
1977 //=======================================================================
1978 void shape(const Vector<double>& s, Shape& psi) const
1979 {
1980 psi[0] = s[0];
1981 psi[1] = s[1];
1982 psi[2] = s[2];
1983 psi[3] = 1.0 - s[0] - s[1] - s[2];
1984 }
1985
1986
1987 //=======================================================================
1988 /// Derivatives of shape functions for specific TElement<3,2>
1989 //=======================================================================
1990 void dshape_local(const Vector<double>& s, Shape& psi, DShape& dpsids) const
1991 {
1992 // ALH: Don't know why object qualifier is needed
1993 this->shape(s, psi);
1994
1995 // Derivatives
1996 dpsids(0, 0) = 1.0;
1997 dpsids(0, 1) = 0.0;
1998 dpsids(0, 2) = 0.0;
1999
2000 dpsids(1, 0) = 0.0;
2001 dpsids(1, 1) = 1.0;
2002 dpsids(1, 2) = 0.0;
2003
2004 dpsids(2, 0) = 0.0;
2005 dpsids(2, 1) = 0.0;
2006 dpsids(2, 2) = 1.0;
2007
2008 dpsids(3, 0) = -1.0;
2009 dpsids(3, 1) = -1.0;
2010 dpsids(3, 2) = -1.0;
2011 }
2012
2013
2014 //=======================================================================
2015 /// Second derivatives of shape functions for specific TElement<3,2>:
2016 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
2017 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
2018 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
2019 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
2020 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
2021 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
2022 //=======================================================================
2023 void d2shape_local(const Vector<double>& s,
2024 Shape& psi,
2025 DShape& dpsids,
2026 DShape& d2psids) const
2027 {
2028 // ALH: Don't know why object qualifier is needed
2029 this->dshape_local(s, psi, dpsids);
2030
2031 for (unsigned i = 0; i < 4; i++)
2032 {
2033 d2psids(i, 0) = 0.0;
2034 d2psids(i, 1) = 0.0;
2035 d2psids(i, 2) = 0.0;
2036 d2psids(i, 3) = 0.0;
2037 d2psids(i, 4) = 0.0;
2038 d2psids(i, 5) = 0.0;
2039 }
2040 }
2041 };
2042
2043
2044 //=======================================================================
2045 /// Return local coordinates of node j
2046 //=======================================================================
2047 template<>
2048 class TElementShape<3, 3>
2049 {
2050 public:
2051 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
2052 {
2053 s.resize(3);
2054
2055 switch (j)
2056 {
2057 case 0:
2058 s[0] = 1.0;
2059 s[1] = 0.0;
2060 s[2] = 0.0;
2061 break;
2062
2063 case 1:
2064 s[0] = 0.0;
2065 s[1] = 1.0;
2066 s[2] = 0.0;
2067 break;
2068
2069 case 2:
2070 s[0] = 0.0;
2071 s[1] = 0.0;
2072 s[2] = 1.0;
2073 break;
2074
2075 case 3:
2076 s[0] = 0.0;
2077 s[1] = 0.0;
2078 s[2] = 0.0;
2079 break;
2080
2081 case 4:
2082 s[0] = 0.5;
2083 s[1] = 0.5;
2084 s[2] = 0.0;
2085 break;
2086
2087 case 5:
2088 s[0] = 0.5;
2089 s[1] = 0.0;
2090 s[2] = 0.5;
2091 break;
2092
2093 case 6:
2094 s[0] = 0.5;
2095 s[1] = 0.0;
2096 s[2] = 0.0;
2097 break;
2098
2099 case 7:
2100 s[0] = 0.0;
2101 s[1] = 0.5;
2102 s[2] = 0.5;
2103 break;
2104
2105 case 8:
2106 s[0] = 0.0;
2107 s[1] = 0.0;
2108 s[2] = 0.5;
2109 break;
2110
2111 case 9:
2112 s[0] = 0.0;
2113 s[1] = 0.5;
2114 s[2] = 0.0;
2115 break;
2116
2117 default:
2118 std::ostringstream error_message;
2119 error_message
2120 << "Element only has ten nodes; called with node number " << j
2121 << std::endl;
2122
2123 throw OomphLibError(error_message.str(),
2124 OOMPH_CURRENT_FUNCTION,
2125 OOMPH_EXCEPTION_LOCATION);
2126 }
2127 }
2128
2129
2130 //=======================================================================
2131 /// Shape function for specific TElement<3,3>
2132 //=======================================================================
2133 void shape(const Vector<double>& s, Shape& psi) const
2134 {
2135 double s3 = 1.0 - s[0] - s[1] - s[2];
2136 psi[0] = (2.0 * s[0] - 1.0) * s[0];
2137 psi[1] = (2.0 * s[1] - 1.0) * s[1];
2138 psi[2] = (2.0 * s[2] - 1.0) * s[2];
2139 psi[3] = (2.0 * s3 - 1.0) * s3;
2140 psi[4] = 4.0 * s[0] * s[1];
2141 psi[5] = 4.0 * s[0] * s[2];
2142 psi[6] = 4.0 * s[0] * s3;
2143 psi[7] = 4.0 * s[1] * s[2];
2144 psi[8] = 4.0 * s[2] * s3;
2145 psi[9] = 4.0 * s[1] * s3;
2146 }
2147
2148
2149 //=======================================================================
2150 /// Derivatives of shape functions for specific TElement<3,3>
2151 //=======================================================================
2152 void dshape_local(const Vector<double>& s, Shape& psi, DShape& dpsids) const
2153 {
2154 // ALH: Don't know why object qualifier is needed
2155 this->shape(s, psi);
2156
2157 // Derivatives
2158 double s3 = 1.0 - s[0] - s[1] - s[2];
2159
2160 dpsids(0, 0) = 4.0 * s[0] - 1.0;
2161 dpsids(0, 1) = 0.0;
2162 dpsids(0, 2) = 0.0;
2163
2164 dpsids(1, 0) = 0.0;
2165 dpsids(1, 1) = 4.0 * s[1] - 1.0;
2166 dpsids(1, 2) = 0.0;
2167
2168 dpsids(2, 0) = 0.0;
2169 dpsids(2, 1) = 0.0;
2170 dpsids(2, 2) = 4.0 * s[2] - 1.0;
2171
2172 dpsids(3, 0) = -4.0 * s3 + 1.0;
2173 dpsids(3, 1) = -4.0 * s3 + 1.0;
2174 dpsids(3, 2) = -4.0 * s3 + 1.0;
2175
2176 dpsids(4, 0) = 4.0 * s[1];
2177 dpsids(4, 1) = 4.0 * s[0];
2178 dpsids(4, 2) = 0.0;
2179
2180 dpsids(5, 0) = 4.0 * s[2];
2181 dpsids(5, 1) = 0.0;
2182 dpsids(5, 2) = 4.0 * s[0];
2183
2184 dpsids(6, 0) = 4.0 * (s3 - s[0]);
2185 dpsids(6, 1) = -4.0 * s[0];
2186 dpsids(6, 2) = -4.0 * s[0];
2187
2188 dpsids(7, 0) = 0.0;
2189 dpsids(7, 1) = 4.0 * s[2];
2190 dpsids(7, 2) = 4.0 * s[1];
2191
2192 dpsids(8, 0) = -4.0 * s[2];
2193 dpsids(8, 1) = -4.0 * s[2];
2194 dpsids(8, 2) = 4.0 * (s3 - s[2]);
2195
2196 dpsids(9, 0) = -4.0 * s[1];
2197 dpsids(9, 1) = 4.0 * (s3 - s[1]);
2198 dpsids(9, 2) = -4.0 * s[1];
2199 }
2200
2201
2202 //=======================================================================
2203 /// Second derivatives of shape functions for specific TElement<3,3>:
2204 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
2205 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
2206 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
2207 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
2208 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
2209 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
2210 //=======================================================================
2211 void d2shape_local(const Vector<double>& s,
2212 Shape& psi,
2213 DShape& dpsids,
2214 DShape& d2psids) const
2215 {
2216 // ALH: Don't know why object qualifier is needed
2217 this->dshape_local(s, psi, dpsids);
2218
2219 //(.,3) for mixed derivative s[0]-s[1]
2220 //(.,4) for mixed derivative s[0]-s[2]
2221 //(.,5) for mixed derivative s[1]-s[2]
2222
2223 d2psids(0, 0) = 4.0;
2224 d2psids(0, 1) = 0.0;
2225 d2psids(0, 2) = 0.0;
2226 d2psids(0, 3) = 0.0;
2227 d2psids(0, 4) = 0.0;
2228 d2psids(0, 5) = 0.0;
2229
2230
2231 d2psids(1, 0) = 0.0;
2232 d2psids(1, 1) = 4.0;
2233 d2psids(1, 2) = 0.0;
2234 d2psids(1, 3) = 0.0;
2235 d2psids(1, 4) = 0.0;
2236 d2psids(1, 5) = 0.0;
2237
2238 d2psids(2, 0) = 0.0;
2239 d2psids(2, 1) = 0.0;
2240 d2psids(2, 2) = 4.0;
2241 d2psids(2, 3) = 0.0;
2242 d2psids(2, 4) = 0.0;
2243 d2psids(2, 5) = 0.0;
2244
2245 d2psids(3, 0) = 4.0;
2246 d2psids(3, 1) = 4.0;
2247 d2psids(3, 2) = 4.0;
2248 d2psids(3, 3) = 4.0;
2249 d2psids(3, 4) = 4.0;
2250 d2psids(3, 5) = 4.0;
2251
2252 d2psids(4, 0) = 0.0;
2253 d2psids(4, 1) = 0.0;
2254 d2psids(4, 2) = 0.0;
2255 d2psids(4, 3) = 4.0;
2256 d2psids(4, 4) = 0.0;
2257 d2psids(4, 5) = 0.0;
2258
2259 d2psids(5, 0) = 0.0;
2260 d2psids(5, 1) = 0.0;
2261 d2psids(5, 2) = 0.0;
2262 d2psids(5, 3) = 0.0;
2263 d2psids(5, 4) = 4.0;
2264 d2psids(5, 5) = 0.0;
2265
2266 d2psids(6, 0) = -8.0;
2267 d2psids(6, 1) = 0.0;
2268 d2psids(6, 2) = 0.0;
2269 d2psids(6, 3) = -4.0;
2270 d2psids(6, 4) = -4.0;
2271 d2psids(6, 5) = 0.0;
2272
2273 d2psids(7, 0) = 0.0;
2274 d2psids(7, 1) = 0.0;
2275 d2psids(7, 2) = 0.0;
2276 d2psids(7, 3) = 0.0;
2277 d2psids(7, 4) = 0.0;
2278 d2psids(7, 5) = 4.0;
2279
2280 d2psids(8, 0) = 0.0;
2281 d2psids(8, 1) = 0.0;
2282 d2psids(8, 2) = -8.0;
2283 d2psids(8, 3) = 0.0;
2284 d2psids(8, 4) = -4.0;
2285 d2psids(8, 5) = -4.0;
2286
2287 d2psids(9, 0) = 0.0;
2288 d2psids(9, 1) = -8.0;
2289 d2psids(9, 2) = 0.0;
2290 d2psids(9, 3) = -4.0;
2291 d2psids(9, 4) = 0.0;
2292 d2psids(9, 5) = -4.0;
2293 }
2294 };
2295
2296 //////////////////////////////////////////////////////////////////////
2297 //////////////////////////////////////////////////////////////////////
2298 //////////////////////////////////////////////////////////////////////
2299
2300 //====================================================================
2301 /// Standard quadratic shape functions enriched by the addition of
2302 /// three cubic "face" bubbles and quartic "volume" bubble,
2303 /// which consists of adding a node at the centroid of
2304 /// each face and a single node at the centroid
2305 /// of the tetrahedron
2306 //=========================================================================
2307
2308 //=======================================================================
2309 /// Return local coordinates of node j
2310 //=======================================================================
2311 template<>
2312 class TBubbleEnrichedElementShape<3, 3>
2313 {
2314 public:
2315 unsigned n_enriched_nodes()
2316 {
2317 return 5;
2318 }
2319
2320 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
2321 {
2322 s.resize(3);
2323
2324 switch (j)
2325 {
2326 case 0:
2327 s[0] = 1.0;
2328 s[1] = 0.0;
2329 s[2] = 0.0;
2330 break;
2331
2332 case 1:
2333 s[0] = 0.0;
2334 s[1] = 1.0;
2335 s[2] = 0.0;
2336 break;
2337
2338 case 2:
2339 s[0] = 0.0;
2340 s[1] = 0.0;
2341 s[2] = 1.0;
2342 break;
2343
2344 case 3:
2345 s[0] = 0.0;
2346 s[1] = 0.0;
2347 s[2] = 0.0;
2348 break;
2349
2350 case 4:
2351 s[0] = 0.5;
2352 s[1] = 0.5;
2353 s[2] = 0.0;
2354 break;
2355
2356 case 5:
2357 s[0] = 0.5;
2358 s[1] = 0.0;
2359 s[2] = 0.5;
2360 break;
2361
2362 case 6:
2363 s[0] = 0.5;
2364 s[1] = 0.0;
2365 s[2] = 0.0;
2366 break;
2367
2368 case 7:
2369 s[0] = 0.0;
2370 s[1] = 0.5;
2371 s[2] = 0.5;
2372 break;
2373
2374 case 8:
2375 s[0] = 0.0;
2376 s[1] = 0.0;
2377 s[2] = 0.5;
2378 break;
2379
2380 case 9:
2381 s[0] = 0.0;
2382 s[1] = 0.5;
2383 s[2] = 0.0;
2384 break;
2385
2386 // Node at centroid of face spanned by nodes 0, 1, 3
2387 case 10:
2388 s[0] = 1.0 / 3.0;
2389 s[1] = 1.0 / 3.0;
2390 s[2] = 0.0;
2391 break;
2392
2393 // Node at centroid of face spanned by nodes 0, 1, 2
2394 case 11:
2395 s[0] = 1.0 / 3.0;
2396 s[1] = 1.0 / 3.0;
2397 s[2] = 1.0 / 3.0;
2398 break;
2399
2400 // Node at centroid of face spanned by nodes 0, 2, 3
2401 case 12:
2402 s[0] = 1.0 / 3.0;
2403 s[1] = 0.0;
2404 s[2] = 1.0 / 3.0;
2405 break;
2406
2407 // Node at centroid of face spannd by nodes 1, 2, 3
2408 case 13:
2409 s[0] = 0.0;
2410 s[1] = 1.0 / 3.0;
2411 s[2] = 1.0 / 3.0;
2412 break;
2413
2414 // Node at centroid of volume
2415 case 14:
2416 s[0] = 0.25;
2417 s[1] = 0.25;
2418 s[2] = 0.25;
2419 break;
2420
2421
2422 default:
2423 std::ostringstream error_message;
2424 error_message
2425 << "Element only has fifteen nodes; called with node number " << j
2426 << std::endl;
2427
2428 throw OomphLibError(error_message.str(),
2429 OOMPH_CURRENT_FUNCTION,
2430 OOMPH_EXCEPTION_LOCATION);
2431 }
2432 }
2433
2434
2435 //=======================================================================
2436 /// Shape function for specific TBubbleEnrichedElement<3,3>
2437 //=======================================================================
2438 void shape(const Vector<double>& s, Shape& psi) const
2439 {
2440 // Constructe the fourth volume coordinate
2441 const double s3 = 1.0 - s[0] - s[1] - s[2];
2442 // calculate the enrichment functions
2443 const double quartic_bubble = s[0] * s[1] * s[2] * s3;
2444 const double cubic_bubble012 = s[0] * s[1] * s[2];
2445 const double cubic_bubble013 = s[0] * s[1] * s3;
2446 const double cubic_bubble023 = s[0] * s[2] * s3;
2447 const double cubic_bubble123 = s[1] * s[2] * s3;
2448
2449 // The appropriate "amount" of cubic and quartic bubble functions are
2450 // added/subtracted
2451 // to each original quadratic shape function to ensure that the new
2452 // shape function is zero at the centroid (0.25,0.25,0.25)
2453 // and at the face centroids
2454 psi[0] = (2.0 * s[0] - 1.0) * s[0] +
2455 3.0 * (cubic_bubble012 + cubic_bubble013 + cubic_bubble023) -
2456 4.0 * quartic_bubble;
2457 psi[1] = (2.0 * s[1] - 1.0) * s[1] +
2458 3.0 * (cubic_bubble012 + cubic_bubble013 + cubic_bubble123) -
2459 4.0 * quartic_bubble;
2460 psi[2] = (2.0 * s[2] - 1.0) * s[2] +
2461 3.0 * (cubic_bubble012 + cubic_bubble023 + cubic_bubble123) -
2462 4.0 * quartic_bubble;
2463 psi[3] = (2.0 * s3 - 1.0) * s3 +
2464 3.0 * (cubic_bubble013 + cubic_bubble023 + cubic_bubble123) -
2465 4.0 * quartic_bubble;
2466 psi[4] = 4.0 * s[0] * s[1] - 12.0 * (cubic_bubble012 + cubic_bubble013) +
2467 32.0 * quartic_bubble;
2468 psi[5] = 4.0 * s[0] * s[2] - 12.0 * (cubic_bubble012 + cubic_bubble023) +
2469 32.0 * quartic_bubble;
2470 psi[6] = 4.0 * s[0] * s3 - 12.0 * (cubic_bubble013 + cubic_bubble023) +
2471 32.0 * quartic_bubble;
2472 psi[7] = 4.0 * s[1] * s[2] - 12.0 * (cubic_bubble012 + cubic_bubble123) +
2473 32.0 * quartic_bubble;
2474 psi[8] = 4.0 * s[2] * s3 - 12.0 * (cubic_bubble023 + cubic_bubble123) +
2475 32.0 * quartic_bubble;
2476 psi[9] = 4.0 * s[1] * s3 - 12.0 * (cubic_bubble013 + cubic_bubble123) +
2477 32.0 * quartic_bubble;
2478 // Add the bubble function on the face spanned by 0 1 3
2479 psi[10] = 27.0 * cubic_bubble013 - 108.0 * quartic_bubble;
2480 // Add the bubble function on the face spanned by 0 1 2
2481 psi[11] = 27.0 * cubic_bubble012 - 108.0 * quartic_bubble;
2482 // Add the bubble function on the face spanned by 0 2 3
2483 psi[12] = 27.0 * cubic_bubble023 - 108.0 * quartic_bubble;
2484 // Add the bubble function on the face spanned by 1 2 3
2485 psi[13] = 27.0 * cubic_bubble123 - 108.0 * quartic_bubble;
2486 // Add the volume bubble function, scaled to have value one
2487 psi[14] = 256.0 * quartic_bubble;
2488 }
2489
2490
2491 //=======================================================================
2492 /// Derivatives of shape functions for specific TElement<3,3>
2493 //=======================================================================
2494 void dshape_local(const Vector<double>& s, Shape& psi, DShape& dpsids) const
2495 {
2496 // ALH: Don't know why object qualifier is needed
2497 this->shape(s, psi);
2498
2499 // Define s3 the fourth volume coordinate
2500 const double s3 = 1.0 - s[0] - s[1] - s[2];
2501
2502 // Calculate derivatives of the bubble function
2503 const double d_quartic_bubble_ds0 =
2504 s[1] * s[2] * (1.0 - s[1] - s[2] - 2.0 * s[0]);
2505 const double d_quartic_bubble_ds1 =
2506 s[0] * s[2] * (1.0 - s[0] - s[2] - 2.0 * s[1]);
2507 const double d_quartic_bubble_ds2 =
2508 s[0] * s[1] * (1.0 - s[0] - s[1] - 2.0 * s[2]);
2509
2510 const double d_cubic_bubble012_ds0 = s[1] * s[2];
2511 const double d_cubic_bubble012_ds1 = s[0] * s[2];
2512 const double d_cubic_bubble012_ds2 = s[0] * s[1];
2513
2514 const double d_cubic_bubble013_ds0 = s[1] * (s3 - s[0]);
2515 const double d_cubic_bubble013_ds1 = s[0] * (s3 - s[1]);
2516 const double d_cubic_bubble013_ds2 = -s[0] * s[1];
2517
2518 const double d_cubic_bubble023_ds0 = s[2] * (s3 - s[0]);
2519 const double d_cubic_bubble023_ds1 = -s[0] * s[2];
2520 const double d_cubic_bubble023_ds2 = s[0] * (s3 - s[2]);
2521
2522 const double d_cubic_bubble123_ds0 = -s[1] * s[2];
2523 const double d_cubic_bubble123_ds1 = s[2] * (s3 - s[1]);
2524 const double d_cubic_bubble123_ds2 = s[1] * (s3 - s[2]);
2525
2526
2527 // Add the appropriate dervatives of the bubble function to the
2528 // shape function derivatives
2529 dpsids(0, 0) = 4.0 * s[0] - 1.0 +
2530 3.0 * (d_cubic_bubble012_ds0 + d_cubic_bubble013_ds0 +
2531 d_cubic_bubble023_ds0) -
2532 4.0 * d_quartic_bubble_ds0;
2533 dpsids(0, 1) = 0.0 +
2534 3.0 * (d_cubic_bubble012_ds1 + d_cubic_bubble013_ds1 +
2535 d_cubic_bubble023_ds1) -
2536 4.0 * d_quartic_bubble_ds1;
2537 dpsids(0, 2) = 0.0 +
2538 3.0 * (d_cubic_bubble012_ds2 + d_cubic_bubble013_ds2 +
2539 d_cubic_bubble023_ds2) -
2540 4.0 * d_quartic_bubble_ds2;
2541
2542 dpsids(1, 0) = 0.0 +
2543 3.0 * (d_cubic_bubble012_ds0 + d_cubic_bubble013_ds0 +
2544 d_cubic_bubble123_ds0) -
2545 4.0 * d_quartic_bubble_ds0;
2546 dpsids(1, 1) = 4.0 * s[1] - 1.0 +
2547 3.0 * (d_cubic_bubble012_ds1 + d_cubic_bubble013_ds1 +
2548 d_cubic_bubble123_ds1) -
2549 4.0 * d_quartic_bubble_ds1;
2550 dpsids(1, 2) = 0.0 +
2551 3.0 * (d_cubic_bubble012_ds2 + d_cubic_bubble013_ds2 +
2552 d_cubic_bubble123_ds2) -
2553 4.0 * d_quartic_bubble_ds2;
2554
2555 dpsids(2, 0) = 0.0 +
2556 3.0 * (d_cubic_bubble012_ds0 + d_cubic_bubble023_ds0 +
2557 d_cubic_bubble123_ds0) -
2558 4.0 * d_quartic_bubble_ds0;
2559 dpsids(2, 1) = 0.0 +
2560 3.0 * (d_cubic_bubble012_ds1 + d_cubic_bubble023_ds1 +
2561 d_cubic_bubble123_ds1) -
2562 4.0 * d_quartic_bubble_ds1;
2563 dpsids(2, 2) = 4.0 * s[2] - 1.0 +
2564 3.0 * (d_cubic_bubble012_ds2 + d_cubic_bubble023_ds2 +
2565 d_cubic_bubble123_ds2) -
2566 4.0 * d_quartic_bubble_ds2;
2567
2568 dpsids(3, 0) = -4.0 * s3 + 1.0 +
2569 3.0 * (d_cubic_bubble013_ds0 + d_cubic_bubble023_ds0 +
2570 d_cubic_bubble123_ds0) -
2571 4.0 * d_quartic_bubble_ds0;
2572 dpsids(3, 1) = -4.0 * s3 + 1.0 +
2573 3.0 * (d_cubic_bubble013_ds1 + d_cubic_bubble023_ds1 +
2574 d_cubic_bubble123_ds1) -
2575 4.0 * d_quartic_bubble_ds1;
2576 dpsids(3, 2) = -4.0 * s3 + 1.0 +
2577 3.0 * (d_cubic_bubble013_ds2 + d_cubic_bubble023_ds2 +
2578 d_cubic_bubble123_ds2) -
2579 4.0 * d_quartic_bubble_ds2;
2580
2581 dpsids(4, 0) = 4.0 * s[1] -
2582 12.0 * (d_cubic_bubble012_ds0 + d_cubic_bubble013_ds0) +
2583 32.0 * d_quartic_bubble_ds0;
2584 dpsids(4, 1) = 4.0 * s[0] -
2585 12.0 * (d_cubic_bubble012_ds1 + d_cubic_bubble013_ds1) +
2586 32.0 * d_quartic_bubble_ds1;
2587 dpsids(4, 2) = 0.0 -
2588 12.0 * (d_cubic_bubble012_ds2 + d_cubic_bubble013_ds2) +
2589 32.0 * d_quartic_bubble_ds2;
2590
2591 dpsids(5, 0) = 4.0 * s[2] -
2592 12.0 * (d_cubic_bubble012_ds0 + d_cubic_bubble023_ds0) +
2593 32.0 * d_quartic_bubble_ds0;
2594 dpsids(5, 1) = 0.0 -
2595 12.0 * (d_cubic_bubble012_ds1 + d_cubic_bubble023_ds1) +
2596 32.0 * d_quartic_bubble_ds1;
2597 dpsids(5, 2) = 4.0 * s[0] -
2598 12.0 * (d_cubic_bubble012_ds2 + d_cubic_bubble023_ds2) +
2599 32.0 * d_quartic_bubble_ds2;
2600
2601 dpsids(6, 0) = 4.0 * (s3 - s[0]) -
2602 12.0 * (d_cubic_bubble013_ds0 + d_cubic_bubble023_ds0) +
2603 32.0 * d_quartic_bubble_ds0;
2604 dpsids(6, 1) = -4.0 * s[0] -
2605 12.0 * (d_cubic_bubble013_ds1 + d_cubic_bubble023_ds1) +
2606 32.0 * d_quartic_bubble_ds1;
2607 dpsids(6, 2) = -4.0 * s[0] -
2608 12.0 * (d_cubic_bubble013_ds2 + d_cubic_bubble023_ds2) +
2609 32.0 * d_quartic_bubble_ds2;
2610
2611 dpsids(7, 0) = 0.0 -
2612 12.0 * (d_cubic_bubble012_ds0 + d_cubic_bubble123_ds0) +
2613 32.0 * d_quartic_bubble_ds0;
2614 dpsids(7, 1) = 4.0 * s[2] -
2615 12.0 * (d_cubic_bubble012_ds1 + d_cubic_bubble123_ds1) +
2616 32.0 * d_quartic_bubble_ds1;
2617 dpsids(7, 2) = 4.0 * s[1] -
2618 12.0 * (d_cubic_bubble012_ds2 + d_cubic_bubble123_ds2) +
2619 32.0 * d_quartic_bubble_ds2;
2620
2621 dpsids(8, 0) = -4.0 * s[2] -
2622 12.0 * (d_cubic_bubble023_ds0 + d_cubic_bubble123_ds0) +
2623 32.0 * d_quartic_bubble_ds0;
2624 dpsids(8, 1) = -4.0 * s[2] -
2625 12.0 * (d_cubic_bubble023_ds1 + d_cubic_bubble123_ds1) +
2626 32.0 * d_quartic_bubble_ds1;
2627 dpsids(8, 2) = 4.0 * (s3 - s[2]) -
2628 12.0 * (d_cubic_bubble023_ds2 + d_cubic_bubble123_ds2) +
2629 32.0 * d_quartic_bubble_ds2;
2630
2631 dpsids(9, 0) = -4.0 * s[1] -
2632 12.0 * (d_cubic_bubble013_ds0 + d_cubic_bubble123_ds0) +
2633 32.0 * d_quartic_bubble_ds0;
2634 dpsids(9, 1) = 4.0 * (s3 - s[1]) -
2635 12.0 * (d_cubic_bubble013_ds1 + d_cubic_bubble123_ds1) +
2636 32.0 * d_quartic_bubble_ds1;
2637 dpsids(9, 2) = -4.0 * s[1] -
2638 12.0 * (d_cubic_bubble013_ds2 + d_cubic_bubble123_ds2) +
2639 32.0 * d_quartic_bubble_ds2;
2640
2641 // Add the bubble function on the face spanned by 0 1 3
2642 dpsids(10, 0) =
2643 27.0 * d_cubic_bubble013_ds0 - 108.0 * d_quartic_bubble_ds0;
2644 dpsids(10, 1) =
2645 27.0 * d_cubic_bubble013_ds1 - 108.0 * d_quartic_bubble_ds1;
2646 dpsids(10, 2) =
2647 27.0 * d_cubic_bubble013_ds2 - 108.0 * d_quartic_bubble_ds2;
2648
2649 // Add the bubble function on the face spanned by 0 1 2
2650 dpsids(11, 0) =
2651 27.0 * d_cubic_bubble012_ds0 - 108.0 * d_quartic_bubble_ds0;
2652 dpsids(11, 1) =
2653 27.0 * d_cubic_bubble012_ds1 - 108.0 * d_quartic_bubble_ds1;
2654 dpsids(11, 2) =
2655 27.0 * d_cubic_bubble012_ds2 - 108.0 * d_quartic_bubble_ds2;
2656
2657 // Add the bubble function on the face spanned by 0 2 3
2658 dpsids(12, 0) =
2659 27.0 * d_cubic_bubble023_ds0 - 108.0 * d_quartic_bubble_ds0;
2660 dpsids(12, 1) =
2661 27.0 * d_cubic_bubble023_ds1 - 108.0 * d_quartic_bubble_ds1;
2662 dpsids(12, 2) =
2663 27.0 * d_cubic_bubble023_ds2 - 108.0 * d_quartic_bubble_ds2;
2664
2665 // Add the bubble function on the face spanned by 1 2 3
2666 dpsids(13, 0) =
2667 27.0 * d_cubic_bubble123_ds0 - 108.0 * d_quartic_bubble_ds0;
2668 dpsids(13, 1) =
2669 27.0 * d_cubic_bubble123_ds1 - 108.0 * d_quartic_bubble_ds1;
2670 dpsids(13, 2) =
2671 27.0 * d_cubic_bubble123_ds2 - 108.0 * d_quartic_bubble_ds2;
2672
2673 // Add the volumetric bubble function derivatives
2674 dpsids(14, 0) = 256.0 * d_quartic_bubble_ds0;
2675 dpsids(14, 1) = 256.0 * d_quartic_bubble_ds1;
2676 dpsids(14, 2) = 256.0 * d_quartic_bubble_ds2;
2677 }
2678
2679
2680 //=======================================================================
2681 /// Second derivatives of shape functions for specific TElement<3,3>:
2682 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
2683 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
2684 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
2685 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
2686 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
2687 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
2688 //=======================================================================
2689 void d2shape_local(const Vector<double>& s,
2690 Shape& psi,
2691 DShape& dpsids,
2692 DShape& d2psids) const
2693 {
2694 // ALH: Don't know why object qualifier is needed
2695 this->dshape_local(s, psi, dpsids);
2696
2697 // Define s3
2698 const double s3 = 1.0 - s[0] - s[1] - s[2];
2699
2700 // Calculate second derivatives of the bubble functions
2701 //(.,3) for mixed derivative s[0]-s[1]
2702 //(.,4) for mixed derivative s[0]-s[2]
2703 //(.,5) for mixed derivative s[1]-s[2]
2704
2705
2706 const double d2_quartic_bubble_ds0 = -2.0 * s[1] * s[2];
2707 const double d2_quartic_bubble_ds1 = -2.0 * s[0] * s[2];
2708 const double d2_quartic_bubble_ds2 = -2.0 * s[0] * s[1];
2709 const double d2_quartic_bubble_ds3 =
2710 s[2] * (1.0 - 2.0 * s[0] - 2.0 * s[1] - s[2]);
2711 const double d2_quartic_bubble_ds4 =
2712 s[1] * (1.0 - 2.0 * s[0] - 2.0 * s[2] - s[1]);
2713 const double d2_quartic_bubble_ds5 =
2714 s[0] * (1.0 - 2.0 * s[1] - 2.0 * s[2] - s[0]);
2715
2716 const double d2_cubic_bubble012_ds0 = 0.0;
2717 const double d2_cubic_bubble012_ds1 = 0.0;
2718 const double d2_cubic_bubble012_ds2 = 0.0;
2719 const double d2_cubic_bubble012_ds3 = s[2];
2720 const double d2_cubic_bubble012_ds4 = s[1];
2721 const double d2_cubic_bubble012_ds5 = s[0];
2722
2723 const double d2_cubic_bubble013_ds0 = -2.0 * s[1];
2724 const double d2_cubic_bubble013_ds1 = -2.0 * s[0];
2725 const double d2_cubic_bubble013_ds2 = 0.0;
2726 const double d2_cubic_bubble013_ds3 = s3 - s[0] - s[1];
2727 const double d2_cubic_bubble013_ds4 = -s[1];
2728 const double d2_cubic_bubble013_ds5 = -s[0];
2729
2730 const double d2_cubic_bubble023_ds0 = -2.0 * s[2];
2731 const double d2_cubic_bubble023_ds1 = 0.0;
2732 const double d2_cubic_bubble023_ds2 = -2.0 * s[0];
2733 const double d2_cubic_bubble023_ds3 = -s[2];
2734 const double d2_cubic_bubble023_ds4 = s3 - s[0] - s[2];
2735 const double d2_cubic_bubble023_ds5 = -s[0];
2736
2737 const double d2_cubic_bubble123_ds0 = 0.0;
2738 const double d2_cubic_bubble123_ds1 = -2.0 * s[2];
2739 const double d2_cubic_bubble123_ds2 = -2.0 * s[1];
2740 const double d2_cubic_bubble123_ds3 = -s[2];
2741 const double d2_cubic_bubble123_ds4 = -s[1];
2742 const double d2_cubic_bubble123_ds5 = s3 - s[1] - s[2];
2743
2744
2745 d2psids(0, 0) = 4.0 +
2746 3.0 * (d2_cubic_bubble012_ds0 + d2_cubic_bubble013_ds0 +
2747 d2_cubic_bubble023_ds0) -
2748 4.0 * d2_quartic_bubble_ds0;
2749 d2psids(0, 1) = 0.0 +
2750 3.0 * (d2_cubic_bubble012_ds1 + d2_cubic_bubble013_ds1 +
2751 d2_cubic_bubble023_ds1) -
2752 4.0 * d2_quartic_bubble_ds1;
2753 d2psids(0, 2) = 0.0 +
2754 3.0 * (d2_cubic_bubble012_ds2 + d2_cubic_bubble013_ds2 +
2755 d2_cubic_bubble023_ds2) -
2756 4.0 * d2_quartic_bubble_ds2;
2757 d2psids(0, 3) = 0.0 +
2758 3.0 * (d2_cubic_bubble012_ds3 + d2_cubic_bubble013_ds3 +
2759 d2_cubic_bubble023_ds3) -
2760 4.0 * d2_quartic_bubble_ds3;
2761 d2psids(0, 4) = 0.0 +
2762 3.0 * (d2_cubic_bubble012_ds4 + d2_cubic_bubble013_ds4 +
2763 d2_cubic_bubble023_ds4) -
2764 4.0 * d2_quartic_bubble_ds4;
2765 d2psids(0, 5) = 0.0 +
2766 3.0 * (d2_cubic_bubble012_ds5 + d2_cubic_bubble013_ds5 +
2767 d2_cubic_bubble023_ds5) -
2768 4.0 * d2_quartic_bubble_ds5;
2769
2770
2771 d2psids(1, 0) = 0.0 +
2772 3.0 * (d2_cubic_bubble012_ds0 + d2_cubic_bubble013_ds0 +
2773 d2_cubic_bubble123_ds0) -
2774 4.0 * d2_quartic_bubble_ds0;
2775 d2psids(1, 1) = 4.0 +
2776 3.0 * (d2_cubic_bubble012_ds1 + d2_cubic_bubble013_ds1 +
2777 d2_cubic_bubble123_ds1) -
2778 4.0 * d2_quartic_bubble_ds1;
2779 d2psids(1, 2) = 0.0 +
2780 3.0 * (d2_cubic_bubble012_ds2 + d2_cubic_bubble013_ds2 +
2781 d2_cubic_bubble123_ds2) -
2782 4.0 * d2_quartic_bubble_ds2;
2783 d2psids(1, 3) = 0.0 +
2784 3.0 * (d2_cubic_bubble012_ds3 + d2_cubic_bubble013_ds3 +
2785 d2_cubic_bubble123_ds3) -
2786 4.0 * d2_quartic_bubble_ds3;
2787 d2psids(1, 4) = 0.0 +
2788 3.0 * (d2_cubic_bubble012_ds4 + d2_cubic_bubble013_ds4 +
2789 d2_cubic_bubble123_ds4) -
2790 4.0 * d2_quartic_bubble_ds4;
2791 d2psids(1, 5) = 0.0 +
2792 3.0 * (d2_cubic_bubble012_ds5 + d2_cubic_bubble013_ds5 +
2793 d2_cubic_bubble123_ds5) -
2794 4.0 * d2_quartic_bubble_ds5;
2795
2796
2797 d2psids(2, 0) = 0.0 +
2798 3.0 * (d2_cubic_bubble012_ds0 + d2_cubic_bubble023_ds0 +
2799 d2_cubic_bubble123_ds0) -
2800 4.0 * d2_quartic_bubble_ds0;
2801 d2psids(2, 1) = 0.0 +
2802 3.0 * (d2_cubic_bubble012_ds1 + d2_cubic_bubble023_ds1 +
2803 d2_cubic_bubble123_ds1) -
2804 4.0 * d2_quartic_bubble_ds1;
2805 d2psids(2, 2) = 4.0 +
2806 3.0 * (d2_cubic_bubble012_ds2 + d2_cubic_bubble023_ds2 +
2807 d2_cubic_bubble123_ds2) -
2808 4.0 * d2_quartic_bubble_ds2;
2809 d2psids(2, 3) = 0.0 +
2810 3.0 * (d2_cubic_bubble012_ds3 + d2_cubic_bubble023_ds3 +
2811 d2_cubic_bubble123_ds3) -
2812 4.0 * d2_quartic_bubble_ds3;
2813 d2psids(2, 4) = 0.0 +
2814 3.0 * (d2_cubic_bubble012_ds4 + d2_cubic_bubble023_ds4 +
2815 d2_cubic_bubble123_ds4) -
2816 4.0 * d2_quartic_bubble_ds4;
2817 d2psids(2, 5) = 0.0 +
2818 3.0 * (d2_cubic_bubble012_ds5 + d2_cubic_bubble023_ds5 +
2819 d2_cubic_bubble123_ds5) -
2820 4.0 * d2_quartic_bubble_ds5;
2821
2822
2823 d2psids(3, 0) = 4.0 +
2824 3.0 * (d2_cubic_bubble013_ds0 + d2_cubic_bubble023_ds0 +
2825 d2_cubic_bubble123_ds0) -
2826 4.0 * d2_quartic_bubble_ds0;
2827 d2psids(3, 1) = 4.0 +
2828 3.0 * (d2_cubic_bubble013_ds1 + d2_cubic_bubble023_ds1 +
2829 d2_cubic_bubble123_ds1) -
2830 4.0 * d2_quartic_bubble_ds1;
2831 d2psids(3, 2) = 4.0 +
2832 3.0 * (d2_cubic_bubble013_ds2 + d2_cubic_bubble023_ds2 +
2833 d2_cubic_bubble123_ds2) -
2834 4.0 * d2_quartic_bubble_ds2;
2835 d2psids(3, 3) = 4.0 +
2836 3.0 * (d2_cubic_bubble013_ds3 + d2_cubic_bubble023_ds3 +
2837 d2_cubic_bubble123_ds3) -
2838 4.0 * d2_quartic_bubble_ds3;
2839 d2psids(3, 4) = 4.0 +
2840 3.0 * (d2_cubic_bubble013_ds4 + d2_cubic_bubble023_ds4 +
2841 d2_cubic_bubble123_ds4) -
2842 4.0 * d2_quartic_bubble_ds4;
2843 d2psids(3, 5) = 4.0 +
2844 3.0 * (d2_cubic_bubble013_ds5 + d2_cubic_bubble023_ds5 +
2845 d2_cubic_bubble123_ds5) -
2846 4.0 * d2_quartic_bubble_ds5;
2847
2848
2849 d2psids(4, 0) = 0.0 -
2850 12.0 * (d2_cubic_bubble012_ds0 + d2_cubic_bubble013_ds0) +
2851 32.0 * d2_quartic_bubble_ds0;
2852 d2psids(4, 1) = 0.0 -
2853 12.0 * (d2_cubic_bubble012_ds1 + d2_cubic_bubble013_ds1) +
2854 32.0 * d2_quartic_bubble_ds1;
2855 d2psids(4, 2) = 0.0 -
2856 12.0 * (d2_cubic_bubble012_ds2 + d2_cubic_bubble013_ds2) +
2857 32.0 * d2_quartic_bubble_ds2;
2858 d2psids(4, 3) = 4.0 -
2859 12.0 * (d2_cubic_bubble012_ds3 + d2_cubic_bubble013_ds3) +
2860 32.0 * d2_quartic_bubble_ds3;
2861 d2psids(4, 4) = 0.0 -
2862 12.0 * (d2_cubic_bubble012_ds4 + d2_cubic_bubble013_ds4) +
2863 32.0 * d2_quartic_bubble_ds4;
2864 d2psids(4, 5) = 0.0 -
2865 12.0 * (d2_cubic_bubble012_ds5 + d2_cubic_bubble013_ds5) +
2866 32.0 * d2_quartic_bubble_ds5;
2867
2868
2869 d2psids(5, 0) = 0.0 -
2870 12.0 * (d2_cubic_bubble012_ds0 + d2_cubic_bubble023_ds0) +
2871 32.0 * d2_quartic_bubble_ds0;
2872 d2psids(5, 1) = 0.0 -
2873 12.0 * (d2_cubic_bubble012_ds1 + d2_cubic_bubble023_ds1) +
2874 32.0 * d2_quartic_bubble_ds1;
2875 d2psids(5, 2) = 0.0 -
2876 12.0 * (d2_cubic_bubble012_ds2 + d2_cubic_bubble023_ds2) +
2877 32.0 * d2_quartic_bubble_ds2;
2878 d2psids(5, 3) = 0.0 -
2879 12.0 * (d2_cubic_bubble012_ds3 + d2_cubic_bubble023_ds3) +
2880 32.0 * d2_quartic_bubble_ds3;
2881 d2psids(5, 4) = 4.0 -
2882 12.0 * (d2_cubic_bubble012_ds4 + d2_cubic_bubble023_ds4) +
2883 32.0 * d2_quartic_bubble_ds4;
2884 d2psids(5, 5) = 0.0 -
2885 12.0 * (d2_cubic_bubble012_ds5 + d2_cubic_bubble023_ds5) +
2886 32.0 * d2_quartic_bubble_ds5;
2887
2888
2889 d2psids(6, 0) = -8.0 -
2890 12.0 * (d2_cubic_bubble013_ds0 + d2_cubic_bubble023_ds0) +
2891 32.0 * d2_quartic_bubble_ds0;
2892 d2psids(6, 1) = 0.0 -
2893 12.0 * (d2_cubic_bubble013_ds1 + d2_cubic_bubble023_ds1) +
2894 32.0 * d2_quartic_bubble_ds1;
2895 d2psids(6, 2) = 0.0 -
2896 12.0 * (d2_cubic_bubble013_ds2 + d2_cubic_bubble023_ds2) +
2897 32.0 * d2_quartic_bubble_ds2;
2898 d2psids(6, 3) = -4.0 -
2899 12.0 * (d2_cubic_bubble013_ds3 + d2_cubic_bubble023_ds3) +
2900 32.0 * d2_quartic_bubble_ds3;
2901 d2psids(6, 4) = -4.0 -
2902 12.0 * (d2_cubic_bubble013_ds4 + d2_cubic_bubble023_ds4) +
2903 32.0 * d2_quartic_bubble_ds4;
2904 d2psids(6, 5) = 0.0 -
2905 12.0 * (d2_cubic_bubble013_ds5 + d2_cubic_bubble023_ds5) +
2906 32.0 * d2_quartic_bubble_ds5;
2907
2908 d2psids(7, 0) = 0.0 -
2909 12.0 * (d2_cubic_bubble012_ds0 + d2_cubic_bubble123_ds0) +
2910 32.0 * d2_quartic_bubble_ds0;
2911 d2psids(7, 1) = 0.0 -
2912 12.0 * (d2_cubic_bubble012_ds1 + d2_cubic_bubble123_ds1) +
2913 32.0 * d2_quartic_bubble_ds1;
2914 d2psids(7, 2) = 0.0 -
2915 12.0 * (d2_cubic_bubble012_ds2 + d2_cubic_bubble123_ds2) +
2916 32.0 * d2_quartic_bubble_ds2;
2917 d2psids(7, 3) = 0.0 -
2918 12.0 * (d2_cubic_bubble012_ds3 + d2_cubic_bubble123_ds3) +
2919 32.0 * d2_quartic_bubble_ds3;
2920 d2psids(7, 4) = 0.0 -
2921 12.0 * (d2_cubic_bubble012_ds4 + d2_cubic_bubble123_ds4) +
2922 32.0 * d2_quartic_bubble_ds4;
2923 d2psids(7, 5) = 4.0 -
2924 12.0 * (d2_cubic_bubble012_ds5 + d2_cubic_bubble123_ds5) +
2925 32.0 * d2_quartic_bubble_ds5;
2926
2927 d2psids(8, 0) = 0.0 -
2928 12.0 * (d2_cubic_bubble023_ds0 + d2_cubic_bubble123_ds0) +
2929 32.0 * d2_quartic_bubble_ds0;
2930 d2psids(8, 1) = 0.0 -
2931 12.0 * (d2_cubic_bubble023_ds1 + d2_cubic_bubble123_ds1) +
2932 32.0 * d2_quartic_bubble_ds1;
2933 d2psids(8, 2) = -8.0 -
2934 12.0 * (d2_cubic_bubble023_ds2 + d2_cubic_bubble123_ds2) +
2935 32.0 * d2_quartic_bubble_ds2;
2936 d2psids(8, 3) = 0.0 -
2937 12.0 * (d2_cubic_bubble023_ds3 + d2_cubic_bubble123_ds3) +
2938 32.0 * d2_quartic_bubble_ds3;
2939 d2psids(8, 4) = -4.0 -
2940 12.0 * (d2_cubic_bubble023_ds4 + d2_cubic_bubble123_ds4) +
2941 32.0 * d2_quartic_bubble_ds4;
2942 d2psids(8, 5) = -4.0 -
2943 12.0 * (d2_cubic_bubble023_ds5 + d2_cubic_bubble123_ds5) +
2944 32.0 * d2_quartic_bubble_ds5;
2945
2946 d2psids(9, 0) = 0.0 -
2947 12.0 * (d2_cubic_bubble013_ds0 + d2_cubic_bubble123_ds0) +
2948 32.0 * d2_quartic_bubble_ds0;
2949 d2psids(9, 1) = -8.0 -
2950 12.0 * (d2_cubic_bubble013_ds1 + d2_cubic_bubble123_ds1) +
2951 32.0 * d2_quartic_bubble_ds1;
2952 d2psids(9, 2) = 0.0 -
2953 12.0 * (d2_cubic_bubble013_ds2 + d2_cubic_bubble123_ds2) +
2954 32.0 * d2_quartic_bubble_ds3;
2955 d2psids(9, 3) = -4.0 -
2956 12.0 * (d2_cubic_bubble013_ds3 + d2_cubic_bubble123_ds3) +
2957 32.0 * d2_quartic_bubble_ds3;
2958 d2psids(9, 4) = 0.0 -
2959 12.0 * (d2_cubic_bubble013_ds4 + d2_cubic_bubble123_ds4) +
2960 32.0 * d2_quartic_bubble_ds4;
2961 d2psids(9, 5) = -4.0 -
2962 12.0 * (d2_cubic_bubble013_ds5 + d2_cubic_bubble123_ds5) +
2963 32.0 * d2_quartic_bubble_ds5;
2964
2965 // Add the bubble function on the face spanned by 0 1 3
2966 d2psids(10, 0) =
2967 27.0 * d2_cubic_bubble013_ds0 - 108.0 * d2_quartic_bubble_ds0;
2968 d2psids(10, 1) =
2969 27.0 * d2_cubic_bubble013_ds1 - 108.0 * d2_quartic_bubble_ds1;
2970 d2psids(10, 2) =
2971 27.0 * d2_cubic_bubble013_ds2 - 108.0 * d2_quartic_bubble_ds2;
2972 d2psids(10, 3) =
2973 27.0 * d2_cubic_bubble013_ds3 - 108.0 * d2_quartic_bubble_ds3;
2974 d2psids(10, 4) =
2975 27.0 * d2_cubic_bubble013_ds4 - 108.0 * d2_quartic_bubble_ds4;
2976 d2psids(10, 5) =
2977 27.0 * d2_cubic_bubble013_ds5 - 108.0 * d2_quartic_bubble_ds5;
2978
2979 // Add the bubble function on the face spanned by 0 1 2
2980 d2psids(11, 0) =
2981 27.0 * d2_cubic_bubble012_ds0 - 108.0 * d2_quartic_bubble_ds0;
2982 d2psids(11, 1) =
2983 27.0 * d2_cubic_bubble012_ds1 - 108.0 * d2_quartic_bubble_ds1;
2984 d2psids(11, 2) =
2985 27.0 * d2_cubic_bubble012_ds2 - 108.0 * d2_quartic_bubble_ds2;
2986 d2psids(11, 3) =
2987 27.0 * d2_cubic_bubble012_ds3 - 108.0 * d2_quartic_bubble_ds3;
2988 d2psids(11, 4) =
2989 27.0 * d2_cubic_bubble012_ds4 - 108.0 * d2_quartic_bubble_ds4;
2990 d2psids(11, 5) =
2991 27.0 * d2_cubic_bubble012_ds5 - 108.0 * d2_quartic_bubble_ds5;
2992
2993 // Add the bubble function on the face spanned by 0 2 3
2994 d2psids(12, 0) =
2995 27.0 * d2_cubic_bubble023_ds0 - 108.0 * d2_quartic_bubble_ds0;
2996 d2psids(12, 1) =
2997 27.0 * d2_cubic_bubble023_ds1 - 108.0 * d2_quartic_bubble_ds1;
2998 d2psids(12, 2) =
2999 27.0 * d2_cubic_bubble023_ds2 - 108.0 * d2_quartic_bubble_ds2;
3000 d2psids(12, 3) =
3001 27.0 * d2_cubic_bubble023_ds3 - 108.0 * d2_quartic_bubble_ds3;
3002 d2psids(12, 4) =
3003 27.0 * d2_cubic_bubble023_ds4 - 108.0 * d2_quartic_bubble_ds4;
3004 d2psids(12, 5) =
3005 27.0 * d2_cubic_bubble023_ds5 - 108.0 * d2_quartic_bubble_ds5;
3006
3007 // Add the bubble function on the face spanned by 1 2 3
3008 d2psids(13, 0) =
3009 27.0 * d2_cubic_bubble123_ds0 - 108.0 * d2_quartic_bubble_ds0;
3010 d2psids(13, 1) =
3011 27.0 * d2_cubic_bubble123_ds1 - 108.0 * d2_quartic_bubble_ds1;
3012 d2psids(13, 2) =
3013 27.0 * d2_cubic_bubble123_ds2 - 108.0 * d2_quartic_bubble_ds2;
3014 d2psids(13, 3) =
3015 27.0 * d2_cubic_bubble123_ds3 - 108.0 * d2_quartic_bubble_ds3;
3016 d2psids(13, 4) =
3017 27.0 * d2_cubic_bubble123_ds4 - 108.0 * d2_quartic_bubble_ds4;
3018 d2psids(13, 5) =
3019 27.0 * d2_cubic_bubble123_ds5 - 108.0 * d2_quartic_bubble_ds5;
3020
3021 // Add the volumetric bubble function derivatives
3022 d2psids(14, 0) = 256.0 * d2_quartic_bubble_ds0;
3023 d2psids(14, 1) = 256.0 * d2_quartic_bubble_ds1;
3024 d2psids(14, 2) = 256.0 * d2_quartic_bubble_ds2;
3025 d2psids(14, 3) = 256.0 * d2_quartic_bubble_ds3;
3026 d2psids(14, 4) = 256.0 * d2_quartic_bubble_ds4;
3027 d2psids(14, 5) = 256.0 * d2_quartic_bubble_ds5;
3028 }
3029 };
3030
3031
3032 //=======================================================================
3033 /// General TElement class specialised to three spatial dimensions (tet)
3034 /// Ordering of nodes inverted from Zienkiewizc sketches: When looking into
3035 /// the tet from vertex node 0. The vertex nodes on the opposite face are 1 -
3036 /// 2 - 3 in anticlockwise direction. Other nodes filled in edge by edge, then
3037 /// the face ones, then the internal ones.
3038 //=======================================================================
3039 template<unsigned NNODE_1D>
3040 class TElement<3, NNODE_1D> : public virtual TElementBase,
3041 public TElementShape<3, NNODE_1D>
3042 {
3043 private:
3044 /// Nodal translation scheme for use when generating face elements
3045 static const unsigned Node_on_face[4][(NNODE_1D * (NNODE_1D + 1)) / 2];
3046
3047 /// Default integration rule: Gaussian integration of same 'order' as
3048 /// the element
3049 // This is sort of optimal, because it means that the integration is exact
3050 // for the shape functions. Can overwrite this in specific element
3051 // defintion.
3052 static TGauss<3, NNODE_1D> Default_integration_scheme;
3053
3054 public:
3055 /// Constructor
3056 TElement()
3057 {
3058 switch (NNODE_1D)
3059 {
3060 case 2:
3061 case 3:
3062 break;
3063
3064 /* case 4: */
3065 /* n_node = 20; */
3066 /* break; */
3067
3068 default:
3069 std::string error_message =
3070 "Tets are currently only implemented for\n";
3071 error_message += "four and ten nodes, i.e. NNODE_1D=2 , 3 \n";
3072
3073 throw OomphLibError(
3074 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3075 }
3076
3077
3078 // Set the number of nodes
3079 unsigned n_node =
3080 (NNODE_1D * (NNODE_1D + 1)) / 2 + 1 + 3 * (NNODE_1D - 2);
3081 this->set_n_node(n_node);
3082
3083 // Set the elemental and nodal dimensions
3084 set_dimension(3);
3085
3086 // Assign default (full) spatial integration scheme
3087 set_integration_scheme(&Default_integration_scheme);
3088 }
3089
3090
3091 /// Broken copy constructor
3092 TElement(const TElement&) = delete;
3093
3094 /// Broken assignment operator
3095 /*void operator=(const TElement&) = delete;*/
3096
3097
3098 /// Destructor
3099 ~TElement() {}
3100
3101 /// Number of nodes along each element edge
3102 unsigned nnode_1d() const
3103 {
3104 return NNODE_1D;
3105 }
3106
3107
3108 /// Number of vertex nodes in the element: One more
3109 /// than spatial dimension
3110 unsigned nvertex_node() const
3111 {
3112 return 4;
3113 }
3114
3115 /// Public access function for Node_on_face.
3116 unsigned get_bulk_node_number(const int& face_index,
3117 const unsigned& i) const
3118 {
3119 return Node_on_face[face_index][i];
3120 }
3121
3122 /// Pointer to the j-th vertex node in the element
3123 Node* vertex_node_pt(const unsigned& j) const
3124 {
3125 // Vertex nodes come first:
3126#ifdef PARANOID
3127 if (j > 3)
3128 {
3129 std::ostringstream error_message;
3130 error_message
3131 << "Element only has four vertex nodes; called with node number " << j
3132 << std::endl;
3133 throw OomphLibError(error_message.str(),
3134 OOMPH_CURRENT_FUNCTION,
3135 OOMPH_EXCEPTION_LOCATION);
3136 }
3137#endif
3138 return node_pt(j);
3139 }
3140
3141 /// Calculate the geometric shape functions at local coordinate s
3142 inline void shape(const Vector<double>& s, Shape& psi) const
3143 {
3144 TElementShape<3, NNODE_1D>::shape(s, psi);
3145 }
3146
3147 /// Compute the geometric shape functions and
3148 /// derivatives w.r.t. local coordinates at local coordinate s
3149 inline void dshape_local(const Vector<double>& s,
3150 Shape& psi,
3151 DShape& dpsids) const
3152 {
3153 TElementShape<3, NNODE_1D>::dshape_local(s, psi, dpsids);
3154 }
3155
3156 /// Compute the geometric shape functions, derivatives and
3157 /// second derivatives w.r.t local coordinates at local coordinate s.
3158 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
3159 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
3160 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
3161 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
3162 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
3163 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
3164 inline void d2shape_local(const Vector<double>& s,
3165 Shape& psi,
3166 DShape& dpsids,
3167 DShape& d2psids) const
3168 {
3169 TElementShape<3, NNODE_1D>::d2shape_local(s, psi, dpsids, d2psids);
3170 }
3171
3172 /// Overload the template-free interface for the calculation of
3173 /// inverse jacobian matrix. This is a three dimensional element, so use
3174 /// the 3D version.
3175 double invert_jacobian_mapping(const DenseMatrix<double>& jacobian,
3176 DenseMatrix<double>& inverse_jacobian) const
3177 {
3178 return FiniteElement::invert_jacobian<3>(jacobian, inverse_jacobian);
3179 }
3180
3181 /// Min. value of local coordinate
3182 double s_min() const
3183 {
3184 return 0.0;
3185 }
3186
3187 /// Max. value of local coordinate
3188 double s_max() const
3189 {
3190 return 1.0;
3191 }
3192
3193 /// Return local coordinates of node j
3194 inline void local_coordinate_of_node(const unsigned& j,
3195 Vector<double>& s) const
3196 {
3197 TElementShape<3, NNODE_1D>::local_coordinate_of_node(j, s);
3198 }
3199
3200 /// Return the number of actual plot points for paraview
3201 /// plot with parameter nplot.
3202 unsigned nplot_points_paraview(const unsigned& nplot) const
3203 {
3204 unsigned node_sum = 0;
3205 for (unsigned j = 1; j <= nplot; j++)
3206 {
3207 for (unsigned i = 1; i <= j; i++)
3208 {
3209 node_sum += i;
3210 }
3211 }
3212 return node_sum;
3213 }
3214
3215 /// Return the number of local sub-elements for paraview plot with
3216 /// parameter nplot.
3217 unsigned nsub_elements_paraview(const unsigned& nplot) const
3218 {
3219 return (nplot - 1) * (nplot - 1) * (nplot - 1);
3220 }
3221
3222 /// Fill in the offset information for paraview plot.
3223 /// Needs to be implemented for each new geometric element type; see
3224 /// http://www.vtk.org/VTK/img/file-formats.pdf
3225 void write_paraview_output_offset_information(std::ofstream& file_out,
3226 const unsigned& nplot,
3227 unsigned& counter) const
3228 {
3229 // Output node lists for sub elements for Paraview (node index
3230 // must start at 0, fixed with magical counter-1)
3231
3232 unsigned paraview_fix = counter - 1;
3233 unsigned nod_count = 1;
3234
3235 for (unsigned i = 0; i < nplot; i++)
3236 {
3237 for (unsigned j = 0; j < nplot - i; j++)
3238 {
3239 for (unsigned k = 0; k < nplot - i - j; k++)
3240 {
3241 if (k < nplot - i - j - 1)
3242 {
3243 file_out << nod_count + paraview_fix << " "
3244 << nod_count + 1 + paraview_fix << " "
3245 << nod_count + nplot - i - j + paraview_fix << " "
3246 << nod_count + nplot - i - j +
3247 ((nplot - 1 - i) * (nplot - i) / 2) + paraview_fix
3248 << std::endl;
3249 if (k < nplot - i - j - 2)
3250 {
3251 file_out << nod_count + 1 + paraview_fix << " "
3252 << nod_count + nplot - i - j + paraview_fix << " "
3253 << nod_count + nplot - i - j +
3254 ((nplot - 1 - i) * (nplot - i) / 2) + paraview_fix
3255 << " "
3256 << nod_count + 2 * (nplot - i - j) - 1 +
3257 ((nplot - 1 - i) * (nplot - i) / 2) + paraview_fix
3258 << std::endl;
3259 file_out << nod_count + 1 + paraview_fix << " "
3260 << nod_count + nplot - i - j + paraview_fix << " "
3261 << nod_count + nplot - i - j + 1 + paraview_fix << " "
3262 << nod_count + 2 * (nplot - i - j) - 1 +
3263 ((nplot - 1 - i) * (nplot - i) / 2) + paraview_fix
3264 << std::endl;
3265 file_out << nod_count + 1 + paraview_fix << " "
3266 << nod_count + nplot - i - j +
3267 ((nplot - 1 - i) * (nplot - i) / 2) + paraview_fix
3268 << " "
3269 << nod_count + nplot - i - j +
3270 ((nplot - 1 - i) * (nplot - i) / 2) + 1 +
3271 paraview_fix
3272 << " "
3273 << nod_count + 2 * (nplot - i - j) - 1 +
3274 ((nplot - 1 - i) * (nplot - i) / 2) + paraview_fix
3275 << std::endl;
3276 file_out << nod_count + 1 + paraview_fix << " "
3277 << nod_count + nplot - i - j + 1 + paraview_fix << " "
3278 << nod_count + nplot - i - j +
3279 ((nplot - 1 - i) * (nplot - i) / 2) + 1 +
3280 paraview_fix
3281 << " "
3282 << nod_count + 2 * (nplot - i - j) - 1 +
3283 ((nplot - 1 - i) * (nplot - i) / 2) + paraview_fix
3284 << std::endl;
3285 }
3286 if (k > 1)
3287 {
3288 file_out << nod_count + nplot - i - j - 1 + paraview_fix << " "
3289 << nod_count + nplot - i - j +
3290 ((nplot - 1 - i) * (nplot - i) / 2) - 1 +
3291 paraview_fix
3292 << " "
3293 << nod_count + 2 * (nplot - i - j - 1) +
3294 ((nplot - 1 - i) * (nplot - i) / 2) - 1 +
3295 paraview_fix
3296 << " "
3297 << nod_count + 2 * (nplot - i - j - 1) +
3298 ((nplot - 1 - i) * (nplot - i) / 2) + paraview_fix
3299 << std::endl;
3300 }
3301 } // end if(k<nplot-i-j-1)
3302 ++nod_count;
3303 }
3304 }
3305 }
3306
3307 // increment the counter to keep track of global connectivity
3308 counter += nplot_points_paraview(nplot);
3309
3310 } // end of write Paraview_element...
3311
3312 /// Return the paraview element type.
3313 /// Needs to be implemented for each new geometric element type; see
3314 /// http://www.vtk.org/VTK/img/file-formats.pdf
3315 void write_paraview_type(std::ofstream& file_out,
3316 const unsigned& nplot) const
3317 {
3318 unsigned local_loop = nsub_elements_paraview(nplot);
3319 for (unsigned i = 0; i < local_loop; i++)
3320 {
3321 file_out << "10" << std::endl;
3322 }
3323 }
3324
3325 /// Return the offsets for the paraview sub-elements. Needs
3326 /// to be implemented for each new geometric element type; see
3327 /// http://www.vtk.org/VTK/img/file-formats.pdf
3328 void write_paraview_offsets(std::ofstream& file_out,
3329 const unsigned& nplot,
3330 unsigned& offset_sum) const
3331 {
3332 unsigned local_loop = nsub_elements_paraview(nplot);
3333 for (unsigned i = 0; i < local_loop; i++)
3334 {
3335 offset_sum += 4;
3336 file_out << offset_sum << std::endl;
3337 }
3338 }
3339
3340 /// Output
3341 void output(std::ostream& output);
3342
3343 /// Output at specified number of plot points
3344 void output(std::ostream& outfile, const unsigned& nplot);
3345
3346 /// C-style output
3347 void output(FILE* file_pt);
3348
3349 /// C_style output at n_plot points
3350 void output(FILE* file_pt, const unsigned& n_plot);
3351
3352 /// Get vector of local coordinates of plot point i (when plotting
3353 /// nplot points in each "coordinate direction).
3354 void get_s_plot(
3355 const unsigned& iplot,
3356 const unsigned& nplot,
3357 Vector<double>& s,
3358 const bool& use_equally_spaced_interior_sample_points = false) const
3359 {
3360 if (nplot > 1)
3361 {
3362 unsigned np = 0;
3363 for (unsigned i = 0; i < nplot; ++i)
3364 {
3365 for (unsigned j = 0; j < nplot - i; ++j)
3366 {
3367 for (unsigned k = 0; k < nplot - i - j; ++k)
3368 {
3369 if (np == iplot)
3370 {
3371 {
3372 s[0] = double(j) / double(nplot - 1);
3373 s[1] = double(i) / double(nplot - 1);
3374 s[2] = double(k) / double(nplot - 1);
3375 if (use_equally_spaced_interior_sample_points)
3376 {
3377 double range = 1.0;
3378 double dx_new = range / double(nplot + 1);
3379 double range_new = double(nplot - 1) * dx_new;
3380 s[0] = 0.5 * dx_new + range_new * s[0] / range;
3381 s[1] = 0.5 * dx_new + range_new * s[1] / range;
3382 s[2] = 0.5 * dx_new + range_new * s[2] / range;
3383 }
3384 return;
3385 }
3386 }
3387 np++;
3388 }
3389 }
3390 }
3391 }
3392 else
3393 {
3394 s[0] = 1.0 / 4.0;
3395 s[1] = 1.0 / 4.0;
3396 s[2] = 1.0 / 4.0;
3397 }
3398 }
3399
3400 /// Return string for tecplot zone header (when plotting
3401 /// nplot points in each "coordinate direction)
3402 std::string tecplot_zone_string(const unsigned& nplot) const
3403 {
3404 std::ostringstream header;
3405 unsigned nel = 0;
3406 nel = (nplot - 1) * (nplot - 1) * (nplot - 1);
3407 header << "ZONE N=" << nplot_points(nplot) << ", E=" << nel
3408 << ", F=FEPOINT, ET=TETRAHEDRON\n";
3409 return header.str();
3410 }
3411
3412 /// Add tecplot zone "footer" to output stream (when plotting
3413 /// nplot points in each "coordinate direction).
3414 /// Empty by default -- can be used, e.g., to add FE connectivity
3415 /// lists to elements that need it.
3416 void write_tecplot_zone_footer(std::ostream& outfile,
3417 const unsigned& nplot) const
3418 {
3419 // Output node lists for sub elements for Tecplot (node index
3420 // must start at 1)
3421 unsigned nod_count = 1;
3422 for (unsigned i = 0; i < nplot; i++)
3423 {
3424 for (unsigned j = 0; j < nplot - i; j++)
3425 {
3426 for (unsigned k = 0; k < nplot - i - j; k++)
3427 {
3428 if (k < nplot - i - j - 1)
3429 {
3430 outfile << nod_count << " " << nod_count + 1 << " "
3431 << nod_count + nplot - i - j << " "
3432 << nod_count + nplot - i - j +
3433 ((nplot - 1 - i) * (nplot - i) / 2)
3434 << std::endl;
3435 if (k < nplot - i - j - 2)
3436 {
3437 outfile << nod_count + 1 << " " << nod_count + nplot - i - j
3438 << " "
3439 << nod_count + nplot - i - j +
3440 ((nplot - 1 - i) * (nplot - i) / 2)
3441 << " "
3442 << nod_count + 2 * (nplot - i - j) - 1 +
3443 ((nplot - 1 - i) * (nplot - i) / 2)
3444 << std::endl;
3445 outfile << nod_count + 1 << " " << nod_count + nplot - i - j
3446 << " " << nod_count + nplot - i - j + 1 << " "
3447 << nod_count + 2 * (nplot - i - j) - 1 +
3448 ((nplot - 1 - i) * (nplot - i) / 2)
3449 << std::endl;
3450 outfile << nod_count + 1 << " "
3451 << nod_count + nplot - i - j +
3452 ((nplot - 1 - i) * (nplot - i) / 2)
3453 << " "
3454 << nod_count + nplot - i - j +
3455 ((nplot - 1 - i) * (nplot - i) / 2) + 1
3456 << " "
3457 << nod_count + 2 * (nplot - i - j) - 1 +
3458 ((nplot - 1 - i) * (nplot - i) / 2)
3459 << std::endl;
3460 outfile << nod_count + 1 << " " << nod_count + nplot - i - j + 1
3461 << " "
3462 << nod_count + nplot - i - j +
3463 ((nplot - 1 - i) * (nplot - i) / 2) + 1
3464 << " "
3465 << nod_count + 2 * (nplot - i - j) - 1 +
3466 ((nplot - 1 - i) * (nplot - i) / 2)
3467 << std::endl;
3468 }
3469 if (k > 1)
3470 {
3471 outfile << nod_count + nplot - i - j - 1 << " "
3472 << nod_count + nplot - i - j +
3473 ((nplot - 1 - i) * (nplot - i) / 2) - 1
3474 << " "
3475 << nod_count + 2 * (nplot - i - j - 1) +
3476 ((nplot - 1 - i) * (nplot - i) / 2) - 1
3477 << " "
3478 << nod_count + 2 * (nplot - i - j - 1) +
3479 ((nplot - 1 - i) * (nplot - i) / 2)
3480 << std::endl;
3481 }
3482 } // end if(k<nplot-i-j-1)
3483 ++nod_count;
3484 }
3485 }
3486 }
3487 } // end of write tecplot...
3488
3489
3490 /// Add tecplot zone "footer" to C-style output. (when plotting
3491 /// nplot points in each "coordinate direction).
3492 /// Empty by default -- can be used, e.g., to add FE connectivity
3493 /// lists to elements that need it.
3494 void write_tecplot_zone_footer(FILE* file_pt, const unsigned& nplot) const
3495 {
3496 // Output node lists for sub elements for Tecplot (node index
3497 // must start at 1)
3498 unsigned nod_count = 1;
3499 for (unsigned i = 0; i < nplot; i++)
3500 {
3501 for (unsigned j = 0; j < nplot - i; j++)
3502 {
3503 for (unsigned k = 0; k < nplot - i - j; k++)
3504 {
3505 if (j < nplot - i - 1)
3506 {
3507 fprintf(file_pt,
3508 "%i %i %i \n",
3509 nod_count,
3510 nod_count + 1,
3511 nod_count + nplot - i);
3512 if (j < nplot - i - 2)
3513 {
3514 fprintf(file_pt,
3515 "%i %i %i \n",
3516 nod_count + 1,
3517 nod_count + nplot - i + 1,
3518 nod_count + nplot - i);
3519 }
3520 }
3521 ++nod_count;
3522 }
3523 }
3524 }
3525 }
3526
3527 /// Return total number of plot points (when plotting
3528 /// nplot points in each "coordinate direction)
3529 unsigned nplot_points(const unsigned& nplot) const
3530 {
3531 unsigned res = 0;
3532 if (nplot > 1)
3533 {
3534 res = 4;
3535 for (unsigned i = 3; i <= nplot; i++)
3536 {
3537 res += (i * (i + 1) / 2);
3538 }
3539 return res;
3540 }
3541 // Otherwise we return 1(?)
3542 return 1;
3543 }
3544
3545 /// Build the lower-dimensional FaceElement (an element of type
3546 /// TElement<2,NNODE_1D>). The face index can take one of four values
3547 /// corresponding to the four possible faces:
3548 /// 0: (left) s[0] = 0.0
3549 /// 1: (bottom) s[1] = 0.0
3550 /// 2: (back) s[2] = 0.0
3551 /// 3: (sloping face) s[0] + s[1] + s[2] = 1
3552 void build_face_element(const int& face_index,
3553 FaceElement* face_element_pt);
3554 };
3555
3556
3557 ///////////////////////////////////////////////////////////////////////
3558 ///////////////////////////////////////////////////////////////////////
3559 ///////////////////////////////////////////////////////////////////////
3560
3561
3562 //=======================================================================
3563 /// TElement class for which the shape functions have been enriched
3564 /// by a single bubble function of the next order
3565 ///
3566 /// Empty, just establishes the template parameters
3567 //=======================================================================
3568 template<unsigned DIM, unsigned NNODE_1D>
3569 class TBubbleEnrichedElement
3570 {
3571 };
3572
3573 //=====================================================================
3574 /// Define integration schemes that are required to exactly integrate
3575 /// the mass matrices of the bubble-enriched elements. The enrichement
3576 /// increases the polynomial order which means that higher-order Gauss
3577 /// rules must be used.
3578 //====================================================================
3579 template<unsigned DIM, unsigned NPTS_1D>
3580 class TBubbleEnrichedGauss
3581 {
3582 };
3583
3584 //====================================================================
3585 /// Specialisation for two-dimensional elements, in which the highest
3586 /// order polynomial is cubic, so we need the integration scheme
3587 /// for the unenriched cubic element
3588 //======================================================================
3589 template<>
3590 class TBubbleEnrichedGauss<2, 3> : public TGauss<2, 4>
3591 {
3592 public:
3593 TBubbleEnrichedGauss() : TGauss<2, 4>() {}
3594 };
3595
3596 //====================================================================
3597 /// Specialisation for three-dimensional elements, in which the highest
3598 /// order polynomial is quartic, so we need the integration scheme
3599 /// for the unenriched quartic element
3600 //======================================================================
3601 template<>
3602 class TBubbleEnrichedGauss<3, 3> : public TGauss<3, 5>
3603 {
3604 public:
3605 TBubbleEnrichedGauss() : TGauss<3, 5>() {}
3606 };
3607
3608
3609 //=======================================================================
3610 /// Enriched TElement class specialised to two spatial dimensions
3611 /// and three nodes per side (quadratic element)
3612 /// Ordering of nodes as in Zienkiwizc sketches: vertex nodes
3613 /// 0 - 1 - 2 anticlockwise. Midside nodes filled in progressing
3614 /// along the consecutive edges. Central node(s) come(s) last.
3615 /// The idea is that we inherit from the existing TElement<2,3>, add
3616 /// the single extra node at the centroid and
3617 /// overload the shape functions to be those corresponding to the
3618 /// enriched element.
3619 //=======================================================================
3620 template<unsigned DIM>
3621 class TBubbleEnrichedElement<DIM, 3>
3622 : public virtual TElement<DIM, 3>,
3623 public TBubbleEnrichedElementShape<DIM, 3>
3624 {
3625 private:
3626 // Static storage for a new integration scheme
3627 static TBubbleEnrichedGauss<DIM, 3> Default_enriched_integration_scheme;
3628
3629 // Static storage for central node
3630 static const unsigned Central_node_on_face[DIM + 1];
3631
3632 public:
3633 /// Constructor
3634 TBubbleEnrichedElement()
3635 : TElement<DIM, 3>(), TBubbleEnrichedElementShape<DIM, 3>()
3636 {
3637 // Add the additional enrichment nodes
3638 unsigned n_node = this->nnode();
3639 this->set_n_node(n_node +
3640 TBubbleEnrichedElementShape<DIM, 3>::n_enriched_nodes());
3641 // Set the new integration scheme
3642 this->set_integration_scheme(&Default_enriched_integration_scheme);
3643 }
3644
3645 /// Broken copy constructor
3646 TBubbleEnrichedElement(const TBubbleEnrichedElement&) = delete;
3647
3648 /// Broken assignment operator
3649 /*void operator=(const TBubbleEnrichedElement&) = delete;*/
3650
3651 /// Destructor
3652 ~TBubbleEnrichedElement() {}
3653
3654 /// Calculate the geometric shape functions at local coordinate s
3655 inline void shape(const Vector<double>& s, Shape& psi) const
3656 {
3657 TBubbleEnrichedElementShape<DIM, 3>::shape(s, psi);
3658 }
3659
3660 /// Compute the geometric shape functions and
3661 /// derivatives w.r.t. local coordinates at local coordinate s
3662 inline void dshape_local(const Vector<double>& s,
3663 Shape& psi,
3664 DShape& dpsids) const
3665 {
3666 TBubbleEnrichedElementShape<DIM, 3>::dshape_local(s, psi, dpsids);
3667 }
3668
3669
3670 /// Compute the geometric shape functions, derivatives and
3671 /// second derivatives w.r.t local coordinates at local coordinate s
3672 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
3673 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
3674 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
3675 inline void d2shape_local(const Vector<double>& s,
3676 Shape& psi,
3677 DShape& dpsids,
3678 DShape& d2psids) const
3679 {
3680 TBubbleEnrichedElementShape<DIM, 3>::d2shape_local(
3681 s, psi, dpsids, d2psids);
3682 }
3683
3684 /// Return local coordinates of node j
3685 inline void local_coordinate_of_node(const unsigned& j,
3686 Vector<double>& s) const
3687 {
3688 TBubbleEnrichedElementShape<DIM, 3>::local_coordinate_of_node(j, s);
3689 }
3690
3691 /// Build the lower-dimensional FaceElement
3692 void build_face_element(const int& face_index,
3693 FaceElement* face_element_pt);
3694 };
3695
3696
3697 //========================================================================
3698 /// Base class for Solid Telements
3699 //========================================================================
3700 class TSolidElementBase : public virtual TElementBase,
3701 public virtual SolidFiniteElement
3702 {
3703 public:
3704 /// Constructor: Empty
3705 TSolidElementBase() {};
3706
3707 /// Broken copy constructor
3708 TSolidElementBase(const TSolidElementBase&) = delete;
3709
3710 /// Broken assignment operator
3711 /*void operator=(const TSolidElementBase&) = delete;*/
3712 };
3713
3714
3715 //////////////////////////////////////////////////////////////////////////
3716 //////////////////////////////////////////////////////////////////////////
3717 //////////////////////////////////////////////////////////////////////////
3718
3719
3720 //=======================================================================
3721 /// SolidTElement elements are triangular/tet elements whose
3722 /// derivatives also include those based upon the lagrangian
3723 /// positions of the nodes.
3724 /// They are the basis for solid mechanics elements.
3725 //=======================================================================
3726 template<unsigned DIM, unsigned NNODE_1D>
3727 class SolidTElement
3728 {
3729 };
3730
3731
3732 //=======================================================================
3733 /// SolidTElement elements, specialised to one spatial dimension
3734 //=======================================================================
3735 template<unsigned NNODE_1D>
3736 class SolidTElement<1, NNODE_1D> : public virtual TElement<1, NNODE_1D>,
3737 public virtual TSolidElementBase
3738 {
3739 public:
3740 /// Constructor
3741 SolidTElement() : TElement<1, NNODE_1D>(), SolidFiniteElement()
3742 {
3743 // Set the Lagrangian dimension of the element
3744 set_lagrangian_dimension(1);
3745 }
3746
3747 /// Broken copy constructor
3748 SolidTElement(const SolidTElement&) = delete;
3749
3750 /// Broken assignment operator
3751 /*void operator=(const SolidTElement&) = delete;*/
3752
3753 /// Build the lower-dimensional FaceElement (an element of type
3754 /// SolidPointElement). The face index takes two values
3755 /// corresponding to the two possible faces:
3756 /// -1 (Left) s[0] = -1.0
3757 /// +1 (Right) s[0] = 1.0
3758 inline void build_face_element(const int& face_index,
3759 FaceElement* face_element_pt);
3760 };
3761
3762
3763 ///////////////////////////////////////////////////////////////////////////
3764 ///////////////////////////////////////////////////////////////////////////
3765 // SolidTElements
3766 ///////////////////////////////////////////////////////////////////////////
3767 ///////////////////////////////////////////////////////////////////////////
3768
3769
3770 ///////////////////////////////////////////////////////////////////////////
3771 // 1D SolidTElements
3772 ///////////////////////////////////////////////////////////////////////////
3773
3774
3775 //===========================================================
3776 /// Function to setup geometrical information for lower-dimensional
3777 /// FaceElements (which are of type SolidTElement<0,1>).
3778 //===========================================================
3779 template<unsigned NNODE_1D>
3780 void SolidTElement<1, NNODE_1D>::build_face_element(
3781 const int& face_index, FaceElement* face_element_pt)
3782 {
3783 // Build the standard non-solid FaceElement
3784 TElement<1, NNODE_1D>::build_face_element(face_index, face_element_pt);
3785
3786 // Set the Lagrangian dimension from the first node of the present element
3787 dynamic_cast<SolidFiniteElement*>(face_element_pt)
3788 ->set_lagrangian_dimension(
3789 static_cast<SolidNode*>(node_pt(0))->nlagrangian());
3790 }
3791
3792
3793 //=======================================================================
3794 /// SolidTElement elements, specialised to two spatial dimensions
3795 //=======================================================================
3796 template<unsigned NNODE_1D>
3797 class SolidTElement<2, NNODE_1D> : public virtual TElement<2, NNODE_1D>,
3798 public virtual TSolidElementBase
3799 {
3800 public:
3801 /// Constructor
3802 SolidTElement()
3803 : TElementBase(),
3804 TElement<2, NNODE_1D>(),
3805 SolidFiniteElement(),
3806 TSolidElementBase()
3807 {
3808 // Set the Lagrangian dimension of the element
3809 set_lagrangian_dimension(2);
3810 }
3811
3812 /// Broken copy constructor
3813 SolidTElement(const SolidTElement&) = delete;
3814
3815 /// Broken assignment operator
3816 /*void operator=(const SolidTElement&) = delete;*/
3817
3818 /// Build the lower-dimensional FaceElement (an element of type
3819 /// SolidTElement<1,NNODE_1D>). The face index takes three possible values:
3820 /// 0 (Left) s[0] = 0.0
3821 /// 1 (Bottom) s[1] = 0.0
3822 /// 2 (Sloping face) s[0] = 1.0 - s[1]
3823 inline void build_face_element(const int& face_index,
3824 FaceElement* face_element_pt);
3825 };
3826
3827
3828 ///////////////////////////////////////////////////////////////////////////
3829 ///////////////////////////////////////////////////////////////////////////
3830 // 2D SolidTElements
3831 ///////////////////////////////////////////////////////////////////////////
3832 ///////////////////////////////////////////////////////////////////////////
3833
3834
3835 //===========================================================
3836 /// Function to setup geometrical information for lower-dimensional
3837 /// FaceElements (which are of type SolidTElement<1,NNODE_1D>).
3838 //===========================================================
3839 template<unsigned NNODE_1D>
3840 void SolidTElement<2, NNODE_1D>::build_face_element(
3841 const int& face_index, FaceElement* face_element_pt)
3842 {
3843 // Build the standard non-solid FaceElement
3844 TElement<2, NNODE_1D>::build_face_element(face_index, face_element_pt);
3845
3846 // Set the Lagrangian dimension from the first node of the present element
3847 dynamic_cast<SolidFiniteElement*>(face_element_pt)
3848 ->set_lagrangian_dimension(
3849 static_cast<SolidNode*>(node_pt(0))->nlagrangian());
3850 }
3851
3852
3853 //=======================================================================
3854 /// SolidTElement elements, specialised to three spatial dimensions
3855 //=======================================================================
3856 template<unsigned NNODE_1D>
3857 class SolidTElement<3, NNODE_1D> : public virtual TElement<3, NNODE_1D>,
3858 public virtual TSolidElementBase
3859 {
3860 public:
3861 /// Constructor
3862 SolidTElement()
3863 : TElementBase(),
3864 TElement<3, NNODE_1D>(),
3865 SolidFiniteElement(),
3866 TSolidElementBase()
3867 {
3868 // Set the Lagrangian dimension of the element
3869 set_lagrangian_dimension(3);
3870 }
3871
3872 /// Broken copy constructor
3873 SolidTElement(const SolidTElement&) = delete;
3874
3875 /// Broken assignment operator
3876 /*void operator=(const SolidTElement&) = delete;*/
3877
3878
3879 /// Build the lower-dimensional FaceElement (an element of type
3880 /// SolidTElement<2,NNODE_1D>). The face index can take one of four values
3881 /// corresponding to the four possible faces:
3882 /// 0: (left) s[0] = 0.0
3883 /// 1: (bottom) s[1] = 0.0
3884 /// 2: (back) s[2] = 0.0
3885 /// 3: (sloping face) s[0] + s[1] + s[2] = 1
3886 inline void build_face_element(const int& face_index,
3887 FaceElement* face_element_pt);
3888 };
3889
3890
3891 ///////////////////////////////////////////////////////////////////////////
3892 // 3D SolidTElements
3893 ///////////////////////////////////////////////////////////////////////////
3894
3895
3896 //===========================================================
3897 /// Function to setup geometrical information for lower-dimensional
3898 /// FaceElements (which are of type SolidTElement<1,NNODE_1D>).
3899 //===========================================================
3900 template<unsigned NNODE_1D>
3901 void SolidTElement<3, NNODE_1D>::build_face_element(
3902 const int& face_index, FaceElement* face_element_pt)
3903 {
3904 // Build the standard non-solid FaceElement
3905 TElement<3, NNODE_1D>::build_face_element(face_index, face_element_pt);
3906
3907 // Set the Lagrangian dimension from the first node of the present element
3908 dynamic_cast<SolidFiniteElement*>(face_element_pt)
3909 ->set_lagrangian_dimension(
3910 static_cast<SolidNode*>(node_pt(0))->nlagrangian());
3911 }
3912
3913
3914 ///////////////////////////////////////////////////////////////////////
3915 ///////////////////////////////////////////////////////////////////////
3916 ///////////////////////////////////////////////////////////////////////
3917
3918 //=======================================================================
3919 /// SolidTBubbleEnrichedElement elements are the enriched version
3920 /// of the SolidTElements. They will simply inherit from the appropriate
3921 /// SolidTElement and TBubblEnrichedElement.
3922 /// They are the basis for solid mechanics elements.
3923 //=======================================================================
3924 template<unsigned DIM, unsigned NNODE_1D>
3925 class SolidTBubbleEnrichedElement
3926 {
3927 };
3928
3929 //===================================================================
3930 /// Specify the SolidTBubbleEnrichedElement corresponding to the
3931 /// quadratic triangle
3932 //===================================================================
3933 template<unsigned DIM>
3934 class SolidTBubbleEnrichedElement<DIM, 3>
3935 : public virtual SolidTElement<DIM, 3>,
3936 public virtual TBubbleEnrichedElement<DIM, 3>
3937 {
3938 public:
3939 /// Constructor
3944
3945 /// Broken copy constructor
3946 SolidTBubbleEnrichedElement(const SolidTBubbleEnrichedElement&) = delete;
3947
3948 /// Broken assignment operator
3949 /*void operator=(const SolidTBubbleEnrichedElement&) = delete;*/
3950
3951 /// Destructor
3952 ~SolidTBubbleEnrichedElement() {}
3953
3954 /// Build the lower-dimensional FaceElement
3955 /// Need to put in a final override here
3956 void build_face_element(const int& face_index,
3957 FaceElement* face_element_pt);
3958 };
3959
3960
3961 //=======================================================================
3962 /// Face geometry for the TElement elements: The spatial
3963 /// dimension of the face elements is one lower than that of the
3964 /// bulk element but they have the same number of points
3965 /// along their 1D edges.
3966 //=======================================================================
3967 template<unsigned DIM, unsigned NNODE_1D>
3968 class FaceGeometry<TElement<DIM, NNODE_1D>>
3969 : public virtual TElement<DIM - 1, NNODE_1D>
3970 {
3971 public:
3972 /// Constructor: Call the constructor for the
3973 /// appropriate lower-dimensional QElement
3974 FaceGeometry() : TElement<DIM - 1, NNODE_1D>() {}
3975 };
3976
3977
3978 //=======================================================================
3979 /// Face geometry for the 1D TElement elements: Point elements
3980 //=======================================================================
3981 template<unsigned NNODE_1D>
3982 class FaceGeometry<TElement<1, NNODE_1D>> : public virtual PointElement
3983 {
3984 public:
3985 /// Constructor: Call the constructor for the
3986 /// appropriate lower-dimensional TElement
3987 FaceGeometry() : PointElement() {}
3988 };
3989
3990
3991 ///////////////////////////////////////////////////////////////////////
3992 ///////////////////////////////////////////////////////////////////////
3993 ///////////////////////////////////////////////////////////////////////
3994
3995
3996 //=======================================================================
3997 /// Face geometry for the 2D TBubbleEnrichedElement elements is exactly
3998 /// the same as for the corresponding TElement. The spatial
3999 /// dimension of the face elements is one lower than that of the
4000 /// bulk element but they have the same number of points
4001 /// along their 1D edges.
4002 //=======================================================================
4003 template<unsigned NNODE_1D>
4004 class FaceGeometry<TBubbleEnrichedElement<2, NNODE_1D>>
4005 : public virtual TElement<1, NNODE_1D>
4006 {
4007 public:
4008 /// Constructor: Call the constructor for the
4009 /// appropriate lower-dimensional QElement
4010 FaceGeometry() : TElement<1, NNODE_1D>() {}
4011 };
4012
4013
4014 //=======================================================================
4015 /// Face geometry for the 3D TBubbleEnrichedElement elements is the
4016 /// 2D TBubbleEnrichedElement. The spatial
4017 /// dimension of the face elements is one lower than that of the
4018 /// bulk element but they have the same number of points
4019 /// along their 1D edges.
4020 //=======================================================================
4021 template<unsigned NNODE_1D>
4022 class FaceGeometry<TBubbleEnrichedElement<3, NNODE_1D>>
4023 : public virtual TBubbleEnrichedElement<2, NNODE_1D>
4024 {
4025 public:
4026 /// Constructor: Call the constructor for the
4027 /// appropriate lower-dimensional QElement
4028 FaceGeometry() : TBubbleEnrichedElement<2, NNODE_1D>() {}
4029 };
4030
4031
4032 ///////////////////////////////////////////////////////////////////////
4033 ///////////////////////////////////////////////////////////////////////
4034 ///////////////////////////////////////////////////////////////////////
4035
4036
4037 //=======================================================================
4038 /// Face geometry for the TElement elements: The spatial
4039 /// dimension of the face elements is one lower than that of the
4040 /// bulk element but they have the same number of points
4041 /// along their 1D edges.
4042 //=======================================================================
4043 template<unsigned DIM, unsigned NNODE_1D>
4044 class FaceGeometry<SolidTElement<DIM, NNODE_1D>>
4045 : public virtual SolidTElement<DIM - 1, NNODE_1D>
4046 {
4047 public:
4048 /// Constructor: Call the constructor for the
4049 /// appropriate lower-dimensional QElement
4050 FaceGeometry() : SolidTElement<DIM - 1, NNODE_1D>() {}
4051 };
4052
4053
4054 //=======================================================================
4055 /// Face geometry for the 1D TElement elements: Point elements
4056 //=======================================================================
4057 template<unsigned NNODE_1D>
4058 class FaceGeometry<SolidTElement<1, NNODE_1D>>
4059 : public virtual SolidPointElement
4060 {
4061 public:
4062 /// Constructor: Call the constructor for the
4063 /// appropriate lower-dimensional TElement
4064 FaceGeometry() : SolidPointElement() {}
4065 };
4066
4067
4068 ///////////////////////////////////////////////////////////////////////
4069 ///////////////////////////////////////////////////////////////////////
4070 ///////////////////////////////////////////////////////////////////////
4071
4072
4073 //=======================================================================
4074 /// Face geometry for the 2D SolidTBubbleEnrichedElement elements is exactly
4075 /// the same as for the corresponding 2D SolidTElement. The spatial
4076 /// dimension of the face elements is one lower than that of the
4077 /// bulk element but they have the same number of points
4078 /// along their 1D edges.
4079 //=======================================================================
4080 template<unsigned NNODE_1D>
4081 class FaceGeometry<SolidTBubbleEnrichedElement<2, NNODE_1D>>
4082 : public virtual SolidTElement<1, NNODE_1D>
4083 {
4084 public:
4085 /// Constructor: Call the constructor for the
4086 /// appropriate lower-dimensional QElement
4087 FaceGeometry() : SolidTElement<1, NNODE_1D>() {}
4088 };
4089
4090
4091 //=======================================================================
4092 /// Face geometry for the 3D SolidTBubbleEnrichedElement elements is
4093 /// the 2D SolidTBubbleEnrichedElement. The spatial
4094 /// dimension of the face elements is one lower than that of the
4095 /// bulk element but they have the same number of points
4096 /// along their 1D edges.
4097 //=======================================================================
4098 template<unsigned NNODE_1D>
4099 class FaceGeometry<SolidTBubbleEnrichedElement<3, NNODE_1D>>
4100 : public virtual SolidTBubbleEnrichedElement<2, NNODE_1D>
4101 {
4102 public:
4103 /// Constructor: Call the constructor for the
4104 /// appropriate lower-dimensional QElement
4105 FaceGeometry() : SolidTBubbleEnrichedElement<2, NNODE_1D>() {}
4106 };
4107
4108} // namespace oomph
4109
4110
4111#endif
s
static char t char * s
Definition cfortran.h:568
i
cstr elem_len * i
Definition cfortran.h:603
oomph::BoundaryNodeBase
A class that contains the information required by Nodes that are located on Mesh boundaries....
Definition nodes.h:1996
oomph::DShape
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition shape.h:278
oomph::FaceElement
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition elements.h:4342
oomph::FaceGeometry< SolidTBubbleEnrichedElement< 2, NNODE_1D > >::FaceGeometry
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
Definition Telements.h:4087
oomph::FaceGeometry< SolidTBubbleEnrichedElement< 3, NNODE_1D > >::FaceGeometry
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
Definition Telements.h:4105
oomph::FaceGeometry< SolidTElement< 1, NNODE_1D > >::FaceGeometry
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional TElement.
Definition Telements.h:4064
oomph::FaceGeometry< SolidTElement< DIM, NNODE_1D > >::FaceGeometry
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
Definition Telements.h:4050
oomph::FaceGeometry< TBubbleEnrichedElement< 2, NNODE_1D > >::FaceGeometry
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
Definition Telements.h:4010
oomph::FaceGeometry< TBubbleEnrichedElement< 3, NNODE_1D > >::FaceGeometry
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
Definition Telements.h:4028
oomph::FaceGeometry< TElement< 1, NNODE_1D > >::FaceGeometry
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional TElement.
Definition Telements.h:3987
oomph::FaceGeometry< TElement< DIM, NNODE_1D > >::FaceGeometry
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
Definition Telements.h:3974
oomph::FaceGeometry
FaceGeometry class definition: This policy class is used to allow construction of face elements that ...
Definition elements.h:5002
oomph::FiniteElement
A general Finite Element class.
Definition elements.h:1317
oomph::FiniteElement::dim
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition elements.h:2615
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
oomph::Node
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition nodes.h:906
oomph::Node::is_on_boundary
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary,...
Definition nodes.h:1373
oomph::Node::get_boundaries_pt
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Return a pointer to set of mesh boundaries that this node occupies; this will be overloaded by Bounda...
Definition nodes.h:1365
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition oomph_definitions.h:229
oomph::PointElement
Point element has just a single node and a single shape function which is identically equal to one.
Definition elements.h:3443
oomph::Shape
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition shape.h:76
oomph::SolidFiniteElement
SolidFiniteElement class.
Definition elements.h:3565
oomph::SolidNode
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition nodes.h:1686
oomph::SolidNode::nlagrangian
unsigned nlagrangian() const
Return number of lagrangian coordinates.
Definition nodes.h:1870
oomph::SolidPointElement
Solid point element.
Definition elements.h:4984
oomph::SolidTBubbleEnrichedElement< DIM, 3 >::SolidTBubbleEnrichedElement
SolidTBubbleEnrichedElement(const SolidTBubbleEnrichedElement &)=delete
Broken copy constructor.
oomph::SolidTBubbleEnrichedElement< DIM, 3 >::build_face_element
void build_face_element(const int &face_index, FaceElement *face_element_pt)
Build the lower-dimensional FaceElement Need to put in a final override here.
oomph::SolidTBubbleEnrichedElement< DIM, 3 >::~SolidTBubbleEnrichedElement
~SolidTBubbleEnrichedElement()
Broken assignment operator.
Definition Telements.h:3952
oomph::SolidTBubbleEnrichedElement
SolidTBubbleEnrichedElement elements are the enriched version of the SolidTElements....
Definition Telements.h:3926
oomph::SolidTElement< 1, NNODE_1D >::SolidTElement
SolidTElement(const SolidTElement &)=delete
Broken copy constructor.
oomph::SolidTElement< 1, NNODE_1D >::SolidTElement
SolidTElement()
Constructor.
Definition Telements.h:3741
oomph::SolidTElement< 2, NNODE_1D >::SolidTElement
SolidTElement()
Constructor.
Definition Telements.h:3802
oomph::SolidTElement< 2, NNODE_1D >::SolidTElement
SolidTElement(const SolidTElement &)=delete
Broken copy constructor.
oomph::SolidTElement< 3, NNODE_1D >::SolidTElement
SolidTElement()
Constructor.
Definition Telements.h:3862
oomph::SolidTElement< 3, NNODE_1D >::SolidTElement
SolidTElement(const SolidTElement &)=delete
Broken copy constructor.
oomph::SolidTElement
SolidTElement elements are triangular/tet elements whose derivatives also include those based upon th...
Definition Telements.h:3728
oomph::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Definition Tadvection_diffusion_reaction_elements.h:66
oomph::TBubbleEnrichedElementShape< 2, 3 >::n_enriched_nodes
unsigned n_enriched_nodes()
Return the number of nodes required for enrichement.
Definition Telements.h:921
oomph::TBubbleEnrichedElementShape< 2, 3 >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Second derivatives of shape functions for specific TBubbleElement<2,3>: d2psids(i,...
Definition Telements.h:1049
oomph::TBubbleEnrichedElementShape< 2, 3 >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TBubbleElement<2,3>
Definition Telements.h:1015
oomph::TBubbleEnrichedElementShape< 2, 3 >::shape
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TBubbleEnrichedElement<2,3>
Definition Telements.h:987
oomph::TBubbleEnrichedElementShape< 2, 3 >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition Telements.h:929
oomph::TBubbleEnrichedElementShape< 3, 3 >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Definition Telements.h:2320
oomph::TBubbleEnrichedElementShape< 3, 3 >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Second derivatives of shape functions for specific TElement<3,3>: d2psids(i,0) = d2psids(i,...
Definition Telements.h:2689
oomph::TBubbleEnrichedElementShape< 3, 3 >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<3,3>
Definition Telements.h:2494
oomph::TBubbleEnrichedElementShape< 3, 3 >::shape
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TBubbleEnrichedElement<3,3>
Definition Telements.h:2438
oomph::TBubbleEnrichedElementShape
A class for those member functions that must be fully specialised for Telements that are enriched by ...
Definition Telements.h:901
oomph::TBubbleEnrichedElement< DIM, 3 >::TBubbleEnrichedElement
TBubbleEnrichedElement(const TBubbleEnrichedElement &)=delete
Broken copy constructor.
oomph::TBubbleEnrichedElement< DIM, 3 >::build_face_element
void build_face_element(const int &face_index, FaceElement *face_element_pt)
Build the lower-dimensional FaceElement.
oomph::TBubbleEnrichedElement< DIM, 3 >::shape
void shape(const Vector< double > &s, Shape &psi) const
Calculate the geometric shape functions at local coordinate s.
Definition Telements.h:3655
oomph::TBubbleEnrichedElement< DIM, 3 >::~TBubbleEnrichedElement
~TBubbleEnrichedElement()
Broken assignment operator.
Definition Telements.h:3652
oomph::TBubbleEnrichedElement< DIM, 3 >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Compute the geometric shape functions and derivatives w.r.t. local coordinates at local coordinate s.
Definition Telements.h:3662
oomph::TBubbleEnrichedElement< DIM, 3 >::Default_enriched_integration_scheme
static TBubbleEnrichedGauss< DIM, 3 > Default_enriched_integration_scheme
Definition Telements.h:3627
oomph::TBubbleEnrichedElement< DIM, 3 >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition Telements.h:3685
oomph::TBubbleEnrichedElement< DIM, 3 >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Compute the geometric shape functions, derivatives and second derivatives w.r.t local coordinates at ...
Definition Telements.h:3675
oomph::TBubbleEnrichedElement
TElement class for which the shape functions have been enriched by a single bubble function of the ne...
Definition Telements.h:3570
oomph::TBubbleEnrichedElement::Central_node_on_face
const unsigned Central_node_on_face[3]
Definition Telements.cc:865
oomph::TBubbleEnrichedGauss
Define integration schemes that are required to exactly integrate the mass matrices of the bubble-enr...
Definition Telements.h:3581
oomph::TElementBase
Empty base class for Telements (created so that we can use dynamic_cast<>() to figure out if a an ele...
Definition Telements.h:1130
oomph::TElementBase::move_local_coord_back_into_element
void move_local_coord_back_into_element(Vector< double > &s) const
Adjust local coordinates so that they're located inside the element.
Definition Telements.h:1175
oomph::TElementBase::TElementBase
TElementBase()
Empty default constructor.
Definition Telements.h:1133
oomph::TElementBase::element_geometry
ElementGeometry::ElementGeometry element_geometry() const
Broken assignment operator.
Definition Telements.h:1142
oomph::TElementBase::TElementBase
TElementBase(const TElementBase &)=delete
Broken copy constructor.
oomph::TElementBase::local_coord_is_valid
bool local_coord_is_valid(const Vector< double > &s)
Check whether the local coordinates are valid or not.
Definition Telements.h:1148
oomph::TElementGeometricBase
Empty base class for Telements (created so that we can use dynamic_cast<>() to figure out if a an ele...
Definition Telements.h:1102
oomph::TElementGeometricBase::TElementGeometricBase
TElementGeometricBase(const TElementGeometricBase &)=delete
Broken copy constructor.
oomph::TElementGeometricBase::TElementGeometricBase
TElementGeometricBase()
Empty default constructor.
Definition Telements.h:1105
oomph::TElementShape< 1, 2 >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition Telements.h:242
oomph::TElementShape< 1, 2 >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<2,2>
Definition Telements.h:281
oomph::TElementShape< 1, 2 >::shape
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<1,2>
Definition Telements.h:271
oomph::TElementShape< 1, 2 >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Second derivatives of shape functions for specific TElement<1,2>: d2psids(i,0) = .
Definition Telements.h:295
oomph::TElementShape< 1, 3 >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition Telements.h:315
oomph::TElementShape< 1, 3 >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<1,3>
Definition Telements.h:359
oomph::TElementShape< 1, 3 >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Second derivatives of shape functions for specific TElement<1,3>: d2psids(i,0) = .
Definition Telements.h:374
oomph::TElementShape< 1, 3 >::shape
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<1,3>
Definition Telements.h:348
oomph::TElementShape< 1, 4 >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition Telements.h:395
oomph::TElementShape< 1, 4 >::shape
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<1,4>
Definition Telements.h:432
oomph::TElementShape< 1, 4 >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Second derivatives of shape functions for specific TElement<2,4>: d2psids(i,0) = .
Definition Telements.h:458
oomph::TElementShape< 1, 4 >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<1,4>
Definition Telements.h:443
oomph::TElementShape< 2, 2 >::shape
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<2,2>
Definition Telements.h:521
oomph::TElementShape< 2, 2 >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition Telements.h:484
oomph::TElementShape< 2, 2 >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Second derivatives of shape functions for specific TElement<2,2>: d2psids(i,0) = d2psids(i,...
Definition Telements.h:552
oomph::TElementShape< 2, 2 >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<2,2>
Definition Telements.h:532
oomph::TElementShape< 2, 3 >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition Telements.h:575
oomph::TElementShape< 2, 3 >::shape
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<2,3>
Definition Telements.h:627
oomph::TElementShape< 2, 3 >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<2,3>
Definition Telements.h:647
oomph::TElementShape< 2, 3 >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Second derivatives of shape functions for specific TElement<2,3>: d2psids(i,0) = d2psids(i,...
Definition Telements.h:673
oomph::TElementShape< 2, 4 >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<2,4>
Definition Telements.h:805
oomph::TElementShape< 2, 4 >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition Telements.h:714
oomph::TElementShape< 2, 4 >::shape
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<2,4>
Definition Telements.h:786
oomph::TElementShape< 2, 4 >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Second derivatives of shape functions for specific TElement<2,4>: d2psids(i,0) = d2psids(i,...
Definition Telements.h:846
oomph::TElementShape< 3, 2 >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Definition Telements.h:1932
oomph::TElementShape< 3, 2 >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<3,2>
Definition Telements.h:1990
oomph::TElementShape< 3, 2 >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Second derivatives of shape functions for specific TElement<3,2>: d2psids(i,0) = d2psids(i,...
Definition Telements.h:2023
oomph::TElementShape< 3, 2 >::shape
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<3,2>
Definition Telements.h:1978
oomph::TElementShape< 3, 3 >::shape
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<3,3>
Definition Telements.h:2133
oomph::TElementShape< 3, 3 >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Definition Telements.h:2051
oomph::TElementShape< 3, 3 >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Second derivatives of shape functions for specific TElement<3,3>: d2psids(i,0) = d2psids(i,...
Definition Telements.h:2211
oomph::TElementShape< 3, 3 >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<3,3>
Definition Telements.h:2152
oomph::TElementShape
A class for those member functions that must be fully specialised for the Telements....
Definition Telements.h:229
oomph::TElement< 1, NNODE_1D >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Compute the geometric shape functions and derivatives w.r.t. local coordinates at local coordinate s.
Definition Telements.h:1322
oomph::TElement< 1, NNODE_1D >::tecplot_zone_string
std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction)
Definition Telements.h:1472
oomph::TElement< 1, NNODE_1D >::nvertex_node
unsigned nvertex_node() const
Number of vertex nodes in the element: One more than spatial dimension.
Definition Telements.h:1282
oomph::TElement< 1, NNODE_1D >::nplot_points
unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction)
Definition Telements.h:1481
oomph::TElement< 1, NNODE_1D >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Compute the geometric shape functions, derivatives and second derivatives w.r.t local coordinates at ...
Definition Telements.h:1334
oomph::TElement< 1, NNODE_1D >::TElement
TElement()
Constructor.
Definition Telements.h:1232
oomph::TElement< 1, NNODE_1D >::get_s_plot
void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &use_equally_spaced_interior_sample_points=false) const
Get vector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition Telements.h:1446
oomph::TElement< 1, NNODE_1D >::write_paraview_type
void write_paraview_type(std::ofstream &file_out, const unsigned &nplot) const
Return the paraview element type. Needs to be implemented for each new geometric element type; see ht...
Definition Telements.h:1405
oomph::TElement< 1, NNODE_1D >::~TElement
~TElement()
Broken assignment operator.
Definition Telements.h:1271
oomph::TElement< 1, NNODE_1D >::s_max
double s_max() const
Max. value of local coordinate.
Definition Telements.h:1358
oomph::TElement< 1, NNODE_1D >::write_paraview_output_offset_information
void write_paraview_output_offset_information(std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const
Fill in the offset information for paraview plot. Needs to be implemented for each new geometric elem...
Definition Telements.h:1387
oomph::TElement< 1, NNODE_1D >::nsub_elements_paraview
unsigned nsub_elements_paraview(const unsigned &nplot) const
Return the number of local sub-elements for paraview plot with parameter nplot.
Definition Telements.h:1379
oomph::TElement< 1, NNODE_1D >::invert_jacobian_mapping
double invert_jacobian_mapping(const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Overload the template-free interface for the calculation of inverse jacobian matrix....
Definition Telements.h:1345
oomph::TElement< 1, NNODE_1D >::TElement
TElement(const TElement &)=delete
Broken copy constructor.
oomph::TElement< 1, NNODE_1D >::nnode_1d
unsigned nnode_1d() const
Number of nodes along each element edge.
Definition Telements.h:1274
oomph::TElement< 1, NNODE_1D >::shape
void shape(const Vector< double > &s, Shape &psi) const
Calculate the geometric shape functions at local coordinate s.
Definition Telements.h:1315
oomph::TElement< 1, NNODE_1D >::Default_integration_scheme
static TGauss< 1, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition Telements.h:1228
oomph::TElement< 1, NNODE_1D >::s_min
double s_min() const
Min. value of local coordinate.
Definition Telements.h:1352
oomph::TElement< 1, NNODE_1D >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition Telements.h:1364
oomph::TElement< 1, NNODE_1D >::vertex_node_pt
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Definition Telements.h:1288
oomph::TElement< 1, NNODE_1D >::write_paraview_offsets
void write_paraview_offsets(std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const
Return the offsets for the paraview sub-elements. Needs to be implemented for each new geometric elem...
Definition Telements.h:1418
oomph::TElement< 1, NNODE_1D >::nplot_points_paraview
unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot.
Definition Telements.h:1372
oomph::TElement< 2, NNODE_1D >::write_paraview_type
void write_paraview_type(std::ofstream &file_out, const unsigned &nplot) const
Return the paraview element type. Needs to be implemented for each new geometric element type; see ht...
Definition Telements.h:1737
oomph::TElement< 2, NNODE_1D >::nvertex_node
unsigned nvertex_node() const
Number of vertex nodes in the element: One more than spatial dimension.
Definition Telements.h:1592
oomph::TElement< 2, NNODE_1D >::s_min
double s_min() const
Min. value of local coordinate.
Definition Telements.h:1661
oomph::TElement< 2, NNODE_1D >::TElement
TElement()
Constructor.
Definition Telements.h:1520
oomph::TElement< 2, NNODE_1D >::write_tecplot_zone_footer
void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction)...
Definition Telements.h:1838
oomph::TElement< 2, NNODE_1D >::shape
void shape(const Vector< double > &s, Shape &psi) const
Calculate the geometric shape functions at local coordinate s.
Definition Telements.h:1624
oomph::TElement< 2, NNODE_1D >::write_paraview_offsets
void write_paraview_offsets(std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const
Return the offsets for the paraview sub-elements. Needs to be implemented for each new geometric elem...
Definition Telements.h:1752
oomph::TElement< 2, NNODE_1D >::tecplot_zone_string
std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction)
Definition Telements.h:1821
oomph::TElement< 2, NNODE_1D >::nplot_points_paraview
unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot.
Definition Telements.h:1681
oomph::TElement< 2, NNODE_1D >::s_max
double s_max() const
Max. value of local coordinate.
Definition Telements.h:1667
oomph::TElement< 2, NNODE_1D >::write_paraview_output_offset_information
void write_paraview_output_offset_information(std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const
Fill in the offset information for paraview plot. Needs to be implemented for each new geometric elem...
Definition Telements.h:1706
oomph::TElement< 2, NNODE_1D >::~TElement
~TElement()
Broken assignment operator.
Definition Telements.h:1582
oomph::TElement< 2, NNODE_1D >::invert_jacobian_mapping
double invert_jacobian_mapping(const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Overload the template-free interface for the calculation of inverse jacobian matrix....
Definition Telements.h:1654
oomph::TElement< 2, NNODE_1D >::get_s_plot
void get_s_plot(const unsigned &iplot, const unsigned &nplot, Vector< double > &s, const bool &use_equally_spaced_interior_sample_points=false) const
Get vector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition Telements.h:1781
oomph::TElement< 2, NNODE_1D >::vertex_node_pt
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Definition Telements.h:1605
oomph::TElement< 2, NNODE_1D >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Compute the geometric shape functions and derivatives w.r.t. local coordinates at local coordinate s.
Definition Telements.h:1631
oomph::TElement< 2, NNODE_1D >::Default_integration_scheme
static TGauss< 2, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition Telements.h:1515
oomph::TElement< 2, NNODE_1D >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Compute the geometric shape functions, derivatives and second derivatives w.r.t local coordinates at ...
Definition Telements.h:1643
oomph::TElement< 2, NNODE_1D >::nsub_elements_paraview
unsigned nsub_elements_paraview(const unsigned &nplot) const
Return the number of local sub-elements for paraview plot with parameter nplot.
Definition Telements.h:1693
oomph::TElement< 2, NNODE_1D >::nnode_1d
unsigned nnode_1d() const
Number of nodes along each element edge.
Definition Telements.h:1585
oomph::TElement< 2, NNODE_1D >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition Telements.h:1673
oomph::TElement< 2, NNODE_1D >::nplot_points
unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction)
Definition Telements.h:1899
oomph::TElement< 2, NNODE_1D >::write_tecplot_zone_footer
void write_tecplot_zone_footer(FILE *file_pt, const unsigned &nplot) const
Add tecplot zone "footer" to C-style output. (when plotting nplot points in each "coordinate directio...
Definition Telements.h:1867
oomph::TElement< 2, NNODE_1D >::get_bulk_node_number
unsigned get_bulk_node_number(const int &face_index, const unsigned &i) const
Public access function for Node_on_face.
Definition Telements.h:1598
oomph::TElement< 2, NNODE_1D >::TElement
TElement(const TElement &)=delete
Broken copy constructor.
oomph::TElement< 2, NNODE_1D >::TElement
TElement(const bool &allow_high_order)
Alternative constructor.
Definition Telements.h:1551
oomph::TElement< 3, NNODE_1D >::get_bulk_node_number
unsigned get_bulk_node_number(const int &face_index, const unsigned &i) const
Public access function for Node_on_face.
Definition Telements.h:3116
oomph::TElement< 3, NNODE_1D >::TElement
TElement(const TElement &)=delete
Broken copy constructor.
oomph::TElement< 3, NNODE_1D >::nnode_1d
unsigned nnode_1d() const
Number of nodes along each element edge.
Definition Telements.h:3102
oomph::TElement< 3, NNODE_1D >::s_max
double s_max() const
Max. value of local coordinate.
Definition Telements.h:3188
oomph::TElement< 3, NNODE_1D >::Default_integration_scheme
static TGauss< 3, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition Telements.h:3052
oomph::TElement< 3, NNODE_1D >::nsub_elements_paraview
unsigned nsub_elements_paraview(const unsigned &nplot) const
Return the number of local sub-elements for paraview plot with parameter nplot.
Definition Telements.h:3217
oomph::TElement< 3, NNODE_1D >::write_tecplot_zone_footer
void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction)...
Definition Telements.h:3416
oomph::TElement< 3, NNODE_1D >::d2shape_local
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Compute the geometric shape functions, derivatives and second derivatives w.r.t local coordinates at ...
Definition Telements.h:3164
oomph::TElement< 3, NNODE_1D >::~TElement
~TElement()
Broken assignment operator.
Definition Telements.h:3099
oomph::TElement< 3, NNODE_1D >::nplot_points_paraview
unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot.
Definition Telements.h:3202
oomph::TElement< 3, NNODE_1D >::tecplot_zone_string
std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction)
Definition Telements.h:3402
oomph::TElement< 3, NNODE_1D >::write_paraview_output_offset_information
void write_paraview_output_offset_information(std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const
Fill in the offset information for paraview plot. Needs to be implemented for each new geometric elem...
Definition Telements.h:3225
oomph::TElement< 3, NNODE_1D >::vertex_node_pt
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Definition Telements.h:3123
oomph::TElement< 3, NNODE_1D >::write_tecplot_zone_footer
void write_tecplot_zone_footer(FILE *file_pt, const unsigned &nplot) const
Add tecplot zone "footer" to C-style output. (when plotting nplot points in each "coordinate directio...
Definition Telements.h:3494
oomph::TElement< 3, NNODE_1D >::shape
void shape(const Vector< double > &s, Shape &psi) const
Calculate the geometric shape functions at local coordinate s.
Definition Telements.h:3142
oomph::TElement< 3, NNODE_1D >::nvertex_node
unsigned nvertex_node() const
Number of vertex nodes in the element: One more than spatial dimension.
Definition Telements.h:3110
oomph::TElement< 3, NNODE_1D >::nplot_points
unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction)
Definition Telements.h:3529
oomph::TElement< 3, NNODE_1D >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition Telements.h:3194
oomph::TElement< 3, NNODE_1D >::write_paraview_type
void write_paraview_type(std::ofstream &file_out, const unsigned &nplot) const
Return the paraview element type. Needs to be implemented for each new geometric element type; see ht...
Definition Telements.h:3315
oomph::TElement< 3, NNODE_1D >::s_min
double s_min() const
Min. value of local coordinate.
Definition Telements.h:3182
oomph::TElement< 3, NNODE_1D >::write_paraview_offsets
void write_paraview_offsets(std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const
Return the offsets for the paraview sub-elements. Needs to be implemented for each new geometric elem...
Definition Telements.h:3328
oomph::TElement< 3, NNODE_1D >::TElement
TElement()
Constructor.
Definition Telements.h:3056
oomph::TElement< 3, NNODE_1D >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Compute the geometric shape functions and derivatives w.r.t. local coordinates at local coordinate s.
Definition Telements.h:3149
oomph::TElement< 3, NNODE_1D >::invert_jacobian_mapping
double invert_jacobian_mapping(const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Overload the template-free interface for the calculation of inverse jacobian matrix....
Definition Telements.h:3175
oomph::TElement< 3, NNODE_1D >::get_s_plot
void get_s_plot(const unsigned &iplot, const unsigned &nplot, Vector< double > &s, const bool &use_equally_spaced_interior_sample_points=false) const
Get vector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition Telements.h:3354
oomph::TElement
General TElement class.
Definition Telements.h:1208
oomph::TElement::Node_on_face
const unsigned Node_on_face[3][2]
Assign the nodal translation schemes.
Definition Telements.cc:272
oomph::TFace
Triangular Face class.
Definition Telements.h:56
oomph::TFace::is_on_boundary
bool is_on_boundary() const
Test whether the face lies on a boundary. Relatively simple test, based on all vertices lying on (som...
Definition Telements.h:158
oomph::TFace::Node2_pt
Node * Node2_pt
Second vertex node.
Definition Telements.h:208
oomph::TFace::node3_pt
Node * node3_pt() const
Access to the third vertex node.
Definition Telements.h:100
oomph::TFace::Node3_pt
Node * Node3_pt
Third vertex node.
Definition Telements.h:211
oomph::TFace::get_boundaries_pt
void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Access to pointer to set of mesh boundaries that this face occupies; NULL if the node is not on any b...
Definition Telements.h:178
oomph::TFace::node1_pt
Node * node1_pt() const
Access to the first vertex node.
Definition Telements.h:88
oomph::TFace::TFace
TFace(Node *node1_pt, Node *node2_pt, Node *node3_pt)
Constructor: Pass in the three vertex nodes.
Definition Telements.h:59
oomph::TFace::Node1_pt
Node * Node1_pt
First vertex node.
Definition Telements.h:205
oomph::TFace::operator<
bool operator<(const TFace &other) const
Less-than operator.
Definition Telements.h:121
oomph::TFace::is_boundary_face
bool is_boundary_face() const
Test whether the face is a boundary face, i.e. does it connnect three boundary nodes?
Definition Telements.h:167
oomph::TFace::node2_pt
Node * node2_pt() const
Access to the second vertex node.
Definition Telements.h:94
oomph::TFace::operator==
bool operator==(const TFace &other) const
Comparison operator.
Definition Telements.h:106
oomph::TGauss
Class for Gaussian integration rules for triangles/tets.
Definition integral.h:626
oomph::TSolidElementBase
Base class for Solid Telements.
Definition Telements.h:3702
oomph::TSolidElementBase::TSolidElementBase
TSolidElementBase()
Constructor: Empty.
Definition Telements.h:3705
oomph::TSolidElementBase::TSolidElementBase
TSolidElementBase(const TSolidElementBase &)=delete
Broken copy constructor.
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition advection_diffusion_elements.cc:30