Telements.cc
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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 //Non-inline member functions for Telements
31 
32 
33 // oomph-lib headers
34 #include "Telements.h"
35 
36 
37 namespace oomph
38 {
39 
40 //=======================================================================
41 /// Assign the static integral
42 //=======================================================================
43 template<unsigned NNODE_1D>
45 template<unsigned NNODE_1D>
47 template<unsigned NNODE_1D>
49 
50 //////////////////////////////////////////////////////////////////
51 // 1D Telements
52 //////////////////////////////////////////////////////////////////
53 
54 //=======================================================================
55 /// The output function for general 1D TElements
56 //=======================================================================
57 template <unsigned NNODE_1D>
58 void TElement<1,NNODE_1D>::output(std::ostream &outfile)
59 {
60  //Tecplot header info
61  outfile << "ZONE I=" << NNODE_1D << std::endl;
62 
63  //Find the dimension of the node
64  unsigned n_dim = this->nodal_dimension();
65 
66  //Loop over element nodes
67  for(unsigned l=0;l<NNODE_1D;l++)
68  {
69  //Loop over the dimensions and output the position
70  for(unsigned i=0;i<n_dim;i++)
71  {
72  outfile << node_pt(l)->x(i) << " ";
73  }
74  //Find out how many data values at the node
75  unsigned initial_nvalue = node_pt(l)->nvalue();
76  //Lopp over the data and output whether pinned or not
77  for(unsigned i=0;i<initial_nvalue;i++)
78  {
79  outfile << node_pt(l)->is_pinned(i) << " ";
80  }
81  outfile << std::endl;
82  }
83  outfile << std::endl;
84 }
85 
86 //=======================================================================
87 /// The output function for n_plot points in each coordinate direction
88 //=======================================================================
89 template <unsigned NNODE_1D>
90 void TElement<1,NNODE_1D>::output(std::ostream &outfile,
91  const unsigned &n_plot)
92 {
93  //Local variables
94  Vector<double> s(1);
95 
96  //Tecplot header info
97  outfile << "ZONE I=" << n_plot << std::endl;
98 
99  //Find the dimension of the nodes
100  unsigned n_dim = this->nodal_dimension();
101 
102  //Loop over plot points
103  for(unsigned l=0;l<n_plot;l++)
104  {
105  s[0] = l*1.0/(n_plot-1);
106  //Output the coordinates
107  for (unsigned i=0;i<n_dim;i++)
108  {
109  outfile << interpolated_x(s,i) << " " ;
110  }
111  outfile << std::endl;
112  }
113  outfile << std::endl;
114 }
115 
116 
117 
118 //=======================================================================
119 /// C style output function for general 1D TElements
120 //=======================================================================
121 template <unsigned NNODE_1D>
122 void TElement<1,NNODE_1D>::output(FILE* file_pt)
123 {
124 
125  //Tecplot header info
126  fprintf(file_pt,"ZONE I=%i\n",NNODE_1D);
127 
128  //Find the dimension of the nodes
129  unsigned n_dim = this->nodal_dimension();
130 
131  //Loop over element nodes
132  for(unsigned l=0;l<NNODE_1D;l++)
133  {
134  //Loop over the dimensions and output the position
135  for(unsigned i=0;i<n_dim;i++)
136  {
137  //outfile << Node_pt[l]->x(i) << " ";
138  fprintf(file_pt,"%g ",node_pt(l)->x(i));
139  }
140  //Find out how many data values at the node
141  unsigned initial_nvalue = node_pt(l)->nvalue();
142  //Lopp over the data and output whether pinned or not
143  for(unsigned i=0;i<initial_nvalue;i++)
144  {
145  //outfile << Node_pt[l]->is_pinned(i) << " ";
146  fprintf(file_pt,"%i ",node_pt(l)->is_pinned(i));
147  }
148  //outfile << std::endl;
149  fprintf(file_pt,"\n");
150  }
151  //outfile << std::endl;
152  fprintf(file_pt,"\n");
153 }
154 
155 //=======================================================================
156 /// C style output function for n_plot points in each coordinate direction
157 //=======================================================================
158 template <unsigned NNODE_1D>
159 void TElement<1,NNODE_1D>::output(FILE* file_pt, const unsigned &n_plot)
160 {
161  //Local variables
162  Vector<double> s(1);
163 
164  //Tecplot header info
165  //outfile << "ZONE I=" << n_plot << std::endl;
166  fprintf(file_pt,"ZONE I=%i\n",n_plot);
167 
168  //Find the dimension of the first node
169  unsigned n_dim = this->nodal_dimension();
170 
171  //Loop over plot points
172  for(unsigned l=0;l<n_plot;l++)
173  {
174  s[0] = l*1.0/(n_plot-1);
175 
176  //Output the coordinates
177  for (unsigned i=0;i<n_dim;i++)
178  {
179  //outfile << interpolated_x(s,i) << " " ;
180  fprintf(file_pt,"%g ",interpolated_x(s,i));
181  }
182  //outfile << std::endl;
183  fprintf(file_pt,"\n");
184  }
185  //outfile << std::endl;
186  fprintf(file_pt,"\n");
187 }
188 
189 //=============================================================
190 ///Namespace for helper functions that return the local
191 ///coordinates in the bulk elements
192 //==============================================================
193 namespace TElement1FaceToBulkCoordinates
194 {
195  ///The translation scheme for the face s0 = 0.0
196  void face0(const Vector<double> &s, Vector<double> &s_bulk)
197  {
198  s_bulk[0] = 0.0;
199  }
200 
201  ///The translation scheme for the face s0 = 1.0
202  void face1(const Vector<double> &s, Vector<double> &s_bulk)
203  {
204  s_bulk[0] = 1.0;
205  }
206 }
207 
208 
209 //=============================================================
210 /// Namespace for helper functions that calculate derivatives
211 /// of the local coordinates in the bulk elements wrt the
212 /// local coordinates in the face element.
213 //=============================================================
214 namespace TElement1BulkCoordinateDerivatives
215 {
216  ///Function for both faces -- the bulk coordinate is fixed on both
217  void faces0(const Vector<double> &s,
218  DenseMatrix<double> &dsbulk_dsface,
219  unsigned &interior_direction)
220  {
221  //Bulk coordinate s[0] does not vary along the face
222  dsbulk_dsface(0,0) = 0.0;
223 
224  //The interior direction is given by s[0]
225  interior_direction=0;
226  }
227 }
228 
229 //===========================================================
230 /// Function to setup geometrical information for lower-dimensional
231 /// FaceElements (which are Point Elements).
232 //===========================================================
233 template<unsigned NNODE_1D>
234 void TElement<1,NNODE_1D>::build_face_element(const int &face_index,
235  FaceElement *face_element_pt)
236 {
237  // Overload the nodal dimension
238  face_element_pt->set_nodal_dimension(nodal_dimension());
239 
240  // Set the pointer to the "bulk" element
241  face_element_pt->bulk_element_pt()=this;
242 
243 #ifdef OOMPH_HAS_MPI
244  // Pass on non-halo proc ID
245  face_element_pt->set_halo(Non_halo_proc_ID);
246 #endif
247 
248  // Resize the storage for the original number of values at the (one and only)
249  // node of the face element.
250  face_element_pt->nbulk_value_resize(1);
251 
252  // Resize the storage for the bulk node number corresponding to the (one
253  // and only) node of the face element
254  face_element_pt->bulk_node_number_resize(1);
255 
256  //Set the face index in the face element
257  face_element_pt->face_index() = face_index;
258 
259  //Now set up the node pointer
260  //The convention is that the "normal", should always point
261  //out of the element
262  switch(face_index)
263  {
264  //Bottom, normal sign is negative (coordinate points into element)
265  case(-1):
266  face_element_pt->node_pt(0) = node_pt(0);
267  face_element_pt->bulk_node_number(0) = 0;
268  face_element_pt->normal_sign() = -1;
269 
270  //Set the pointer to the function that determines the bulk coordinates
271  //in the face element
272  face_element_pt->face_to_bulk_coordinate_fct_pt() =
274 
275  //Set the pointer to the function that determines the mapping of
276  //derivatives
277  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
279 
280  //Set the number of values stored when the node is part of the "bulk"
281  //element. The required_nvalue() must be used, rather than nvalue(),
282  //because otherwise nodes on boundaries will be resized multiple
283  //times. If you want any other behaviour, you MUST set nbulk_value()
284  //manually after construction of your specific face element.
285  face_element_pt->nbulk_value(0) = required_nvalue(0);
286  break;
287 
288  //Top, normal sign is positive (coordinate points out of element)
289  case(1):
290  face_element_pt->node_pt(0) = node_pt(NNODE_1D-1);
291  face_element_pt->bulk_node_number(0) = NNODE_1D-1;
292  face_element_pt->normal_sign() = +1;
293 
294  //Set the pointer to the function that determines the bulk coordinates
295  //in the face element
296  face_element_pt->face_to_bulk_coordinate_fct_pt() =
298 
299  //Set the pointer to the function that determines the mapping of
300  //derivatives
301  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
303 
304 
305  //Set the number of values stored when the node is part of the "bulk"
306  //element.
307  face_element_pt->nbulk_value(0) = required_nvalue(NNODE_1D-1);
308  break;
309 
310  //All other cases throw an error
311  default:
312  std::ostringstream error_message;
313  error_message << "Face_index should only take "
314  << "the values +/-1, not " << face_index << std::endl;
315 
316  throw OomphLibError(error_message.str(),
317  OOMPH_CURRENT_FUNCTION,
318  OOMPH_EXCEPTION_LOCATION);
319  }
320 }
321 
322 
323 
324 ////////////////////////////////////////////////////////////////
325 // 2D Telements
326 ////////////////////////////////////////////////////////////////
327 
328 /// Assign the nodal translation schemes
329 template<>
330 const unsigned TElement<2,2>::Node_on_face[3][2] ={{2,1},{2,0},{0,1}};
331 
332 template<>
333 const unsigned TElement<2,3>::Node_on_face[3][3] ={{2,4,1},{2,5,0},{0,3,1}};
334 
335 template<>
336 const unsigned TElement<2,4>::Node_on_face[3][4] =
337 {{2,6,5,1},{2,7,8,0},{0,3,4,1}};
338 
339 
340 //===================================================================
341 /// Namespace for the functions that translate local face coordinates
342 /// to the coordinates in the bulk element
343 //==================================================================
344 namespace TElement2FaceToBulkCoordinates
345 {
346  ///The translation scheme for the face s0 = 0
347  void face0(const Vector<double> &s, Vector<double> &s_bulk)
348  {
349  s_bulk[0] = 0.0;
350  s_bulk[1] = s[0];
351  }
352 
353  ///The translation scheme for the face s1 = 0
354  void face1(const Vector<double> &s, Vector<double> &s_bulk)
355  {
356  s_bulk[0] = s[0];
357  s_bulk[1] = 0.0;
358  }
359 
360  ///The translation scheme for the face s2 = 0
361  void face2(const Vector<double> &s, Vector<double> &s_bulk)
362  {
363  s_bulk[0] = 1.0 - s[0];
364  s_bulk[1] = s[0];
365  }
366 }
367 
368 
369 //=============================================================
370 /// Namespace for helper functions that calculate derivatives
371 /// of the local coordinates in the bulk elements wrt the
372 /// local coordinates in the face element.
373 //=============================================================
374 namespace TElement2BulkCoordinateDerivatives
375 {
376  ///Function for the "left" face along which s0 is fixed
377  void face0(const Vector<double> &s,
378  DenseMatrix<double> &dsbulk_dsface,
379  unsigned &interior_direction)
380  {
381  //Bulk coordinate s[0] does not vary along the face
382  dsbulk_dsface(0,0) = 0.0;
383  //Bulk coordinate s[1] is the face coordinate
384  dsbulk_dsface(1,0) = 1.0;
385 
386  //The interior direction is given by s[0]
387  interior_direction=0;
388  }
389 
390 
391  ///Function for the "bottom" face along which s1 is fixed
392  void face1(const Vector<double> &s,
393  DenseMatrix<double> &dsbulk_dsface,
394  unsigned &interior_direction)
395  {
396  //Bulk coordinate s[0] is the face coordinate
397  dsbulk_dsface(0,0) = 1.0;
398  //Bulk coordinate s[1] does not vary along the face
399  dsbulk_dsface(1,0) = 0.0;
400 
401  //The interior direction is given by s[1]
402  interior_direction=1;
403  }
404 
405  ///Function for the sloping face
406  void face2(const Vector<double> &s,
407  DenseMatrix<double> &dsbulk_dsface,
408  unsigned &interior_direction)
409  {
410  //Bulk coordinate s[0] decreases along the face
411  dsbulk_dsface(0,0) = -1.0;
412  //Bulk coordinate s[1] increases along the face
413  dsbulk_dsface(1,0) = 1.0;
414 
415  //The interior direction is given by s[0] (or s[1])
416  interior_direction=0;
417  }
418 
419 }
420 
421 
422 //=======================================================================
423 /// Function to setup geometrical information for lower-dimensional
424 /// FaceElements (which are of type TElement<2,NNODE_1D>).
425 //=======================================================================
426  template<unsigned NNODE_1D>
427  void TElement<2,NNODE_1D>::build_face_element(const int &face_index,
428  FaceElement* face_element_pt)
429 {
430  //Set the nodal dimension from the first node
431  face_element_pt->set_nodal_dimension(nodal_dimension());
432 
433  //Set the pointer to the orginal "bulk" element
434  face_element_pt->bulk_element_pt()=this;
435 
436 #ifdef OOMPH_HAS_MPI
437  // Pass on non-halo proc ID
438  face_element_pt->set_halo(Non_halo_proc_ID);
439 #endif
440 
441  //Calculate the number of nodes in the face element
442  const unsigned n_face_nodes = NNODE_1D;
443 
444  // Resize storage for the number of values originally stored
445  // at the face element's nodes.
446  face_element_pt->nbulk_value_resize(n_face_nodes);
447 
448  // Resize storage for the bulk node numbers corresponding to
449  // the nodes of the face
450  face_element_pt->bulk_node_number_resize(n_face_nodes);
451 
452  //Set the face index in the face element
453  face_element_pt->face_index() = face_index;
454 
455  //So the faces are
456  // 0 : s_0 fixed
457  // 1 : s_1 fixed
458  // 2 : sloping face
459 
460  // Copy nodes across to the face
461  for (unsigned i=0;i<n_face_nodes;i++)
462  {
463  //The number is just offset by one
464  unsigned bulk_number = Node_on_face[face_index][i];
465  face_element_pt->node_pt(i)=node_pt(bulk_number);
466  face_element_pt->bulk_node_number(i) = bulk_number;
467  //set the number of values originally stored at this node
468  face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
469  }
470 
471  //Now set up the node pointers and the normal vectors
472  switch(face_index)
473  {
474  //
475  //-----The face s0 is constant
476  case 0:
477 
478  //Set the pointer to the function that determines the bulk
479  //coordinates in the face element
480  face_element_pt->face_to_bulk_coordinate_fct_pt() =
482 
483  //Set the pointer to the function that determines the mapping of
484  //derivatives
485  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
487 
488  // Outer unit normal is the positive cross product of two in plane
489  // tangent vectors
490  face_element_pt->normal_sign()=+1;
491 
492  break;
493 
494  //-----The face s1 is constant
495  case 1:
496 
497  //Set the pointer to the function that determines the bulk
498  //coordinates in the face element
499  face_element_pt->face_to_bulk_coordinate_fct_pt() =
501 
502  //Set the pointer to the function that determines the mapping of
503  //derivatives
504  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
506 
507  // Outer unit normal is the positive cross product of two in plane
508  // tangent vectors
509  face_element_pt->normal_sign()=+1;
510 
511  break;
512 
513  //
514  //-----The face s2 is constant
515  case 2:
516 
517  //Set the pointer to the function that determines the bulk
518  //coordinates in the face element
519  face_element_pt->face_to_bulk_coordinate_fct_pt() =
521 
522  //Set the pointer to the function that determines the mapping of
523  //derivatives
524  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
526 
527  // Outer unit normal is the negative cross product of two in plane
528  // tangent vectors
529  face_element_pt->normal_sign()=-1;
530 
531  break;
532 
533  default:
534 
535  std::ostringstream error_message;
536  error_message << "Face_index should only take "
537  << "the values 0, 1 or 2 not " << face_index << std::endl;
538 
539  throw OomphLibError(error_message.str(),
540  OOMPH_CURRENT_FUNCTION,
541  OOMPH_EXCEPTION_LOCATION);
542  } //end switch
543 
544 }
545 
546 
547 
548 
549 //=======================================================================
550 /// The output function for TElement<2,NNODE_1D>
551 //=======================================================================
552 template<unsigned NNODE_1D>
553 void TElement<2,NNODE_1D>::output(std::ostream &outfile)
554 {
555  unsigned n_node = nnode();
556  outfile << "ZONE N=" << n_node << ", F=FEPOINT, ET=TRIANGLE" << std::endl;
557 
558  //Find the dimensions of the nodes
559  unsigned n_dim = this->nodal_dimension();
560 
561  //Loop over local nodes
562  for (unsigned n=0;n<n_node;++n)
563  {
564  //Loop over the dimensions and output the position
565  for(unsigned i=0;i<n_dim;i++)
566  {
567  outfile << node_pt(n)->x(i) << " ";
568  }
569  //Find out how many data values at the node
570  unsigned initial_nvalue = node_pt(n)->nvalue();
571  //Loop over the data and output whether pinned or not
572  for(unsigned i=0;i<initial_nvalue;i++)
573  {
574  outfile << node_pt(n)->is_pinned(i) << " ";
575  }
576  outfile << std::endl;
577  }
578 
579  //Output local node connectivity list (for tecplot)
580  for (unsigned n=0;n<n_node;++n)
581  {
582  outfile << n+1 << " ";
583  }
584  outfile << std::endl;
585 
586 }
587 
588 
589 
590 //=======================================================================
591 /// The output function for TElement<2,NNODE_1D>
592 //=======================================================================
593 template<unsigned NNODE_1D>
594 void TElement<2,NNODE_1D>::output(std::ostream &outfile,const unsigned &nplot)
595 {
596 
597  //Vector of local coordinates
598  Vector<double> s(2);
599 
600  //Get the dimension of the node
601  unsigned n_dim = this->nodal_dimension();
602 
603  // Tecplot header info
604  outfile << tecplot_zone_string(nplot);
605 
606  // Loop over plot points
607  unsigned num_plot_points=nplot_points(nplot);
608  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
609  {
610  // Get local coordinates of plot point
611  get_s_plot(iplot,nplot,s);
612 
613  for(unsigned i=0;i<n_dim;i++)
614  {
615  outfile << interpolated_x(s,i) << " ";
616  }
617  outfile << std::endl;
618  }
619 
620  // Write tecplot footer (e.g. FE connectivity lists)
621  write_tecplot_zone_footer(outfile,nplot);
622 
623 }
624 
625 //=======================================================================
626 /// The C-style output function for TElement<2,NNODE_1D>
627 //=======================================================================
628 template<unsigned NNODE_1D>
629 void TElement<2,NNODE_1D>::output(FILE* file_pt)
630 {
631  unsigned n_node = nnode();
632  fprintf(file_pt,"ZONE N=%i, F=FEPOINT, ET=TRIANGLE\n",n_node);
633  //outfile << "ZONE N=" << n_node << ", F=FEPOINT, ET=TRIANGLE" << std::endl;
634 
635  //Find the dimension of the nodes
636  unsigned n_dim = this->nodal_dimension();
637 
638  //Loop over local nodes
639  for (unsigned n=0;n<n_node;++n)
640  {
641  //Loop over the dimensions and output the position
642  for(unsigned i=0;i<n_dim;i++)
643  {
644  fprintf(file_pt,"%g ", node_pt(n)->x(i));
645  //outfile << node_pt(n)->x(i) << " ";
646  }
647  //Find out how many data values at the node
648  unsigned initial_nvalue = node_pt(n)->nvalue();
649  //Loop over the data and output whether pinned or not
650  for(unsigned i=0;i<initial_nvalue;i++)
651  {
652  fprintf(file_pt,"%i ",node_pt(n)->is_pinned(i));
653  //outfile << node_pt(n)->is_pinned(i) << " ";
654  }
655  fprintf(file_pt,"\n");
656  //outfile << std::endl;
657  }
658 
659  //Output local node connectivity list (for tecplot)
660  for (unsigned n=0;n<n_node;++n)
661  {
662  fprintf(file_pt,"%i ",n+1);
663  //outfile << n+1 << " ";
664  }
665  fprintf(file_pt,"\n");
666  //outfile << std::endl;
667 }
668 
669 
670 
671 //=======================================================================
672 /// The C-style output function for TElement<2,NNODE_1D>
673 //=======================================================================
674 template<unsigned NNODE_1D>
675 void TElement<2,NNODE_1D>::output(FILE* file_pt,const unsigned &nplot)
676 {
677 
678  //Vector of local coordinates
679  Vector<double> s(2);
680 
681  //Find the dimensions of the nodes
682  unsigned n_dim = this->nodal_dimension();
683 
684  // Tecplot header info
685  fprintf(file_pt,"%s \n",tecplot_zone_string(nplot).c_str());
686 
687  // Loop over plot points
688  unsigned num_plot_points=nplot_points(nplot);
689  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
690  {
691  // Get local coordinates of plot point
692  get_s_plot(iplot,nplot,s);
693 
694  for(unsigned i=0;i<n_dim;i++)
695  {
696  fprintf(file_pt,"%g ", interpolated_x(s,i));
697  //outfile << interpolated_x(s,i) << " ";
698  }
699  fprintf(file_pt,"\n");
700  //outfile << std::endl;
701  }
702 
703  // Write tecplot footer (e.g. FE connectivity lists)
704  write_tecplot_zone_footer(file_pt,nplot);
705 
706 }
707 
708 
709 ///////////////////////////////////////////////////////////////////////
710 ///////////////////////////////////////////////////////////////////////
711 ///////////////////////////////////////////////////////////////////////
712 
713 
714 
715 //=======================================================================
716 /// The output function for TElement<3,NNODE_1D>
717 //=======================================================================
718 template<unsigned NNODE_1D>
719 void TElement<3,NNODE_1D>::output(std::ostream &outfile)
720 {
721  unsigned n_node = nnode();
722  outfile << "ZONE N=" << n_node << ", F=FEPOINT, ET=TETRAHEDRON" << std::endl;
723 
724  //Find the dimension of the nodes
725  unsigned n_dim = this->nodal_dimension();
726 
727  //Loop over local nodes
728  for (unsigned n=0;n<n_node;++n)
729  {
730  //Loop over the dimensions and output the position
731  for(unsigned i=0;i<n_dim;i++)
732  {
733  outfile << node_pt(n)->x(i) << " ";
734  }
735 // //Find out how many data values at the node
736 // unsigned initial_nvalue = node_pt(n)->nvalue();
737 // //Loop over the data and output whether pinned or not
738 // for(unsigned i=0;i<initial_nvalue;i++)
739 // {
740 // outfile << node_pt(n)->is_pinned(i) << " ";
741 // }
742  outfile << std::endl;
743  }
744 
745  // Write tecplot footer (e.g. FE connectivity lists)
746  write_tecplot_zone_footer(outfile,NNODE_1D);
747 
748 }
749 
750 
751 
752 //=======================================================================
753 /// The output function for TElement<3,NNODE_1D>
754 //=======================================================================
755 template<unsigned NNODE_1D>
756 void TElement<3,NNODE_1D>::output(std::ostream &outfile,const unsigned &nplot)
757 {
758 
759  //Vector of local coordinates
760  Vector<double> s(3);
761 
762  //Find the dimension of the nodes
763  unsigned n_dim = this->nodal_dimension();
764 
765  // Tecplot header info
766  outfile << tecplot_zone_string(nplot);
767 
768  // Loop over plot points
769  unsigned num_plot_points=nplot_points(nplot);
770  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
771  {
772  // Get local coordinates of plot point
773  get_s_plot(iplot,nplot,s);
774 
775  for(unsigned i=0;i<n_dim;i++)
776  {
777  outfile << interpolated_x(s,i) << " ";
778  }
779  outfile << "\n";
780  }
781 
782  // Write tecplot footer (e.g. FE connectivity lists)
783  write_tecplot_zone_footer(outfile,nplot);
784 
785 }
786 
787 
788 //=======================================================================
789 /// The C-style output function for TElement<3,NNODE_1D>
790 //=======================================================================
791 template<unsigned NNODE_1D>
792 void TElement<3,NNODE_1D>::output(FILE* file_pt)
793 {
794  unsigned n_node = nnode();
795  fprintf(file_pt,"ZONE N=%i, F=FEPOINT, ET=TETRAHEDRON\n",n_node);
796  //outfile << "ZONE N=" << n_node << ", F=FEPOINT, ET=TRIANGLE" << std::endl;
797 
798  //Find the dimension of the nodes
799  unsigned n_dim = this->nodal_dimension();
800 
801  //Loop over local nodes
802  for (unsigned n=0;n<n_node;++n)
803  {
804  //Loop over the dimensions and output the position
805  for(unsigned i=0;i<n_dim;i++)
806  {
807  fprintf(file_pt,"%g ", node_pt(n)->x(i));
808  //outfile << node_pt(n)->x(i) << " ";
809  }
810  //Find out how many data values at the node
811  unsigned initial_nvalue = node_pt(n)->nvalue();
812  //Loop over the data and output whether pinned or not
813  for(unsigned i=0;i<initial_nvalue;i++)
814  {
815  fprintf(file_pt,"%i ",node_pt(n)->is_pinned(i));
816  //outfile << node_pt(n)->is_pinned(i) << " ";
817  }
818  fprintf(file_pt,"\n");
819  //outfile << std::endl;
820  }
821 
822  //Output local node connectivity list (for tecplot)
823  for (unsigned n=0;n<n_node;++n)
824  {
825  fprintf(file_pt,"%i ",n+1);
826  //outfile << n+1 << " ";
827  }
828  fprintf(file_pt,"\n");
829  //outfile << std::endl;
830 }
831 
832 
833 //=======================================================================
834 /// The C-style output function for TElement<3,NNODE_1D>
835 //=======================================================================
836 template<unsigned NNODE_1D>
837 void TElement<3,NNODE_1D>::output(FILE* file_pt,const unsigned &nplot)
838 {
839 
840  //Vector of local coordinates
841  Vector<double> s(3);
842 
843  //Find the dimension of the nodes
844  unsigned n_dim = this->nodal_dimension();
845 
846  // Tecplot header info
847  fprintf(file_pt,"%s \n",tecplot_zone_string(nplot).c_str());
848 
849  // Loop over plot points
850  unsigned num_plot_points=nplot_points(nplot);
851  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
852  {
853  // Get local coordinates of plot point
854  get_s_plot(iplot,nplot,s);
855 
856  for(unsigned i=0;i<n_dim;i++)
857  {
858  fprintf(file_pt,"%g ", interpolated_x(s,i));
859  //outfile << interpolated_x(s,i) << " ";
860  }
861  fprintf(file_pt,"\n");
862  //outfile << std::endl;
863  }
864 
865  // Write tecplot footer (e.g. FE connectivity lists)
866  write_tecplot_zone_footer(file_pt,nplot);
867 
868 }
869 
870 //===================================================================
871 /// Namespace for the functions that translate local face coordinates
872 /// to the coordinates in the bulk element
873 //==================================================================
874 namespace TElement3FaceToBulkCoordinates
875 {
876  ///The translation scheme for the face s0 = 0
877  void face0(const Vector<double> &s, Vector<double> &s_bulk)
878  {
879  s_bulk[0] = 0.0;
880  s_bulk[1] = s[0];
881  s_bulk[2] = s[1];
882  }
883 
884  ///The translation scheme for the face s1 = 0
885  void face1(const Vector<double> &s, Vector<double> &s_bulk)
886  {
887  s_bulk[0] = s[0];
888  s_bulk[1] = 0.0;
889  s_bulk[2] = s[1];
890  }
891 
892  ///The translation scheme for the face s2 = 0
893  void face2(const Vector<double> &s, Vector<double> &s_bulk)
894  {
895  s_bulk[0] = s[0];
896  s_bulk[1] = s[1];
897  s_bulk[2] = 0.0;
898  }
899 
900  ///The translation scheme for the sloping face
901  void face3(const Vector<double> &s, Vector<double> &s_bulk)
902  {
903  s_bulk[0] = 1 - s[0] - s[1];
904  s_bulk[1] = s[0];
905  s_bulk[2] = s[1];
906  }
907 
908 }
909 
910 ///Assign the nodal translation scheme for linear elements
911 template<>
912 const unsigned TElement<3,2>::Node_on_face[4][3] =
913 {{1,2,3},{0,2,3},{0,1,3},{1,2,0}};
914 
915 ///Assign the nodal translation scheme for quadratic elements
916 template<>
917 const unsigned TElement<3,3>::Node_on_face[4][6] =
918 {{1,2,3,7,8,9},{0,2,3,5,8,6},{0,1,3,4,9,6},{1,2,0,7,5,4}};
919 
920 
921 //=======================================================================
922 /// Function to setup geometrical information for lower-dimensional
923 /// FaceElements (which are of type TElement<2,NNODE_1D>).
924 //=======================================================================
925  template<unsigned NNODE_1D>
926  void TElement<3,NNODE_1D>::build_face_element(const int &face_index,
927  FaceElement* face_element_pt)
928 {
929  //Set the nodal dimension
930  face_element_pt->set_nodal_dimension(nodal_dimension());
931 
932  //Set the pointer to the orginal "bulk" element
933  face_element_pt->bulk_element_pt()=this;
934 
935 #ifdef OOMPH_HAS_MPI
936  // Pass on non-halo proc ID
937  face_element_pt->set_halo(Non_halo_proc_ID);
938 #endif
939 
940  //Calculate the number of nodes in the face element
941  const unsigned n_face_nodes = (NNODE_1D*(NNODE_1D+1))/2;
942 
943  // Resize storage for the number of values originally stored
944  // at the face element's nodes.
945  face_element_pt->nbulk_value_resize(n_face_nodes);
946 
947  // Resize storage for the bulk node numbers corresponding to
948  // the nodes of the face
949  face_element_pt->bulk_node_number_resize(n_face_nodes);
950 
951  //Set the face index in the face element
952  face_element_pt->face_index() = face_index;
953 
954  //So the faces are
955  // 0 : s_0 fixed
956  // 1 : s_1 fixed
957  // 2 : s_2 fixed
958  // 3 : sloping face
959 
960  // Copy nodes across to the face
961  for (unsigned i=0;i<n_face_nodes;i++)
962  {
963  //The number is just offset by one
964  unsigned bulk_number = Node_on_face[face_index][i];
965  face_element_pt->node_pt(i)=node_pt(bulk_number);
966  face_element_pt->bulk_node_number(i) = bulk_number;
967  //set the number of values originally stored at this node
968  face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
969  }
970 
971  //Now set up the node pointers and the normal vectors
972  switch(face_index)
973  {
974  //
975  //-----The face s0 is constant
976  case 0:
977 
978  //Set the pointer to the function that determines the bulk
979  //coordinates in the face element
980  face_element_pt->face_to_bulk_coordinate_fct_pt() =
982 
983  // Outer unit normal is the negative cross product of two in plane
984  // tangent vectors
985  face_element_pt->normal_sign()=-1;
986 
987  break;
988 
989  //-----The face s1 is constant
990  case 1:
991 
992  //Set the pointer to the function that determines the bulk
993  //coordinates in the face element
994  face_element_pt->face_to_bulk_coordinate_fct_pt() =
996 
997  // Outer unit normal is the positive cross product of two in plane
998  // tangent vectors
999  face_element_pt->normal_sign()=+1;
1000 
1001  break;
1002 
1003  //
1004  //-----The face s2 is constant
1005  case 2:
1006 
1007  //Set the pointer to the function that determines the bulk
1008  //coordinates in the face element
1009  face_element_pt->face_to_bulk_coordinate_fct_pt() =
1011 
1012  // Outer unit normal is the negative cross product of two in plane
1013  // tangent vectors
1014  face_element_pt->normal_sign()=-1;
1015 
1016  break;
1017 
1018  //-----The sloping face of the tetrahedron
1019  case 3:
1020 
1021  //Set the pointer to the function that determines the bulk
1022  //coordinates in the face element
1023  face_element_pt->face_to_bulk_coordinate_fct_pt() =
1025 
1026  // Outer unit normal is the positive cross product of two in plane
1027  // tangent vectors
1028  face_element_pt->normal_sign()=+1;
1029 
1030  break;
1031 
1032 
1033 
1034  default:
1035 
1036  std::ostringstream error_message;
1037  error_message << "Face_index should only take "
1038  << "the values 0, 1, 2 or 3, not "
1039  << face_index << std::endl;
1040 
1041  throw OomphLibError(error_message.str(),
1042  OOMPH_CURRENT_FUNCTION,
1043  OOMPH_EXCEPTION_LOCATION);
1044  } //end switch
1045 
1046 }
1047 
1048 
1049 //==================================================================
1050 //Default integration scheme for the TBubbleEnrichedElement<2,3>
1051 //===================================================================
1052 template<unsigned DIM>
1055 
1056 //===================================================================
1057 //Central node on the face of the TBubbleEnrichedelement<2,3>
1058 //===================================================================
1059 template<>
1061 
1062 
1063 //================================================================
1064 ///The face element for is the same in the two-dimesional case
1065 //================================================================
1066 template<>
1068  const int &face_index,FaceElement* face_element_pt)
1069 {TElement<2,3>::build_face_element(face_index,face_element_pt);}
1070 
1071 
1072 //===================================================================
1073 //Central node on the face of the TBubbleEnrichedelement<3,3>
1074 //===================================================================
1075 template<>
1077 {13,12,10,11};
1078 
1079 
1080 
1081 //=======================================================================
1082 /// Function to setup geometrical information for lower-dimensional
1083 /// FaceElements (which are of type TBubbleEnrichedElement<2,3>).
1084 //=======================================================================
1085  template<>
1087  const int &face_index,
1088  FaceElement* face_element_pt)
1089 {
1090  //Call the standard unenriched build function
1091  TElement<3,3>::build_face_element(face_index,face_element_pt);
1092 
1093  //Set the enriched number of total face nodes
1094  const unsigned n_face_nodes = 7;
1095 
1096  // Resize storage for the number of values originally stored
1097  // at the face element's nodes.
1098  face_element_pt->nbulk_value_resize(n_face_nodes);
1099 
1100  // Resize storage for the bulk node numbers corresponding to
1101  // the nodes of the face
1102  face_element_pt->bulk_node_number_resize(n_face_nodes);
1103 
1104  //So the faces are
1105  // 0 : s_0 fixed
1106  // 1 : s_1 fixed
1107  // 2 : s_2 fixed
1108  // 3 : sloping face
1109 
1110  //Copy central node across
1111  unsigned bulk_number = Central_node_on_face[face_index];
1112  face_element_pt->node_pt(n_face_nodes-1)=node_pt(bulk_number);
1113  face_element_pt->bulk_node_number(n_face_nodes-1) = bulk_number;
1114  //set the number of values originally stored at this node
1115  face_element_pt->nbulk_value(n_face_nodes-1) = required_nvalue(bulk_number);
1116 }
1117 
1118 /////////////////////////////////////////////////////////////////////////
1119 /////////////////////////////////////////////////////////////////////////
1120 /////////////////////////////////////////////////////////////////////////
1121 
1122 //===========================================================
1123 /// Final override for
1124 /// function to setup geometrical information for lower-dimensional
1125 /// FaceElements (which are of type SolidTBubbleEnrichedElement<1,3>).
1126 //===========================================================
1127 template<>
1129 build_face_element(const int &face_index, FaceElement *face_element_pt)
1130 {
1131  //Build the standard non-solid FaceElement
1133  build_face_element(face_index,face_element_pt);
1134 
1135  //Set the Lagrangian dimension from the first node of the present element
1136  dynamic_cast<SolidFiniteElement*>(face_element_pt)->
1137  set_lagrangian_dimension(static_cast<SolidNode*>(node_pt(0))->nlagrangian());
1138 }
1139 
1140 
1141 //===========================================================
1142 /// Final override for
1143 /// function to setup geometrical information for lower-dimensional
1144 /// FaceElements (which are of type SolidTBubbleEnrichedElement<2,3>).
1145 //===========================================================
1146 template<>
1148 build_face_element(const int &face_index, FaceElement *face_element_pt)
1149 {
1150  //Build the standard non-solid FaceElement
1152  build_face_element(face_index,face_element_pt);
1153 
1154  //Set the Lagrangian dimension from the first node of the present element
1155  dynamic_cast<SolidFiniteElement*>(face_element_pt)->
1156  set_lagrangian_dimension(static_cast<SolidNode*>(node_pt(0))->nlagrangian());
1157 }
1158 
1159 
1160 
1161 
1162 //===================================================================
1163 // Build required templates
1164 //===================================================================
1165 template class TElement<1,2>;
1166 template class TElement<1,3>;
1167 template class TElement<1,4>;
1168 template class TElement<2,2>;
1169 template class TElement<2,3>;
1170 template class TElement<2,4>;
1171 template class TElement<3,2>;
1172 template class TElement<3,3>;
1173 template class TBubbleEnrichedElement<2,3>;
1174 template class TBubbleEnrichedElement<3,3>;
1175 template class SolidTBubbleEnrichedElement<2,3>;
1176 template class SolidTBubbleEnrichedElement<3,3>;
1177 template class TBubbleEnrichedGauss<2,3>;
1178 template class TBubbleEnrichedGauss<3,3>;
1179 }
cstr elem_len * i
Definition: cfortran.h:607
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = 0.
Definition: Telements.cc:877
void face3(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the sloping face.
Definition: Telements.cc:901
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
BulkCoordinateDerivativesFctPt & bulk_coordinate_derivatives_fct_pt()
Return the pointer to the function that returns the derivatives of the bulk coordinates wrt the face ...
Definition: elements.h:4495
void faces0(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for both faces – the bulk coordinate is fixed on both.
Definition: Telements.cc:217
unsigned & bulk_node_number(const unsigned &n)
Return the bulk node number that corresponds to the n-th local node number.
Definition: elements.h:4539
int & normal_sign()
Sign of outer unit normal (relative to cross-products of tangent vectors in the corresponding "bulk" ...
Definition: elements.h:4359
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = 0.0.
Definition: Telements.cc:196
static TGauss< 1, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Telements.h:1287
void set_halo(const unsigned &non_halo_proc_ID)
Label the element as halo and specify processor that holds non-halo counterpart.
Definition: elements.h:1132
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s1 = 0.
Definition: Telements.cc:885
static char t char * s
Definition: cfortran.h:572
static TGauss< 2, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Telements.h:1546
void nbulk_value_resize(const unsigned &i)
Resize the storage for the number of values originally stored at the local nodes to i entries...
Definition: elements.h:4562
const unsigned Node_on_face[3][2]
Definition: Telements.cc:330
void face2(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the sloping face.
Definition: Telements.cc:406
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void output(std::ostream &outfile)
Output with default number of plot points.
static TGauss< 3, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Telements.h:3009
void face0(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the "left" face along which s0 is fixed.
Definition: Telements.cc:377
void face2(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s2 = 0.
Definition: Telements.cc:361
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s1 = 0.
Definition: Telements.cc:354
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4470
void face2(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s2 = 0.
Definition: Telements.cc:893
void set_nodal_dimension(const unsigned &nodal_dim)
Set the dimension of the nodes in the element. This will typically only be required when constructing...
Definition: elements.h:1348
CoordinateMappingFctPt & face_to_bulk_coordinate_fct_pt()
Return the pointer to the function that maps the face coordinate to the bulk coordinate.
Definition: elements.h:4484
void bulk_node_number_resize(const unsigned &i)
Resize the storage for the bulk node numbers.
Definition: elements.h:4535
unsigned & nbulk_value(const unsigned &n)
Return the number of values originally stored at local node n (before the FaceElement added additiona...
Definition: elements.h:4552
void build_face_element(const int &face_index, FaceElement *face_element_pt)
Definition: Telements.cc:1067
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = 1.0.
Definition: Telements.cc:202
SolidFiniteElement class.
Definition: elements.h:3320
void face1(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the "bottom" face along which s1 is fixed.
Definition: Telements.cc:392
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = 0.
Definition: Telements.cc:347