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: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 2016) $
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
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() &&
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
208 
209  /// Second vertex node
211 
212  /// Third vertex node
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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>
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<>
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 //=======================================================================
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 //=======================================================================
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
1122 
1123  /// Broken copy constructor
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
1159 
1160  /// Broken copy constructor
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!
1174  {
1175  return ElementGeometry::T;
1176  }
1177 
1178  ///Check whether the local coordinates are valid or not
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 Check whether the local coordinates are valid or not and
1207  /// allow for rounding tolerance. If the local coordinate specifed
1208  /// by s is "slightly" outside the element (as specified by
1209  /// rounding_tolerance) we move the point back into the element.
1211  const double &rounding_tolerance)
1212  {
1213 
1214  // Check coordinates
1215  unsigned ncoord=dim();
1216  double sum=0.0;
1217  for (unsigned i=0;i<ncoord;i++)
1218  {
1219  //We allow a slight rounding tolerance
1220  //Each coordinate must be positive
1221  if (s[i]<0.0)
1222  {
1223  if (fabs(s[i])<rounding_tolerance)
1224  {
1225  s[i]=0.0;
1226  }
1227  else
1228  {
1229  return false;
1230  }
1231  }
1232  sum+=s[i];
1233  }
1234 
1235 
1236  //Sum must be less than 1
1237  double excess=sum-1.0;
1238  if (sum<=1.0)
1239  {
1240  return true;
1241  }
1242  else if ( (excess >= 0.0) && ( excess < rounding_tolerance) )
1243  {
1244  //In this case, we subtract 10% more than excess to refit
1245  double sub = (1.1*excess)/3.0 ;
1246  for (unsigned i=0;i<ncoord;i++)
1247  {
1248  s[i]-=sub;
1249  }
1250  return true;
1251  }
1252 
1253  // We're genuinely outside the element -- bail out.
1254  return false;
1255 
1256  }
1257 
1258 };
1259 
1260 //=======================================================================
1261 ///General TElement class
1262 ///
1263 /// Empty, just establishes the template parameters
1264 //=======================================================================
1265 template<unsigned DIM, unsigned NNODE_1D>
1266 class TElement
1267 {
1268 };
1269 
1270 
1271 //=======================================================================
1272 /// General TElement class specialised to one spatial dimensions
1273 /// Ordering of nodes is 0 at local coordinate s[0] = 0, 1 at local
1274 /// coordinate s[0] = 1 and then filling in the intermediate values
1275 /// from s[0]=0 to 1.
1276 //=======================================================================
1277 template<unsigned NNODE_1D>
1278 class TElement<1,NNODE_1D> : public virtual TElementBase,
1279  public TElementShape<1,NNODE_1D>
1280 {
1281  private:
1282 
1283  /// \short Default integration rule: Gaussian integration of same 'order' as
1284  /// the element
1285  //This is sort of optimal, because it means that the integration is exact
1286  //for the shape functions. Can overwrite this in specific element defintion.
1288 
1289 public:
1290 
1291  /// Constructor
1293  {
1294  // Number of nodes
1295  switch (NNODE_1D)
1296  {
1297  case 2:
1298  case 3:
1299  case 4:
1300  break;
1301 
1302  default:
1303  std::string error_message =
1304  "One-dimensional TElements are currently only implemented for\n";
1305  error_message +=
1306  "three and six nodes, i.e. NNODE_1D=2 , 3 , 4\n";
1307 
1308  throw OomphLibError(error_message,
1309  OOMPH_CURRENT_FUNCTION,
1310  OOMPH_EXCEPTION_LOCATION);
1311  }
1312 
1313  // Set the number of nodes
1314  unsigned n_node = NNODE_1D;
1315  this->set_n_node(n_node);
1316 
1317  // Set the elemental and nodal dimension
1318  set_dimension(1);
1319 
1320  //Assign default (full) spatial integration scheme
1321  set_integration_scheme(&Default_integration_scheme);
1322  }
1323 
1324 
1325  /// Broken copy constructor
1327  {
1328  BrokenCopy::broken_copy("TElement");
1329  }
1330 
1331  /// Broken assignment operator
1332  /*void operator=(const TElement&)
1333  {
1334  BrokenCopy::broken_assign("TElement");
1335  }*/
1336 
1337 
1338  /// Destructor
1340 
1341  /// Number of nodes along each element edge
1342  unsigned nnode_1d() const {return NNODE_1D;}
1343 
1344 
1345  /// \short Number of vertex nodes in the element: One more
1346  /// than spatial dimension
1347  unsigned nvertex_node() const {return 2;}
1348 
1349  /// \short Pointer to the j-th vertex node in the element
1350  Node* vertex_node_pt(const unsigned& j) const
1351  {
1352  switch (j)
1353  {
1354  case 0:
1355 
1356  return node_pt(0);
1357  break;
1358 
1359  case 1:
1360 
1361  return node_pt(NNODE_1D-1);
1362  break;
1363 
1364  default:
1365 
1366  std::ostringstream error_message;
1367  error_message
1368  << "Element only has two vertex nodes; called with node number "
1369  << j << std::endl;
1370  throw OomphLibError(error_message.str(),
1371  OOMPH_CURRENT_FUNCTION,
1372  OOMPH_EXCEPTION_LOCATION);
1373  }
1374  }
1375 
1376  /// Calculate the geometric shape functions at local coordinate s
1377  inline void shape(const Vector<double> &s, Shape &psi) const
1379 
1380  /// \short Compute the geometric shape functions and
1381  /// derivatives w.r.t. local coordinates at local coordinate s
1382  inline void dshape_local(const Vector<double> &s, Shape &psi,
1383  DShape &dpsids) const
1385 
1386  /// \short Compute the geometric shape functions, derivatives and
1387  /// second derivatives w.r.t local coordinates at local coordinate s
1388  /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
1389  /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
1390  /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
1391  inline void d2shape_local(const Vector<double> &s, Shape &psi,
1392  DShape &dpsids, DShape &d2psids) const
1393  {TElementShape<1,NNODE_1D>::d2shape_local(s,psi,dpsids,d2psids);}
1394 
1395  /// \short Overload the template-free interface for the calculation of
1396  /// inverse jacobian matrix. This is a one dimensional element, so use
1397  /// the 1D version.
1399  DenseMatrix<double>&inverse_jacobian) const
1400  {return FiniteElement::invert_jacobian<1>(jacobian,inverse_jacobian);}
1401 
1402  /// Min. value of local coordinate
1403  double s_min() const {return 0.0;}
1404 
1405  /// Max. value of local coordinate
1406  double s_max() const {return 1.0;}
1407 
1408  /// Return local coordinates of node j
1409  inline void local_coordinate_of_node(const unsigned& j,Vector<double>& s) const
1411 
1412  /// \short Return the number of actual plot points for paraview
1413  /// plot with parameter nplot.
1414  unsigned nplot_points_paraview(const unsigned& nplot) const
1415  {
1416  return nplot;
1417  }
1418 
1419  /// \short Return the number of local sub-elements for paraview plot with
1420  /// parameter nplot.
1421  unsigned nsub_elements_paraview(const unsigned& nplot) const
1422  {
1423  return (nplot-1);
1424  }
1425 
1426  /// \short Fill in the offset information for paraview plot.
1427  /// Needs to be implemented for each new geometric element type; see
1428  /// http://www.vtk.org/VTK/img/file-formats.pdf
1429  void write_paraview_output_offset_information(std::ofstream& file_out,
1430  const unsigned& nplot,
1431  unsigned& counter) const
1432  {
1433  // Number of local elements we want to plot over
1434  unsigned plot=nsub_elements_paraview(nplot);
1435 
1436  // loops over the i-th local element in parent element
1437  for(unsigned i=0;i<plot;i++)
1438  {
1439  file_out << i+counter << " "
1440  << i+1+counter
1441  << std::endl;
1442  }
1443  counter+=nplot_points_paraview(nplot);
1444  }
1445 
1446  /// \short Return the paraview element type.
1447  /// Needs to be implemented for each new geometric element type; see
1448  /// http://www.vtk.org/VTK/img/file-formats.pdf
1449  void write_paraview_type(std::ofstream& file_out,
1450  const unsigned& nplot) const
1451  {
1452  unsigned local_loop=nsub_elements_paraview(nplot);
1453  for(unsigned i=0;i<local_loop;i++)
1454  {
1455  file_out << "3" << std::endl;
1456  }
1457  }
1458 
1459  /// \short Return the offsets for the paraview sub-elements. Needs
1460  /// to be implemented for each new geometric element type; see
1461  /// http://www.vtk.org/VTK/img/file-formats.pdf
1462  void write_paraview_offsets(std::ofstream& file_out,
1463  const unsigned& nplot,
1464  unsigned& offset_sum) const
1465  {
1466  // Loop over all local elements and add its offset to the overall offset_sum
1467  unsigned local_loop=nsub_elements_paraview(nplot);
1468  for(unsigned i=0;i<local_loop;i++)
1469  {
1470  offset_sum+=2;
1471  file_out << offset_sum << std::endl;
1472  }
1473  }
1474 
1475  /// Output
1476  void output(std::ostream &output);
1477 
1478  /// Output at specified number of plot points
1479  void output(std::ostream &outfile, const unsigned &nplot);
1480 
1481  /// C-style output
1482  void output(FILE* file_pt);
1483 
1484  /// C_style output at n_plot points
1485  void output(FILE* file_pt, const unsigned &n_plot);
1486 
1487  /// \short Get vector of local coordinates of plot point i (when plotting
1488  /// nplot points in each "coordinate direction").
1489  void get_s_plot(const unsigned& i, const unsigned& nplot,
1490  Vector<double>& s) const
1491  {
1492  if (nplot>1)
1493  {
1494  s[0] = double(i)/double(nplot-1);
1495  }
1496  else
1497  {
1498  s[0]=0.5;
1499  }
1500  }
1501 
1502  /// \short Return string for tecplot zone header (when plotting
1503  /// nplot points in each "coordinate direction)
1504  std::string tecplot_zone_string(const unsigned& nplot) const
1505  {
1506  std::ostringstream header;
1507  header << "ZONE I=" << nplot << "\n";
1508  return header.str();
1509  }
1510 
1511  /// Return total number of plot points (when plotting
1512  /// nplot points in each "coordinate direction)
1513  unsigned nplot_points(const unsigned& nplot) const
1514  {return nplot;}
1515 
1516  /// \short Build the lower-dimensional FaceElement (an element of type
1517  /// PointElement). The face index takes two values
1518  /// corresponding to the two possible faces:
1519  /// -1 (Left) s[0] = -1.0
1520  /// +1 (Right) s[0] = 1.0
1521  void build_face_element(const int &face_index,
1522  FaceElement* face_element_pt);
1523 
1524 };
1525 
1526 
1527 //=======================================================================
1528 /// General TElement class specialised to two spatial dimensions
1529 /// Ordering of nodes as in Zienkiwizc sketches: vertex nodes
1530 /// 0 - 1 - 2 anticlockwise. Midside nodes filled in progressing
1531 /// along the consecutive edges. Central node(s) come(s) last.
1532 //=======================================================================
1533 template<unsigned NNODE_1D>
1534 class TElement<2,NNODE_1D> : public virtual TElementBase,
1535  public TElementShape<2,NNODE_1D>
1536 {
1537  private:
1538 
1539  /// Nodal translation scheme for use when generating face elements
1540  static const unsigned Node_on_face[3][NNODE_1D];
1541 
1542  /// \short Default integration rule: Gaussian integration of same 'order' as
1543  /// the element
1544  //This is sort of optimal, because it means that the integration is exact
1545  //for the shape functions. Can overwrite this in specific element defintion.
1547 
1548 
1549 public:
1550 
1551  /// Constructor
1553  {
1554  // Number of nodes
1555  switch (NNODE_1D)
1556  {
1557  case 2:
1558  case 3:
1559  case 4:
1560  break;
1561 
1562  default:
1563  std::string error_message =
1564  "Triangles are currently only implemented for\n";
1565  error_message +=
1566  "three and six nodes, i.e. NNODE_1D=2 , 3 , 4\n";
1567 
1568  throw OomphLibError(error_message,
1569  OOMPH_CURRENT_FUNCTION,
1570  OOMPH_EXCEPTION_LOCATION);
1571  }
1572 
1573  // Set the number of nodes
1574  unsigned n_node = (NNODE_1D*(NNODE_1D+1))/2;
1575  this->set_n_node(n_node);
1576 
1577  // Set the elemental and nodal dimension
1578  set_dimension(2);
1579 
1580  //Assign default (full) spatial integration scheme
1581  set_integration_scheme(&Default_integration_scheme);
1582  }
1583 
1584  /// Alternative constructor
1585  TElement(const bool &allow_high_order)
1586  {
1587  // Check if we are overriding the restriction on NNODE_1D
1588  if(!allow_high_order)
1589  {
1590  // Call the default constructor
1592  }
1593  else
1594  {
1595  // Set the number of nodes
1596  unsigned n_node = (NNODE_1D*(NNODE_1D+1))/2;
1597  this->set_n_node(n_node);
1598 
1599  // Set the elemental and nodal dimension
1600  set_dimension(2);
1601 
1602  //Assign default (full) spatial integration scheme
1603  set_integration_scheme(&Default_integration_scheme);
1604  }
1605  }
1606 
1607 
1608  /// Broken copy constructor
1610  {
1611  BrokenCopy::broken_copy("TElement");
1612  }
1613 
1614  /// Broken assignment operator
1615  /*void operator=(const TElement&)
1616  {
1617  BrokenCopy::broken_assign("TElement");
1618  }*/
1619 
1620 
1621  /// Destructor
1623 
1624  /// Number of nodes along each element edge
1625  unsigned nnode_1d() const {return NNODE_1D;}
1626 
1627  /// \short Number of vertex nodes in the element: One more
1628  /// than spatial dimension
1629  unsigned nvertex_node() const {return 3;}
1630 
1631  /// \short Public access function for Node_on_face.
1632  unsigned get_bulk_node_number(const int& face_index,
1633  const unsigned& i) const
1634  {
1635  return Node_on_face[face_index][i];
1636  }
1637 
1638  /// \short Pointer to the j-th vertex node in the element
1639  Node* vertex_node_pt(const unsigned& j) const
1640  {
1641  // Vertex nodes come first:
1642 #ifdef PARANOID
1643  if (j>2)
1644  {
1645  std::ostringstream error_message;
1646  error_message
1647  << "Element only has three vertex nodes; called with node number "
1648  << j << std::endl;
1649  throw OomphLibError(error_message.str(),
1650  OOMPH_CURRENT_FUNCTION,
1651  OOMPH_EXCEPTION_LOCATION);
1652  }
1653 #endif
1654  return node_pt(j);
1655  }
1656 
1657  /// Calculate the geometric shape functions at local coordinate s
1658  inline void shape(const Vector<double> &s, Shape &psi) const
1660 
1661  /// \short Compute the geometric shape functions and
1662  /// derivatives w.r.t. local coordinates at local coordinate s
1663  inline void dshape_local(const Vector<double> &s, Shape &psi,
1664  DShape &dpsids) const
1666 
1667  /// \short Compute the geometric shape functions, derivatives and
1668  /// second derivatives w.r.t local coordinates at local coordinate s
1669  /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
1670  /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
1671  /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
1672  inline void d2shape_local(const Vector<double> &s, Shape &psi,
1673  DShape &dpsids, DShape &d2psids) const
1674  {TElementShape<2,NNODE_1D>::d2shape_local(s,psi,dpsids,d2psids);}
1675 
1676  /// \short Overload the template-free interface for the calculation of
1677  /// inverse jacobian matrix. This is a two dimensional element, so use
1678  /// the 2D version.
1680  DenseMatrix<double>&inverse_jacobian) const
1681  {return FiniteElement::invert_jacobian<2>(jacobian,inverse_jacobian);}
1682 
1683  /// Min. value of local coordinate
1684  double s_min() const {return 0.0;}
1685 
1686  /// Max. value of local coordinate
1687  double s_max() const {return 1.0;}
1688 
1689  /// Return local coordinates of node j
1690  inline void local_coordinate_of_node(const unsigned& j,Vector<double>& s) const
1692 
1693  /// \short Return the number of actual plot points for paraview
1694  /// plot with parameter nplot.
1695  unsigned nplot_points_paraview(const unsigned& nplot) const
1696  {
1697  unsigned node_sum=0;
1698  for(unsigned i=1;i<=nplot;i++) {node_sum+=i;}
1699  return node_sum;
1700  }
1701 
1702  /// \short Return the number of local sub-elements for paraview plot with
1703  /// parameter nplot.
1704  unsigned nsub_elements_paraview(const unsigned& nplot) const
1705  {
1706  unsigned local_sum=0;
1707  for(unsigned i=1;i<nplot;i++) {local_sum+=2*(nplot-i-1)+1;}
1708  return local_sum;
1709  }
1710 
1711  /// \short Fill in the offset information for paraview plot.
1712  /// Needs to be implemented for each new geometric element type; see
1713  /// http://www.vtk.org/VTK/img/file-formats.pdf
1714  void write_paraview_output_offset_information(std::ofstream& file_out,
1715  const unsigned& nplot,
1716  unsigned& counter) const
1717  {
1718 
1719  // Outputs list of connectivity of Paraview elements,
1720  // whilst remembering the overall ordering
1721 
1722  unsigned node_count=0;
1723  for(unsigned i=0;i<nplot-1;i++)
1724  {
1725  for(unsigned j=0;j<nplot-i-1;j++)
1726  {
1727  file_out << j+node_count+counter<< " "
1728  << j+node_count+1+counter << " "
1729  << j+nplot+node_count-i+counter << std::endl;
1730 
1731  if(j<nplot-i-2)
1732  {
1733  file_out << j+node_count+1+counter << " "
1734  << j+nplot+node_count-i+1+counter << " "
1735  << j+nplot+node_count-i+counter << std::endl;
1736  }
1737  }
1738  node_count+=(nplot-i);
1739  }
1740  counter+=nplot_points_paraview(nplot);
1741  }
1742 
1743  /// \short Return the paraview element type.
1744  /// Needs to be implemented for each new geometric element type; see
1745  /// http://www.vtk.org/VTK/img/file-formats.pdf
1746  void write_paraview_type(std::ofstream& file_out,
1747  const unsigned& nplot) const
1748  {
1749  unsigned local_loop=nsub_elements_paraview(nplot);
1750 
1751  // Loop over all the local elements and print its paraview type
1752  for(unsigned i=0;i<local_loop;i++)
1753  {
1754  file_out << "5" << std::endl;
1755  }
1756  }
1757 
1758  /// \short Return the offsets for the paraview sub-elements. Needs
1759  /// to be implemented for each new geometric element type; see
1760  /// http://www.vtk.org/VTK/img/file-formats.pdf
1761  void write_paraview_offsets(std::ofstream& file_out,
1762  const unsigned& nplot,
1763  unsigned& offset_sum) const
1764  {
1765  unsigned local_loop=nsub_elements_paraview(nplot);
1766 
1767  // Loop over all local elements and add its offset to the overall offset_sum
1768  for(unsigned i=0;i<local_loop;i++)
1769  {
1770  offset_sum+=3;
1771  file_out << offset_sum << std::endl;
1772  }
1773  }
1774 
1775  /// Output
1776  void output(std::ostream &output);
1777 
1778  /// Output at specified number of plot points
1779  void output(std::ostream &outfile, const unsigned &nplot);
1780 
1781  /// C-style output
1782  void output(FILE* file_pt);
1783 
1784  /// C_style output at n_plot points
1785  void output(FILE* file_pt, const unsigned &n_plot);
1786 
1787  /// \short Get vector of local coordinates of plot point i (when plotting
1788  /// nplot points in each "coordinate direction").
1789  void get_s_plot(const unsigned& iplot, const unsigned& nplot,
1790  Vector<double>& s) const
1791  {
1792  if (nplot>1)
1793  {
1794  unsigned np=0,i,j;
1795  for(i=0;i<nplot;++i)
1796  {
1797  for(j=0;j<nplot-i;++j)
1798  {
1799  if(np==iplot)
1800  {
1801  s[0] = double(j)/double(nplot-1);
1802  s[1] = double(i)/double(nplot-1);
1803  return;
1804  }
1805  ++np;
1806  }
1807  }
1808  }
1809  else
1810  {
1811  s[0] = 1.0/3.0;
1812  s[1] = 1.0/3.0;
1813  }
1814  }
1815 
1816  /// \short Return string for tecplot zone header (when plotting
1817  /// nplot points in each "coordinate direction)
1818  std::string tecplot_zone_string(const unsigned& nplot) const
1819  {
1820  std::ostringstream header;
1821  unsigned nel=0;
1822  for (unsigned i=1;i<nplot;i++) {nel+=2*i-1;}
1823  header << "ZONE N=" << nplot_points(nplot) << ", E="
1824  << nel << ", F=FEPOINT, ET=TRIANGLE\n";
1825  return header.str();
1826  }
1827 
1828  /// \short Add tecplot zone "footer" to output stream (when plotting
1829  /// nplot points in each "coordinate direction).
1830  /// Empty by default -- can be used, e.g., to add FE connectivity
1831  /// lists to elements that need it.
1832  void write_tecplot_zone_footer(std::ostream& outfile,
1833  const unsigned& nplot) const
1834  {
1835  //Output node lists for sub elements for Tecplot (node index
1836  //must start at 1)
1837  unsigned nod_count=1;
1838  for(unsigned i=0;i<nplot;i++)
1839  {
1840  for(unsigned j=0;j<nplot-i;j++)
1841  {
1842  if(j<nplot-i-1)
1843  {
1844  outfile << nod_count << " " << nod_count+1
1845  << " " << nod_count+nplot-i << std::endl;
1846  if(j<nplot-i-2)
1847  {
1848  outfile << nod_count+1 << " "
1849  << nod_count+nplot-i+1 << " "
1850  << nod_count+nplot-i << std::endl;
1851  }
1852  }
1853  ++nod_count;
1854  }
1855  }
1856  }
1857 
1858  /// \short Add tecplot zone "footer" to C-style output. (when plotting
1859  /// nplot points in each "coordinate direction).
1860  /// Empty by default -- can be used, e.g., to add FE connectivity
1861  /// lists to elements that need it.
1862  void write_tecplot_zone_footer(FILE* file_pt,
1863  const unsigned& nplot) const
1864  {
1865  //Output node lists for sub elements for Tecplot (node index
1866  //must start at 1)
1867  unsigned nod_count=1;
1868  for(unsigned i=0;i<nplot;i++)
1869  {
1870  for(unsigned j=0;j<nplot-i;j++)
1871  {
1872  if(j<nplot-i-1)
1873  {
1874  fprintf(file_pt,"%i %i %i \n",nod_count,nod_count+1,
1875  nod_count+nplot-i);
1876  if(j<nplot-i-2)
1877  {
1878  fprintf(file_pt,"%i %i %i \n",nod_count+1,nod_count+nplot-i+1,
1879  nod_count+nplot-i);
1880  }
1881  }
1882  ++nod_count;
1883  }
1884  }
1885  }
1886 
1887  /// Return total number of plot points (when plotting
1888  /// nplot points in each "coordinate direction)
1889  unsigned nplot_points(const unsigned& nplot) const
1890  {
1891  unsigned np=0;
1892  for (unsigned i=1;i<=nplot;i++) {np+=i;}
1893  return np;
1894  }
1895 
1896 
1897  /// \short Build the lower-dimensional FaceElement (an element of type
1898  /// TElement<1,NNODE_1D>). The face index takes three possible values:
1899  /// 0 (Left) s[0] = 0.0
1900  /// 1 (Bottom) s[1] = 0.0
1901  /// 2 (Sloping face) s[0] = 1.0 - s[1]
1902  void build_face_element(const int &face_index,
1903  FaceElement* face_element_pt);
1904 
1905 
1906 };
1907 
1908 
1909 ///////////////////////////////////////////////////////////////////////
1910 ///////////////////////////////////////////////////////////////////////
1911 ///////////////////////////////////////////////////////////////////////
1912 
1913 
1914 //=======================================================================
1915 /// Return local coordinates of node j
1916 //=======================================================================
1917 template<>
1918 class TElementShape<3,2>
1919 {
1920  public:
1921  void local_coordinate_of_node(const unsigned& j,
1922  Vector<double>& s) const
1923  {
1924  s.resize(3);
1925 
1926  switch (j)
1927  {
1928  case 0:
1929  s[0]=1.0;
1930  s[1]=0.0;
1931  s[2]=0.0;
1932  break;
1933 
1934  case 1:
1935  s[0]=0.0;
1936  s[1]=1.0;
1937  s[2]=0.0;
1938  break;
1939 
1940  case 2:
1941  s[0]=0.0;
1942  s[1]=0.0;
1943  s[2]=1.0;
1944  break;
1945 
1946  case 3:
1947  s[0]=0.0;
1948  s[1]=0.0;
1949  s[2]=0.0;
1950  break;
1951 
1952  default:
1953  std::ostringstream error_message;
1954  error_message << "Element only has four nodes; called with node number "
1955  << j << std::endl;
1956 
1957  throw OomphLibError(error_message.str(),
1958  OOMPH_CURRENT_FUNCTION,
1959  OOMPH_EXCEPTION_LOCATION);
1960  }
1961  }
1962 
1963 
1964 
1965 //=======================================================================
1966 /// Shape function for specific TElement<3,2>
1967 //=======================================================================
1968  void shape(const Vector<double> &s, Shape &psi) const
1969  {
1970  psi[0] = s[0];
1971  psi[1] = s[1];
1972  psi[2] = s[2];
1973  psi[3] = 1.0-s[0]-s[1]-s[2];
1974  }
1975 
1976 
1977 //=======================================================================
1978 /// Derivatives of shape functions for specific TElement<3,2>
1979 //=======================================================================
1981  Shape &psi, DShape &dpsids) const
1982  {
1983  //ALH: Don't know why object qualifier is needed
1984  this->shape(s, psi);
1985 
1986  // Derivatives
1987  dpsids(0,0) = 1.0;
1988  dpsids(0,1) = 0.0;
1989  dpsids(0,2) = 0.0;
1990 
1991  dpsids(1,0) = 0.0;
1992  dpsids(1,1) = 1.0;
1993  dpsids(1,2) = 0.0;
1994 
1995  dpsids(2,0) = 0.0;
1996  dpsids(2,1) = 0.0;
1997  dpsids(2,2) = 1.0;
1998 
1999  dpsids(3,0) = -1.0;
2000  dpsids(3,1) = -1.0;
2001  dpsids(3,2) = -1.0;
2002 }
2003 
2004 
2005 
2006 //=======================================================================
2007 /// Second derivatives of shape functions for specific TElement<3,2>:
2008 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
2009 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
2010 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
2011 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
2012 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
2013 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
2014 //=======================================================================
2016  Shape &psi,
2017  DShape &dpsids,
2018  DShape &d2psids) const
2019 {
2020  //ALH: Don't know why object qualifier is needed
2021  this->dshape_local(s, psi,dpsids);
2022 
2023  for(unsigned i=0;i<4;i++)
2024  {
2025  d2psids(i,0) = 0.0;
2026  d2psids(i,1) = 0.0;
2027  d2psids(i,2) = 0.0;
2028  d2psids(i,3) = 0.0;
2029  d2psids(i,4) = 0.0;
2030  d2psids(i,5) = 0.0;
2031  }
2032 }
2033 
2034 };
2035 
2036 
2037 
2038 //=======================================================================
2039 /// Return local coordinates of node j
2040 //=======================================================================
2041 template<>
2042  class TElementShape<3,3>
2043  {
2044  public:
2045  void local_coordinate_of_node(const unsigned& j,
2046  Vector<double>& s) const
2047 {
2048  s.resize(3);
2049 
2050  switch (j)
2051  {
2052  case 0:
2053  s[0]=1.0;
2054  s[1]=0.0;
2055  s[2]=0.0;
2056  break;
2057 
2058  case 1:
2059  s[0]=0.0;
2060  s[1]=1.0;
2061  s[2]=0.0;
2062  break;
2063 
2064  case 2:
2065  s[0]=0.0;
2066  s[1]=0.0;
2067  s[2]=1.0;
2068  break;
2069 
2070  case 3:
2071  s[0]=0.0;
2072  s[1]=0.0;
2073  s[2]=0.0;
2074  break;
2075 
2076  case 4:
2077  s[0]=0.5;
2078  s[1]=0.5;
2079  s[2]=0.0;
2080  break;
2081 
2082  case 5:
2083  s[0]=0.5;
2084  s[1]=0.0;
2085  s[2]=0.5;
2086  break;
2087 
2088  case 6:
2089  s[0]=0.5;
2090  s[1]=0.0;
2091  s[2]=0.0;
2092  break;
2093 
2094  case 7:
2095  s[0]=0.0;
2096  s[1]=0.5;
2097  s[2]=0.5;
2098  break;
2099 
2100  case 8:
2101  s[0]=0.0;
2102  s[1]=0.0;
2103  s[2]=0.5;
2104  break;
2105 
2106  case 9:
2107  s[0]=0.0;
2108  s[1]=0.5;
2109  s[2]=0.0;
2110  break;
2111 
2112  default:
2113  std::ostringstream error_message;
2114  error_message << "Element only has ten nodes; called with node number "
2115  << j << std::endl;
2116 
2117  throw OomphLibError(error_message.str(),
2118  OOMPH_CURRENT_FUNCTION,
2119  OOMPH_EXCEPTION_LOCATION);
2120  }
2121 }
2122 
2123 
2124 
2125 //=======================================================================
2126 /// Shape function for specific TElement<3,3>
2127 //=======================================================================
2128 void shape(const Vector<double> &s, Shape &psi) const
2129 {
2130  double s3=1.0-s[0]-s[1]-s[2];
2131  psi[0] = (2.0*s[0]-1.0)*s[0];
2132  psi[1] = (2.0*s[1]-1.0)*s[1];
2133  psi[2] = (2.0*s[2]-1.0)*s[2];
2134  psi[3] = (2.0*s3-1.0)*s3;
2135  psi[4] = 4.0*s[0]*s[1];
2136  psi[5] = 4.0*s[0]*s[2];
2137  psi[6] = 4.0*s[0]*s3;
2138  psi[7] = 4.0*s[1]*s[2];
2139  psi[8] = 4.0*s[2]*s3;
2140  psi[9] = 4.0*s[1]*s3;
2141 }
2142 
2143 
2144 //=======================================================================
2145 /// Derivatives of shape functions for specific TElement<3,3>
2146 //=======================================================================
2148  Shape &psi, DShape &dpsids) const
2149 {
2150  //ALH: Don't know why object qualifier is needed
2151  this->shape(s, psi);
2152 
2153  // Derivatives
2154  double s3=1.0-s[0]-s[1]-s[2];
2155 
2156  dpsids(0,0) = 4.0*s[0]-1.0;
2157  dpsids(0,1) = 0.0;
2158  dpsids(0,2) = 0.0;
2159 
2160  dpsids(1,0) = 0.0;
2161  dpsids(1,1) = 4.0*s[1]-1.0;
2162  dpsids(1,2) = 0.0;
2163 
2164  dpsids(2,0) = 0.0;
2165  dpsids(2,1) = 0.0;
2166  dpsids(2,2) = 4.0*s[2]-1.0;
2167 
2168  dpsids(3,0) = -4.0*s3+1.0;
2169  dpsids(3,1) = -4.0*s3+1.0;
2170  dpsids(3,2) = -4.0*s3+1.0;
2171 
2172  dpsids(4,0) = 4.0*s[1];
2173  dpsids(4,1) = 4.0*s[0];
2174  dpsids(4,2) = 0.0;
2175 
2176  dpsids(5,0) = 4.0*s[2];
2177  dpsids(5,1) = 0.0;
2178  dpsids(5,2) = 4.0*s[0];
2179 
2180  dpsids(6,0) = 4.0*(s3-s[0]);
2181  dpsids(6,1) = -4.0*s[0];
2182  dpsids(6,2) = -4.0*s[0];
2183 
2184  dpsids(7,0) = 0.0;
2185  dpsids(7,1) = 4.0*s[2];
2186  dpsids(7,2) = 4.0*s[1];
2187 
2188  dpsids(8,0) = -4.0*s[2];
2189  dpsids(8,1) = -4.0*s[2];
2190  dpsids(8,2) = 4.0*(s3-s[2]);
2191 
2192  dpsids(9,0) = -4.0*s[1];
2193  dpsids(9,1) = 4.0*(s3-s[1]);
2194  dpsids(9,2) = -4.0*s[1];
2195 
2196 
2197 }
2198 
2199 
2200 //=======================================================================
2201 /// Second derivatives of shape functions for specific TElement<3,3>:
2202 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
2203 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
2204 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
2205 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
2206 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
2207 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
2208 //=======================================================================
2210  Shape &psi,
2211  DShape &dpsids,
2212  DShape &d2psids) const
2213 {
2214  //ALH: Don't know why object qualifier is needed
2215  this->dshape_local(s, psi,dpsids);
2216 
2217  //(.,3) for mixed derivative s[0]-s[1]
2218  //(.,4) for mixed derivative s[0]-s[2]
2219  //(.,5) for mixed derivative s[1]-s[2]
2220 
2221  d2psids(0,0) = 4.0;
2222  d2psids(0,1) = 0.0;
2223  d2psids(0,2) = 0.0;
2224  d2psids(0,3) = 0.0;
2225  d2psids(0,4) = 0.0;
2226  d2psids(0,5) = 0.0;
2227 
2228 
2229  d2psids(1,0) = 0.0;
2230  d2psids(1,1) = 4.0;
2231  d2psids(1,2) = 0.0;
2232  d2psids(1,3) = 0.0;
2233  d2psids(1,4) = 0.0;
2234  d2psids(1,5) = 0.0;
2235 
2236  d2psids(2,0) = 0.0;
2237  d2psids(2,1) = 0.0;
2238  d2psids(2,2) = 4.0;
2239  d2psids(2,3) = 0.0;
2240  d2psids(2,4) = 0.0;
2241  d2psids(2,5) = 0.0;
2242 
2243  d2psids(3,0) = 4.0;
2244  d2psids(3,1) = 4.0;
2245  d2psids(3,2) = 4.0;
2246  d2psids(3,3) = 4.0;
2247  d2psids(3,4) = 4.0;
2248  d2psids(3,5) = 4.0;
2249 
2250  d2psids(4,0) = 0.0;
2251  d2psids(4,1) = 0.0;
2252  d2psids(4,2) = 0.0;
2253  d2psids(4,3) = 4.0;
2254  d2psids(4,4) = 0.0;
2255  d2psids(4,5) = 0.0;
2256 
2257  d2psids(5,0) = 0.0;
2258  d2psids(5,1) = 0.0;
2259  d2psids(5,2) = 0.0;
2260  d2psids(5,3) = 0.0;
2261  d2psids(5,4) = 4.0;
2262  d2psids(5,5) = 0.0;
2263 
2264  d2psids(6,0) =-8.0;
2265  d2psids(6,1) = 0.0;
2266  d2psids(6,2) = 0.0;
2267  d2psids(6,3) = -4.0;
2268  d2psids(6,4) = -4.0;
2269  d2psids(6,5) = 0.0;
2270 
2271  d2psids(7,0) = 0.0;
2272  d2psids(7,1) = 0.0;
2273  d2psids(7,2) = 0.0;
2274  d2psids(7,3) = 0.0;
2275  d2psids(7,4) = 0.0;
2276  d2psids(7,5) = 4.0;
2277 
2278  d2psids(8,0) = 0.0;
2279  d2psids(8,1) = 0.0;
2280  d2psids(8,2) = -8.0;
2281  d2psids(8,3) = 0.0;
2282  d2psids(8,4) = -4.0;
2283  d2psids(8,5) = -4.0;
2284 
2285  d2psids(9,0) = 0.0;
2286  d2psids(9,1) = -8.0;
2287  d2psids(9,2) = 0.0;
2288  d2psids(9,3) = -4.0;
2289  d2psids(9,4) = 0.0;
2290  d2psids(9,5) = -4.0;
2291 
2292 }
2293 
2294 
2295 };
2296 
2297 //////////////////////////////////////////////////////////////////////
2298 //////////////////////////////////////////////////////////////////////
2299 //////////////////////////////////////////////////////////////////////
2300 
2301 //====================================================================
2302 ///Standard quadratic shape functions enriched by the addition of
2303 ///three cubic "face" bubbles and quartic "volume" bubble,
2304 ///which consists of adding a node at the centroid of
2305 ///each face and a single node at the centroid
2306 ///of the tetrahedron
2307 //=========================================================================
2308 
2309 //=======================================================================
2310 /// Return local coordinates of node j
2311 //=======================================================================
2312 template<>
2314  {
2315  public:
2316 
2317  unsigned n_enriched_nodes() {return 5;}
2318 
2319  void local_coordinate_of_node(const unsigned& j,
2320  Vector<double>& s) const
2321 {
2322  s.resize(3);
2323 
2324  switch (j)
2325  {
2326  case 0:
2327  s[0]=1.0;
2328  s[1]=0.0;
2329  s[2]=0.0;
2330  break;
2331 
2332  case 1:
2333  s[0]=0.0;
2334  s[1]=1.0;
2335  s[2]=0.0;
2336  break;
2337 
2338  case 2:
2339  s[0]=0.0;
2340  s[1]=0.0;
2341  s[2]=1.0;
2342  break;
2343 
2344  case 3:
2345  s[0]=0.0;
2346  s[1]=0.0;
2347  s[2]=0.0;
2348  break;
2349 
2350  case 4:
2351  s[0]=0.5;
2352  s[1]=0.5;
2353  s[2]=0.0;
2354  break;
2355 
2356  case 5:
2357  s[0]=0.5;
2358  s[1]=0.0;
2359  s[2]=0.5;
2360  break;
2361 
2362  case 6:
2363  s[0]=0.5;
2364  s[1]=0.0;
2365  s[2]=0.0;
2366  break;
2367 
2368  case 7:
2369  s[0]=0.0;
2370  s[1]=0.5;
2371  s[2]=0.5;
2372  break;
2373 
2374  case 8:
2375  s[0]=0.0;
2376  s[1]=0.0;
2377  s[2]=0.5;
2378  break;
2379 
2380  case 9:
2381  s[0]=0.0;
2382  s[1]=0.5;
2383  s[2]=0.0;
2384  break;
2385 
2386  //Node at centroid of face spanned by nodes 0, 1, 3
2387  case 10:
2388  s[0]=1.0/3.0;
2389  s[1]=1.0/3.0;
2390  s[2]=0.0;
2391  break;
2392 
2393  //Node at centroid of face spanned by nodes 0, 1, 2
2394  case 11:
2395  s[0]=1.0/3.0;
2396  s[1]=1.0/3.0;
2397  s[2]=1.0/3.0;
2398  break;
2399 
2400  //Node at centroid of face spanned by nodes 0, 2, 3
2401  case 12:
2402  s[0]=1.0/3.0;
2403  s[1]=0.0;
2404  s[2]=1.0/3.0;
2405  break;
2406 
2407  //Node at centroid of face spannd by nodes 1, 2, 3
2408  case 13:
2409  s[0]=0.0;
2410  s[1]=1.0/3.0;
2411  s[2]=1.0/3.0;
2412  break;
2413 
2414  //Node at centroid of volume
2415  case 14:
2416  s[0] = 0.25;
2417  s[1] = 0.25;
2418  s[2] = 0.25;
2419  break;
2420 
2421 
2422  default:
2423  std::ostringstream error_message;
2424  error_message << "Element only has fifteen nodes; called with node number "
2425  << j << std::endl;
2426 
2427  throw OomphLibError(error_message.str(),
2428  OOMPH_CURRENT_FUNCTION,
2429  OOMPH_EXCEPTION_LOCATION);
2430  }
2431 }
2432 
2433 
2434 
2435 //=======================================================================
2436 /// Shape function for specific TBubbleEnrichedElement<3,3>
2437 //=======================================================================
2438 void shape(const Vector<double> &s, Shape &psi) const
2439  {
2440  //Constructe the fourth volume coordinate
2441  const double s3=1.0-s[0]-s[1]-s[2];
2442  //calculate the enrichment functions
2443  const double quartic_bubble = s[0]*s[1]*s[2]*s3;
2444  const double cubic_bubble012 = s[0]*s[1]*s[2];
2445  const double cubic_bubble013 = s[0]*s[1]*s3;
2446  const double cubic_bubble023 = s[0]*s[2]*s3;
2447  const double cubic_bubble123 = s[1]*s[2]*s3;
2448 
2449  //The appropriate "amount" of cubic and quartic bubble functions are
2450  //added/subtracted
2451  //to each original quadratic shape function to ensure that the new
2452  //shape function is zero at the centroid (0.25,0.25,0.25)
2453  //and at the face centroids
2454  psi[0] = (2.0*s[0]-1.0)*s[0]
2455  + 3.0*(cubic_bubble012 + cubic_bubble013 + cubic_bubble023)
2456  - 4.0*quartic_bubble;
2457  psi[1] = (2.0*s[1]-1.0)*s[1]
2458  + 3.0*(cubic_bubble012 + cubic_bubble013 + cubic_bubble123)
2459  - 4.0*quartic_bubble;
2460  psi[2] = (2.0*s[2]-1.0)*s[2]
2461  + 3.0*(cubic_bubble012 + cubic_bubble023 + cubic_bubble123)
2462  - 4.0*quartic_bubble;
2463  psi[3] = (2.0*s3-1.0)*s3
2464  + 3.0*(cubic_bubble013 + cubic_bubble023 + cubic_bubble123)
2465  -4.0*quartic_bubble;
2466  psi[4] = 4.0*s[0]*s[1]
2467  - 12.0*(cubic_bubble012 + cubic_bubble013)
2468  + 32.0*quartic_bubble;
2469  psi[5] = 4.0*s[0]*s[2]
2470  - 12.0*(cubic_bubble012 + cubic_bubble023)
2471  + 32.0*quartic_bubble;
2472  psi[6] = 4.0*s[0]*s3
2473  - 12.0*(cubic_bubble013 + cubic_bubble023)
2474  + 32.0*quartic_bubble;
2475  psi[7] = 4.0*s[1]*s[2]
2476  - 12.0*(cubic_bubble012 + cubic_bubble123)
2477  + 32.0*quartic_bubble;
2478  psi[8] = 4.0*s[2]*s3
2479  - 12.0*(cubic_bubble023 + cubic_bubble123)
2480  + 32.0*quartic_bubble;
2481  psi[9] = 4.0*s[1]*s3
2482  - 12.0*(cubic_bubble013 + cubic_bubble123)
2483  + 32.0*quartic_bubble;
2484  //Add the bubble function on the face spanned by 0 1 3
2485  psi[10] = 27.0*cubic_bubble013 - 108.0*quartic_bubble;
2486  //Add the bubble function on the face spanned by 0 1 2
2487  psi[11] = 27.0*cubic_bubble012 - 108.0*quartic_bubble;
2488  //Add the bubble function on the face spanned by 0 2 3
2489  psi[12] = 27.0*cubic_bubble023 - 108.0*quartic_bubble;
2490  //Add the bubble function on the face spanned by 1 2 3
2491  psi[13] = 27.0*cubic_bubble123 - 108.0*quartic_bubble;
2492  //Add the volume bubble function, scaled to have value one
2493  psi[14] = 256.0*quartic_bubble;
2494 }
2495 
2496 
2497 //=======================================================================
2498 /// Derivatives of shape functions for specific TElement<3,3>
2499 //=======================================================================
2501  Shape &psi, DShape &dpsids) const
2502 {
2503  //ALH: Don't know why object qualifier is needed
2504  this->shape(s, psi);
2505 
2506  //Define s3 the fourth volume coordinate
2507  const double s3=1.0-s[0]-s[1]-s[2];
2508 
2509  //Calculate derivatives of the bubble function
2510  const double d_quartic_bubble_ds0 = s[1]*s[2]*(1.0 - s[1] - s[2] - 2.0*s[0]);
2511  const double d_quartic_bubble_ds1 = s[0]*s[2]*(1.0 - s[0] - s[2] - 2.0*s[1]);
2512  const double d_quartic_bubble_ds2 = s[0]*s[1]*(1.0 - s[0] - s[1] - 2.0*s[2]);
2513 
2514  const double d_cubic_bubble012_ds0 = s[1]*s[2];
2515  const double d_cubic_bubble012_ds1 = s[0]*s[2];
2516  const double d_cubic_bubble012_ds2 = s[0]*s[1];
2517 
2518  const double d_cubic_bubble013_ds0 = s[1]*(s3 - s[0]);
2519  const double d_cubic_bubble013_ds1 = s[0]*(s3 - s[1]);
2520  const double d_cubic_bubble013_ds2 = -s[0]*s[1];
2521 
2522  const double d_cubic_bubble023_ds0 = s[2]*(s3 - s[0]);
2523  const double d_cubic_bubble023_ds1 = -s[0]*s[2];
2524  const double d_cubic_bubble023_ds2 = s[0]*(s3 - s[2]);
2525 
2526  const double d_cubic_bubble123_ds0 = -s[1]*s[2];
2527  const double d_cubic_bubble123_ds1 = s[2]*(s3 - s[1]);
2528  const double d_cubic_bubble123_ds2 = s[1]*(s3 - s[2]);
2529 
2530 
2531  //Add the appropriate dervatives of the bubble function to the
2532  //shape function derivatives
2533  dpsids(0,0) = 4.0*s[0]-1.0
2534  + 3.0*(d_cubic_bubble012_ds0 + d_cubic_bubble013_ds0 + d_cubic_bubble023_ds0)
2535  - 4.0*d_quartic_bubble_ds0;
2536  dpsids(0,1) = 0.0
2537  + 3.0*(d_cubic_bubble012_ds1 + d_cubic_bubble013_ds1 + d_cubic_bubble023_ds1)
2538  - 4.0*d_quartic_bubble_ds1;
2539  dpsids(0,2) = 0.0
2540  + 3.0*(d_cubic_bubble012_ds2 + d_cubic_bubble013_ds2 + d_cubic_bubble023_ds2)
2541  - 4.0*d_quartic_bubble_ds2;
2542 
2543  dpsids(1,0) = 0.0
2544  + 3.0*(d_cubic_bubble012_ds0 + d_cubic_bubble013_ds0 + d_cubic_bubble123_ds0)
2545  - 4.0*d_quartic_bubble_ds0;
2546  dpsids(1,1) = 4.0*s[1]-1.0
2547  + 3.0*(d_cubic_bubble012_ds1 + d_cubic_bubble013_ds1 + d_cubic_bubble123_ds1)
2548  - 4.0*d_quartic_bubble_ds1;
2549  dpsids(1,2) = 0.0
2550  + 3.0*(d_cubic_bubble012_ds2 + d_cubic_bubble013_ds2 + d_cubic_bubble123_ds2)
2551  - 4.0*d_quartic_bubble_ds2;
2552 
2553  dpsids(2,0) = 0.0
2554  + 3.0*(d_cubic_bubble012_ds0 + d_cubic_bubble023_ds0 + d_cubic_bubble123_ds0)
2555  - 4.0*d_quartic_bubble_ds0;
2556  dpsids(2,1) = 0.0
2557  + 3.0*(d_cubic_bubble012_ds1 + d_cubic_bubble023_ds1 + d_cubic_bubble123_ds1)
2558  - 4.0*d_quartic_bubble_ds1;
2559  dpsids(2,2) = 4.0*s[2]-1.0
2560  + 3.0*(d_cubic_bubble012_ds2 + d_cubic_bubble023_ds2 + d_cubic_bubble123_ds2)
2561  - 4.0*d_quartic_bubble_ds2;
2562 
2563  dpsids(3,0) = -4.0*s3+1.0
2564  + 3.0*(d_cubic_bubble013_ds0 + d_cubic_bubble023_ds0 + d_cubic_bubble123_ds0)
2565  -4.0*d_quartic_bubble_ds0;
2566  dpsids(3,1) = -4.0*s3+1.0
2567  + 3.0*(d_cubic_bubble013_ds1 + d_cubic_bubble023_ds1 + d_cubic_bubble123_ds1)
2568  -4.0*d_quartic_bubble_ds1;
2569  dpsids(3,2) = -4.0*s3+1.0
2570  + 3.0*(d_cubic_bubble013_ds2 + d_cubic_bubble023_ds2 + d_cubic_bubble123_ds2)
2571  -4.0*d_quartic_bubble_ds2;
2572 
2573  dpsids(4,0) = 4.0*s[1]
2574  - 12.0*(d_cubic_bubble012_ds0 + d_cubic_bubble013_ds0)
2575  + 32.0*d_quartic_bubble_ds0;
2576  dpsids(4,1) = 4.0*s[0]
2577  - 12.0*(d_cubic_bubble012_ds1 + d_cubic_bubble013_ds1)
2578  + 32.0*d_quartic_bubble_ds1;
2579  dpsids(4,2) = 0.0
2580  - 12.0*(d_cubic_bubble012_ds2 + d_cubic_bubble013_ds2)
2581  + 32.0*d_quartic_bubble_ds2;
2582 
2583  dpsids(5,0) = 4.0*s[2]
2584  - 12.0*(d_cubic_bubble012_ds0 + d_cubic_bubble023_ds0)
2585  + 32.0*d_quartic_bubble_ds0;
2586  dpsids(5,1) = 0.0
2587  - 12.0*(d_cubic_bubble012_ds1 + d_cubic_bubble023_ds1)
2588  + 32.0*d_quartic_bubble_ds1;
2589  dpsids(5,2) = 4.0*s[0]
2590  - 12.0*(d_cubic_bubble012_ds2 + d_cubic_bubble023_ds2)
2591  + 32.0*d_quartic_bubble_ds2;
2592 
2593  dpsids(6,0) = 4.0*(s3-s[0])
2594  - 12.0*(d_cubic_bubble013_ds0 + d_cubic_bubble023_ds0)
2595  + 32.0*d_quartic_bubble_ds0;
2596  dpsids(6,1) = -4.0*s[0]
2597  - 12.0*(d_cubic_bubble013_ds1 + d_cubic_bubble023_ds1)
2598  + 32.0*d_quartic_bubble_ds1;
2599  dpsids(6,2) = -4.0*s[0]
2600  - 12.0*(d_cubic_bubble013_ds2 + d_cubic_bubble023_ds2)
2601  + 32.0*d_quartic_bubble_ds2;
2602 
2603  dpsids(7,0) = 0.0
2604  - 12.0*(d_cubic_bubble012_ds0 + d_cubic_bubble123_ds0)
2605  + 32.0*d_quartic_bubble_ds0;
2606  dpsids(7,1) = 4.0*s[2]
2607  - 12.0*(d_cubic_bubble012_ds1 + d_cubic_bubble123_ds1)
2608  + 32.0*d_quartic_bubble_ds1;
2609  dpsids(7,2) = 4.0*s[1]
2610  - 12.0*(d_cubic_bubble012_ds2 + d_cubic_bubble123_ds2)
2611  + 32.0*d_quartic_bubble_ds2;
2612 
2613  dpsids(8,0) = -4.0*s[2]
2614  - 12.0*(d_cubic_bubble023_ds0 + d_cubic_bubble123_ds0)
2615  + 32.0*d_quartic_bubble_ds0;
2616  dpsids(8,1) = -4.0*s[2]
2617  - 12.0*(d_cubic_bubble023_ds1 + d_cubic_bubble123_ds1)
2618  + 32.0*d_quartic_bubble_ds1;
2619  dpsids(8,2) = 4.0*(s3-s[2])
2620  - 12.0*(d_cubic_bubble023_ds2 + d_cubic_bubble123_ds2)
2621  + 32.0*d_quartic_bubble_ds2;
2622 
2623  dpsids(9,0) = -4.0*s[1]
2624  - 12.0*(d_cubic_bubble013_ds0 + d_cubic_bubble123_ds0)
2625  + 32.0*d_quartic_bubble_ds0;
2626  dpsids(9,1) = 4.0*(s3-s[1])
2627  - 12.0*(d_cubic_bubble013_ds1 + d_cubic_bubble123_ds1)
2628  + 32.0*d_quartic_bubble_ds1;
2629  dpsids(9,2) = -4.0*s[1]
2630  - 12.0*(d_cubic_bubble013_ds2 + d_cubic_bubble123_ds2)
2631  + 32.0*d_quartic_bubble_ds2;
2632 
2633  //Add the bubble function on the face spanned by 0 1 3
2634  dpsids(10,0) = 27.0*d_cubic_bubble013_ds0 - 108.0*d_quartic_bubble_ds0;
2635  dpsids(10,1) = 27.0*d_cubic_bubble013_ds1 - 108.0*d_quartic_bubble_ds1;
2636  dpsids(10,2) = 27.0*d_cubic_bubble013_ds2 - 108.0*d_quartic_bubble_ds2;
2637 
2638  //Add the bubble function on the face spanned by 0 1 2
2639  dpsids(11,0) = 27.0*d_cubic_bubble012_ds0 - 108.0*d_quartic_bubble_ds0;
2640  dpsids(11,1) = 27.0*d_cubic_bubble012_ds1 - 108.0*d_quartic_bubble_ds1;
2641  dpsids(11,2) = 27.0*d_cubic_bubble012_ds2 - 108.0*d_quartic_bubble_ds2;
2642 
2643  //Add the bubble function on the face spanned by 0 2 3
2644  dpsids(12,0) = 27.0*d_cubic_bubble023_ds0 - 108.0*d_quartic_bubble_ds0;
2645  dpsids(12,1) = 27.0*d_cubic_bubble023_ds1 - 108.0*d_quartic_bubble_ds1;
2646  dpsids(12,2) = 27.0*d_cubic_bubble023_ds2 - 108.0*d_quartic_bubble_ds2;
2647 
2648  //Add the bubble function on the face spanned by 1 2 3
2649  dpsids(13,0) = 27.0*d_cubic_bubble123_ds0 - 108.0*d_quartic_bubble_ds0;
2650  dpsids(13,1) = 27.0*d_cubic_bubble123_ds1 - 108.0*d_quartic_bubble_ds1;
2651  dpsids(13,2) = 27.0*d_cubic_bubble123_ds2 - 108.0*d_quartic_bubble_ds2;
2652 
2653  //Add the volumetric bubble function derivatives
2654  dpsids(14,0) = 256.0*d_quartic_bubble_ds0;
2655  dpsids(14,1) = 256.0*d_quartic_bubble_ds1;
2656  dpsids(14,2) = 256.0*d_quartic_bubble_ds2;
2657 }
2658 
2659 
2660 //=======================================================================
2661 /// Second derivatives of shape functions for specific TElement<3,3>:
2662 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
2663 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
2664 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
2665 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
2666 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
2667 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
2668 //=======================================================================
2670  Shape &psi,
2671  DShape &dpsids,
2672  DShape &d2psids) const
2673 {
2674  //ALH: Don't know why object qualifier is needed
2675  this->dshape_local(s, psi,dpsids);
2676 
2677  // Define s3
2678  const double s3=1.0-s[0]-s[1]-s[2];
2679 
2680  //Calculate second derivatives of the bubble functions
2681  //(.,3) for mixed derivative s[0]-s[1]
2682  //(.,4) for mixed derivative s[0]-s[2]
2683  //(.,5) for mixed derivative s[1]-s[2]
2684 
2685 
2686  const double d2_quartic_bubble_ds0 = -2.0*s[1]*s[2];
2687  const double d2_quartic_bubble_ds1 = -2.0*s[0]*s[2];
2688  const double d2_quartic_bubble_ds2 = -2.0*s[0]*s[1];
2689  const double d2_quartic_bubble_ds3 = s[2]*(1.0 - 2.0*s[0] - 2.0*s[1] - s[2]);
2690  const double d2_quartic_bubble_ds4 = s[1]*(1.0 - 2.0*s[0] - 2.0*s[2] - s[1]);
2691  const double d2_quartic_bubble_ds5 = s[0]*(1.0 - 2.0*s[1] - 2.0*s[2] - s[0]);
2692 
2693  const double d2_cubic_bubble012_ds0 = 0.0;
2694  const double d2_cubic_bubble012_ds1 = 0.0;
2695  const double d2_cubic_bubble012_ds2 = 0.0;
2696  const double d2_cubic_bubble012_ds3 = s[2];
2697  const double d2_cubic_bubble012_ds4 = s[1];
2698  const double d2_cubic_bubble012_ds5 = s[0];
2699 
2700  const double d2_cubic_bubble013_ds0 = -2.0*s[1];
2701  const double d2_cubic_bubble013_ds1 = -2.0*s[0];
2702  const double d2_cubic_bubble013_ds2 = 0.0;
2703  const double d2_cubic_bubble013_ds3 = s3 - s[0] - s[1];
2704  const double d2_cubic_bubble013_ds4 = -s[1];
2705  const double d2_cubic_bubble013_ds5 = -s[0];
2706 
2707  const double d2_cubic_bubble023_ds0 = -2.0*s[2];
2708  const double d2_cubic_bubble023_ds1 = 0.0;
2709  const double d2_cubic_bubble023_ds2 = -2.0*s[0];
2710  const double d2_cubic_bubble023_ds3 = -s[2];
2711  const double d2_cubic_bubble023_ds4 = s3 - s[0] - s[2];
2712  const double d2_cubic_bubble023_ds5 = -s[0];
2713 
2714  const double d2_cubic_bubble123_ds0 = 0.0;
2715  const double d2_cubic_bubble123_ds1 = -2.0*s[2];
2716  const double d2_cubic_bubble123_ds2 = -2.0*s[1];
2717  const double d2_cubic_bubble123_ds3 = -s[2];
2718  const double d2_cubic_bubble123_ds4 = -s[1];
2719  const double d2_cubic_bubble123_ds5 = s3 - s[1] - s[2];
2720 
2721 
2722  d2psids(0,0) = 4.0
2723  + 3.0*(d2_cubic_bubble012_ds0 + d2_cubic_bubble013_ds0
2724  + d2_cubic_bubble023_ds0)
2725  - 4.0*d2_quartic_bubble_ds0;
2726  d2psids(0,1) = 0.0
2727  + 3.0*(d2_cubic_bubble012_ds1 + d2_cubic_bubble013_ds1
2728  + d2_cubic_bubble023_ds1)
2729  - 4.0*d2_quartic_bubble_ds1;
2730  d2psids(0,2) = 0.0
2731  + 3.0*(d2_cubic_bubble012_ds2 + d2_cubic_bubble013_ds2
2732  + d2_cubic_bubble023_ds2)
2733  - 4.0*d2_quartic_bubble_ds2;
2734  d2psids(0,3) = 0.0
2735  + 3.0*(d2_cubic_bubble012_ds3 + d2_cubic_bubble013_ds3
2736  + d2_cubic_bubble023_ds3)
2737  - 4.0*d2_quartic_bubble_ds3;
2738  d2psids(0,4) = 0.0
2739  + 3.0*(d2_cubic_bubble012_ds4 + d2_cubic_bubble013_ds4
2740  + d2_cubic_bubble023_ds4)
2741  - 4.0*d2_quartic_bubble_ds4;
2742  d2psids(0,5) = 0.0
2743  + 3.0*(d2_cubic_bubble012_ds5 + d2_cubic_bubble013_ds5
2744  + d2_cubic_bubble023_ds5)
2745  - 4.0*d2_quartic_bubble_ds5;
2746 
2747 
2748  d2psids(1,0) = 0.0
2749  + 3.0*(d2_cubic_bubble012_ds0 + d2_cubic_bubble013_ds0
2750  + d2_cubic_bubble123_ds0)
2751  - 4.0*d2_quartic_bubble_ds0;
2752  d2psids(1,1) = 4.0
2753  + 3.0*(d2_cubic_bubble012_ds1 + d2_cubic_bubble013_ds1
2754  + d2_cubic_bubble123_ds1)
2755  - 4.0*d2_quartic_bubble_ds1;
2756  d2psids(1,2) = 0.0
2757  + 3.0*(d2_cubic_bubble012_ds2 + d2_cubic_bubble013_ds2
2758  + d2_cubic_bubble123_ds2)
2759  - 4.0*d2_quartic_bubble_ds2;
2760  d2psids(1,3) = 0.0
2761  + 3.0*(d2_cubic_bubble012_ds3 + d2_cubic_bubble013_ds3
2762  + d2_cubic_bubble123_ds3)
2763  - 4.0*d2_quartic_bubble_ds3;
2764  d2psids(1,4) = 0.0
2765  + 3.0*(d2_cubic_bubble012_ds4 + d2_cubic_bubble013_ds4
2766  + d2_cubic_bubble123_ds4)
2767  - 4.0*d2_quartic_bubble_ds4;
2768  d2psids(1,5) = 0.0
2769  + 3.0*(d2_cubic_bubble012_ds5 + d2_cubic_bubble013_ds5
2770  + d2_cubic_bubble123_ds5)
2771  - 4.0*d2_quartic_bubble_ds5;
2772 
2773 
2774  d2psids(2,0) = 0.0
2775  + 3.0*(d2_cubic_bubble012_ds0 + d2_cubic_bubble023_ds0
2776  + d2_cubic_bubble123_ds0)
2777  - 4.0*d2_quartic_bubble_ds0;
2778  d2psids(2,1) = 0.0
2779  + 3.0*(d2_cubic_bubble012_ds1 + d2_cubic_bubble023_ds1
2780  + d2_cubic_bubble123_ds1)
2781  - 4.0*d2_quartic_bubble_ds1;
2782  d2psids(2,2) = 4.0
2783  + 3.0*(d2_cubic_bubble012_ds2 + d2_cubic_bubble023_ds2
2784  + d2_cubic_bubble123_ds2)
2785  - 4.0*d2_quartic_bubble_ds2;
2786  d2psids(2,3) = 0.0
2787  + 3.0*(d2_cubic_bubble012_ds3 + d2_cubic_bubble023_ds3
2788  + d2_cubic_bubble123_ds3)
2789  - 4.0*d2_quartic_bubble_ds3;
2790  d2psids(2,4) = 0.0
2791  + 3.0*(d2_cubic_bubble012_ds4 + d2_cubic_bubble023_ds4
2792  + d2_cubic_bubble123_ds4)
2793  - 4.0*d2_quartic_bubble_ds4;
2794  d2psids(2,5) = 0.0
2795  + 3.0*(d2_cubic_bubble012_ds5 + d2_cubic_bubble023_ds5
2796  + d2_cubic_bubble123_ds5)
2797  - 4.0*d2_quartic_bubble_ds5;
2798 
2799 
2800  d2psids(3,0) = 4.0
2801  + 3.0*(d2_cubic_bubble013_ds0 + d2_cubic_bubble023_ds0
2802  + d2_cubic_bubble123_ds0)
2803  -4.0*d2_quartic_bubble_ds0;
2804  d2psids(3,1) = 4.0
2805  + 3.0*(d2_cubic_bubble013_ds1 + d2_cubic_bubble023_ds1
2806  + d2_cubic_bubble123_ds1)
2807  -4.0*d2_quartic_bubble_ds1;
2808  d2psids(3,2) = 4.0
2809  + 3.0*(d2_cubic_bubble013_ds2 + d2_cubic_bubble023_ds2
2810  + d2_cubic_bubble123_ds2)
2811  -4.0*d2_quartic_bubble_ds2;
2812  d2psids(3,3) = 4.0
2813  + 3.0*(d2_cubic_bubble013_ds3 + d2_cubic_bubble023_ds3
2814  + d2_cubic_bubble123_ds3)
2815  -4.0*d2_quartic_bubble_ds3;
2816  d2psids(3,4) = 4.0
2817  + 3.0*(d2_cubic_bubble013_ds4 + d2_cubic_bubble023_ds4
2818  + d2_cubic_bubble123_ds4)
2819  -4.0*d2_quartic_bubble_ds4;
2820  d2psids(3,5) = 4.0
2821  + 3.0*(d2_cubic_bubble013_ds5 + d2_cubic_bubble023_ds5
2822  + d2_cubic_bubble123_ds5)
2823  -4.0*d2_quartic_bubble_ds5;
2824 
2825 
2826  d2psids(4,0) = 0.0
2827  - 12.0*(d2_cubic_bubble012_ds0 + d2_cubic_bubble013_ds0)
2828  + 32.0*d2_quartic_bubble_ds0;
2829  d2psids(4,1) = 0.0
2830  - 12.0*(d2_cubic_bubble012_ds1 + d2_cubic_bubble013_ds1)
2831  + 32.0*d2_quartic_bubble_ds1;
2832  d2psids(4,2) = 0.0
2833  - 12.0*(d2_cubic_bubble012_ds2 + d2_cubic_bubble013_ds2)
2834  + 32.0*d2_quartic_bubble_ds2;
2835  d2psids(4,3) = 4.0
2836  - 12.0*(d2_cubic_bubble012_ds3 + d2_cubic_bubble013_ds3)
2837  + 32.0*d2_quartic_bubble_ds3;
2838  d2psids(4,4) = 0.0
2839  - 12.0*(d2_cubic_bubble012_ds4 + d2_cubic_bubble013_ds4)
2840  + 32.0*d2_quartic_bubble_ds4;
2841  d2psids(4,5) = 0.0
2842  - 12.0*(d2_cubic_bubble012_ds5 + d2_cubic_bubble013_ds5)
2843  + 32.0*d2_quartic_bubble_ds5;
2844 
2845 
2846  d2psids(5,0) = 0.0
2847  - 12.0*(d2_cubic_bubble012_ds0 + d2_cubic_bubble023_ds0)
2848  + 32.0*d2_quartic_bubble_ds0;
2849  d2psids(5,1) = 0.0
2850  - 12.0*(d2_cubic_bubble012_ds1 + d2_cubic_bubble023_ds1)
2851  + 32.0*d2_quartic_bubble_ds1;
2852  d2psids(5,2) = 0.0
2853  - 12.0*(d2_cubic_bubble012_ds2 + d2_cubic_bubble023_ds2)
2854  + 32.0*d2_quartic_bubble_ds2;
2855  d2psids(5,3) = 0.0
2856  - 12.0*(d2_cubic_bubble012_ds3 + d2_cubic_bubble023_ds3)
2857  + 32.0*d2_quartic_bubble_ds3;
2858  d2psids(5,4) = 4.0
2859  - 12.0*(d2_cubic_bubble012_ds4 + d2_cubic_bubble023_ds4)
2860  + 32.0*d2_quartic_bubble_ds4;
2861  d2psids(5,5) = 0.0
2862  - 12.0*(d2_cubic_bubble012_ds5 + d2_cubic_bubble023_ds5)
2863  + 32.0*d2_quartic_bubble_ds5;
2864 
2865 
2866  d2psids(6,0) =-8.0
2867  - 12.0*(d2_cubic_bubble013_ds0 + d2_cubic_bubble023_ds0)
2868  + 32.0*d2_quartic_bubble_ds0;
2869  d2psids(6,1) = 0.0
2870  - 12.0*(d2_cubic_bubble013_ds1 + d2_cubic_bubble023_ds1)
2871  + 32.0*d2_quartic_bubble_ds1;
2872  d2psids(6,2) = 0.0
2873  - 12.0*(d2_cubic_bubble013_ds2 + d2_cubic_bubble023_ds2)
2874  + 32.0*d2_quartic_bubble_ds2;
2875  d2psids(6,3) = -4.0
2876  - 12.0*(d2_cubic_bubble013_ds3 + d2_cubic_bubble023_ds3)
2877  + 32.0*d2_quartic_bubble_ds3;
2878  d2psids(6,4) = -4.0
2879  - 12.0*(d2_cubic_bubble013_ds4 + d2_cubic_bubble023_ds4)
2880  + 32.0*d2_quartic_bubble_ds4;
2881  d2psids(6,5) = 0.0
2882  - 12.0*(d2_cubic_bubble013_ds5 + d2_cubic_bubble023_ds5)
2883  + 32.0*d2_quartic_bubble_ds5;
2884 
2885  d2psids(7,0) = 0.0
2886  - 12.0*(d2_cubic_bubble012_ds0 + d2_cubic_bubble123_ds0)
2887  + 32.0*d2_quartic_bubble_ds0;
2888  d2psids(7,1) = 0.0
2889  - 12.0*(d2_cubic_bubble012_ds1 + d2_cubic_bubble123_ds1)
2890  + 32.0*d2_quartic_bubble_ds1;
2891  d2psids(7,2) = 0.0
2892  - 12.0*(d2_cubic_bubble012_ds2 + d2_cubic_bubble123_ds2)
2893  + 32.0*d2_quartic_bubble_ds2;
2894  d2psids(7,3) = 0.0
2895  - 12.0*(d2_cubic_bubble012_ds3 + d2_cubic_bubble123_ds3)
2896  + 32.0*d2_quartic_bubble_ds3;
2897  d2psids(7,4) = 0.0
2898  - 12.0*(d2_cubic_bubble012_ds4 + d2_cubic_bubble123_ds4)
2899  + 32.0*d2_quartic_bubble_ds4;
2900  d2psids(7,5) = 4.0
2901  - 12.0*(d2_cubic_bubble012_ds5 + d2_cubic_bubble123_ds5)
2902  + 32.0*d2_quartic_bubble_ds5;
2903 
2904  d2psids(8,0) = 0.0
2905  - 12.0*(d2_cubic_bubble023_ds0 + d2_cubic_bubble123_ds0)
2906  + 32.0*d2_quartic_bubble_ds0;
2907  d2psids(8,1) = 0.0
2908  - 12.0*(d2_cubic_bubble023_ds1 + d2_cubic_bubble123_ds1)
2909  + 32.0*d2_quartic_bubble_ds1;
2910  d2psids(8,2) = -8.0
2911  - 12.0*(d2_cubic_bubble023_ds2 + d2_cubic_bubble123_ds2)
2912  + 32.0*d2_quartic_bubble_ds2;
2913  d2psids(8,3) = 0.0
2914  - 12.0*(d2_cubic_bubble023_ds3 + d2_cubic_bubble123_ds3)
2915  + 32.0*d2_quartic_bubble_ds3;
2916  d2psids(8,4) = -4.0
2917  - 12.0*(d2_cubic_bubble023_ds4 + d2_cubic_bubble123_ds4)
2918  + 32.0*d2_quartic_bubble_ds4;
2919  d2psids(8,5) = -4.0
2920  - 12.0*(d2_cubic_bubble023_ds5 + d2_cubic_bubble123_ds5)
2921  + 32.0*d2_quartic_bubble_ds5;
2922 
2923  d2psids(9,0) = 0.0
2924  - 12.0*(d2_cubic_bubble013_ds0 + d2_cubic_bubble123_ds0)
2925  + 32.0*d2_quartic_bubble_ds0;
2926  d2psids(9,1) = -8.0
2927  - 12.0*(d2_cubic_bubble013_ds1 + d2_cubic_bubble123_ds1)
2928  + 32.0*d2_quartic_bubble_ds1;
2929  d2psids(9,2) = 0.0
2930  - 12.0*(d2_cubic_bubble013_ds2 + d2_cubic_bubble123_ds2)
2931  + 32.0*d2_quartic_bubble_ds3;
2932  d2psids(9,3) = -4.0
2933  - 12.0*(d2_cubic_bubble013_ds3 + d2_cubic_bubble123_ds3)
2934  + 32.0*d2_quartic_bubble_ds3;
2935  d2psids(9,4) = 0.0
2936  - 12.0*(d2_cubic_bubble013_ds4 + d2_cubic_bubble123_ds4)
2937  + 32.0*d2_quartic_bubble_ds4;
2938  d2psids(9,5) = -4.0
2939  - 12.0*(d2_cubic_bubble013_ds5 + d2_cubic_bubble123_ds5)
2940  + 32.0*d2_quartic_bubble_ds5;
2941 
2942  //Add the bubble function on the face spanned by 0 1 3
2943  d2psids(10,0) = 27.0*d2_cubic_bubble013_ds0 - 108.0*d2_quartic_bubble_ds0;
2944  d2psids(10,1) = 27.0*d2_cubic_bubble013_ds1 - 108.0*d2_quartic_bubble_ds1;
2945  d2psids(10,2) = 27.0*d2_cubic_bubble013_ds2 - 108.0*d2_quartic_bubble_ds2;
2946  d2psids(10,3) = 27.0*d2_cubic_bubble013_ds3 - 108.0*d2_quartic_bubble_ds3;
2947  d2psids(10,4) = 27.0*d2_cubic_bubble013_ds4 - 108.0*d2_quartic_bubble_ds4;
2948  d2psids(10,5) = 27.0*d2_cubic_bubble013_ds5 - 108.0*d2_quartic_bubble_ds5;
2949 
2950  //Add the bubble function on the face spanned by 0 1 2
2951  d2psids(11,0) = 27.0*d2_cubic_bubble012_ds0 - 108.0*d2_quartic_bubble_ds0;
2952  d2psids(11,1) = 27.0*d2_cubic_bubble012_ds1 - 108.0*d2_quartic_bubble_ds1;
2953  d2psids(11,2) = 27.0*d2_cubic_bubble012_ds2 - 108.0*d2_quartic_bubble_ds2;
2954  d2psids(11,3) = 27.0*d2_cubic_bubble012_ds3 - 108.0*d2_quartic_bubble_ds3;
2955  d2psids(11,4) = 27.0*d2_cubic_bubble012_ds4 - 108.0*d2_quartic_bubble_ds4;
2956  d2psids(11,5) = 27.0*d2_cubic_bubble012_ds5 - 108.0*d2_quartic_bubble_ds5;
2957 
2958  //Add the bubble function on the face spanned by 0 2 3
2959  d2psids(12,0) = 27.0*d2_cubic_bubble023_ds0 - 108.0*d2_quartic_bubble_ds0;
2960  d2psids(12,1) = 27.0*d2_cubic_bubble023_ds1 - 108.0*d2_quartic_bubble_ds1;
2961  d2psids(12,2) = 27.0*d2_cubic_bubble023_ds2 - 108.0*d2_quartic_bubble_ds2;
2962  d2psids(12,3) = 27.0*d2_cubic_bubble023_ds3 - 108.0*d2_quartic_bubble_ds3;
2963  d2psids(12,4) = 27.0*d2_cubic_bubble023_ds4 - 108.0*d2_quartic_bubble_ds4;
2964  d2psids(12,5) = 27.0*d2_cubic_bubble023_ds5 - 108.0*d2_quartic_bubble_ds5;
2965 
2966  //Add the bubble function on the face spanned by 1 2 3
2967  d2psids(13,0) = 27.0*d2_cubic_bubble123_ds0 - 108.0*d2_quartic_bubble_ds0;
2968  d2psids(13,1) = 27.0*d2_cubic_bubble123_ds1 - 108.0*d2_quartic_bubble_ds1;
2969  d2psids(13,2) = 27.0*d2_cubic_bubble123_ds2 - 108.0*d2_quartic_bubble_ds2;
2970  d2psids(13,3) = 27.0*d2_cubic_bubble123_ds3 - 108.0*d2_quartic_bubble_ds3;
2971  d2psids(13,4) = 27.0*d2_cubic_bubble123_ds4 - 108.0*d2_quartic_bubble_ds4;
2972  d2psids(13,5) = 27.0*d2_cubic_bubble123_ds5 - 108.0*d2_quartic_bubble_ds5;
2973 
2974  //Add the volumetric bubble function derivatives
2975  d2psids(14,0) = 256.0*d2_quartic_bubble_ds0;
2976  d2psids(14,1) = 256.0*d2_quartic_bubble_ds1;
2977  d2psids(14,2) = 256.0*d2_quartic_bubble_ds2;
2978  d2psids(14,3) = 256.0*d2_quartic_bubble_ds3;
2979  d2psids(14,4) = 256.0*d2_quartic_bubble_ds4;
2980  d2psids(14,5) = 256.0*d2_quartic_bubble_ds5;
2981 }
2982 
2983 
2984 };
2985 
2986 
2987 
2988 
2989 //=======================================================================
2990 /// General TElement class specialised to three spatial dimensions (tet)
2991 /// Ordering of nodes inverted from Zienkiewizc sketches: When looking into the
2992 /// tet from vertex node 0. The vertex nodes on the opposite face are
2993 /// 1 - 2 - 3 in anticlockwise direction. Other nodes filled in edge by
2994 /// edge, then the face ones, then the internal ones.
2995 //=======================================================================
2996 template<unsigned NNODE_1D>
2997 class TElement<3,NNODE_1D> : public virtual TElementBase,
2998  public TElementShape<3,NNODE_1D>
2999 {
3000  private:
3001 
3002  /// Nodal translation scheme for use when generating face elements
3003  static const unsigned Node_on_face[4][(NNODE_1D*(NNODE_1D+1))/2];
3004 
3005  /// \short Default integration rule: Gaussian integration of same 'order' as
3006  /// the element
3007  //This is sort of optimal, because it means that the integration is exact
3008  //for the shape functions. Can overwrite this in specific element defintion.
3010 
3011 public:
3012 
3013  /// Constructor
3015  {
3016  switch (NNODE_1D)
3017  {
3018  case 2:
3019  case 3:
3020  break;
3021 
3022 /* case 4: */
3023 /* n_node = 20; */
3024 /* break; */
3025 
3026  default:
3027  std::string error_message =
3028  "Tets are currently only implemented for\n";
3029  error_message +=
3030  "four and ten nodes, i.e. NNODE_1D=2 , 3 \n";
3031 
3032  throw OomphLibError(error_message,
3033  OOMPH_CURRENT_FUNCTION,
3034  OOMPH_EXCEPTION_LOCATION);
3035  }
3036 
3037 
3038  // Set the number of nodes
3039  unsigned n_node = (NNODE_1D*(NNODE_1D+1))/2 + 1 + 3*(NNODE_1D-2);
3040  this->set_n_node(n_node);
3041 
3042  // Set the elemental and nodal dimensions
3043  set_dimension(3);
3044 
3045  //Assign default (full) spatial integration scheme
3046  set_integration_scheme(&Default_integration_scheme);
3047  }
3048 
3049 
3050 
3051  /// Broken copy constructor
3053  {
3054  BrokenCopy::broken_copy("TElement");
3055  }
3056 
3057  /// Broken assignment operator
3058  /*void operator=(const TElement&)
3059  {
3060  BrokenCopy::broken_assign("TElement");
3061  }*/
3062 
3063 
3064  /// Destructor
3066 
3067  /// Number of nodes along each element edge
3068  unsigned nnode_1d() const {return NNODE_1D;}
3069 
3070 
3071  /// \short Number of vertex nodes in the element: One more
3072  /// than spatial dimension
3073  unsigned nvertex_node() const {return 4;}
3074 
3075  /// \short Public access function for Node_on_face.
3076  unsigned get_bulk_node_number(const int& face_index,
3077  const unsigned& i) const
3078  {
3079  return Node_on_face[face_index][i];
3080  }
3081 
3082  /// \short Pointer to the j-th vertex node in the element
3083  Node* vertex_node_pt(const unsigned& j) const
3084  {
3085  // Vertex nodes come first:
3086 #ifdef PARANOID
3087  if (j>3)
3088  {
3089  std::ostringstream error_message;
3090  error_message
3091  << "Element only has four vertex nodes; called with node number "
3092  << j << std::endl;
3093  throw OomphLibError(error_message.str(),
3094  OOMPH_CURRENT_FUNCTION,
3095  OOMPH_EXCEPTION_LOCATION);
3096  }
3097 #endif
3098  return node_pt(j);
3099  }
3100 
3101  /// Calculate the geometric shape functions at local coordinate s
3102  inline void shape(const Vector<double> &s, Shape &psi) const
3104 
3105  /// \short Compute the geometric shape functions and
3106  /// derivatives w.r.t. local coordinates at local coordinate s
3107  inline void dshape_local(const Vector<double> &s, Shape &psi,
3108  DShape &dpsids) const
3110 
3111  /// \short Compute the geometric shape functions, derivatives and
3112  /// second derivatives w.r.t local coordinates at local coordinate s.
3113  /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
3114  /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
3115  /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
3116  /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
3117  /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
3118  /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
3119  inline void d2shape_local(const Vector<double> &s, Shape &psi,
3120  DShape &dpsids, DShape &d2psids) const
3121  {TElementShape<3,NNODE_1D>::d2shape_local(s,psi,dpsids,d2psids);}
3122 
3123  /// \short Overload the template-free interface for the calculation of
3124  /// inverse jacobian matrix. This is a three dimensional element, so use
3125  /// the 3D version.
3127  DenseMatrix<double>&inverse_jacobian) const
3128  {return FiniteElement::invert_jacobian<3>(jacobian,inverse_jacobian);}
3129 
3130  /// Min. value of local coordinate
3131  double s_min() const {return 0.0;}
3132 
3133  /// Max. value of local coordinate
3134  double s_max() const {return 1.0;}
3135 
3136  /// Return local coordinates of node j
3137  inline void local_coordinate_of_node(const unsigned& j,Vector<double>& s) const
3139 
3140  /// \short Return the number of actual plot points for paraview
3141  /// plot with parameter nplot.
3142  unsigned nplot_points_paraview(const unsigned& nplot) const
3143  {
3144  unsigned node_sum=0;
3145  for(unsigned j=1;j<=nplot;j++)
3146  {
3147  for(unsigned i=1;i<=j;i++)
3148  {
3149  node_sum+=i;
3150  }
3151  }
3152  return node_sum;
3153  }
3154 
3155  /// \short Return the number of local sub-elements for paraview plot with
3156  /// parameter nplot.
3157  unsigned nsub_elements_paraview(const unsigned& nplot) const
3158  {
3159  return (nplot-1)*(nplot-1)*(nplot-1);
3160  }
3161 
3162  /// \short Fill in the offset information for paraview plot.
3163  /// Needs to be implemented for each new geometric element type; see
3164  /// http://www.vtk.org/VTK/img/file-formats.pdf
3165  void write_paraview_output_offset_information(std::ofstream& file_out,
3166  const unsigned& nplot,
3167  unsigned& counter) const
3168  {
3169  //Output node lists for sub elements for Paraview (node index
3170  //must start at 0, fixed with magical counter-1)
3171 
3172  unsigned paraview_fix=counter-1;
3173  unsigned nod_count=1;
3174 
3175  for(unsigned i=0;i<nplot;i++)
3176  {
3177  for(unsigned j=0;j<nplot-i;j++)
3178  {
3179  for(unsigned k=0;k<nplot-i-j;k++)
3180  {
3181  if(k<nplot-i-j-1)
3182  {
3183  file_out<< nod_count+paraview_fix << " "
3184  << nod_count+1+paraview_fix << " "
3185  << nod_count+nplot-i-j+paraview_fix <<" "
3186  <<nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2)
3187  +paraview_fix
3188  << std::endl;
3189  if(k<nplot-i-j-2)
3190  {
3191  file_out << nod_count+1+paraview_fix << " "
3192  << nod_count+nplot-i-j+paraview_fix << " "
3193  << nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2)
3194  +paraview_fix << " "
3195  << nod_count+2*(nplot-i-j)-1+((nplot-1-i)*(nplot-i)/2)
3196  +paraview_fix
3197  << std::endl;
3198  file_out << nod_count+1+paraview_fix << " "
3199  << nod_count+nplot-i-j+paraview_fix << " "
3200  << nod_count+nplot-i-j+1+paraview_fix << " "
3201  << nod_count+2*(nplot-i-j)-1+((nplot-1-i)*(nplot-i)/2)
3202  +paraview_fix
3203  << std::endl;
3204  file_out << nod_count+1+paraview_fix << " "
3205  << nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2)
3206  +paraview_fix << " "
3207  << nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2)+1
3208  +paraview_fix << " "
3209  << nod_count+2*(nplot-i-j)-1+((nplot-1-i)*(nplot-i)/2)
3210  +paraview_fix
3211  << std::endl;
3212  file_out << nod_count+1+paraview_fix << " "
3213  << nod_count+nplot-i-j+1+paraview_fix << " "
3214  << nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2)+1
3215  +paraview_fix
3216  << " "
3217  << nod_count+2*(nplot-i-j)-1+((nplot-1-i)*(nplot-i)/2)
3218  +paraview_fix
3219  << std::endl;
3220  }
3221  if(k>1)
3222  {
3223  file_out << nod_count+nplot-i-j-1+paraview_fix << " "
3224  << nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2)-1
3225  +paraview_fix
3226  << " "
3227  << nod_count+2*(nplot-i-j-1)+((nplot-1-i)*(nplot-i)/2)-1
3228  +paraview_fix
3229  <<" "
3230  << nod_count+2*(nplot-i-j-1)+((nplot-1-i)*(nplot-i)/2)
3231  +paraview_fix
3232  << std::endl;
3233  }
3234  }//end if(k<nplot-i-j-1)
3235  ++nod_count;
3236  }
3237  }
3238  }
3239 
3240  //increment the counter to keep track of global connectivity
3241  counter+=nplot_points_paraview(nplot);
3242 
3243  } //end of write Paraview_element...
3244 
3245  /// \short Return the paraview element type.
3246  /// Needs to be implemented for each new geometric element type; see
3247  /// http://www.vtk.org/VTK/img/file-formats.pdf
3248  void write_paraview_type(std::ofstream& file_out,
3249  const unsigned& nplot) const
3250  {
3251  unsigned local_loop=nsub_elements_paraview(nplot);
3252  for(unsigned i=0;i<local_loop;i++)
3253  {
3254  file_out << "10" << std::endl;
3255  }
3256  }
3257 
3258  /// \short Return the offsets for the paraview sub-elements. Needs
3259  /// to be implemented for each new geometric element type; see
3260  /// http://www.vtk.org/VTK/img/file-formats.pdf
3261  void write_paraview_offsets(std::ofstream& file_out,
3262  const unsigned& nplot,
3263  unsigned& offset_sum) const
3264  {
3265  unsigned local_loop=nsub_elements_paraview(nplot);
3266  for(unsigned i=0;i<local_loop;i++)
3267  {
3268  offset_sum+=4;
3269  file_out << offset_sum << std::endl;
3270  }
3271  }
3272 
3273  /// Output
3274  void output(std::ostream &output);
3275 
3276  /// Output at specified number of plot points
3277  void output(std::ostream &outfile, const unsigned &nplot);
3278 
3279  /// C-style output
3280  void output(FILE* file_pt);
3281 
3282  /// C_style output at n_plot points
3283  void output(FILE* file_pt, const unsigned &n_plot);
3284 
3285  /// \short Get vector of local coordinates of plot point i (when plotting
3286  /// nplot points in each "coordinate direction).
3287  void get_s_plot(const unsigned& iplot, const unsigned& nplot,
3288  Vector<double>& s) const
3289  {
3290  if (nplot>1)
3291  {
3292  unsigned np=0;
3293  for(unsigned i=0;i<nplot;++i)
3294  {
3295  for(unsigned j=0;j<nplot-i;++j)
3296  {
3297  for(unsigned k=0;k<nplot-i-j;++k)
3298  {
3299  if(np==iplot)
3300  {
3301  {
3302  s[0] = double(j)/double(nplot-1);
3303  s[1] = double(i)/double(nplot-1);
3304  s[2] = double(k)/double(nplot-1);
3305  return;
3306  }
3307  }
3308  np++;
3309  }
3310  }
3311  }
3312  }
3313  else
3314  {
3315  s[0] = 1.0/4.0;
3316  s[1] = 1.0/4.0;
3317  s[2] = 1.0/4.0;
3318  }
3319  }
3320 
3321  /// \short Return string for tecplot zone header (when plotting
3322  /// nplot points in each "coordinate direction)
3323  std::string tecplot_zone_string(const unsigned& nplot) const
3324  {
3325  std::ostringstream header;
3326  unsigned nel=0;
3327  nel=(nplot-1)*(nplot-1)*(nplot-1);
3328  header << "ZONE N=" << nplot_points(nplot) << ", E="
3329  << nel << ", F=FEPOINT, ET=TETRAHEDRON\n";
3330  return header.str();
3331  }
3332 
3333  /// \short Add tecplot zone "footer" to output stream (when plotting
3334  /// nplot points in each "coordinate direction).
3335  /// Empty by default -- can be used, e.g., to add FE connectivity
3336  /// lists to elements that need it.
3337  void write_tecplot_zone_footer(std::ostream& outfile,
3338  const unsigned& nplot) const
3339  {
3340 
3341  //Output node lists for sub elements for Tecplot (node index
3342  //must start at 1)
3343  unsigned nod_count=1;
3344  for(unsigned i=0;i<nplot;i++)
3345  {
3346  for(unsigned j=0;j<nplot-i;j++)
3347  {
3348  for(unsigned k=0;k<nplot-i-j;k++)
3349  {
3350  if(k<nplot-i-j-1)
3351  {
3352  outfile<< nod_count << " "
3353  << nod_count+1 << " "
3354  << nod_count+nplot-i-j <<" "
3355  <<nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2)
3356  << std::endl;
3357  if(k<nplot-i-j-2)
3358  {
3359  outfile << nod_count+1<< " "
3360  << nod_count+nplot-i-j << " "
3361  << nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2) << " "
3362  << nod_count+2*(nplot-i-j)-1+((nplot-1-i)*(nplot-i)/2)
3363  << std::endl;
3364  outfile << nod_count+1 << " "
3365  << nod_count+nplot-i-j << " "
3366  << nod_count+nplot-i-j+1<< " "
3367  << nod_count+2*(nplot-i-j)-1+((nplot-1-i)*(nplot-i)/2)
3368  << std::endl;
3369  outfile << nod_count+1<< " "
3370  << nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2)<< " "
3371  << nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2)+1<< " "
3372  << nod_count+2*(nplot-i-j)-1+((nplot-1-i)*(nplot-i)/2)
3373  << std::endl;
3374  outfile << nod_count+1<< " "
3375  << nod_count+nplot-i-j+1<< " "
3376  << nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2)+1
3377  << " "
3378  << nod_count+2*(nplot-i-j)-1+((nplot-1-i)*(nplot-i)/2)
3379  << std::endl;
3380  }
3381  if(k>1)
3382  {
3383  outfile << nod_count+nplot-i-j-1<< " "
3384  << nod_count+nplot-i-j+((nplot-1-i)*(nplot-i)/2)-1
3385  << " "
3386  << nod_count+2*(nplot-i-j-1)+((nplot-1-i)*(nplot-i)/2)-1
3387  <<" "
3388  << nod_count+2*(nplot-i-j-1)+((nplot-1-i)*(nplot-i)/2)
3389  << std::endl;
3390  }
3391  }//end if(k<nplot-i-j-1)
3392  ++nod_count;
3393  }
3394  }
3395  }
3396  } //end of write tecplot...
3397 
3398 
3399  /// \short Add tecplot zone "footer" to C-style output. (when plotting
3400  /// nplot points in each "coordinate direction).
3401  /// Empty by default -- can be used, e.g., to add FE connectivity
3402  /// lists to elements that need it.
3403  void write_tecplot_zone_footer(FILE* file_pt,
3404  const unsigned& nplot) const
3405  {
3406  //Output node lists for sub elements for Tecplot (node index
3407  //must start at 1)
3408  unsigned nod_count=1;
3409  for(unsigned i=0;i<nplot;i++)
3410  {
3411  for(unsigned j=0;j<nplot-i;j++)
3412  {
3413  for(unsigned k=0;k<nplot-i-j;k++)
3414  {
3415  if(j<nplot-i-1)
3416  {
3417  fprintf(file_pt,"%i %i %i \n",nod_count,
3418  nod_count+1,nod_count+nplot-i);
3419  if(j<nplot-i-2)
3420  {
3421  fprintf(file_pt,"%i %i %i \n",nod_count+1,nod_count+nplot-i+1,
3422  nod_count+nplot-i);
3423  }
3424  }
3425  ++nod_count;
3426  }
3427  }
3428  }
3429  }
3430 
3431  /// \short Return total number of plot points (when plotting
3432  /// nplot points in each "coordinate direction)
3433  unsigned nplot_points(const unsigned& nplot) const
3434  {
3435  unsigned res=0;
3436  if(nplot>1)
3437  {
3438  res=4;
3439  for(unsigned i=3;i<=nplot;i++)
3440  {
3441  res+=(i*(i+1)/2);
3442  }
3443  return res;
3444  }
3445  //Otherwise we return 1(?)
3446  return 1;
3447  }
3448 
3449  /// \short Build the lower-dimensional FaceElement (an element of type
3450  /// TElement<2,NNODE_1D>). The face index can take one of four values
3451  /// corresponding to the four possible faces:
3452  /// 0: (left) s[0] = 0.0
3453  /// 1: (bottom) s[1] = 0.0
3454  /// 2: (back) s[2] = 0.0
3455  /// 3: (sloping face) s[0] + s[1] + s[2] = 1
3456  void build_face_element(const int &face_index,
3457  FaceElement* face_element_pt);
3458 
3459 };
3460 
3461 
3462 ///////////////////////////////////////////////////////////////////////
3463 ///////////////////////////////////////////////////////////////////////
3464 ///////////////////////////////////////////////////////////////////////
3465 
3466 
3467 
3468 //=======================================================================
3469 /// TElement class for which the shape functions have been enriched
3470 /// by a single bubble function of the next order
3471 ///
3472 /// Empty, just establishes the template parameters
3473 //=======================================================================
3474 template<unsigned DIM, unsigned NNODE_1D>
3476 {
3477 };
3478 
3479 //=====================================================================
3480 /// Define integration schemes that are required to exactly integrate
3481 /// the mass matrices of the bubble-enriched elements. The enrichement
3482 /// increases the polynomial order which means that higher-order Gauss
3483 /// rules must be used.
3484 //====================================================================
3485 template<unsigned DIM, unsigned NPTS_1D>
3487 {
3488 };
3489 
3490 //====================================================================
3491 ///Specialisation for two-dimensional elements, in which the highest
3492 ///order polynomial is cubic, so we need the integration scheme
3493 ///for the unenriched cubic element
3494 //======================================================================
3495 template<>
3496 class TBubbleEnrichedGauss<2,3> : public TGauss<2,4>
3497 {
3498  public:
3500 };
3501 
3502 //====================================================================
3503 ///Specialisation for three-dimensional elements, in which the highest
3504 ///order polynomial is quartic, so we need the integration scheme
3505 ///for the unenriched quartic element
3506 //======================================================================
3507 template<>
3508 class TBubbleEnrichedGauss<3,3> : public TGauss<3,5>
3509 {
3510  public:
3512 };
3513 
3514 
3515 
3516 //=======================================================================
3517 /// Enriched TElement class specialised to two spatial dimensions
3518 /// and three nodes per side (quadratic element)
3519 /// Ordering of nodes as in Zienkiwizc sketches: vertex nodes
3520 /// 0 - 1 - 2 anticlockwise. Midside nodes filled in progressing
3521 /// along the consecutive edges. Central node(s) come(s) last.
3522 /// The idea is that we inherit from the existing TElement<2,3>, add
3523 /// the single extra node at the centroid and
3524 /// overload the shape functions to be those corresponding to the
3525 /// enriched element.
3526 //=======================================================================
3527 template<unsigned DIM>
3528 class TBubbleEnrichedElement<DIM,3> : public virtual TElement<DIM,3>,
3529  public TBubbleEnrichedElementShape<DIM,3>
3530 {
3531  private:
3532 
3533  //Static storage for a new integration scheme
3535 
3536  //Static storage for central node
3537  static const unsigned Central_node_on_face[DIM+1];
3538 
3539 public:
3540 
3541  ///Constructor
3544  {
3545  //Add the additional enrichment nodes
3546  unsigned n_node = this->nnode();
3547  this->set_n_node(n_node+
3549  //Set the new integration scheme
3550  this->set_integration_scheme(&Default_enriched_integration_scheme);
3551  }
3552 
3553  /// Broken copy constructor
3555  {
3556  BrokenCopy::broken_copy("TBubbleEnrichedElement");
3557  }
3558 
3559  /// Broken assignment operator
3560  /*void operator=(const TBubbleEnrichedElement&)
3561  {
3562  BrokenCopy::broken_assign("TBubbleEnrichedElement");
3563  }*/
3564 
3565  /// Destructor
3567 
3568  /// Calculate the geometric shape functions at local coordinate s
3569  inline void shape(const Vector<double> &s, Shape &psi) const
3571 
3572  /// \short Compute the geometric shape functions and
3573  /// derivatives w.r.t. local coordinates at local coordinate s
3574  inline void dshape_local(const Vector<double> &s, Shape &psi,
3575  DShape &dpsids) const
3577 
3578 
3579  /// \short Compute the geometric shape functions, derivatives and
3580  /// second derivatives w.r.t local coordinates at local coordinate s
3581  /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
3582  /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
3583  /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
3584  inline void d2shape_local(const Vector<double> &s, Shape &psi,
3585  DShape &dpsids, DShape &d2psids) const
3586  {TBubbleEnrichedElementShape<DIM,3>::d2shape_local(s,psi,dpsids,d2psids);}
3587 
3588  /// Return local coordinates of node j
3589  inline void local_coordinate_of_node(const unsigned& j,Vector<double>& s) const
3591 
3592  /// \short Build the lower-dimensional FaceElement
3593  void build_face_element(const int &face_index,
3594  FaceElement* face_element_pt);
3595 
3596 };
3597 
3598 
3599 //========================================================================
3600 /// Base class for Solid Telements
3601 //========================================================================
3602 class TSolidElementBase : public virtual TElementBase,
3603  public virtual SolidFiniteElement
3604 {
3605 
3606 
3607  public:
3608 
3609  /// Constructor: Empty
3611 
3612  /// Broken copy constructor
3614  {
3615  BrokenCopy::broken_copy("TSolidElementBase");
3616  }
3617 
3618  /// Broken assignment operator
3619  /*void operator=(const TSolidElementBase&)
3620  {
3621  BrokenCopy::broken_assign("TSolidElementBase");
3622  }*/
3623 
3624 };
3625 
3626 
3627 
3628 
3629 //////////////////////////////////////////////////////////////////////////
3630 //////////////////////////////////////////////////////////////////////////
3631 //////////////////////////////////////////////////////////////////////////
3632 
3633 
3634 //=======================================================================
3635 /// SolidTElement elements are triangular/tet elements whose
3636 /// derivatives also include those based upon the lagrangian
3637 /// positions of the nodes.
3638 /// They are the basis for solid mechanics elements.
3639 //=======================================================================
3640 template <unsigned DIM, unsigned NNODE_1D>
3642 {};
3643 
3644 
3645 //=======================================================================
3646 ///SolidTElement elements, specialised to one spatial dimension
3647 //=======================================================================
3648 template <unsigned NNODE_1D>
3649 class SolidTElement<1,NNODE_1D> : public virtual TElement<1,NNODE_1D>,
3650  public virtual TSolidElementBase
3651 {
3652  public:
3653 
3654  /// Constructor
3656  {
3657  //Set the Lagrangian dimension of the element
3658  set_lagrangian_dimension(1);
3659  }
3660 
3661  /// Broken copy constructor
3663  {
3664  BrokenCopy::broken_copy("SolidTElement");
3665  }
3666 
3667  /// Broken assignment operator
3668  /*void operator=(const SolidTElement&)
3669  {
3670  BrokenCopy::broken_assign("SolidTElement");
3671  }*/
3672 
3673  /// \short Build the lower-dimensional FaceElement (an element of type
3674  /// SolidPointElement). The face index takes two values
3675  /// corresponding to the two possible faces:
3676  /// -1 (Left) s[0] = -1.0
3677  /// +1 (Right) s[0] = 1.0
3678  inline void build_face_element(const int &face_index,
3679  FaceElement* face_element_pt);
3680 
3681 
3682 };
3683 
3684 
3685 ///////////////////////////////////////////////////////////////////////////
3686 ///////////////////////////////////////////////////////////////////////////
3687 // SolidTElements
3688 ///////////////////////////////////////////////////////////////////////////
3689 ///////////////////////////////////////////////////////////////////////////
3690 
3691 
3692 ///////////////////////////////////////////////////////////////////////////
3693 // 1D SolidTElements
3694 ///////////////////////////////////////////////////////////////////////////
3695 
3696 
3697 //===========================================================
3698 /// Function to setup geometrical information for lower-dimensional
3699 /// FaceElements (which are of type SolidTElement<0,1>).
3700 //===========================================================
3701 template<unsigned NNODE_1D>
3703 build_face_element(const int &face_index,
3704  FaceElement *face_element_pt)
3705 {
3706  //Build the standard non-solid FaceElement
3707  TElement<1,NNODE_1D>::build_face_element(face_index,face_element_pt);
3708 
3709  //Set the Lagrangian dimension from the first node of the present element
3710  dynamic_cast<SolidFiniteElement*>(face_element_pt)->
3711  set_lagrangian_dimension(static_cast<SolidNode*>(node_pt(0))->nlagrangian());
3712 }
3713 
3714 
3715 
3716 //=======================================================================
3717 /// SolidTElement elements, specialised to two spatial dimensions
3718 //=======================================================================
3719 template <unsigned NNODE_1D>
3720 class SolidTElement<2,NNODE_1D> : public virtual TElement<2,NNODE_1D>,
3721  public virtual TSolidElementBase
3722 {
3723  public:
3724 
3725  /// Constructor
3726  SolidTElement() : TElementBase(), TElement<2,NNODE_1D>(),
3728  {
3729  //Set the Lagrangian dimension of the element
3730  set_lagrangian_dimension(2);
3731  }
3732 
3733  /// Broken copy constructor
3735  {
3736  BrokenCopy::broken_copy("SolidTElement");
3737  }
3738 
3739  /// Broken assignment operator
3740  /*void operator=(const SolidTElement&)
3741  {
3742  BrokenCopy::broken_assign("SolidTElement");
3743  }*/
3744 
3745  /// \short Build the lower-dimensional FaceElement (an element of type
3746  /// SolidTElement<1,NNODE_1D>). The face index takes three possible values:
3747  /// 0 (Left) s[0] = 0.0
3748  /// 1 (Bottom) s[1] = 0.0
3749  /// 2 (Sloping face) s[0] = 1.0 - s[1]
3750  inline void build_face_element(const int &face_index,
3751  FaceElement* face_element_pt);
3752 
3753 };
3754 
3755 
3756 ///////////////////////////////////////////////////////////////////////////
3757 ///////////////////////////////////////////////////////////////////////////
3758 // 2D SolidTElements
3759 ///////////////////////////////////////////////////////////////////////////
3760 ///////////////////////////////////////////////////////////////////////////
3761 
3762 
3763 //===========================================================
3764 /// Function to setup geometrical information for lower-dimensional
3765 /// FaceElements (which are of type SolidTElement<1,NNODE_1D>).
3766 //===========================================================
3767 template<unsigned NNODE_1D>
3769 build_face_element(const int &face_index, FaceElement *face_element_pt)
3770 {
3771  //Build the standard non-solid FaceElement
3772  TElement<2,NNODE_1D>::build_face_element(face_index,face_element_pt);
3773 
3774  //Set the Lagrangian dimension from the first node of the present element
3775  dynamic_cast<SolidFiniteElement*>(face_element_pt)->
3776  set_lagrangian_dimension(static_cast<SolidNode*>(node_pt(0))->nlagrangian());
3777 }
3778 
3779 
3780 
3781 //=======================================================================
3782 /// SolidTElement elements, specialised to three spatial dimensions
3783 //=======================================================================
3784 template <unsigned NNODE_1D>
3785 class SolidTElement<3,NNODE_1D> : public virtual TElement<3,NNODE_1D>,
3786  public virtual TSolidElementBase
3787 {
3788  public:
3789 
3790  /// Constructor
3791  SolidTElement() : TElementBase(), TElement<3,NNODE_1D>(),
3793  {
3794  //Set the Lagrangian dimension of the element
3795  set_lagrangian_dimension(3);
3796  }
3797 
3798  /// Broken copy constructor
3800  {
3801  BrokenCopy::broken_copy("SolidTElement");
3802  }
3803 
3804  /// Broken assignment operator
3805  /*void operator=(const SolidTElement&)
3806  {
3807  BrokenCopy::broken_assign("SolidTElement");
3808  }*/
3809 
3810 
3811  /// \short Build the lower-dimensional FaceElement (an element of type
3812  /// SolidTElement<2,NNODE_1D>). The face index can take one of four values
3813  /// corresponding to the four possible faces:
3814  /// 0: (left) s[0] = 0.0
3815  /// 1: (bottom) s[1] = 0.0
3816  /// 2: (back) s[2] = 0.0
3817  /// 3: (sloping face) s[0] + s[1] + s[2] = 1
3818  inline void build_face_element(const int &face_index,
3819  FaceElement* face_element_pt);
3820 
3821 };
3822 
3823 
3824 
3825 ///////////////////////////////////////////////////////////////////////////
3826 // 3D SolidTElements
3827 ///////////////////////////////////////////////////////////////////////////
3828 
3829 
3830 //===========================================================
3831 /// Function to setup geometrical information for lower-dimensional
3832 /// FaceElements (which are of type SolidTElement<1,NNODE_1D>).
3833 //===========================================================
3834 template<unsigned NNODE_1D>
3836  build_face_element(const int &face_index,
3837  FaceElement *face_element_pt)
3838 {
3839  //Build the standard non-solid FaceElement
3840  TElement<3,NNODE_1D>::build_face_element(face_index,face_element_pt);
3841 
3842  //Set the Lagrangian dimension from the first node of the present element
3843  dynamic_cast<SolidFiniteElement*>(face_element_pt)->
3844  set_lagrangian_dimension(static_cast<SolidNode*>(node_pt(0))->nlagrangian());
3845 }
3846 
3847 
3848 
3849 
3850 ///////////////////////////////////////////////////////////////////////
3851 ///////////////////////////////////////////////////////////////////////
3852 ///////////////////////////////////////////////////////////////////////
3853 
3854 //=======================================================================
3855 /// SolidTBubbleEnrichedElement elements are the enriched version
3856 /// of the SolidTElements. They will simply inherit from the appropriate
3857 /// SolidTElement and TBubblEnrichedElement.
3858 /// They are the basis for solid mechanics elements.
3859 //=======================================================================
3860 template <unsigned DIM, unsigned NNODE_1D>
3862 {};
3863 
3864 //===================================================================
3865 ///Specify the SolidTBubbleEnrichedElement corresponding to the
3866 ///quadratic triangle
3867 //===================================================================
3868 template<unsigned DIM>
3870 public virtual SolidTElement<DIM,3>,
3871 public virtual TBubbleEnrichedElement<DIM,3>
3872 {
3873 
3874 public:
3875 
3876  ///Constructor
3878  TBubbleEnrichedElement<DIM,3>() {}
3879 
3880  /// Broken copy constructor
3882  {
3883  BrokenCopy::broken_copy("SolidTBubbleEnrichedElement");
3884  }
3885 
3886  /// Broken assignment operator
3887  /*void operator=(const SolidTBubbleEnrichedElement&)
3888  {
3889  BrokenCopy::broken_assign("SolidTBubbleEnrichedElement");
3890  }*/
3891 
3892  /// Destructor
3894 
3895  /// \short Build the lower-dimensional FaceElement
3896  /// Need to put in a final override here
3897  void build_face_element(const int &face_index,
3898  FaceElement* face_element_pt);
3899 
3900 };
3901 
3902 
3903 //=======================================================================
3904 /// Face geometry for the TElement elements: The spatial
3905 /// dimension of the face elements is one lower than that of the
3906 /// bulk element but they have the same number of points
3907 /// along their 1D edges.
3908 //=======================================================================
3909 template<unsigned DIM, unsigned NNODE_1D>
3910 class FaceGeometry<TElement<DIM,NNODE_1D> >:
3911  public virtual TElement<DIM-1,NNODE_1D>
3912 {
3913 
3914  public:
3915 
3916  /// \short Constructor: Call the constructor for the
3917  /// appropriate lower-dimensional QElement
3918  FaceGeometry() : TElement<DIM-1,NNODE_1D>() {}
3919 
3920 };
3921 
3922 
3923 
3924 //=======================================================================
3925 /// Face geometry for the 1D TElement elements: Point elements
3926 //=======================================================================
3927 template<unsigned NNODE_1D>
3928 class FaceGeometry<TElement<1,NNODE_1D> >:
3929  public virtual PointElement
3930 {
3931 
3932  public:
3933 
3934  /// \short Constructor: Call the constructor for the
3935  /// appropriate lower-dimensional TElement
3937 
3938 };
3939 
3940 
3941 
3942 
3943 ///////////////////////////////////////////////////////////////////////
3944 ///////////////////////////////////////////////////////////////////////
3945 ///////////////////////////////////////////////////////////////////////
3946 
3947 
3948 //=======================================================================
3949 /// Face geometry for the 2D TBubbleEnrichedElement elements is exactly
3950 /// the same as for the corresponding TElement. The spatial
3951 /// dimension of the face elements is one lower than that of the
3952 /// bulk element but they have the same number of points
3953 /// along their 1D edges.
3954 //=======================================================================
3955 template<unsigned NNODE_1D>
3957  public virtual TElement<1,NNODE_1D>
3958 {
3959 
3960  public:
3961 
3962  /// \short Constructor: Call the constructor for the
3963  /// appropriate lower-dimensional QElement
3964  FaceGeometry() : TElement<1,NNODE_1D>() {}
3965 
3966 };
3967 
3968 
3969 //=======================================================================
3970 /// Face geometry for the 3D TBubbleEnrichedElement elements is the
3971 /// 2D TBubbleEnrichedElement. The spatial
3972 /// dimension of the face elements is one lower than that of the
3973 /// bulk element but they have the same number of points
3974 /// along their 1D edges.
3975 //=======================================================================
3976 template<unsigned NNODE_1D>
3978  public virtual TBubbleEnrichedElement<2,NNODE_1D>
3979 {
3980 
3981  public:
3982 
3983  /// \short Constructor: Call the constructor for the
3984  /// appropriate lower-dimensional QElement
3986 
3987 };
3988 
3989 
3990 ///////////////////////////////////////////////////////////////////////
3991 ///////////////////////////////////////////////////////////////////////
3992 ///////////////////////////////////////////////////////////////////////
3993 
3994 
3995 
3996 
3997 //=======================================================================
3998 /// Face geometry for the TElement elements: The spatial
3999 /// dimension of the face elements is one lower than that of the
4000 /// bulk element but they have the same number of points
4001 /// along their 1D edges.
4002 //=======================================================================
4003 template<unsigned DIM, unsigned NNODE_1D>
4004 class FaceGeometry<SolidTElement<DIM,NNODE_1D> >:
4005  public virtual SolidTElement<DIM-1,NNODE_1D>
4006 {
4007 
4008  public:
4009 
4010  /// \short Constructor: Call the constructor for the
4011  /// appropriate lower-dimensional QElement
4012  FaceGeometry() : SolidTElement<DIM-1,NNODE_1D>() {}
4013 
4014 };
4015 
4016 
4017 
4018 //=======================================================================
4019 /// Face geometry for the 1D TElement elements: Point elements
4020 //=======================================================================
4021 template<unsigned NNODE_1D>
4022 class FaceGeometry<SolidTElement<1,NNODE_1D> >:
4023  public virtual SolidPointElement
4024 {
4025 
4026  public:
4027 
4028  /// \short Constructor: Call the constructor for the
4029  /// appropriate lower-dimensional TElement
4031 
4032 };
4033 
4034 
4035 ///////////////////////////////////////////////////////////////////////
4036 ///////////////////////////////////////////////////////////////////////
4037 ///////////////////////////////////////////////////////////////////////
4038 
4039 
4040 //=======================================================================
4041 /// Face geometry for the 2D SolidTBubbleEnrichedElement elements is exactly
4042 /// the same as for the corresponding 2D SolidTElement. The spatial
4043 /// dimension of the face elements is one lower than that of the
4044 /// bulk element but they have the same number of points
4045 /// along their 1D edges.
4046 //=======================================================================
4047 template<unsigned NNODE_1D>
4049  public virtual SolidTElement<1,NNODE_1D>
4050 {
4051 
4052  public:
4053 
4054  /// \short Constructor: Call the constructor for the
4055  /// appropriate lower-dimensional QElement
4056  FaceGeometry() : SolidTElement<1,NNODE_1D>() {}
4057 
4058 };
4059 
4060 
4061 
4062 //=======================================================================
4063 /// Face geometry for the 3D SolidTBubbleEnrichedElement elements is
4064 /// the 2D SolidTBubbleEnrichedElement. The spatial
4065 /// dimension of the face elements is one lower than that of the
4066 /// bulk element but they have the same number of points
4067 /// along their 1D edges.
4068 //=======================================================================
4069 template<unsigned NNODE_1D>
4071  public virtual SolidTBubbleEnrichedElement<2,NNODE_1D>
4072 {
4073 
4074  public:
4075 
4076  /// \short Constructor: Call the constructor for the
4077  /// appropriate lower-dimensional QElement
4079 
4080 };
4081 
4082 }
4083 
4084 
4085 
4086 
4087 #endif
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:1504
Triangular Face class.
Definition: Telements.h:63
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition: Telements.h:404
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:1746
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition: Telements.h:3589
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<3,2>
Definition: Telements.h:1980
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
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:1714
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Definition: Telements.h:382
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:1679
Node * Node1_pt
First vertex node.
Definition: Telements.h:207
void shape(const Vector< double > &s, Shape &psi) const
Calculate the geometric shape functions at local coordinate s.
Definition: Telements.h:3102
static TBubbleEnrichedGauss< DIM, 3 > Default_enriched_integration_scheme
Definition: Telements.h:3534
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
Definition: Telements.h:4078
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:3584
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:1391
~TElement()
Broken assignment operator.
Definition: Telements.h:1622
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition: Telements.h:1409
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
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:1672
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:1462
Base class for Solid Telements.
Definition: Telements.h:3602
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:3165
TElementGeometricBase()
Empty default constructor.
Definition: Telements.h:1121
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition: Telements.h:321
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:3403
cstr elem_len * i
Definition: cfortran.h:607
unsigned nnode_1d() const
Number of nodes along each element edge.
Definition: Telements.h:1625
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:3323
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Definition: Telements.h:468
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Definition: Telements.h:2669
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:3337
bool operator<(const TFace &other) const
Less-than operator.
Definition: Telements.h:125
double s_max() const
Max. value of local coordinate.
Definition: Telements.h:1687
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Definition: Telements.h:1921
~SolidTBubbleEnrichedElement()
Broken assignment operator.
Definition: Telements.h:3893
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
void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get vector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: Telements.h:1489
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition: Telements.h:940
TSolidElementBase()
Constructor: Empty.
Definition: Telements.h:3610
double s_max() const
Max. value of local coordinate.
Definition: Telements.h:3134
TElement(const TElement &)
Broken copy constructor.
Definition: Telements.h:3052
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition: Telements.h:247
unsigned nvertex_node() const
Number of vertex nodes in the element: One more than spatial dimension.
Definition: Telements.h:3073
unsigned get_bulk_node_number(const int &face_index, const unsigned &i) const
Public access function for Node_on_face.
Definition: Telements.h:1632
A general Finite Element class.
Definition: elements.h:1271
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Definition: Telements.h:1062
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Definition: Telements.h:2319
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:1293
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:3126
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
SolidTElement(const SolidTElement &)
Broken copy constructor.
Definition: Telements.h:3662
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:3574
const unsigned Central_node_on_face[3]
Definition: Telements.cc:1060
bool local_coord_is_valid(Vector< double > &s, const double &rounding_tolerance)
Check whether the local coordinates are valid or not and allow for rounding tolerance. If the local coordinate specifed by s is "slightly" outside the element (as specified by rounding_tolerance) we move the point back into the element.
Definition: Telements.h:1210
TElementBase(const TElementBase &)
Broken copy constructor.
Definition: Telements.h:1161
TFace(Node *node1_pt, Node *node2_pt, Node *node3_pt)
Constructor: Pass in the three vertex nodes.
Definition: Telements.h:69
Node * node3_pt() const
Access to the third vertex node.
Definition: Telements.h:105
void shape(const Vector< double > &s, Shape &psi) const
Calculate the geometric shape functions at local coordinate s.
Definition: Telements.h:1377
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional TElement. ...
Definition: Telements.h:4030
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Definition: Telements.h:688
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:1818
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Definition: Telements.h:3083
TElement(const bool &allow_high_order)
Alternative constructor.
Definition: Telements.h:1585
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Definition: Telements.h:2045
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Definition: Telements.h:2015
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
unsigned nnode_1d() const
Number of nodes along each element edge.
Definition: Telements.h:3068
unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction) ...
Definition: Telements.h:3433
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<3,3>
Definition: Telements.h:2128
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:3248
static TGauss< 1, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Telements.h:1287
void get_s_plot(const unsigned &iplot, const unsigned &nplot, Vector< double > &s) const
Get vector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: Telements.h:3287
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:1449
unsigned nplot_points(const unsigned &nplot) const
Definition: Telements.h:1513
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition: Telements.h:497
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<3,2>
Definition: Telements.h:1968
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
Definition: Telements.h:3985
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:1287
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:3261
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
Definition: Telements.h:3964
double s_min() const
Min. value of local coordinate.
Definition: Telements.h:1403
TElementBase()
Empty default constructor.
Definition: Telements.h:1158
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<1,2>
Definition: Telements.h:276
unsigned nsub_elements_paraview(const unsigned &nplot) const
Return the number of local sub-elements for paraview plot with.
Definition: Telements.h:3157
unsigned nsub_elements_paraview(const unsigned &nplot) const
Return the number of local sub-elements for paraview plot with parameter nplot.
Definition: Telements.h:1704
void shape(const Vector< double > &s, Shape &psi) const
Calculate the geometric shape functions at local coordinate s.
Definition: Telements.h:1658
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition: Telements.h:730
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:1398
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
Definition: Telements.h:3918
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
unsigned nsub_elements_paraview(const unsigned &nplot) const
Return the number of local sub-elements for paraview plot with parameter nplot.
Definition: Telements.h:1421
TBubbleEnrichedElement(const TBubbleEnrichedElement &)
Broken copy constructor.
Definition: Telements.h:3554
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
Definition: Telements.h:4056
unsigned nplot_points(const unsigned &nplot) const
Definition: Telements.h:1889
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<1,3>
Definition: Telements.h:355
void shape(const Vector< double > &s, Shape &psi) const
Calculate the geometric shape functions at local coordinate s.
Definition: Telements.h:3569
Node * node2_pt() const
Access to the second vertex node.
Definition: Telements.h:102
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
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition: Telements.h:589
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
static char t char * s
Definition: cfortran.h:572
void build_face_element(const int &face_index, FaceElement *face_element_pt)
Definition: Telements.cc:1129
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:3119
SolidTElement(const SolidTElement &)
Broken copy constructor.
Definition: Telements.h:3734
TSolidElementBase(const TSolidElementBase &)
Broken copy constructor.
Definition: Telements.h:3613
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<3,3>
Definition: Telements.h:2500
static TGauss< 2, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Telements.h:1546
unsigned n_enriched_nodes()
Return the number of nodes required for enrichement.
Definition: Telements.h:935
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:1862
const unsigned Node_on_face[3][2]
Definition: Telements.cc:330
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition: Telements.h:3137
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:1382
void get_s_plot(const unsigned &iplot, const unsigned &nplot, Vector< double > &s) const
Get vector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: Telements.h:1789
double s_min() const
Min. value of local coordinate.
Definition: Telements.h:1684
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
Definition: Telements.h:4012
double s_min() const
Min. value of local coordinate.
Definition: Telements.h:3131
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
void output(std::ostream &outfile)
Output with default number of plot points.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional TElement. ...
Definition: Telements.h:3936
static TGauss< 3, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Telements.h:3009
bool local_coord_is_valid(const Vector< double > &s)
Check whether the local coordinates are valid or not.
Definition: Telements.h:1179
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Definition: Telements.h:301
TElement(const TElement &)
Broken copy constructor.
Definition: Telements.h:1609
unsigned nnode_1d() const
Number of nodes along each element edge.
Definition: Telements.h:1342
bool operator==(const TFace &other) const
Comparison operator.
Definition: Telements.h:108
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Definition: Telements.h:2209
~TBubbleEnrichedElement()
Broken assignment operator.
Definition: Telements.h:3566
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Definition: Telements.h:1639
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Return local coordinates of node j.
Definition: Telements.h:1690
Node * node1_pt() const
Access to the first vertex node.
Definition: Telements.h:99
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<2,3>
Definition: Telements.h:641
unsigned nvertex_node() const
Number of vertex nodes in the element: One more than spatial dimension.
Definition: Telements.h:1347
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TBubbleEnrichedElement<2,3>
Definition: Telements.h:999
ElementGeometry::ElementGeometry element_geometry() const
Broken assignment operator.
Definition: Telements.h:1173
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
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Definition: Telements.h:1350
~TElement()
Broken assignment operator.
Definition: Telements.h:1339
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot.
Definition: Telements.h:3142
unsigned get_bulk_node_number(const int &face_index, const unsigned &i) const
Public access function for Node_on_face.
Definition: Telements.h:3076
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TBubbleEnrichedElement<3,3>
Definition: Telements.h:2438
TElementGeometricBase(const TElementGeometricBase &)
Broken copy constructor.
Definition: Telements.h:1124
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:1761
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<2,2>
Definition: Telements.h:534
Node * Node3_pt
Third vertex node.
Definition: Telements.h:213
unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot.
Definition: Telements.h:1695
unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot.
Definition: Telements.h:1414
TElement(const TElement &)
Broken copy constructor.
Definition: Telements.h:1326
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Definition: Telements.h:859
unsigned nvertex_node() const
Number of vertex nodes in the element: One more than spatial dimension.
Definition: Telements.h:1629
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Derivatives of shape functions for specific TElement<3,3>
Definition: Telements.h:2147
void build_face_element(const int &face_index, FaceElement *face_element_pt)
Definition: Telements.cc:1067
SolidTBubbleEnrichedElement(const SolidTBubbleEnrichedElement &)
Broken copy constructor.
Definition: Telements.h:3881
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
void build_face_element(const int &face_index, FaceElement *face_element_pt)
Broken assignment operator.
Definition: Telements.h:3836
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:1663
Node * Node2_pt
Second vertex node.
Definition: Telements.h:210
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<1,4>
Definition: Telements.h:441
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
SolidFiniteElement class.
Definition: elements.h:3320
Solid point element.
Definition: elements.h:4685
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:1429
~TElement()
Broken assignment operator.
Definition: Telements.h:3065
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:3107
void shape(const Vector< double > &s, Shape &psi) const
Shape function for specific TElement<2,4>
Definition: Telements.h:802
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Definition: Telements.h:566
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:1832
void build_face_element(const int &face_index, FaceElement *face_element_pt)
Broken assignment operator.
Definition: Telements.h:3769
void build_face_element(const int &face_index, FaceElement *face_element_pt)
Broken assignment operator.
Definition: Telements.h:3703
double s_max() const
Max. value of local coordinate.
Definition: Telements.h:1406
SolidTElement(const SolidTElement &)
Broken copy constructor.
Definition: Telements.h:3799