Qspectral_elements.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: 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 //Header functions for classes that define Qelements
31 
32 //Include guards to prevent multiple inclusion of the header
33 #ifndef OOMPH_QSPECTRAL_ELEMENT_HEADER
34 #define OOMPH_QSPECTRAL_ELEMENT_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38  #include <oomph-lib-config.h>
39 #endif
40 
41 //oomph-lib headers
42 #include "Qelements.h"
43 
44 namespace oomph
45 {
46 
47 //=====================================================================
48 /// Class that returns the shape functions associated with legendre
49 //=====================================================================
51 {
52  public:
53  static std::map<unsigned,Vector<double> > z;
54 
55  /// Static function used to populate the stored positions
56  static inline void calculate_nodal_positions(const unsigned &order)
57  {
58  //If we haven't already calculated these
59  if(z.find(order)==z.end())
60  {
61  Orthpoly::gll_nodes(order,z[order]);
62  }
63  }
64 
65  static inline double nodal_position(const unsigned &order, const unsigned &n)
66  {return z[order][n];}
67 
68  /// Constructor
69  OneDLegendreShapeParam(const unsigned &order,const double &s) :
70  Shape(order)
71  {
72  using namespace Orthpoly;
73 
74  unsigned p = order-1;
75  //Now populate the shape function
76  for(unsigned i=0;i<order;i++)
77  {
78  //If we're at one of the nodes, the value must be 1.0
79  if(std::fabs(s-z[order][i]) < Orthpoly::eps)
80  {
81  (*this)[i] = 1.0;
82  }
83  //Otherwise use the lagrangian interpolant
84  else
85  {
86  (*this)[i] = (1.0 - s*s)*dlegendre(p,s)/
87  (p*(p+1)*legendre(p,z[order][i])*(z[order][i] - s));
88  }
89  }
90  }
91 };
92 
93 
95 {
96  public:
97  // Constructor
98  OneDLegendreDShapeParam(const unsigned &order, const double &s) :
99  Shape(order)
100  {
101  unsigned p = order - 1;
103 
104  bool root=false;
105  for(unsigned i=0;i<order;i++)
106  {
107  unsigned rootnum=0;
108  for(unsigned j=0;j<order;j++)
109  { // Loop over roots to check if
110  if(std::fabs(s-z[j])<10.0*Orthpoly::eps)
111  { // s happens to be a root.
112  root=true;
113  break;
114  }
115  rootnum+=1;
116  }
117  if(root==true)
118  {
119  if(i==rootnum && i==0)
120  {
121  (*this)[i]=-(1.0+p)*p/4.0;
122  }
123  else if(i==rootnum && i==p)
124  {
125  (*this)[i]=(1.0+p)*p/4.0;
126  }
127  else if(i==rootnum)
128  {
129  (*this)[i]=0.0;
130  }
131  else
132  {
133  (*this)[i]=Orthpoly::legendre(p,z[rootnum])
134  /Orthpoly::legendre(p,z[i])/(z[rootnum]-z[i]);
135  }
136  }
137  else
138  {
139  (*this)[i]=((1+s*(s-2*z[i]))/(s-z[i])* Orthpoly::dlegendre(p,s)
140  -(1-s*s)* Orthpoly::ddlegendre(p,s))
141  /p/(p+1.0)/Orthpoly::legendre(p,z[i])/(s-z[i]);
142  }
143  root = false;
144  }
145 
146 
147  }
148 
149 };
150 
151 
152 
153 //========================================================================
154 // A Base class for Spectral elements
155 //========================================================================
156 class SpectralElement : public virtual FiniteElement
157 {
158  protected:
159 
160  /// Additional Storage for shared spectral data
162 
163  /// Vector that represents the spectral order in each dimension
165 
166  /// Vector that represents the nodal spectral order
168 
169  /// Local equation numbers for the spectral degrees of freedom
171 
172  public:
173 
174  SpectralElement() : FiniteElement(), Spectral_data_pt(0) {}
175 
177  {
178  if(Spectral_data_pt!=0)
179  {
180  delete Spectral_data_pt;
181  Spectral_data_pt=0;
182  }
183  }
184 
185 
186  ///\short Return the i-th data object associated with the polynomials
187  ///of order p. Note that i <= p.
188  Data* spectral_data_pt(const unsigned &i) const
189  {return (*Spectral_data_pt)[i];}
190 
191  /// \short Function to describe the local dofs of the element. The ostream
192  /// specifies the output stream to which the description
193  /// is written; the string stores the currently
194  /// assembled output that is ultimately written to the
195  /// output stream by Data::describe_dofs(...); it is typically
196  /// built up incrementally as we descend through the
197  /// call hierarchy of this function when called from
198  /// Problem::describe_dofs(...)
199  virtual void describe_local_dofs(std::ostream& out,
200  const std::string& current_string) const
201  {
202  //Do the standard finite element stuff
203  FiniteElement::describe_local_dofs(out,current_string);
204  //Now get the number of spectral data.
205  unsigned n_spectral = nspectral();
206 
207  //Now loop over the nodes again and assign local equation numbers
208  for(unsigned n=0;n<n_spectral;n++)
209  {
210  //Can we upcast to a node
211  Node* cast_node_pt = dynamic_cast<Node*>(spectral_data_pt(n));
212  if(cast_node_pt)
213  {
214  std::stringstream conversion;
215  conversion <<" of Node "<<n<<current_string;
216  std::string in(conversion.str());
217  cast_node_pt->describe_dofs(out,in);
218  }
219  // Otherwise it is data.
220  else
221  {
222  Data* data_pt = spectral_data_pt(n);
223  std::stringstream conversion;
224  conversion <<" of Data "<<n<<current_string;
225  std::string in(conversion.str());
226  data_pt->describe_dofs(out,in);
227  }
228  }
229  }
230 
231  /// \short Assign the local equation numbers. If the boolean argument is
232  /// true then store degrees of freedom at Dof_pt
234  const bool &store_local_dof_pt)
235  {
236  //Do the standard finite element stuff
238 
239  //Now need to loop over the spectral data
240  unsigned n_spectral = nspectral();
241  if(n_spectral > 0)
242  {
243  //Find the number of local equations assigned so far
244  unsigned local_eqn_number = ndof();
245 
246  //Find number of values stored at the first node
247  unsigned max_n_value = spectral_data_pt(0)->nvalue();
248  //Loop over the other nodes and find the maximum number of values stored
249  for(unsigned n=1;n<n_spectral;n++)
250  {
251  unsigned n_value = spectral_data_pt(n)->nvalue();
252  if(n_value > max_n_value) {max_n_value = n_value;}
253  }
254 
255  //Resize the storage for the nodal local equation numbers
256  //initially all local equations are unclassified
257  Spectral_local_eqn.resize(n_spectral,max_n_value,Data::Is_unclassified);
258 
259  //Construct a set of pointers to the nodes
260  std::set<Data*> set_of_node_pt;
261  unsigned n_node = nnode();
262  for(unsigned n=0;n<n_node;n++) {set_of_node_pt.insert(node_pt(n));}
263 
264  //A local queue to store the global equation numbers
265  std::deque<unsigned long> global_eqn_number_queue;
266 
267  //Now loop over the nodes again and assign local equation numbers
268  for(unsigned n=0;n<n_spectral;n++)
269  {
270  //Need to find whether the data is, in fact a node
271  //(and hence already assigned)
272 
273  //Can we upcast to a node
274  Node* cast_node_pt = dynamic_cast<Node*>(spectral_data_pt(n));
275  if((cast_node_pt) &&
276  (set_of_node_pt.find(cast_node_pt)!=set_of_node_pt.end()))
277  {
278  //Find the number of values
279  unsigned n_value = cast_node_pt->nvalue();
280  //Copy them across
281  for(unsigned j=0;j<n_value;j++)
282  {
283  Spectral_local_eqn(n,j) =
284  nodal_local_eqn(get_node_number(cast_node_pt),j);
285  }
286  //Remove from the set
287  set_of_node_pt.erase(cast_node_pt);
288  }
289  //Otherwise it's just data
290  else
291  {
292  Data* const data_pt = spectral_data_pt(n);
293  unsigned n_value = data_pt->nvalue();
294  //Loop over the number of values
295  for(unsigned j=0;j<n_value;j++)
296  {
297  //Get the GLOBAL equation number
298  long eqn_number = data_pt->eqn_number(j);
299  //If the GLOBAL equation number is positive (a free variable)
300  if(eqn_number >= 0)
301  {
302  //Add the GLOBAL equation number to the
303  //local-to-global translation
304  //scheme
305  global_eqn_number_queue.push_back(eqn_number);
306  //Add pointer to the dof to the queue if required
307  if(store_local_dof_pt)
308  {
310  data_pt->value_pt(j));
311  }
312  //Add the local equation number to the local scheme
313  Spectral_local_eqn(n,j) = local_eqn_number;
314  //Increase the local number
315  local_eqn_number++;
316  }
317  else
318  {
319  //Set the local scheme to be pinned
320  Spectral_local_eqn(n,j) = Data::Is_pinned;
321  }
322  }
323  } //End of case when it's not a nodal dof
324  }
325 
326  //Now add our global equations numbers to the internal element storage
327  add_global_eqn_numbers(global_eqn_number_queue,
329  //Clear the memory used in the deque
330  if(store_local_dof_pt)
331  {std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);}
332 
333  } //End of case when there are spectral degrees of freedom
334  }
335 
336  unsigned nspectral() const
337  {
338  if(Spectral_data_pt==0) {return 0;}
339  else
340  {
341  return Spectral_data_pt->size();
342  }
343  }
344 
345 };
346 
347 
348 
349 //=======================================================================
350 ///General QLegendreElement class
351 ///
352 /// Empty, just establishes the template parameters
353 //=======================================================================
354 template<unsigned DIM, unsigned NNODE_1D>
356 {
357 };
358 
359 
360 //=======================================================================
361 ///General QSpectralElement class specialised to one spatial dimension
362 //=======================================================================
363 template<unsigned NNODE_1D>
364 class QSpectralElement<1,NNODE_1D> : public virtual SpectralElement,
365  public virtual LineElementBase
366 {
367  private:
368 
369  /// \short Default integration rule: Gaussian integration of same 'order'
370  /// as the element
371  //This is sort of optimal, because it means that the integration is exact
372  //for the shape functions. Can overwrite this in specific element defintion.
374 
375 public:
376 
377  /// Constructor
379  {
380  //There are NNODE_1D nodes in this element
381  this->set_n_node(NNODE_1D);
382  Spectral_order.resize(1,NNODE_1D);
383  Nodal_spectral_order.resize(1,NNODE_1D);
384  //Set the elemental and nodal dimensions
385  this->set_dimension(1);
386  //Assign integral point
387  this->set_integration_scheme(&integral);
388  //Calculate the nodal positions for the shape functions
390  //OneDLegendreShapeParam::calculate_nodal_positions(NNODE_1D);
391  }
392 
393  /// Min. value of local coordinate
394  double s_min() const {return -1.0;}
395 
396  /// Max. value of local coordinate
397  double s_max() const {return 1.0;}
398 
399  /// Number of vertex nodes in the element
400  unsigned nvertex_node() const
401  { return 2; }
402 
403  /// Pointer to the j-th vertex node in the element
404  Node* vertex_node_pt(const unsigned &j) const
405  {
406  unsigned n_node_1d = nnode_1d();
407  Node* nod_pt;
408  switch(j)
409  {
410  case 0:
411  nod_pt = this->node_pt(0);
412  break;
413  case 1:
414  nod_pt = this->node_pt(n_node_1d-1);
415  break;
416  default:
417  std::ostringstream error_message;
418  error_message << "Vertex node number is " << j <<
419  " but must be from 0 to 1\n";
420 
421  throw OomphLibError(error_message.str(),
422  OOMPH_CURRENT_FUNCTION,
423  OOMPH_EXCEPTION_LOCATION);
424  }
425  return nod_pt;
426  }
427 
428  /// Get local coordinates of node j in the element; vector sets its own size
429  void local_coordinate_of_node(const unsigned& n, Vector<double>& s) const
430  {
431  s.resize(1);
433  }
434 
435  /// Get the local fractino of node j in the element
436  void local_fraction_of_node(const unsigned &n, Vector<double> &s_fraction)
437  {
438  this->local_coordinate_of_node(n,s_fraction);
439  //Resize
440  s_fraction[0] = 0.5*(s_fraction[0] + 1.0);
441  }
442 
443  /// The local one-d fraction is the same
444  double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
445  {
446  return
448  }
449 
450  /// Calculate the geometric shape functions at local coordinate s
451  inline void shape(const Vector<double> &s, Shape &psi) const;
452 
453  /// \short Compute the geometric shape functions and
454  /// derivatives w.r.t. local coordinates at local coordinate s
455  inline void dshape_local(const Vector<double> &s, Shape &psi, DShape &dpsids)
456  const;
457 
458  /// \short Compute the geometric shape functions, derivatives and
459  /// second derivatives w.r.t. local coordinates at local coordinate s
460  /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
461  inline void d2shape_local(const Vector<double> &s, Shape &psi, DShape &dpsids,
462  DShape &d2psids) const;
463 
464  /// \short Overload the template-free interface for the calculation of
465  /// inverse jacobian matrix. This is a one-dimensional element, so
466  /// use the 1D version.
468  DenseMatrix<double> &inverse_jacobian) const
469  {return FiniteElement::invert_jacobian<1>(jacobian,inverse_jacobian);}
470 
471 
472  /// Number of nodes along each element edge
473  unsigned nnode_1d() const {return NNODE_1D;}
474 
475  /// C-style output
476  void output(FILE* file_pt)
477  {FiniteElement::output(file_pt);}
478 
479  /// C_style output at n_plot points
480  void output(FILE* file_pt, const unsigned &n_plot)
481  {FiniteElement::output(file_pt,n_plot);}
482 
483  /// Output
484  void output(std::ostream &outfile);
485 
486  /// Output at n_plot points
487  void output(std::ostream &outfile, const unsigned &n_plot);
488 
489 /// \short Get cector of local coordinates of plot point i (when plotting
490  /// nplot points in each "coordinate direction).
491  void get_s_plot(const unsigned& i, const unsigned& nplot,
492  Vector<double>& s) const
493  {
494  if (nplot>1)
495  {
496  s[0]=-1.0+2.0*double(i)/double(nplot-1);
497  }
498  else
499  {
500  s[0]=0.0;
501  }
502  }
503 
504  /// \short Return string for tecplot zone header (when plotting
505  /// nplot points in each "coordinate direction)
506  std::string tecplot_zone_string(const unsigned& nplot) const
507  {
508  std::ostringstream header;
509  header << "ZONE I=" << nplot << "\n";
510  return header.str();
511  }
512 
513  /// \short Return total number of plot points (when plotting
514  /// nplot points in each "coordinate direction)
515  unsigned nplot_points(const unsigned& nplot) const
516  {
517  unsigned DIM=1;
518  unsigned np=1;
519  for (unsigned i=0;i<DIM;i++) {np*=nplot;}
520  return np;
521  }
522 
523  /// \short Build the lower-dimensional FaceElement (an element of type
524  /// QSpectralElement<0,NNODE_1D>). The face index takes two values
525  /// corresponding to the two possible faces:
526  /// -1 (Left) s[0] = -1.0
527  /// +1 (Right) s[0] = 1.0
528  void build_face_element(const int &face_index,
529  FaceElement* face_element_pt);
530 
531 };
532 
533 
534 //=======================================================================
535 ///Shape function for specific QSpectralElement<1,..>
536 //=======================================================================
537 template <unsigned NNODE_1D>
539  const
540 {
541  //Call the OneDimensional Shape functions
543  //OneDLegendreShapeParam psi1(NNODE_1D,s[0]);
544 
545  //Now let's loop over the nodal points in the element
546  //and copy the values back in
547  for(unsigned i=0;i<NNODE_1D;i++) {psi(i) = psi1[i];}
548 }
549 
550 //=======================================================================
551 ///Derivatives of shape functions for specific QSpectralElement<1,..>
552 //=======================================================================
553 template <unsigned NNODE_1D>
555  Shape &psi,
556  DShape &dpsids) const
557 {
558  //Call the shape functions and derivatives
561 
562  //OneDLegendreShapeParam psi1(NNODE_1D,s[0]);
563  //OneDLegendreDShapeParam dpsi1ds(NNODE_1D,s[0]);
564 
565  //Loop over shape functions in element and assign to psi
566  for(unsigned l=0;l<NNODE_1D;l++)
567  {
568  psi(l) = psi1[l];
569  dpsids(l,0) = dpsi1ds[l];
570  }
571 }
572 
573 //=======================================================================
574 /// Second derivatives of shape functions for specific QSpectralElement<1,..>
575 /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
576 //=======================================================================
577 template <unsigned NNODE_1D>
579  DShape &dpsids, DShape &d2psids) const
580 {
581  std::ostringstream error_message;
582  error_message <<"\nd2shpe_local currently not implemented for this element\n";
583  throw OomphLibError(error_message.str(),
584  OOMPH_CURRENT_FUNCTION,
585  OOMPH_EXCEPTION_LOCATION);
586 
587 /* //Call the shape functions and derivatives */
588 /* OneDimensionalShape<NNODE_1D> psi1(s[0]); */
589 /* OneDimensionalDShape<NNODE_1D> dpsi1ds(s[0]); */
590 /* OneDimensionalD2Shape<NNODE_1D> d2psi1ds(s[0]); */
591 
592 /* //Loop over shape functions in element and assign to psi */
593 /* for(unsigned l=0;l<NNODE_1D;l++) */
594 /* { */
595 /* psi[l] = psi1[l]; */
596 /* dpsids[l][0] = dpsi1ds[l]; */
597 /* d2psids[l][0] = d2psi1ds[l]; */
598 /* } */
599 }
600 
601 
602 //=======================================================================
603 ///General QSpectralElement class specialised to two spatial dimensions
604 //=======================================================================
605 template<unsigned NNODE_1D>
606 class QSpectralElement<2,NNODE_1D> : public virtual SpectralElement,
607  public virtual QuadElementBase
608 {
609  private:
610 
611  /// \short Default integration rule: Gaussian integration of same 'order'
612  /// as the element
613  //This is sort of optimal, because it means that the integration is exact
614  //for the shape functions. Can overwrite this in specific element defintion.
616 
617 public:
618 
619  /// Constructor
621  {
622  //There are NNODE_1D^2 nodes in this element
623  this->set_n_node(NNODE_1D*NNODE_1D);
624  Spectral_order.resize(2,NNODE_1D);
625  Nodal_spectral_order.resize(2,NNODE_1D);
626  //Set the elemental and nodal dimensions
627  this->set_dimension(2);
628  //Assign integral pointer
629  this->set_integration_scheme(&integral);
630  //Calculate the nodal positions for the shape functions
632  //OneDLegendreShapeParam::calculate_nodal_positions(NNODE_1D);
633  }
634 
635  /// Min. value of local coordinate
636  double s_min() const {return -1.0;}
637 
638  /// Max. value of local coordinate
639  double s_max() const {return 1.0;}
640 
641  /// Number of vertex nodes in the element
642  unsigned nvertex_node() const
643  { return 4;}
644 
645  /// Pointer to the j-th vertex node in the element
646  Node* vertex_node_pt(const unsigned &j) const
647  {
648  unsigned n_node_1d = nnode_1d();
649  Node* nod_pt;
650  switch(j)
651  {
652  case 0:
653  nod_pt = this->node_pt(0);
654  break;
655  case 1:
656  nod_pt = this->node_pt(n_node_1d-1);
657  break;
658  case 2:
659  nod_pt = this->node_pt(n_node_1d*(n_node_1d-1));
660  break;
661  case 3:
662  nod_pt = this->node_pt(n_node_1d*n_node_1d-1);
663  break;
664  default:
665  std::ostringstream error_message;
666  error_message << "Vertex node number is " << j <<
667  " but must be from 0 to 3\n";
668 
669  throw OomphLibError(error_message.str(),
670 OOMPH_CURRENT_FUNCTION,
671  OOMPH_EXCEPTION_LOCATION);
672  }
673  return nod_pt;
674  }
675 
676 
677  /// Get local coordinates of node j in the element; vector sets its own size
678  void local_coordinate_of_node(const unsigned& n, Vector<double>& s) const
679  {
680  s.resize(2);
681  unsigned n0 = n%NNODE_1D;
682  unsigned n1 = unsigned(double(n)/double(NNODE_1D));
685  }
686 
687  /// Get the local fractino of node j in the element
688  void local_fraction_of_node(const unsigned &n, Vector<double> &s_fraction)
689  {
690  this->local_coordinate_of_node(n,s_fraction);
691  //Resize
692  s_fraction[0] = 0.5*(s_fraction[0] + 1.0);
693  s_fraction[1] = 0.5*(s_fraction[1] + 1.0);
694  }
695 
696  /// The local one-d fraction is the same
697  double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
698  {
699  return
701  }
702 
703  /// Calculate the geometric shape functions at local coordinate s
704  inline void shape(const Vector<double> &s, Shape &psi) const;
705 
706  /// \short Compute the geometric shape functions and
707  /// derivatives w.r.t. local coordinates at local coordinate s
708  inline void dshape_local(const Vector<double> &s, Shape &psi, DShape &dpsids)
709  const;
710 
711  /// \short Compute the geometric shape functions, derivatives and
712  /// second derivatives w.r.t. local coordinates at local coordinate s
713  /// d2psids(i,0) = \f$ \partial ^2 \psi_j / \partial s_0^2 \f$
714  /// d2psids(i,1) = \f$ \partial ^2 \psi_j / \partial s_1^2 \f$
715  /// d2psids(i,2) = \f$ \partial ^2 \psi_j / \partial s_0 \partial s_1 \f$
716  inline void d2shape_local(const Vector<double> &s, Shape &psi, DShape &dpsids,
717  DShape &d2psids) const;
718 
719  /// \short Overload the template-free interface for the calculation of
720  /// inverse jacobian matrix. This is a one-dimensional element, so
721  /// use the 1D version.
723  DenseMatrix<double> &inverse_jacobian) const
724  {return FiniteElement::invert_jacobian<2>(jacobian,inverse_jacobian);}
725 
726 
727  /// Number of nodes along each element edge
728  unsigned nnode_1d() const {return NNODE_1D;}
729 
730  /// C-style output
731  void output(FILE* file_pt)
732  {FiniteElement::output(file_pt);}
733 
734  /// C_style output at n_plot points
735  void output(FILE* file_pt, const unsigned &n_plot)
736  {FiniteElement::output(file_pt,n_plot);}
737 
738  /// Output
739  void output(std::ostream &outfile);
740 
741  /// Output at n_plot points
742  void output(std::ostream &outfile, const unsigned &n_plot);
743 
744 /// \short Get cector of local coordinates of plot point i (when plotting
745  /// nplot points in each "coordinate direction).
746  void get_s_plot(const unsigned& i, const unsigned& nplot,
747  Vector<double>& s) const
748  {
749  if (nplot>1)
750  {
751  unsigned i0=i%nplot;
752  unsigned i1=(i-i0)/nplot;
753 
754  s[0]=-1.0+2.0*double(i0)/double(nplot-1);
755  s[1]=-1.0+2.0*double(i1)/double(nplot-1);
756  }
757  else
758  {
759  s[0]=0.0;
760  s[1]=0.0;
761  }
762  }
763 
764 
765  /// \short Return string for tecplot zone header (when plotting
766  /// nplot points in each "coordinate direction)
767  std::string tecplot_zone_string(const unsigned& nplot) const
768  {
769  std::ostringstream header;
770  header << "ZONE I=" << nplot << ", J=" << nplot << "\n";
771  return header.str();
772  }
773 
774  /// \short Return total number of plot points (when plotting
775  /// nplot points in each "coordinate direction)
776  unsigned nplot_points(const unsigned& nplot) const
777  {
778  unsigned DIM=2;
779  unsigned np=1;
780  for (unsigned i=0;i<DIM;i++) {np*=nplot;}
781  return np;
782  }
783 
784  /// \short Build the lower-dimensional FaceElement (an element of type
785  /// QSpectralElement<1,NNODE_1D>). The face index takes four values
786  /// corresponding to the four possible faces:
787  /// -1 (West) s[0] = -1.0
788  /// +1 (East) s[0] = 1.0
789  /// -2 (South) s[1] = -1.0
790  /// +2 (North) s[1] = 1.0
791  void build_face_element(const int &face_index,
792  FaceElement* face_element_pt);
793 };
794 
795 
796 //=======================================================================
797 ///Shape function for specific QSpectralElement<2,..>
798 //=======================================================================
799 template <unsigned NNODE_1D>
801  const
802 {
803  //Call the OneDimensional Shape functions
804  OneDimensionalLegendreShape<NNODE_1D> psi1(s[0]), psi2(s[1]);
805  //OneDLegendreShapeParam psi1(NNODE_1D,s[0]), psi2(NNODE_1D,s[1]);
806 
807  //Now let's loop over the nodal points in the element
808  //and copy the values back in
809  for(unsigned i=0;i<NNODE_1D;i++)
810  {
811  for(unsigned j=0;j<NNODE_1D;j++)
812  {
813  psi(NNODE_1D*i + j) = psi2[i]*psi1[j];
814  }
815  }
816 }
817 
818 //=======================================================================
819 ///Derivatives of shape functions for specific QSpectralElement<2,..>
820 //=======================================================================
821 template <unsigned NNODE_1D>
823  Shape &psi,
824  DShape &dpsids) const
825 {
826  //Call the shape functions and derivatives
827  OneDimensionalLegendreShape<NNODE_1D> psi1(s[0]), psi2(s[1]);
828  OneDimensionalLegendreDShape<NNODE_1D> dpsi1ds(s[0]), dpsi2ds(s[1]);
829 
830  //OneDLegendreShapeParam psi1(NNODE_1D,s[0]), psi2(NNODE_1D,s[1]);
831  //OneDLegendreDShapeParam dpsi1ds(NNODE_1D,s[0]), dpsi2ds(NNODE_1D,s[1]);
832 
833  //Index for the shape functions
834  unsigned index=0;
835  //Loop over shape functions in element
836  for(unsigned i=0;i<NNODE_1D;i++)
837  {
838  for(unsigned j=0;j<NNODE_1D;j++)
839  {
840  //Assign the values
841  dpsids(index,0) = psi2[i]*dpsi1ds[j];
842  dpsids(index,1) = dpsi2ds[i]*psi1[j];
843  psi[index] = psi2[i]*psi1[j];
844  //Increment the index
845  ++index;
846  }
847  }
848 }
849 
850 
851 //=======================================================================
852 ///Second derivatives of shape functions for specific QSpectralElement<2,..>
853 /// d2psids(i,0) = \f$ \partial ^2 \psi_j / \partial s_0^2 \f$
854 /// d2psids(i,1) = \f$ \partial ^2 \psi_j / \partial s_1^2 \f$
855 /// d2psids(i,2) = \f$ \partial ^2 \psi_j / \partial s_0 \partial s_1 \f$
856 //=======================================================================
857 template <unsigned NNODE_1D>
859  DShape &dpsids, DShape &d2psids) const
860 {
861  std::ostringstream error_message;
862  error_message <<"\nd2shpe_local currently not implemented for this element\n";
863  throw OomphLibError(error_message.str(),
864  OOMPH_CURRENT_FUNCTION,
865  OOMPH_EXCEPTION_LOCATION);
866 
867 /* //Call the shape functions and derivatives */
868 /* OneDimensionalShape<NNODE_1D> psi1(s[0]); */
869 /* OneDimensionalDShape<NNODE_1D> dpsi1ds(s[0]); */
870 /* OneDimensionalD2Shape<NNODE_1D> d2psi1ds(s[0]); */
871 
872 /* //Loop over shape functions in element and assign to psi */
873 /* for(unsigned l=0;l<NNODE_1D;l++) */
874 /* { */
875 /* psi[l] = psi1[l]; */
876 /* dpsids[l][0] = dpsi1ds[l]; */
877 /* d2psids[l][0] = d2psi1ds[l]; */
878 /* } */
879 }
880 
881 
882 //=======================================================================
883 ///General QSpectralElement class specialised to three spatial dimensions
884 //=======================================================================
885 template<unsigned NNODE_1D>
886 class QSpectralElement<3,NNODE_1D> : public virtual SpectralElement,
887  public virtual BrickElementBase
888 {
889  private:
890 
891  /// \short Default integration rule: Gaussian integration of same 'order'
892  /// as the element
893  //This is sort of optimal, because it means that the integration is exact
894  //for the shape functions. Can overwrite this in specific element defintion.
896 
897 public:
898 
899  /// Constructor
901  {
902  //There are NNODE_1D^3 nodes in this element
903  this->set_n_node(NNODE_1D*NNODE_1D*NNODE_1D);
904  Spectral_order.resize(3,NNODE_1D);
905  Nodal_spectral_order.resize(3,NNODE_1D);
906  //Set the elemental and nodal dimensions
907  this->set_dimension(3);
908  //Assign integral point
909  this->set_integration_scheme(&integral);
910  //Calculate the nodal positions for the shape functions
912  //OneDLegendreShapeParam::calculate_nodal_positions(NNODE_1D);
913  }
914 
915  /// Min. value of local coordinate
916  double s_min() const {return -1.0;}
917 
918  /// Max. value of local coordinate
919  double s_max() const {return 1.0;}
920 
921  /// Number of vertex nodes in the element
922  unsigned nvertex_node() const
923  { return 8;}
924 
925  /// Pointer to the j-th vertex node in the element
926  Node* vertex_node_pt(const unsigned &j) const
927  {
928  unsigned n_node_1d = nnode_1d();
929  Node* nod_pt;
930  switch(j)
931  {
932  case 0:
933  nod_pt = this->node_pt(0);
934  break;
935  case 1:
936  nod_pt=this->node_pt(n_node_1d-1);
937  break;
938  case 2:
939  nod_pt=this->node_pt(n_node_1d*(n_node_1d-1));
940  break;
941  case 3:
942  nod_pt=this->node_pt(n_node_1d*n_node_1d-1);
943  break;
944  case 4:
945  nod_pt=this->node_pt(n_node_1d*n_node_1d*(n_node_1d-1));
946  break;
947  case 5:
948  nod_pt=this->node_pt(n_node_1d*n_node_1d*(n_node_1d-1)+(n_node_1d-1));
949  break;
950  case 6:
951  nod_pt=this->node_pt(n_node_1d*n_node_1d*n_node_1d-n_node_1d);
952  break;
953  case 7:
954  nod_pt=this->node_pt(n_node_1d*n_node_1d*n_node_1d-1);
955  break;
956  default:
957  std::ostringstream error_message;
958  error_message << "Vertex node number is " << j <<
959  " but must be from 0 to 7\n";
960 
961  throw OomphLibError(error_message.str(),
962 OOMPH_CURRENT_FUNCTION,
963  OOMPH_EXCEPTION_LOCATION);
964  }
965  return nod_pt;
966  }
967 
968 
969  /// Get local coordinates of node j in the element; vector sets its own size
970  void local_coordinate_of_node(const unsigned& n, Vector<double>& s) const
971  {
972  s.resize(3);
973  unsigned n0 = n%NNODE_1D;
974  unsigned n1 = unsigned(double(n)/double(NNODE_1D))%NNODE_1D;
975  unsigned n2 = unsigned(double(n)/double(NNODE_1D*NNODE_1D));
979  }
980 
981  /// Get the local fractino of node j in the element
982  void local_fraction_of_node(const unsigned &n, Vector<double> &s_fraction)
983  {
984  this->local_coordinate_of_node(n,s_fraction);
985  //Resize
986  s_fraction[0] = 0.5*(s_fraction[0] + 1.0);
987  s_fraction[1] = 0.5*(s_fraction[1] + 1.0);
988  s_fraction[2] = 0.5*(s_fraction[2] + 1.0);
989  }
990 
991  /// The local one-d fraction is the same
992  double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
993  {
994  return
996  }
997 
998  /// Calculate the geometric shape functions at local coordinate s
999  inline void shape(const Vector<double> &s, Shape &psi) const;
1000 
1001  /// \short Compute the geometric shape functions and
1002  /// derivatives w.r.t. local coordinates at local coordinate s
1003  inline void dshape_local(const Vector<double> &s, Shape &psi, DShape &dpsids)
1004  const;
1005 
1006  /// \short Compute the geometric shape functions, derivatives and
1007  /// second derivatives w.r.t. local coordinates at local coordinate s
1008  /// d2psids(i,0) = \f$ \partial ^2 \psi_j / \partial s_0^2 \f$
1009  /// d2psids(i,1) = \f$ \partial ^2 \psi_j / \partial s_1^2 \f$
1010  /// d2psids(i,2) = \f$ \partial ^2 \psi_j / \partial s_2^2 \f$
1011  /// d2psids(i,3) = \f$ \partial ^2 \psi_j / \partial s_0 \partial s_1 \f$
1012  /// d2psids(i,4) = \f$ \partial ^2 \psi_j / \partial s_0 \partial s_2 \f$
1013  /// d2psids(i,5) = \f$ \partial ^2 \psi_j / \partial s_1 \partial s_2 \f$
1014  inline void d2shape_local(const Vector<double> &s, Shape &psi, DShape &dpsids,
1015  DShape &d2psids) const;
1016 
1017 
1018  /// \short Overload the template-free interface for the calculation of
1019  /// inverse jacobian matrix. This is a one-dimensional element, so
1020  /// use the 3D version.
1022  DenseMatrix<double> &inverse_jacobian) const
1023  {return FiniteElement::invert_jacobian<3>(jacobian,inverse_jacobian);}
1024 
1025  /// Number of nodes along each element edge
1026  unsigned nnode_1d() const {return NNODE_1D;}
1027 
1028  /// C-style output
1029  void output(FILE* file_pt)
1030  {FiniteElement::output(file_pt);}
1031 
1032  /// C_style output at n_plot points
1033  void output(FILE* file_pt, const unsigned &n_plot)
1034  {FiniteElement::output(file_pt,n_plot);}
1035 
1036  /// Output
1037  void output(std::ostream &outfile);
1038 
1039  /// Output at nplot points
1040  void output(std::ostream &outfile, const unsigned& nplot);
1041 
1042 /// \short Get cector of local coordinates of plot point i (when plotting
1043  /// nplot points in each "coordinate direction).
1044  void get_s_plot(const unsigned& i, const unsigned& nplot,
1045  Vector<double>& s) const
1046  {
1047  if (nplot>1)
1048  {
1049  unsigned i01=i%(nplot*nplot);
1050  unsigned i0=i01%nplot;
1051  unsigned i1=(i01-i0)/nplot;
1052  unsigned i2=(i-i01)/(nplot*nplot);
1053 
1054  s[0]=-1.0+2.0*double(i0)/double(nplot-1);
1055  s[1]=-1.0+2.0*double(i1)/double(nplot-1);
1056  s[2]=-1.0+2.0*double(i2)/double(nplot-1);
1057  }
1058  else
1059  {
1060  s[0]=0.0;
1061  s[1]=0.0;
1062  s[2]=0.0;
1063  }
1064  }
1065 
1066 
1067  /// \short Return string for tecplot zone header (when plotting
1068  /// nplot points in each "coordinate direction)
1069  std::string tecplot_zone_string(const unsigned& nplot) const
1070  {
1071  std::ostringstream header;
1072  header << "ZONE I=" << nplot << ", J=" << nplot << ", K="
1073  << nplot << "\n";
1074  return header.str();
1075  }
1076 
1077  /// \short Return total number of plot points (when plotting
1078  /// nplot points in each "coordinate direction)
1079  unsigned nplot_points(const unsigned& nplot) const
1080  {
1081  unsigned DIM=3;
1082  unsigned np=1;
1083  for (unsigned i=0;i<DIM;i++) {np*=nplot;}
1084  return np;
1085  }
1086 
1087  /// \short Build the lower-dimensional FaceElement (an element of type
1088  /// QSpectralElement<2,NNODE_1D>). The face index takes six values
1089  /// corresponding to the six possible faces:
1090  /// -1 (Left) s[0] = -1.0
1091  /// +1 (Right) s[0] = 1.0
1092  /// -2 (Down) s[1] = -1.0
1093  /// +2 (Up) s[1] = 1.0
1094  /// -3 (Back) s[2] = -1.0
1095  /// +3 (Front) s[2] = 1.0
1096  void build_face_element(const int &face_index,
1097  FaceElement* face_element_pt);
1098 };
1099 
1100 
1101 //=======================================================================
1102 ///Shape function for specific QSpectralElement<3,..>
1103 //=======================================================================
1104 template <unsigned NNODE_1D>
1106  const
1107 {
1108  //Call the OneDimensional Shape functions
1109  OneDimensionalLegendreShape<NNODE_1D> psi1(s[0]), psi2(s[1]), psi3(s[2]);
1110  //OneDLegendreShapeParam psi1(NNODE_1D,s[0]), psi2(NNODE_1D,s[1]),
1111  // psi3(NNODE_1D,s[2]);
1112 
1113  //Now let's loop over the nodal points in the element
1114  //and copy the values back in
1115  for(unsigned i=0;i<NNODE_1D;i++)
1116  {
1117  for(unsigned j=0;j<NNODE_1D;j++)
1118  {
1119  for(unsigned k=0;k<NNODE_1D;k++)
1120  {
1121  psi(NNODE_1D*NNODE_1D*i + NNODE_1D*j + k) = psi3[i]*psi2[j]*psi1[k];
1122  }
1123  }
1124  }
1125 }
1126 
1127 //=======================================================================
1128 ///Derivatives of shape functions for specific QSpectralElement<3,..>
1129 //=======================================================================
1130 template <unsigned NNODE_1D>
1132  Shape &psi,
1133  DShape &dpsids) const
1134 {
1135  //Call the shape functions and derivatives
1136  //OneDLegendreShapeParam psi1(NNODE_1D,s[0]), psi2(NNODE_1D,s[1]),
1137  // psi3(NNODE_1D,s[2]);
1138  //OneDLegendreDShapeParam dpsi1ds(NNODE_1D,s[0]), dpsi2ds(NNODE_1D,s[1]),
1139  // dpsi3ds(NNODE_1D,s[2]);
1140  OneDimensionalLegendreShape<NNODE_1D> psi1(s[0]), psi2(s[1]), psi3(s[2]);
1141  OneDimensionalLegendreDShape<NNODE_1D> dpsi1ds(s[0]), dpsi2ds(s[1]),
1142  dpsi3ds(s[2]);
1143 
1144 
1145  //Index for the shape functions
1146  unsigned index=0;
1147 
1148  //Loop over shape functions in element
1149  for(unsigned i=0;i<NNODE_1D;i++)
1150  {
1151  for(unsigned j=0;j<NNODE_1D;j++)
1152  {
1153  for(unsigned k=0;k<NNODE_1D;k++)
1154  {
1155  //Assign the values
1156  dpsids(index,0) = psi3[i]*psi2[j]*dpsi1ds[k];
1157  dpsids(index,1) = psi3[i]*dpsi2ds[j]*psi1[k];
1158  dpsids(index,2) = dpsi3ds[i]*psi2[j]*psi1[k];
1159  psi[index] = psi3[i]*psi2[j]*psi1[k];
1160  //Increment the index
1161  ++index;
1162  }
1163  }
1164  }
1165 }
1166 
1167 
1168 //=======================================================================
1169 ///Second derivatives of shape functions for specific QSpectralElement<3,..>
1170 /// d2psids(i,0) = \f$ \partial ^2 \psi_j / \partial s_0^2 \f$
1171 /// d2psids(i,1) = \f$ \partial ^2 \psi_j / \partial s_1^2 \f$
1172 /// d2psids(i,2) = \f$ \partial ^2 \psi_j / \partial s_2^2 \f$
1173 /// d2psids(i,3) = \f$ \partial ^2 \psi_j / \partial s_0 \partial s_1 \f$
1174 /// d2psids(i,4) = \f$ \partial ^2 \psi_j / \partial s_0 \partial s_2 \f$
1175 /// d2psids(i,5) = \f$ \partial ^2 \psi_j / \partial s_1 \partial s_2 \f$
1176 //=======================================================================
1177 template <unsigned NNODE_1D>
1179  DShape &dpsids, DShape &d2psids) const
1180 {
1181  std::ostringstream error_message;
1182  error_message <<"\nd2shpe_local currently not implemented for this element\n";
1183  throw OomphLibError(error_message.str(),
1184  OOMPH_CURRENT_FUNCTION,
1185  OOMPH_EXCEPTION_LOCATION);
1186 
1187 /* //Call the shape functions and derivatives */
1188 /* OneDimensionalShape<NNODE_1D> psi1(s[0]); */
1189 /* OneDimensionalDShape<NNODE_1D> dpsi1ds(s[0]); */
1190 /* OneDimensionalD2Shape<NNODE_1D> d2psi1ds(s[0]); */
1191 
1192 /* //Loop over shape functions in element and assign to psi */
1193 /* for(unsigned l=0;l<NNODE_1D;l++) */
1194 /* { */
1195 /* psi[l] = psi1[l]; */
1196 /* dpsids[l][0] = dpsi1ds[l]; */
1197 /* d2psids[l][0] = d2psi1ds[l]; */
1198 /* } */
1199 }
1200 
1201 //==============================================================
1202 /// A class that is used to template the refineable Q spectral elements
1203 /// by dimension. It's really nothing more than a policy class
1204 //=============================================================
1205 template<unsigned DIM>
1207 {
1208  public:
1209 
1210  /// Empty constuctor
1212 };
1213 
1214 }
1215 
1216 #endif
double s_min() const
Min. value of local coordinate.
void local_fraction_of_node(const unsigned &n, Vector< double > &s_fraction)
Get the local fractino of node j in the element.
void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting.
double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
The local one-d fraction is the same.
static GaussLobattoLegendre< 2, NNODE_1D > integral
Default integration rule: Gaussian integration of same 'order' as the element.
Base class for all line elements.
Definition: Qelements.h:504
virtual void assign_all_generic_local_eqn_numbers(const bool &store_local_dof_pt)
Overloaded version of the calculation of the local equation numbers. If the boolean argument is true ...
Definition: elements.h:2087
double dlegendre(const unsigned &p, const double &x)
Calculates first derivative of Legendre polynomial of degree p at x using three term recursive formul...
Definition: orthpoly.h:118
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
unsigned nvertex_node() const
Number of vertex nodes in the element.
DenseMatrix< int > Spectral_local_eqn
Local equation numbers for the spectral degrees of freedom.
cstr elem_len * i
Definition: cfortran.h:607
unsigned nnode_1d() const
Number of nodes along each element edge.
Base class for all brick elements.
Definition: Qelements.h:1175
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
static GaussLobattoLegendre< 3, NNODE_1D > integral
Default integration rule: Gaussian integration of same 'order' as the element.
double s_max() const
Max. value of local coordinate.
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 3D version.
void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting.
A general Finite Element class.
Definition: elements.h:1271
static std::map< unsigned, Vector< double > > z
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
Definition: nodes.h:322
double legendre(const unsigned &p, const double &x)
Calculates Legendre polynomial of degree p at x using the three term recurrence relation ...
Definition: orthpoly.h:62
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
Data * spectral_data_pt(const unsigned &i) const
Return the i-th data object associated with the polynomials of order p. Note that i <= p...
static double nodal_position(const unsigned &n)
Definition: shape.h:768
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
OneDLegendreShapeParam(const unsigned &order, const double &s)
Constructor.
RefineableQSpectralElement()
Empty constuctor.
std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction) ...
unsigned nnode_1d() const
Number of nodes along each element edge.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
virtual void assign_all_generic_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers. If the boolean argument is true then store degrees of freedom at D...
double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
The local one-d fraction is the same.
void local_fraction_of_node(const unsigned &n, Vector< double > &s_fraction)
Get the local fractino of node j in the element.
std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction) ...
void local_coordinate_of_node(const unsigned &n, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size.
Base class for all quad elements.
Definition: Qelements.h:823
void gll_nodes(const unsigned &Nnode, Vector< double > &x)
Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1.
Definition: orthpoly.cc:38
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Class that returns the shape functions associated with legendre.
double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
The local one-d fraction is the same.
double ddlegendre(const unsigned &p, const double &x)
Calculates second derivative of Legendre polynomial of degree p at x using three term recursive formu...
Definition: orthpoly.h:137
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned nvertex_node() const
Number of vertex nodes in the element.
void output(FILE *file_pt, const unsigned &n_plot)
C_style output at n_plot points.
static long Is_pinned
Static "Magic number" used in place of the equation number to indicate that the value is pinned...
Definition: nodes.h:192
virtual void describe_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the dofs of the Node. The ostream specifies the output stream to which the descr...
Definition: nodes.cc:905
static char t char * s
Definition: cfortran.h:572
Vector< unsigned > Spectral_order
Vector that represents the spectral order in each dimension.
virtual void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
OneDLegendreDShapeParam(const unsigned &order, const double &s)
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
double s_max() const
Max. value of local coordinate.
Vector< unsigned > Nodal_spectral_order
Vector that represents the nodal spectral order.
static GaussLobattoLegendre< 1, NNODE_1D > integral
Default integration rule: Gaussian integration of same 'order' as the element.
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.
unsigned nspectral() const
void output(FILE *file_pt, const unsigned &n_plot)
C_style output at n_plot points.
std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction) ...
void output(FILE *file_pt)
C-style output.
double s_max() const
Max. value of local coordinate.
void output(FILE *file_pt)
C-style output.
void local_coordinate_of_node(const unsigned &n, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size.
static void calculate_nodal_positions(const unsigned &order)
Static function used to populate the stored positions.
void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting.
unsigned nnode_1d() const
Number of nodes along each element edge.
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.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
static long Is_unclassified
Static "Magic number" used in place of the equation number to denote a value that hasn't been class...
Definition: nodes.h:201
static void calculate_nodal_positions()
Static function used to populate the stored positions.
Definition: shape.h:759
static std::deque< double * > Dof_pt_deque
Static storage for deque used to add_global_equation_numbers when pointers to the dofs in each elemen...
Definition: elements.h:232
double s_min() const
Min. value of local coordinate.
unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction) ...
void local_fraction_of_node(const unsigned &n, Vector< double > &s_fraction)
Get the local fractino of node j in the element.
void output(FILE *file_pt)
C-style output.
void output(FILE *file_pt, const unsigned &n_plot)
C_style output at n_plot points.
Vector< Data * > * Spectral_data_pt
Additional Storage for shared spectral data.
virtual void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element[s]. The ostream specifies the output stream to whi...
Definition: elements.cc:1680
unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction) ...
unsigned nvertex_node() const
Number of vertex nodes in the element.
const double eps
Definition: orthpoly.h:57
void local_coordinate_of_node(const unsigned &n, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size.
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.
double s_min() const
Min. value of local coordinate.
Class that returns the shape functions associated with legendre.
Definition: shape.h:751
static double nodal_position(const unsigned &order, const unsigned &n)
unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction) ...