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-2022 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Header functions for classes that define Telements
27 #ifndef OOMPH_TELEMENT_HEADER
28 #define OOMPH_TELEMENT_HEADER
29 
30 // Config header generated by autoconfig
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 
45 namespace oomph
46 {
50 
51 
52  //===================================================
54  //===================================================
55  class TFace
56  {
57  public:
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 
88  Node* node1_pt() const
89  {
90  return Node1_pt;
91  }
92 
94  Node* node2_pt() const
95  {
96  return Node2_pt;
97  }
98 
100  Node* node3_pt() const
101  {
102  return Node3_pt;
103  }
104 
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 
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 
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 
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 
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:
205  Node* Node1_pt;
206 
208  Node* Node2_pt;
209 
211  Node* Node3_pt;
212  };
213 
214 
218 
219 
220  //========================================================================
226  //========================================================================
227  template<unsigned DIM, unsigned NNODE_1D>
228  class TElementShape
229  {
230  };
231 
235  template<>
236  class TElementShape<1, 2>
237  {
238  public:
239  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //========================================================================
898  //========================================================================
899  template<unsigned DIM, unsigned NNODE_1D>
900  class TBubbleEnrichedElementShape
901  {
902  };
903 
904 
908 
909  //===============================================================
913  //=============================================================
914  template<>
915  class TBubbleEnrichedElementShape<2, 3>
916  {
917  public:
918  //=====================================================================
920  //====================================================================
921  unsigned n_enriched_nodes()
922  {
923  return 1;
924  }
925 
926  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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 
1095 
1096  //========================================================================
1100  //========================================================================
1101  class TElementGeometricBase : public virtual FiniteElement
1102  {
1103  public:
1105  TElementGeometricBase() {}
1106 
1108  TElementGeometricBase(const TElementGeometricBase&) = delete;
1109 
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 
1123 
1124  //========================================================================
1128  //========================================================================
1129  class TElementBase : public virtual TElementGeometricBase
1130  {
1131  public:
1133  TElementBase() {}
1134 
1136  TElementBase(const TElementBase&) = delete;
1137 
1139  /*void operator=(const TElementBase&) = delete;*/
1140 
1142  ElementGeometry::ElementGeometry element_geometry() const
1143  {
1144  return ElementGeometry::T;
1145  }
1146 
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 
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  //=======================================================================
1205  //=======================================================================
1206  template<unsigned DIM, unsigned NNODE_1D>
1207  class TElement
1208  {
1209  };
1210 
1211 
1212  //=======================================================================
1217  //=======================================================================
1218  template<unsigned NNODE_1D>
1219  class TElement<1, NNODE_1D> : public virtual TElementBase,
1220  public TElementShape<1, NNODE_1D>
1221  {
1222  private:
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:
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 
1264  TElement(const TElement&) = delete;
1265 
1267  /*void operator=(const TElement&) = delete;*/
1268 
1269 
1271  ~TElement() {}
1272 
1274  unsigned nnode_1d() const
1275  {
1276  return NNODE_1D;
1277  }
1278 
1279 
1282  unsigned nvertex_node() const
1283  {
1284  return 2;
1285  }
1286 
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 
1315  inline void shape(const Vector<double>& s, Shape& psi) const
1316  {
1317  TElementShape<1, NNODE_1D>::shape(s, psi);
1318  }
1319 
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 
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 
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 
1352  double s_min() const
1353  {
1354  return 0.0;
1355  }
1356 
1358  double s_max() const
1359  {
1360  return 1.0;
1361  }
1362 
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 
1372  unsigned nplot_points_paraview(const unsigned& nplot) const
1373  {
1374  return nplot;
1375  }
1376 
1379  unsigned nsub_elements_paraview(const unsigned& nplot) const
1380  {
1381  return (nplot - 1);
1382  }
1383 
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 
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 
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 
1433  void output(std::ostream& output);
1434 
1436  void output(std::ostream& outfile, const unsigned& nplot);
1437 
1439  void output(FILE* file_pt);
1440 
1442  void output(FILE* file_pt, const unsigned& n_plot);
1443 
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 
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 
1481  unsigned nplot_points(const unsigned& nplot) const
1482  {
1483  return nplot;
1484  }
1485 
1491  void build_face_element(const int& face_index,
1492  FaceElement* face_element_pt);
1493  };
1494 
1495 
1496  //=======================================================================
1501  //=======================================================================
1502  template<unsigned NNODE_1D>
1503  class TElement<2, NNODE_1D> : public virtual TElementBase,
1504  public TElementShape<2, NNODE_1D>
1505  {
1506  private:
1508  static const unsigned Node_on_face[3][NNODE_1D];
1509 
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:
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 
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 
1575  TElement(const TElement&) = delete;
1576 
1578  /*void operator=(const TElement&) = delete;*/
1579 
1580 
1582  ~TElement() {}
1583 
1585  unsigned nnode_1d() const
1586  {
1587  return NNODE_1D;
1588  }
1589 
1592  unsigned nvertex_node() const
1593  {
1594  return 3;
1595  }
1596 
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 
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 
1624  inline void shape(const Vector<double>& s, Shape& psi) const
1625  {
1626  TElementShape<2, NNODE_1D>::shape(s, psi);
1627  }
1628 
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 
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 
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 
1661  double s_min() const
1662  {
1663  return 0.0;
1664  }
1665 
1667  double s_max() const
1668  {
1669  return 1.0;
1670  }
1671 
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 
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 
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 
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 
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 
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 
1768  void output(std::ostream& output);
1769 
1771  void output(std::ostream& outfile, const unsigned& nplot);
1772 
1774  void output(FILE* file_pt);
1775 
1777  void output(FILE* file_pt, const unsigned& n_plot);
1778 
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 
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 
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 
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 
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 
1915  void build_face_element(const int& face_index,
1916  FaceElement* face_element_pt);
1917  };
1918 
1919 
1923 
1924 
1925  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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 
2299 
2300  //====================================================================
2306  //=========================================================================
2307 
2308  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
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  //=======================================================================
3038  //=======================================================================
3039  template<unsigned NNODE_1D>
3040  class TElement<3, NNODE_1D> : public virtual TElementBase,
3041  public TElementShape<3, NNODE_1D>
3042  {
3043  private:
3045  static const unsigned Node_on_face[4][(NNODE_1D * (NNODE_1D + 1)) / 2];
3046 
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:
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 
3092  TElement(const TElement&) = delete;
3093 
3095  /*void operator=(const TElement&) = delete;*/
3096 
3097 
3099  ~TElement() {}
3100 
3102  unsigned nnode_1d() const
3103  {
3104  return NNODE_1D;
3105  }
3106 
3107 
3110  unsigned nvertex_node() const
3111  {
3112  return 4;
3113  }
3114 
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 
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 
3142  inline void shape(const Vector<double>& s, Shape& psi) const
3143  {
3144  TElementShape<3, NNODE_1D>::shape(s, psi);
3145  }
3146 
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 
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 
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 
3182  double s_min() const
3183  {
3184  return 0.0;
3185  }
3186 
3188  double s_max() const
3189  {
3190  return 1.0;
3191  }
3192 
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 
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 
3217  unsigned nsub_elements_paraview(const unsigned& nplot) const
3218  {
3219  return (nplot - 1) * (nplot - 1) * (nplot - 1);
3220  }
3221 
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 
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 
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 
3341  void output(std::ostream& output);
3342 
3344  void output(std::ostream& outfile, const unsigned& nplot);
3345 
3347  void output(FILE* file_pt);
3348 
3350  void output(FILE* file_pt, const unsigned& n_plot);
3351 
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 
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 
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 
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 
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 
3552  void build_face_element(const int& face_index,
3553  FaceElement* face_element_pt);
3554  };
3555 
3556 
3560 
3561 
3562  //=======================================================================
3567  //=======================================================================
3568  template<unsigned DIM, unsigned NNODE_1D>
3569  class TBubbleEnrichedElement
3570  {
3571  };
3572 
3573  //=====================================================================
3578  //====================================================================
3579  template<unsigned DIM, unsigned NPTS_1D>
3580  class TBubbleEnrichedGauss
3581  {
3582  };
3583 
3584  //====================================================================
3588  //======================================================================
3589  template<>
3590  class TBubbleEnrichedGauss<2, 3> : public TGauss<2, 4>
3591  {
3592  public:
3593  TBubbleEnrichedGauss() : TGauss<2, 4>() {}
3594  };
3595 
3596  //====================================================================
3600  //======================================================================
3601  template<>
3602  class TBubbleEnrichedGauss<3, 3> : public TGauss<3, 5>
3603  {
3604  public:
3605  TBubbleEnrichedGauss() : TGauss<3, 5>() {}
3606  };
3607 
3608 
3609  //=======================================================================
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:
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 
3646  TBubbleEnrichedElement(const TBubbleEnrichedElement&) = delete;
3647 
3649  /*void operator=(const TBubbleEnrichedElement&) = delete;*/
3650 
3652  ~TBubbleEnrichedElement() {}
3653 
3655  inline void shape(const Vector<double>& s, Shape& psi) const
3656  {
3657  TBubbleEnrichedElementShape<DIM, 3>::shape(s, psi);
3658  }
3659 
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 
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 
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 
3692  void build_face_element(const int& face_index,
3693  FaceElement* face_element_pt);
3694  };
3695 
3696 
3697  //========================================================================
3699  //========================================================================
3700  class TSolidElementBase : public virtual TElementBase,
3701  public virtual SolidFiniteElement
3702  {
3703  public:
3705  TSolidElementBase(){};
3706 
3708  TSolidElementBase(const TSolidElementBase&) = delete;
3709 
3711  /*void operator=(const TSolidElementBase&) = delete;*/
3712  };
3713 
3714 
3718 
3719 
3720  //=======================================================================
3725  //=======================================================================
3726  template<unsigned DIM, unsigned NNODE_1D>
3727  class SolidTElement
3728  {
3729  };
3730 
3731 
3732  //=======================================================================
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:
3741  SolidTElement() : TElement<1, NNODE_1D>(), SolidFiniteElement()
3742  {
3743  // Set the Lagrangian dimension of the element
3744  set_lagrangian_dimension(1);
3745  }
3746 
3748  SolidTElement(const SolidTElement&) = delete;
3749 
3751  /*void operator=(const SolidTElement&) = delete;*/
3752 
3758  inline void build_face_element(const int& face_index,
3759  FaceElement* face_element_pt);
3760  };
3761 
3762 
3765  // SolidTElements
3768 
3769 
3771  // 1D SolidTElements
3773 
3774 
3775  //===========================================================
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  //=======================================================================
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:
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 
3813  SolidTElement(const SolidTElement&) = delete;
3814 
3816  /*void operator=(const SolidTElement&) = delete;*/
3817 
3823  inline void build_face_element(const int& face_index,
3824  FaceElement* face_element_pt);
3825  };
3826 
3827 
3830  // 2D SolidTElements
3833 
3834 
3835  //===========================================================
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  //=======================================================================
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:
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 
3873  SolidTElement(const SolidTElement&) = delete;
3874 
3876  /*void operator=(const SolidTElement&) = delete;*/
3877 
3878 
3886  inline void build_face_element(const int& face_index,
3887  FaceElement* face_element_pt);
3888  };
3889 
3890 
3892  // 3D SolidTElements
3894 
3895 
3896  //===========================================================
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 
3917 
3918  //=======================================================================
3923  //=======================================================================
3924  template<unsigned DIM, unsigned NNODE_1D>
3925  class SolidTBubbleEnrichedElement
3926  {
3927  };
3928 
3929  //===================================================================
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:
3940  SolidTBubbleEnrichedElement()
3941  : SolidTElement<DIM, 3>(), TBubbleEnrichedElement<DIM, 3>()
3942  {
3943  }
3944 
3946  SolidTBubbleEnrichedElement(const SolidTBubbleEnrichedElement&) = delete;
3947 
3949  /*void operator=(const SolidTBubbleEnrichedElement&) = delete;*/
3950 
3952  ~SolidTBubbleEnrichedElement() {}
3953 
3956  void build_face_element(const int& face_index,
3957  FaceElement* face_element_pt);
3958  };
3959 
3960 
3961  //=======================================================================
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:
3974  FaceGeometry() : TElement<DIM - 1, NNODE_1D>() {}
3975  };
3976 
3977 
3978  //=======================================================================
3980  //=======================================================================
3981  template<unsigned NNODE_1D>
3982  class FaceGeometry<TElement<1, NNODE_1D>> : public virtual PointElement
3983  {
3984  public:
3987  FaceGeometry() : PointElement() {}
3988  };
3989 
3990 
3994 
3995 
3996  //=======================================================================
4002  //=======================================================================
4003  template<unsigned NNODE_1D>
4004  class FaceGeometry<TBubbleEnrichedElement<2, NNODE_1D>>
4005  : public virtual TElement<1, NNODE_1D>
4006  {
4007  public:
4010  FaceGeometry() : TElement<1, NNODE_1D>() {}
4011  };
4012 
4013 
4014  //=======================================================================
4020  //=======================================================================
4021  template<unsigned NNODE_1D>
4022  class FaceGeometry<TBubbleEnrichedElement<3, NNODE_1D>>
4023  : public virtual TBubbleEnrichedElement<2, NNODE_1D>
4024  {
4025  public:
4028  FaceGeometry() : TBubbleEnrichedElement<2, NNODE_1D>() {}
4029  };
4030 
4031 
4035 
4036 
4037  //=======================================================================
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:
4050  FaceGeometry() : SolidTElement<DIM - 1, NNODE_1D>() {}
4051  };
4052 
4053 
4054  //=======================================================================
4056  //=======================================================================
4057  template<unsigned NNODE_1D>
4058  class FaceGeometry<SolidTElement<1, NNODE_1D>>
4059  : public virtual SolidPointElement
4060  {
4061  public:
4064  FaceGeometry() : SolidPointElement() {}
4065  };
4066 
4067 
4071 
4072 
4073  //=======================================================================
4079  //=======================================================================
4080  template<unsigned NNODE_1D>
4081  class FaceGeometry<SolidTBubbleEnrichedElement<2, NNODE_1D>>
4082  : public virtual SolidTElement<1, NNODE_1D>
4083  {
4084  public:
4087  FaceGeometry() : SolidTElement<1, NNODE_1D>() {}
4088  };
4089 
4090 
4091  //=======================================================================
4097  //=======================================================================
4098  template<unsigned NNODE_1D>
4099  class FaceGeometry<SolidTBubbleEnrichedElement<3, NNODE_1D>>
4100  : public virtual SolidTBubbleEnrichedElement<2, NNODE_1D>
4101  {
4102  public:
4105  FaceGeometry() : SolidTBubbleEnrichedElement<2, NNODE_1D>() {}
4106  };
4107 
4108 } // namespace oomph
4109 
4110 
4111 #endif
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
oomph::BoundaryNodeBase
Definition: nodes.h:1996
oomph::DShape
Definition: shape.h:278
oomph::FaceElement
Definition: elements.h:4338
oomph::FaceGeometry< SolidTBubbleEnrichedElement< 2, NNODE_1D > >::FaceGeometry
FaceGeometry()
Definition: Telements.h:4087
oomph::FaceGeometry< SolidTBubbleEnrichedElement< 3, NNODE_1D > >::FaceGeometry
FaceGeometry()
Definition: Telements.h:4105
oomph::FaceGeometry< SolidTElement< 1, NNODE_1D > >::FaceGeometry
FaceGeometry()
Definition: Telements.h:4064
oomph::FaceGeometry< SolidTElement< DIM, NNODE_1D > >::FaceGeometry
FaceGeometry()
Definition: Telements.h:4050
oomph::FaceGeometry< TBubbleEnrichedElement< 2, NNODE_1D > >::FaceGeometry
FaceGeometry()
Definition: Telements.h:4010
oomph::FaceGeometry< TBubbleEnrichedElement< 3, NNODE_1D > >::FaceGeometry
FaceGeometry()
Definition: Telements.h:4028
oomph::FaceGeometry< TElement< 1, NNODE_1D > >::FaceGeometry
FaceGeometry()
Definition: Telements.h:3987
oomph::FaceGeometry< TElement< DIM, NNODE_1D > >::FaceGeometry
FaceGeometry()
Definition: Telements.h:3974
oomph::FaceGeometry
Definition: elements.h:4998
oomph::FiniteElement
Definition: elements.h:1313
oomph::FiniteElement::dim
unsigned dim() const
Definition: elements.h:2611
oomph::Node
Definition: nodes.h:906
oomph::Node::is_on_boundary
virtual bool is_on_boundary() const
Definition: nodes.h:1373
oomph::Node::get_boundaries_pt
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Definition: nodes.h:1365
oomph::OomphLibError
Definition: oomph_definitions.h:222
oomph::PointElement
Definition: elements.h:3439
oomph::Shape
Definition: shape.h:76
oomph::SolidFiniteElement
Definition: elements.h:3561
oomph::SolidNode
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:4980
oomph::SolidTBubbleEnrichedElement< DIM, 3 >::SolidTBubbleEnrichedElement
SolidTBubbleEnrichedElement(const SolidTBubbleEnrichedElement &)=delete
Broken copy constructor.
oomph::SolidTBubbleEnrichedElement< DIM, 3 >::SolidTBubbleEnrichedElement
SolidTBubbleEnrichedElement()
Constructor.
Definition: Telements.h:3940
oomph::SolidTBubbleEnrichedElement< DIM, 3 >::build_face_element
void build_face_element(const int &face_index, FaceElement *face_element_pt)
oomph::SolidTBubbleEnrichedElement< DIM, 3 >::~SolidTBubbleEnrichedElement
~SolidTBubbleEnrichedElement()
Broken assignment operator.
Definition: Telements.h:3952
oomph::SolidTBubbleEnrichedElement
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
Definition: Telements.h:3728
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
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
Definition: Telements.h:2689
oomph::TBubbleEnrichedElementShape< 3, 3 >::n_enriched_nodes
unsigned n_enriched_nodes()
Definition: Telements.h:2315
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
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
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 >::TBubbleEnrichedElement
TBubbleEnrichedElement()
Constructor.
Definition: Telements.h:3634
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
Definition: Telements.h:3675
oomph::TBubbleEnrichedElement
Definition: Telements.h:3570
oomph::TBubbleEnrichedElement::Central_node_on_face
const unsigned Central_node_on_face[3]
Definition: Telements.cc:865
oomph::TBubbleEnrichedGauss< 2, 3 >::TBubbleEnrichedGauss
TBubbleEnrichedGauss()
Definition: Telements.h:3593
oomph::TBubbleEnrichedGauss< 3, 3 >::TBubbleEnrichedGauss
TBubbleEnrichedGauss()
Definition: Telements.h:3605
oomph::TBubbleEnrichedGauss
Definition: Telements.h:3581
oomph::TElementBase
Definition: Telements.h:1130
oomph::TElementBase::move_local_coord_back_into_element
void move_local_coord_back_into_element(Vector< double > &s) const
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
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
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
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
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
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
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
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
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
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
Definition: Telements.h:229
oomph::TElement< 1, NNODE_1D >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Definition: Telements.h:1322
oomph::TElement< 1, NNODE_1D >::tecplot_zone_string
std::string tecplot_zone_string(const unsigned &nplot) const
Definition: Telements.h:1472
oomph::TElement< 1, NNODE_1D >::nvertex_node
unsigned nvertex_node() const
Definition: Telements.h:1282
oomph::TElement< 1, NNODE_1D >::nplot_points
unsigned nplot_points(const unsigned &nplot) const
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
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
Definition: Telements.h:1446
oomph::TElement< 1, NNODE_1D >::write_paraview_type
void write_paraview_type(std::ofstream &file_out, const unsigned &nplot) const
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
Definition: Telements.h:1387
oomph::TElement< 1, NNODE_1D >::nsub_elements_paraview
unsigned nsub_elements_paraview(const unsigned &nplot) const
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
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
Assign the static integral.
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
Definition: Telements.h:1418
oomph::TElement< 1, NNODE_1D >::nplot_points_paraview
unsigned nplot_points_paraview(const unsigned &nplot) const
Definition: Telements.h:1372
oomph::TElement< 2, NNODE_1D >
Definition: Telements.h:1505
oomph::TElement< 2, NNODE_1D >::write_paraview_type
void write_paraview_type(std::ofstream &file_out, const unsigned &nplot) const
Definition: Telements.h:1737
oomph::TElement< 2, NNODE_1D >::nvertex_node
unsigned nvertex_node() const
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
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
Definition: Telements.h:1752
oomph::TElement< 2, NNODE_1D >::tecplot_zone_string
std::string tecplot_zone_string(const unsigned &nplot) const
Definition: Telements.h:1821
oomph::TElement< 2, NNODE_1D >::nplot_points_paraview
unsigned nplot_points_paraview(const unsigned &nplot) const
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
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
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
Definition: Telements.h:1781
oomph::TElement< 2, NNODE_1D >::dshape_local
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Definition: Telements.h:1631
oomph::TElement< 2, NNODE_1D >::Default_integration_scheme
static TGauss< 2, NNODE_1D > Default_integration_scheme
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
Definition: Telements.h:1643
oomph::TElement< 2, NNODE_1D >::nsub_elements_paraview
unsigned nsub_elements_paraview(const unsigned &nplot) const
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
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
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< 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< 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
Definition: Telements.h:3052
oomph::TElement< 3, NNODE_1D >::nsub_elements_paraview
unsigned nsub_elements_paraview(const unsigned &nplot) const
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
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
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
Definition: Telements.h:3202
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 >::tecplot_zone_string
std::string tecplot_zone_string(const unsigned &nplot) const
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
Definition: Telements.h:3225
oomph::TElement< 3, NNODE_1D >::write_tecplot_zone_footer
void write_tecplot_zone_footer(FILE *file_pt, const unsigned &nplot) const
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
Definition: Telements.h:3110
oomph::TElement< 3, NNODE_1D >::nplot_points
unsigned nplot_points(const unsigned &nplot) const
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
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
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
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
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
Definition: Telements.h:3354
oomph::TElement
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
Definition: Telements.h:158
oomph::TFace::Node2_pt
Node * Node2_pt
Second vertex node.
Definition: Telements.h:208
oomph::TFace::node2_pt
Node * node2_pt() const
Access to the second vertex node.
Definition: Telements.h:94
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)
Definition: Telements.h:178
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() const
Access to the first vertex node.
Definition: Telements.h:88
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
Definition: Telements.h:167
oomph::TFace::node3_pt
Node * node3_pt() const
Access to the third vertex node.
Definition: Telements.h:100
oomph::TFace::operator==
bool operator==(const TFace &other) const
Comparison operator.
Definition: Telements.h:106
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.
s
RealScalar s
Definition: level1_cplx_impl.h:130
k
char char char int int * k
Definition: level2_impl.h:374
DIM
#define DIM
Definition: linearised_navier_stokes_elements.h:44
PlanarWave::plot
void plot()
Plot.
Definition: sphere_scattering.cc:180
oomph::ElementGeometry::ElementGeometry
ElementGeometry
Definition: elements.h:1240
oomph::ElementGeometry::T
@ T
Definition: elements.h:1242
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
oomph::OneDimLagrange::shape
void shape(const double &s, double *Psi)
Definition: shape.h:564
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition: AnisotropicHookean.h:10
oomph::output
void output(std::ostream &outfile)
Output with default number of plot points.
Definition: gen_axisym_advection_diffusion_elements.h:161
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
sub
EIGEN_DONT_INLINE T sub(T a, T b)
Definition: svd_common.h:238
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2