Qspectral_elements.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 #include "Qspectral_elements.h"
32 
33 
34 namespace oomph
35 {
36 
37 
38 std::map<unsigned,Vector<double> > OneDLegendreShapeParam::z;
39 
40 
41 //////////////////////////////////////////////////////////////////////////
42 /// 1D QLegendreElements
43 //////////////////////////////////////////////////////////////////////////
44 
45 
46 //=======================================================================
47 /// Assign the static integral
48 //=======================================================================
49 template<unsigned NNODE_1D>
50 GaussLobattoLegendre<1,NNODE_1D> QSpectralElement<1,NNODE_1D>::integral;
51 
52 //=======================================================================
53 /// The output function for general 1D QSpectralElements
54 //=======================================================================
55 template <unsigned NNODE_1D>
56 void QSpectralElement<1,NNODE_1D>::output(std::ostream &outfile)
57 {
58  //Tecplot header info
59  outfile << "ZONE I=" << NNODE_1D << std::endl;
60 
61  //Find the dimension of the nodes
62  unsigned n_dim = this->nodal_dimension();
63 
64  //Loop over element nodes
65  for(unsigned l=0;l<NNODE_1D;l++)
66  {
67  //Loop over the dimensions and output the position
68  for(unsigned i=0;i<n_dim;i++)
69  {
70  outfile << this->node_pt(l)->x(i) << " ";
71  }
72  //Find out how many data values at the node
73  unsigned initial_nvalue = this->node_pt(l)->nvalue();
74  //Lopp over the data and output whether pinned or not
75  for(unsigned i=0;i<initial_nvalue;i++)
76  {
77  outfile << this->node_pt(l)->is_pinned(i) << " ";
78  }
79  outfile << std::endl;
80  }
81  outfile << std::endl;
82 }
83 
84 //=======================================================================
85 /// The output function for nplot points in each coordinate direction
86 //=======================================================================
87 template <unsigned NNODE_1D>
88 void QSpectralElement<1,NNODE_1D>::output(std::ostream &outfile,
89  const unsigned& nplot)
90 {
91  //Local variables
92  Vector<double> s(1);
93  //Shape functions
94  Shape psi(NNODE_1D);
95 
96  //Find the dimension of the nodes
97  unsigned n_dim = this->nodal_dimension();
98 
99  //Tecplot header info
100  outfile << "ZONE I=" << nplot << std::endl;
101  //Loop over element nodes
102  for(unsigned l=0;l<nplot;l++)
103  {
104  s[0] = -1.0 + l*2.0/(nplot-1);
105  shape(s,psi);
106  for(unsigned i=0;i<n_dim;i++)
107  {
108  //Output the x and y positions
109  outfile << this->interpolated_x(s,i) << " ";
110  }
111  for(unsigned i=0;i<NNODE_1D;i++)
112  {
113  outfile << psi(i) << " ";
114  }
115  outfile << std::endl;
116  }
117  outfile << std::endl;
118 }
119 
120 //===========================================================
121 /// Function to setup geometrical information for lower-dimensional
122 /// FaceElements (which are of type QSpectralElement<0,1>).
123 //===========================================================
124 template<unsigned NNODE_1D>
126 build_face_element(const int &face_index,
127  FaceElement *face_element_pt)
128 {
129  /*throw OomphLibError("Untested",
130  OOMPH_CURRENT_FUNCTION,
131  OOMPH_EXCEPTION_LOCATION);*/
132 
133  // Overload the nodal dimension by reading out the value from the node
134  face_element_pt->set_nodal_dimension(this->node_pt(0)->ndim());
135 
136  // Set the pointer to the "bulk" element
137  face_element_pt->bulk_element_pt()=this;
138 
139 #ifdef OOMPH_HAS_MPI
140  // If the bulk element is halo then the face element must be too
141  //if (this->is_halo())
142  {
143  face_element_pt->set_halo(Non_halo_proc_ID);
144  }
145 #endif
146 
147  // Resize the storage for the original number of values at the (one and only)
148  // node of the face element.
149  face_element_pt->nbulk_value_resize(1);
150 
151  // Resize the storage for the bulk node number corresponding to the (one
152  // and only) node of the face element
153  face_element_pt->bulk_node_number_resize(1);
154 
155  //Set the face index in the face element
156  face_element_pt->face_index() = face_index;
157 
158  //Now set up the node pointer
159  //The convention is that the "normal", should always point
160  //out of the element
161  switch(face_index)
162  {
163  //Bottom, normal sign is negative (coordinate points into element)
164  case(-1):
165  face_element_pt->node_pt(0) = this->node_pt(0);
166  face_element_pt->bulk_node_number(0) = 0;
167  face_element_pt->normal_sign() = -1;
168 
169  //Set the pointer to the function that determines the bulk coordinates
170  //in the face element
171  face_element_pt->face_to_bulk_coordinate_fct_pt() =
173 
174  //Set the pointer to the function that determines the mapping of
175  //derivatives
176  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
178 
179  //Set the number of values stored when the node is part of the "bulk"
180  //element. The required_nvalue() must be used, rather than nvalue(),
181  //because otherwise nodes on boundaries will be resized multiple
182  //times. If you want any other behaviour, you MUST set nbulk_value()
183  //manually after construction of your specific face element.
184  face_element_pt->nbulk_value(0) = this->required_nvalue(0);
185  break;
186 
187  //Top, normal sign is positive (coordinate points out of element)
188  case(1):
189  face_element_pt->node_pt(0) = this->node_pt(NNODE_1D-1);
190  face_element_pt->bulk_node_number(0) = NNODE_1D-1;
191  face_element_pt->normal_sign() = +1;
192 
193 
194  //Set the pointer to the function that determines the bulk coordinates
195  //in the face element
196  face_element_pt->face_to_bulk_coordinate_fct_pt() =
198 
199  //Set the pointer to the function that determines the mapping of
200  //derivatives
201  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
203 
204 
205  //Set the number of values stored when the node is part of the "bulk"
206  //element.
207  face_element_pt->nbulk_value(0) = this->required_nvalue(NNODE_1D-1);
208  break;
209 
210  //Other cases
211  default:
212  std::ostringstream error_message;
213  error_message << "Face_index should only take "
214  << "the values +/-1, not " << face_index << std::endl;
215 
216  throw OomphLibError(error_message.str(),
217  OOMPH_CURRENT_FUNCTION,
218  OOMPH_EXCEPTION_LOCATION);
219  }
220 }
221 
222 
223 
224 //////////////////////////////////////////////////////////////////////////
225 /// 2D QLegendreElements
226 //////////////////////////////////////////////////////////////////////////
227 
228 //=======================================================================
229 /// Assign the static integral
230 //=======================================================================
231 template<unsigned NNODE_1D>
233 
234 //=======================================================================
235 /// The output function for general 1D QSpectralElements
236 //=======================================================================
237 template <unsigned NNODE_1D>
238 void QSpectralElement<2,NNODE_1D>::output(std::ostream &outfile)
239 {
240  //Tecplot header info
241  outfile << "ZONE I=" << NNODE_1D << ", J=" << NNODE_1D << std::endl;
242 
243  //Find the dimension of the nodes
244  unsigned n_dim = this->nodal_dimension();
245 
246  //Loop over element nodes
247  for(unsigned l2=0;l2<NNODE_1D;l2++)
248  {
249  for(unsigned l1=0;l1<NNODE_1D;l1++)
250  {
251  unsigned l = l2*NNODE_1D + l1;
252 
253  //Loop over the dimensions and output the position
254  for(unsigned i=0;i<n_dim;i++)
255  {
256  outfile << this->node_pt(l)->x(i) << " ";
257  }
258  //Find out how many data values at the node
259  unsigned initial_nvalue = this->node_pt(l)->nvalue();
260  //Loop over the data and output whether pinned or not
261  for(unsigned i=0;i<initial_nvalue;i++)
262  {
263  outfile << this->node_pt(l)->is_pinned(i) << " ";
264  }
265  outfile << std::endl;
266  }
267  }
268  outfile << std::endl;
269 }
270 
271 
272 //=======================================================================
273 /// The output function for n+plot points in each coordinate direction
274 //=======================================================================
275 template <unsigned NNODE_1D>
276 void QSpectralElement<2,NNODE_1D>::output(std::ostream &outfile,
277  const unsigned &n_plot)
278 {
279  //Local variables
280  Vector<double> s(2);
281 
282  //Tecplot header info
283  outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
284 
285  //Find the dimension of the first node
286  unsigned n_dim = this->nodal_dimension();
287 
288  //Loop over plot points
289  for(unsigned l2=0;l2<n_plot;l2++)
290  {
291  s[1] = -1.0 + l2*2.0/(n_plot-1);
292  for(unsigned l1=0;l1<n_plot;l1++)
293  {
294  s[0] = -1.0 + l1*2.0/(n_plot-1);
295 
296  //Output the coordinates
297  for (unsigned i=0;i<n_dim;i++)
298  {
299  outfile << this->interpolated_x(s,i) << " " ;
300  }
301  outfile << std::endl;
302  }
303  }
304  outfile << std::endl;
305 }
306 
307 
308 
309 
310 //===========================================================
311 /// Function to setup geometrical information for lower-dimensional
312 /// FaceElements (which are of type QSpectralElement<0,1>).
313 //===========================================================
314 template<unsigned NNODE_1D>
316  const int &face_index,
317  FaceElement *face_element_pt)
318 {
319  // Set the nodal dimension from the "bulk"
320  face_element_pt->set_nodal_dimension(this->node_pt(0)->ndim());
321 
322  // Set the pointer to the "bulk" element
323  face_element_pt->bulk_element_pt()=this;
324 
325 #ifdef OOMPH_HAS_MPI
326  // If the bulk element is halo then the face element must be too
327  //if (this->is_halo())
328  {
329  face_element_pt->set_halo(Non_halo_proc_ID);
330  }
331 #endif
332 
333  // Resize the storage for the original number of values at
334  // NNODE_1D nodes of the FaceElement
335  face_element_pt->nbulk_value_resize(NNODE_1D);
336 
337  // Resize the storage for the bulk node numbers corresponding
338  // to the NNODE_1D nodes of the FaceElement
339  face_element_pt->bulk_node_number_resize(NNODE_1D);
340 
341  // Set the face index in the face element
342  //The faces are
343  // +1 East
344  // -1 West
345  // +2 North
346  // -2 South
347 
348  //Set the face index in the face element
349  face_element_pt->face_index() = face_index;
350 
351  //Now set up the node pointers
352  // =================================
353  //The convention here is that interior_tangent X tangent X tangent
354  //is the OUTWARD normal
355  switch(face_index)
356  {
357  unsigned bulk_number;
358  //West face, normal sign is positive
359  case(-1):
360  //Set the pointer to the bulk coordinate translation scheme
361  face_element_pt->face_to_bulk_coordinate_fct_pt() =
363 
364  //Set the pointer to the derivative mappings
365  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
367 
368  for(unsigned i=0;i<NNODE_1D;i++)
369  {
370  bulk_number = i*NNODE_1D;
371  face_element_pt->node_pt(i) = this->node_pt(bulk_number);
372  face_element_pt->bulk_node_number(i) = bulk_number;
373  face_element_pt->normal_sign() = 1;
374  //Set the number of values originally stored at this node
375  face_element_pt->nbulk_value(i) = this->required_nvalue(bulk_number);
376  }
377  break;
378  //South face, normal sign is positive
379  case(-2):
380  //Set the pointer to the bulk coordinate translation scheme
381  face_element_pt->face_to_bulk_coordinate_fct_pt() =
383 
384  //Set the pointer to the derivative mappings
385  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
387 
388  for(unsigned i=0;i<NNODE_1D;i++)
389  {
390  bulk_number = i;
391  face_element_pt->node_pt(i) = this->node_pt(bulk_number);
392  face_element_pt->bulk_node_number(i) = bulk_number;
393  face_element_pt->normal_sign() = 1;
394  //Set the number of values originally stored at this node
395  face_element_pt->nbulk_value(i) = this->required_nvalue(bulk_number);
396  }
397  break;
398  //East face, normal sign is negative
399  case(1):
400  //Set the pointer to the bulk coordinate translation scheme
401  face_element_pt->face_to_bulk_coordinate_fct_pt() =
403 
404  //Set the pointer to the derivative mappings
405  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
407 
408  for(unsigned i=0;i<NNODE_1D;i++)
409  {
410  bulk_number = NNODE_1D*i + NNODE_1D-1;
411  face_element_pt->node_pt(i) = this->node_pt(bulk_number);
412  face_element_pt->bulk_node_number(i) = bulk_number;
413  face_element_pt->normal_sign() = -1;
414  //Set the number of values originally stored at this node
415  face_element_pt->nbulk_value(i)=this->required_nvalue(bulk_number);
416  }
417  break;
418  //North face, normal sign is negative
419  case(2):
420  //Set the pointer to the bulk coordinate translation scheme
421  face_element_pt->face_to_bulk_coordinate_fct_pt() =
423 
424  //Set the pointer to the derivative mappings
425  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
427 
428  for(unsigned i=0;i<NNODE_1D;i++)
429  {
430  bulk_number = NNODE_1D*(NNODE_1D-1) + i;
431  face_element_pt->node_pt(i) = this->node_pt(bulk_number);
432  face_element_pt->bulk_node_number(i) = bulk_number;
433  face_element_pt->normal_sign() = -1;
434  //Set the number of values originally stored at this node
435  face_element_pt->nbulk_value(i)=this->required_nvalue(bulk_number);
436  }
437  break;
438 
439  //Now cover the other cases
440  default:
441  std::ostringstream error_message;
442  error_message << "Face index should only take the values +/- 1 or +/- 2,"
443  << " not " << face_index << std::endl;
444  throw OomphLibError(error_message.str(),
445  OOMPH_CURRENT_FUNCTION,
446  OOMPH_EXCEPTION_LOCATION);
447  }
448 }
449 
450 
451 //////////////////////////////////////////////////////////////////////////
452 /// 3D QLegendreElements
453 //////////////////////////////////////////////////////////////////////////
454 
455 //=======================================================================
456 /// Assign the static integral
457 //=======================================================================
458 template<unsigned NNODE_1D>
460 
461 //=======================================================================
462 /// The output function for general 1D QSpectralElements
463 //=======================================================================
464 template <unsigned NNODE_1D>
465 void QSpectralElement<3,NNODE_1D>::output(std::ostream &outfile)
466 {
467  //Tecplot header info
468  outfile << "ZONE I=" << NNODE_1D << ", J=" << NNODE_1D
469  << ", K=" << NNODE_1D<< std::endl;
470 
471  //Find the dimension of the nodes
472  unsigned n_dim = this->nodal_dimension();
473 
474  //Loop over element nodes
475  for(unsigned l3=0;l3<NNODE_1D;l3++)
476  {
477  for(unsigned l2=0;l2<NNODE_1D;l2++)
478  {
479  for(unsigned l1=0;l1<NNODE_1D;l1++)
480  {
481  unsigned l = l3*NNODE_1D*NNODE_1D+l2*NNODE_1D + l1;
482 
483  //Loop over the dimensions and output the position
484  for(unsigned i=0;i<n_dim;i++)
485  {
486  outfile << this->node_pt(l)->x(i) << " ";
487  }
488  //Find out how many data values at the node
489  unsigned initial_nvalue = this->node_pt(l)->nvalue();
490  //Loop over the data and output whether pinned or not
491  for(unsigned i=0;i<initial_nvalue;i++)
492  {
493  outfile << this->node_pt(l)->is_pinned(i) << " ";
494  }
495  outfile << std::endl;
496  }
497  }
498  }
499  outfile << std::endl;
500 }
501 
502 //=======================================================================
503 /// The output function for n_plot points in each coordinate direction
504 //=======================================================================
505 template <unsigned NNODE_1D>
507 output(std::ostream &outfile, const unsigned &n_plot)
508 {
509  //Local variables
510  Vector<double> s(3);
511 
512  //Tecplot header info
513  outfile << "ZONE I=" << n_plot << ", J=" << n_plot
514  << ", K=" << n_plot << std::endl;
515 
516  //Find the dimension of the first node
517  unsigned n_dim = this->nodal_dimension();
518 
519  //Loop over element nodes
520  for(unsigned l3=0;l3<n_plot;l3++)
521  {
522  s[2] = -1.0 + l3*2.0/(n_plot-1);
523  for(unsigned l2=0;l2<n_plot;l2++)
524  {
525  s[1] = -1.0 + l2*2.0/(n_plot-1);
526  for(unsigned l1=0;l1<n_plot;l1++)
527  {
528  s[0] = -1.0 + l1*2.0/(n_plot-1);
529 
530  //Output the coordinates
531  for (unsigned i=0;i<n_dim;i++)
532  {
533  outfile << this->interpolated_x(s,i) << " " ;
534  }
535  outfile << std::endl;
536  }
537  }
538  }
539  outfile << std::endl;
540 }
541 
542 
543 //=======================================================================
544 /// Function to setup geometrical information for lower-dimensional
545 /// FaceElements (which are of type QElement<3,NNODE_1D>).
546 //=======================================================================
547  template<unsigned NNODE_1D>
549  build_face_element(const int &face_index,
550  FaceElement* face_element_pt)
551 {
552  oomph_info <<" WARNING UNTESTED CODE" << std::endl;
553 
554  // Set the nodal dimension from the "bulk"
555  face_element_pt->set_nodal_dimension(this->node_pt(0)->ndim());
556 
557  //Set the pointer to the orginal "bulk" element
558  face_element_pt->bulk_element_pt()=this;
559 
560 #ifdef OOMPH_HAS_MPI
561  // If the bulk element is halo then the face element must be too
562  //if (this->is_halo())
563  {
564  face_element_pt->set_halo(Non_halo_proc_ID);
565  }
566 #endif
567 
568  // Resize storage for the number of values originally stored
569  // at the face element's NNODE_1D*NNODE_1D nodes.
570  face_element_pt->nbulk_value_resize(NNODE_1D*NNODE_1D);
571 
572  // Set the face index in the element
573  // The faces are
574  // -3 : BACK (OLD: Bottom
575  // -2 : DOWN (OLD: Front
576  // -1 : LEFT (OLD: Left Side
577  // 1 : RIGHT (OLD: Right Side
578  // 2 : UP (OLD: Back
579  // 3 : FRONT (OLD: Top
580 
581  face_element_pt->face_index() = face_index;
582 
583  //Now set up the node pointers and the normal vectors
584  switch(face_index)
585  {
586 
587  // BACK
588  //-----
589  case -3:
590  //Set the pointer to the bulk coordinate translation scheme
591  face_element_pt->face_to_bulk_coordinate_fct_pt() =
593 
594  //Set the pointer to the derivative mappings
595  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
597 
598  // Copy nodes
599  for (unsigned i=0;i<(NNODE_1D*NNODE_1D);i++)
600  {
601  face_element_pt->node_pt(i)= this->node_pt(i);
602  }
603  // Outer unit normal is negative of cross product of two in plane
604  // tangent vectors
605  face_element_pt->normal_sign()=-1;
606 
607  break;
608 
609  // FRONT
610  //------
611  case 3:
612 
613  //Set the pointer to the bulk coordinate translation scheme
614  face_element_pt->face_to_bulk_coordinate_fct_pt() =
616 
617  //Set the pointer to the derivative mappings
618  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
620 
621  // Copy nodes
622  for (unsigned i=0;i<(NNODE_1D*NNODE_1D);i++)
623  {
624  face_element_pt->node_pt(i)
625  = this->node_pt(i+(NNODE_1D*NNODE_1D)*(NNODE_1D-1));
626  }
627  // Outer unit normal is cross product of two in plane
628  // tangent vectors
629  face_element_pt->normal_sign()=1;
630 
631 
632  break;
633 
634  // DOWN:
635  //------
636  case -2:
637 
638  {
639  //Set the pointer to the bulk coordinate translation scheme
640  face_element_pt->face_to_bulk_coordinate_fct_pt() =
642 
643  //Set the pointer to the derivative mappings
644  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
646 
647  // Copy nodes
648  unsigned count=0;
649  for (unsigned i=0;i<NNODE_1D;i++)
650  {
651  for (unsigned j=0;j<NNODE_1D;j++)
652  {
653  face_element_pt->node_pt(count)
654  = this->node_pt(j+i*(NNODE_1D*NNODE_1D));
655  count++;
656  }
657  }
658 
659  // Outer unit normal is cross product of two in plane
660  // tangent vectors
661  face_element_pt->normal_sign()=1;
662  }
663  break;
664 
665 
666  // UP:
667  //----
668  case 2:
669 
670  {
671  //Set the pointer to the bulk coordinate translation scheme
672  face_element_pt->face_to_bulk_coordinate_fct_pt() =
674 
675  //Set the pointer to the derivative mappings
676  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
678 
679  // Copy nodes
680  unsigned count=0;
681  for (unsigned i=0;i<NNODE_1D;i++)
682  {
683  for (unsigned j=0;j<NNODE_1D;j++)
684  {
685  face_element_pt->node_pt(count)=
686  this->node_pt(j+i*(NNODE_1D*NNODE_1D)+(NNODE_1D*(NNODE_1D-1)));
687  count++;
688  }
689  }
690 
691  // Outer unit normal is negative of cross product of two in plane
692  // tangent vectors
693  face_element_pt->normal_sign()=-1;
694  }
695  break;
696 
697  // LEFT:
698  //------
699  case -1:
700 
701  {
702  //Set the pointer to the bulk coordinate translation scheme
703  face_element_pt->face_to_bulk_coordinate_fct_pt() =
705 
706  //Set the pointer to the derivative mappings
707  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
709 
710  // Copy nodes
711  unsigned count=0;
712  for (unsigned i=0;i<NNODE_1D;i++)
713  {
714  for (unsigned j=0;j<NNODE_1D;j++)
715  {
716  unsigned jj=j*NNODE_1D+i*(NNODE_1D*NNODE_1D);
717  face_element_pt->node_pt(count) = this->node_pt(jj);
718  count++;
719  }
720  }
721 
722  // Outer unit normal is negative of cross product of two in plane
723  // tangent vectors
724  face_element_pt->normal_sign()=-1;
725  }
726  break;
727 
728 
729  // RIGHT:
730  //-------
731  case 1:
732 
733  {
734  //Set the pointer to the bulk coordinate translation scheme
735  face_element_pt->face_to_bulk_coordinate_fct_pt() =
737 
738  //Set the pointer to the derivative mappings
739  face_element_pt->bulk_coordinate_derivatives_fct_pt() =
741 
742  // Copy nodes
743  unsigned count=0;
744  for (unsigned i=0;i<NNODE_1D;i++)
745  {
746  for (unsigned j=0;j<NNODE_1D;j++)
747  {
748  unsigned jj=j*NNODE_1D+i*(NNODE_1D*NNODE_1D)+(NNODE_1D-1);
749  face_element_pt->node_pt(count) = this->node_pt(jj);
750  count++;
751  }
752  }
753 
754  // Outer unit normal is cross product of two in plane
755  // tangent vectors
756  face_element_pt->normal_sign()=1;
757  }
758  break;
759 
760 
761  //Cover all other cases
762  default:
763  std::ostringstream error_message;
764  error_message
765  << "Face index should only take the values +/- 1, +/- 2 or +/- 3,"
766  << " not " << face_index << std::endl;
767  throw OomphLibError(error_message.str(),
768  OOMPH_CURRENT_FUNCTION,
769  OOMPH_EXCEPTION_LOCATION);
770  } //end switch
771 }
772 
773 
774 template class QSpectralElement<1,2>;
775 template class QSpectralElement<1,3>;
776 template class QSpectralElement<1,4>;
777 template class QSpectralElement<1,5>;
778 template class QSpectralElement<1,6>;
779 template class QSpectralElement<1,7>;
780 template class QSpectralElement<1,8>;
781 template class QSpectralElement<1,9>;
782 template class QSpectralElement<1,10>;
783 template class QSpectralElement<1,11>;
784 template class QSpectralElement<1,12>;
785 template class QSpectralElement<1,13>;
786 template class QSpectralElement<1,14>;
787 
788 template class QSpectralElement<2,2>;
789 template class QSpectralElement<2,3>;
790 template class QSpectralElement<2,4>;
791 template class QSpectralElement<2,5>;
792 template class QSpectralElement<2,6>;
793 template class QSpectralElement<2,7>;
794 template class QSpectralElement<2,8>;
795 
796 template class QSpectralElement<3,2>;
797 template class QSpectralElement<3,3>;
798 template class QSpectralElement<3,4>;
799 template class QSpectralElement<3,5>;
800 template class QSpectralElement<3,6>;
801 template class QSpectralElement<3,7>;
802 template class QSpectralElement<3,8>;
803 
804 }
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the south face (s1 = -1.0)
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = -1.0.
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = 1.0.
cstr elem_len * i
Definition: cfortran.h:607
void face2(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the east face (s0 = 1.0)
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
void face2(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the back face (s2 = -1.0)
void faces0(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the east and west faces, along which s0 is fixed.
void faces1(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the up and down faces, along which s1 is fixed.
void face4(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the up face (s1 = 1.0)
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
static std::map< unsigned, Vector< double > > z
OomphInfo oomph_info
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 faces1(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the north and south faces, along which s1 is fixed.
void faces0(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the back and front faces, along which s0 is fixed.
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
static char t char * s
Definition: cfortran.h:572
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
static GaussLobattoLegendre< 1, NNODE_1D > integral
Default integration rule: Gaussian integration of same 'order' as the element.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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.
void faces0(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for both faces – the bulk coordinate is fixed on both.
void faces2(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the left and right faces, along which s2 is fixed.
void face3(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the right face (s0 = 1.0)
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the west face (s0 = -1.0)
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4470
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
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the down face (s1 = -1.0)
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the left face (s0 = -1.0)
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 face5(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the front face (s2 = 1.0)
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 face3(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the north face (s1 = 1.0)