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