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