scalar_advection_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: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for scalar advection elements
31 
32 #ifndef OOMPH_SCALAR_ADVECTION_ELEMENTS_HEADER
33 #define OOMPH_SCALAR_ADVECTION_ELEMENTS_HEADER
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37 #include <oomph-lib-config.h>
38 #endif
39 
41 #include "../generic/Qelements.h"
42 #include "../generic/Qspectral_elements.h"
43 #include "../generic/dg_elements.h"
44 
45 namespace oomph
46 {
47 
48 //==============================================================
49 ///Base class for advection equations
50 //=============================================================
51 template<unsigned DIM>
53 {
54  /// \short Typedef for a wind function as a possible function of position
55  typedef void (*ScalarAdvectionWindFctPt)(const Vector<double>& x,
56  Vector<double>& wind);
57 
58  /// \short Function pointer to the wind function
60 
61 protected:
62 
63  /// \short A single flux is interpolated
64  inline unsigned nflux() const {return 1;}
65 
66  /// \short Return the flux as a function of the unknown
67  void flux(const Vector<double> &u, DenseMatrix<double> &f);
68 
69  /// \short Return the flux derivatives as a function of the unknowns
70  void dflux_du(const Vector<double> &u, RankThreeTensor<double> &df_du);
71 
72  ///\short Return the wind at a given position
73  inline virtual void get_wind_scalar_adv(const unsigned& ipt,
74  const Vector<double> &s,
75  const Vector<double>& x,
76  Vector<double>& wind) const
77  {
78  //If no wind function has been set, return zero
79  if(Wind_fct_pt==0)
80  {
81  for(unsigned i=0;i<DIM;i++) {wind[i]= 0.0;}
82  }
83  else
84  {
85  // Get wind
86  (*Wind_fct_pt)(x,wind);
87  }
88  }
89 
90 public:
91 
92  ///Constructor
94  Wind_fct_pt(0){ }
95 
96  /// Access function: Pointer to wind function
98 
99  /// Access function: Pointer to wind function. Const version
101 
102  ///\short The number of unknowns at each node is the number of values
103  unsigned required_nvalue(const unsigned &n) const {return 1;}
104 
105  ///Compute the error and norm of solution integrated over the element
106  ///Does not plot the error in the outfile
107  void compute_error(std::ostream &outfile,
109  initial_condition_pt, const double &t,
110  Vector<double> &error, Vector<double> &norm)
111  {
112  //Find the number of fluxes
113  const unsigned n_flux = this->nflux();
114  //Find the number of nodes
115  const unsigned n_node = this->nnode();
116  //Storage for the shape function and derivatives of shape function
117  Shape psi(n_node);
118  DShape dpsidx(n_node,DIM);
119 
120  //Find the number of integration points
121  unsigned n_intpt = this->integral_pt()->nweight();
122 
123  error.resize(n_flux); norm.resize(n_flux);
124  for(unsigned i=0;i<n_flux;i++) {error[i] = 0.0; norm[i] = 0.0;}
125 
126  Vector<double> s(DIM);
127 
128  //Loop over the integration points
129  for(unsigned ipt=0;ipt<n_intpt;ipt++)
130  {
131  //Get the shape functions at the knot
132  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
133 
134  //Get the integral weight
135  double W = this->integral_pt()->weight(ipt)*J;
136 
137  //Storage for the local functions
139  Vector<double> interpolated_u(n_flux,0.0);
140 
141  //Loop over the shape functions
142  for(unsigned l=0;l<n_node;l++)
143  {
144  //Locally cache the shape function
145  const double psi_ = psi(l);
146  for(unsigned i=0;i<DIM;i++)
147  {
148  interpolated_x[i] += this->nodal_position(l,i)*psi_;
149  }
150 
151  for(unsigned i=0;i<n_flux;i++)
152  {
153  //Calculate the velocity and tangent vector
154  interpolated_u[i] += this->nodal_value(l,i)*psi_;
155  }
156  }
157 
158  //Get the global wind
159  Vector<double> wind(DIM);
160  this->get_wind_scalar_adv(ipt,s,interpolated_x,wind);
161 
162  //Rescale the position
163  for(unsigned i=0;i<DIM;i++)
164  {
165  interpolated_x[i] -= wind[i]*t;
166  }
167 
168  //Now get the initial condition at this value of x
169  Vector<double> exact_u(n_flux,0.0);
170  (*initial_condition_pt)(0.0,interpolated_x,exact_u);
171 
172  //Loop over the unknowns
173  for(unsigned i=0;i<n_flux;i++)
174  {
175  error[i] +=
176  pow((interpolated_u[i] - exact_u[i]),2.0)*W;
177  norm[i] += interpolated_u[i]*interpolated_u[i]*W;
178  }
179  }
180  }
181 
182 
183 };
184 
185 
186 template <unsigned DIM, unsigned NNODE_1D>
188  public virtual QSpectralElement<DIM,NNODE_1D>,
189  public virtual ScalarAdvectionEquations<DIM>
190  {
191  public:
192 
193 
194  ///\short Constructor: Call constructors for QElement and
195  /// Advection Diffusion equations
198  { }
199 
200  /// Broken copy constructor
203  {
204  BrokenCopy::broken_copy("QSpectralScalarAdvectionElement");
205  }
206 
207  /// Broken assignment operator
208 //Commented out broken assignment operator because this can lead to a conflict warning
209 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
210 //realise that two separate implementations of the broken function are the same and so,
211 //quite rightly, it shouts.
212  /*void operator=(
213  const QSpectralScalarAdvectionElement<DIM,NNODE_1D>&)
214  {
215  BrokenCopy::broken_assign("QSpectralScalarAdvectionElement");
216  }*/
217 
218  /// \short Required # of `values' (pinned or dofs)
219  /// at node n
220  inline unsigned required_nvalue(const unsigned &n) const {return 1;}
221 
222  /// \short Output function:
223  /// x,y,u or x,y,z,u
224  void output(std::ostream &outfile)
226 
227  /// \short Output function:
228  /// x,y,u or x,y,z,u at n_plot^DIM plot points
229  void output(std::ostream &outfile, const unsigned &n_plot)
230  {ScalarAdvectionEquations<DIM>::output(outfile,n_plot);}
231 
232 
233  /*/// \short C-style output function:
234  /// x,y,u or x,y,z,u
235  void output(FILE* file_pt)
236  {
237  ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt);
238  }
239 
240  /// \short C-style output function:
241  /// x,y,u or x,y,z,u at n_plot^DIM plot points
242  void output(FILE* file_pt, const unsigned &n_plot)
243  {
244  ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt,n_plot);
245  }
246 
247  /// \short Output function for an exact solution:
248  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
249  void output_fct(std::ostream &outfile, const unsigned &n_plot,
250  FiniteElement::SteadyExactSolutionFctPt
251  exact_soln_pt)
252  {
253  ScalarAdvectionEquations<NFLUX,DIM>::
254  output_fct(outfile,n_plot,exact_soln_pt);}
255 
256 
257  /// \short Output function for a time-dependent exact solution.
258  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
259  /// (Calls the steady version)
260  void output_fct(std::ostream &outfile, const unsigned &n_plot,
261  const double& time,
262  FiniteElement::UnsteadyExactSolutionFctPt
263  exact_soln_pt)
264  {
265  ScalarAdvectionEquations<NFLUX,DIM>::
266  output_fct(outfile,n_plot,time,exact_soln_pt);
267  }*/
268 
269 
270 protected:
271 
272  /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
274  const Vector<double> &s,
275  Shape &psi,
276  DShape &dpsidx,
277  Shape &test,
278  DShape &dtestdx) const;
279 
280  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
281  /// integration point ipt. Return Jacobian.
283  const unsigned& ipt,
284  Shape &psi,
285  DShape &dpsidx,
286  Shape &test,
287  DShape &dtestdx)
288  const;
289 
290 };
291 
292 //Inline functions:
293 
294 
295 //======================================================================
296 /// \short Define the shape functions and test functions and derivatives
297 /// w.r.t. global coordinates and return Jacobian of mapping.
298 ///
299 /// Galerkin: Test functions = shape functions
300 //======================================================================
301 template<unsigned DIM, unsigned NNODE_1D>
304  Shape &psi,
305  DShape &dpsidx,
306  Shape &test,
307  DShape &dtestdx) const
308 {
309  //Call the geometrical shape functions and derivatives
310  double J = this->dshape_eulerian(s,psi,dpsidx);
311 
312  //Loop over the test functions and derivatives and set them equal to the
313  //shape functions
314  for(unsigned i=0;i<NNODE_1D;i++)
315  {
316  test[i] = psi[i];
317  for(unsigned j=0;j<DIM;j++)
318  {
319  dtestdx(i,j) = dpsidx(i,j);
320  }
321  }
322 
323  //Return the jacobian
324  return J;
325 }
326 
327 
328 
329 //======================================================================
330 /// Define the shape functions and test functions and derivatives
331 /// w.r.t. global coordinates and return Jacobian of mapping.
332 ///
333 /// Galerkin: Test functions = shape functions
334 //======================================================================
335 template<unsigned DIM, unsigned NNODE_1D>
338  const unsigned &ipt,
339  Shape &psi,
340  DShape &dpsidx,
341  Shape &test,
342  DShape &dtestdx) const
343 {
344  //Call the geometrical shape functions and derivatives
345  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
346 
347  //Set the test functions equal to the shape functions (pointer copy)
348  test = psi;
349  dtestdx = dpsidx;
350 
351  //Return the jacobian
352  return J;
353 }
354 
355 
356 ////////////////////////////////////////////////////////////////////////
357 ////////////////////////////////////////////////////////////////////////
358 ////////////////////////////////////////////////////////////////////////
359 
360 
361 
362 //=======================================================================
363 /// \short Face geometry for the QScalarAdvectionElement elements:
364 /// The spatial dimension of the face elements is one lower than that
365 /// of the bulk element but they have the same number of points along
366 /// their 1D edges.
367 //=======================================================================
368 template<unsigned DIM, unsigned NNODE_1D>
370  QSpectralScalarAdvectionElement<DIM,NNODE_1D> >:
371 public virtual QSpectralElement<DIM-1,NNODE_1D>
372 {
373 
374  public:
375 
376  /// \short Constructor: Call the constructor for the
377  /// appropriate lower-dimensional QElement
378  FaceGeometry() : QSpectralElement<DIM-1,NNODE_1D>() {}
379 
380 };
381 
382 
383 
384 
385  //======================================================================
386  /// FaceElement for Discontinuous Galerkin Problems
387  //======================================================================
388  template<class ELEMENT>
389  class DGScalarAdvectionFaceElement : public virtual FaceGeometry<ELEMENT>,
390  public virtual DGFaceElement
391  {
392  public:
393 
394  /// Constructor
396  const int &face_index) :
397  FaceGeometry<ELEMENT>(), DGFaceElement()
398  {
399  //Attach geometric information to the element
400  //N.B. This function also assigns nbulk_value from required_nvalue
401  //of the bulk element.
402  element_pt->build_face_element(face_index,this);
403  }
404 
405 
406  //There is a single required n_flux
407  unsigned required_nflux() {return 1;}
408 
409  ///Calculate the normal flux, so this is the dot product of the
410  ///numerical flux with n_in
411  void numerical_flux(const Vector<double> &n_out,
412  const Vector<double> &u_int,
413  const Vector<double> &u_ext,
414  Vector<double> &flux)
415  {
416  const unsigned dim = this->nodal_dimension();
417  Vector<double> Wind(dim);
418  Vector<double> s, x;
419  // Dummy integration point
420  unsigned ipt=0;
421  dynamic_cast<ELEMENT*>(this->bulk_element_pt())->
422  get_wind_scalar_adv(ipt,s,x,Wind);
423 
424  //Now we can work this out for standard upwind advection
425  double dot=0.0;
426  for(unsigned i=0;i<dim;i++) {dot += Wind[i]*n_out[i];}
427 
428  const unsigned n_value = this->required_nflux();
429  if(dot >= 0.0)
430  {
431  for(unsigned n=0;n<n_value;n++) {flux[n] = dot*u_int[n];}
432  }
433  else
434  {
435  for(unsigned n=0;n<n_value;n++) {flux[n] = dot*u_ext[n];}
436  }
437 
438  //flux[0] = 0.5*(u_int[0]+u_ext[0])*n_out[0];
439  }
440 
441 
442  };
443 
444 //=================================================================
445 ///General DGScalarAdvectionClass. Establish the template parameters
446 //===================================================================
447 template<unsigned DIM, unsigned NNODE_1D>
449 {
450 
451 };
452 
453 
454 //==================================================================
455 //Specialization for 1D DG Advection element
456 //==================================================================
457 template<unsigned NNODE_1D>
459  public QSpectralScalarAdvectionElement<1,NNODE_1D>,
460  public DGElement
461 {
462  friend class
464 
465  Vector<double> Inverse_mass_diagonal;
466 
467 public:
468 
469  //There is a single required n_flux
470  unsigned required_nflux() {return 1;}
471 
472  //Calculate averages
473  void calculate_element_averages(double* &average_value)
475 
476  //Constructor
478  QSpectralScalarAdvectionElement<1,NNODE_1D>(),
479  DGElement() {}
480 
482 
483  void build_all_faces()
484  {
485  //Make the two faces
486  Face_element_pt.resize(2);
487  //Make the face on the left
488  Face_element_pt[0] =
491  (this,-1);
492  //Make the face on the right
493  Face_element_pt[1] =
495  DGSpectralScalarAdvectionElement<1,NNODE_1D> >
496  (this,+1);
497  }
498 
499 
500  ///\short Compute the residuals for the Navier--Stokes equations;
501  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
502  void fill_in_generic_residual_contribution_flux_transport(
503  Vector<double> &residuals, DenseMatrix<double> &jacobian,
504  DenseMatrix<double> &mass_matrix, unsigned flag)
505  {
508  jacobian,
509  mass_matrix,
510  flag);
511 
512  this->add_flux_contributions_to_residuals(residuals,jacobian,flag);
513  }
514 
515 
516 //============================================================================
517 ///Function that returns the current value of the residuals
518 ///multiplied by the inverse mass matrix (virtual so that it can be overloaded
519 ///specific elements in which time and memory saving tricks can be applied)
520 //============================================================================
521 void get_inverse_mass_matrix_times_residuals(Vector<double> &minv_res)
522 {
523  //If there are external data this is not going to work
524  if(nexternal_data() > 0)
525  {
526  std::ostringstream error_stream;
527  error_stream <<
528  "Cannot use a discontinuous formulation for the mass matrix when\n"
529  <<
530  "there are external data.\n " <<
531  "Do not call Problem::enable_discontinuous_formulation()\n";
532 
533  throw OomphLibError(error_stream.str(),
534  OOMPH_CURRENT_FUNCTION,
535  OOMPH_EXCEPTION_LOCATION);
536  }
537 
538 
539  //Now let's assemble stuff
540  const unsigned n_dof = this->ndof();
541 
542  //Resize and initialise the vector that will holds the residuals
543  minv_res.resize(n_dof);
544  for(unsigned n=0;n<n_dof;n++) {minv_res[n] = 0.0;}
545 
546  //If we are recycling the mass matrix
547  if(Mass_matrix_reuse_is_enabled && Mass_matrix_has_been_computed)
548  {
549  //Get the residuals
550  this->fill_in_contribution_to_residuals(minv_res);
551  }
552  //Otherwise
553  else
554  {
555  //Temporary mass matrix
556  DenseDoubleMatrix M(n_dof,n_dof,0.0);
557 
558  //Get the local mass matrix and residuals
559  this->fill_in_contribution_to_mass_matrix(minv_res,M);
560 
561  //Store the diagonal entries
562  Inverse_mass_diagonal.clear();
563  for(unsigned n=0;n<n_dof;n++) {Inverse_mass_diagonal.push_back(1.0/M(n,n));}
564 
565  //The mass matrix has been computed
566  Mass_matrix_has_been_computed=true;
567  }
568 
569  for(unsigned n=0;n<n_dof;n++) {minv_res[n] *= Inverse_mass_diagonal[n];}
570 
571 }
572 
573 
574 
575 };
576 
577 
578 //=======================================================================
579 /// Face geometry of the 1D DG elements
580 //=======================================================================
581 template<unsigned NNODE_1D>
583 public virtual PointElement
584 {
585  public:
587 };
588 
589 
590 
591 //==================================================================
592 ///Specialisation for 2D DG Elements
593 //==================================================================
594 template<unsigned NNODE_1D>
596  public QSpectralScalarAdvectionElement<2,NNODE_1D>,
597  public DGElement
598 {
599  friend class
601 
602 public:
603 
604  //Calculate averages
605  void calculate_element_averages(double* &average_value)
607 
608  //There is a single required n_flux
609  unsigned required_nflux() {return 1;}
610 
611  //Constructor
613  QSpectralScalarAdvectionElement<2,NNODE_1D>(),
614  DGElement() {}
615 
617 
618  void build_all_faces()
619  {
620  Face_element_pt.resize(4);
621  Face_element_pt[0] =
624  (this,2);
625  Face_element_pt[1] =
627  DGSpectralScalarAdvectionElement<2,NNODE_1D> >
628  (this,1);
629  Face_element_pt[2] =
631  DGSpectralScalarAdvectionElement<2,NNODE_1D> >
632  (this,-2);
633  Face_element_pt[3] =
635  DGSpectralScalarAdvectionElement<2,NNODE_1D> >
636  (this,-1);
637  }
638 
639 
640  ///\short Compute the residuals for the Navier--Stokes equations;
641  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
642  void fill_in_generic_residual_contribution_flux_transport(
643  Vector<double> &residuals, DenseMatrix<double> &jacobian,
644  DenseMatrix<double> &mass_matrix, unsigned flag)
645  {
648  jacobian,
649  mass_matrix,
650  flag);
651 
652  this->add_flux_contributions_to_residuals(residuals,jacobian,flag);
653  }
654 
655 
656 };
657 
658 
659 //=======================================================================
660 /// Face geometry of the DG elements
661 //=======================================================================
662 template<unsigned NNODE_1D>
664  public virtual QSpectralElement<1,NNODE_1D>
665 {
666  public:
667  FaceGeometry() : QSpectralElement<1,NNODE_1D>() {}
668 };
669 
670 
671 
672 //=============================================================
673 ///Non-spectral version of the classes
674 //============================================================
675 template <unsigned DIM, unsigned NNODE_1D>
677  public virtual QElement<DIM,NNODE_1D>,
678  public virtual ScalarAdvectionEquations<DIM>
679  {
680  public:
681 
682 
683  ///\short Constructor: Call constructors for QElement and
684  /// Advection Diffusion equations
685  QScalarAdvectionElement() : QElement<DIM,NNODE_1D>(),
687  { }
688 
689  /// Broken copy constructor
692  {
693  BrokenCopy::broken_copy("QScalarAdvectionElement");
694  }
695 
696  /// Broken assignment operator
697  /*void operator=(
698  const QScalarAdvectionElement<DIM,NNODE_1D>&)
699  {
700  BrokenCopy::broken_assign("QScalarAdvectionElement");
701  }*/
702 
703  /// \short Required # of `values' (pinned or dofs)
704  /// at node n
705  inline unsigned required_nvalue(const unsigned &n) const {return 1;}
706 
707  /// \short Output function:
708  /// x,y,u or x,y,z,u
709  void output(std::ostream &outfile)
711 
712  /// \short Output function:
713  /// x,y,u or x,y,z,u at n_plot^DIM plot points
714  void output(std::ostream &outfile, const unsigned &n_plot)
715  {ScalarAdvectionEquations<DIM>::output(outfile,n_plot);}
716 
717 
718  /*/// \short C-style output function:
719  /// x,y,u or x,y,z,u
720  void output(FILE* file_pt)
721  {
722  ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt);
723  }
724 
725  /// \short C-style output function:
726  /// x,y,u or x,y,z,u at n_plot^DIM plot points
727  void output(FILE* file_pt, const unsigned &n_plot)
728  {
729  ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt,n_plot);
730  }
731 
732  /// \short Output function for an exact solution:
733  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
734  void output_fct(std::ostream &outfile, const unsigned &n_plot,
735  FiniteElement::SteadyExactSolutionFctPt
736  exact_soln_pt)
737  {
738  ScalarAdvectionEquations<NFLUX,DIM>::
739  output_fct(outfile,n_plot,exact_soln_pt);}
740 
741 
742  /// \short Output function for a time-dependent exact solution.
743  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
744  /// (Calls the steady version)
745  void output_fct(std::ostream &outfile, const unsigned &n_plot,
746  const double& time,
747  FiniteElement::UnsteadyExactSolutionFctPt
748  exact_soln_pt)
749  {
750  ScalarAdvectionEquations<NFLUX,DIM>::
751  output_fct(outfile,n_plot,time,exact_soln_pt);
752  }*/
753 
754 
755 protected:
756 
757  /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
759  const Vector<double> &s,
760  Shape &psi,
761  DShape &dpsidx,
762  Shape &test,
763  DShape &dtestdx) const;
764 
765  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
766  /// integration point ipt. Return Jacobian.
768  const unsigned& ipt,
769  Shape &psi,
770  DShape &dpsidx,
771  Shape &test,
772  DShape &dtestdx)
773  const;
774 
775 };
776 
777 //Inline functions:
778 
779 
780 //======================================================================
781 /// \short Define the shape functions and test functions and derivatives
782 /// w.r.t. global coordinates and return Jacobian of mapping.
783 ///
784 /// Galerkin: Test functions = shape functions
785 //======================================================================
786 template<unsigned DIM, unsigned NNODE_1D>
789  Shape &psi,
790  DShape &dpsidx,
791  Shape &test,
792  DShape &dtestdx) const
793 {
794  //Call the geometrical shape functions and derivatives
795  double J = this->dshape_eulerian(s,psi,dpsidx);
796 
797  //Loop over the test functions and derivatives and set them equal to the
798  //shape functions
799  for(unsigned i=0;i<NNODE_1D;i++)
800  {
801  test[i] = psi[i];
802  for(unsigned j=0;j<DIM;j++)
803  {
804  dtestdx(i,j) = dpsidx(i,j);
805  }
806  }
807 
808  //Return the jacobian
809  return J;
810 }
811 
812 
813 
814 //======================================================================
815 /// Define the shape functions and test functions and derivatives
816 /// w.r.t. global coordinates and return Jacobian of mapping.
817 ///
818 /// Galerkin: Test functions = shape functions
819 //======================================================================
820 template<unsigned DIM, unsigned NNODE_1D>
823  const unsigned &ipt,
824  Shape &psi,
825  DShape &dpsidx,
826  Shape &test,
827  DShape &dtestdx) const
828 {
829  //Call the geometrical shape functions and derivatives
830  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
831 
832  //Set the test functions equal to the shape functions (pointer copy)
833  test = psi;
834  dtestdx = dpsidx;
835 
836  //Return the jacobian
837  return J;
838 }
839 
840 
841 ////////////////////////////////////////////////////////////////////////
842 ////////////////////////////////////////////////////////////////////////
843 ////////////////////////////////////////////////////////////////////////
844 
845 
846 
847 //=======================================================================
848 /// \short Face geometry for the QScalarAdvectionElement elements:
849 /// The spatial dimension of the face elements is one lower than that
850 /// of the bulk element but they have the same number of points along
851 /// their 1D edges.
852 //=======================================================================
853 template<unsigned DIM, unsigned NNODE_1D>
855  QScalarAdvectionElement<DIM,NNODE_1D> >:
856 public virtual QElement<DIM-1,NNODE_1D>
857 {
858 
859  public:
860 
861  /// \short Constructor: Call the constructor for the
862  /// appropriate lower-dimensional QElement
863  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
864 
865 };
866 
867 
868 
869 //=================================================================
870 ///General DGScalarAdvectionClass. Establish the template parameters
871 //===================================================================
872 template<unsigned DIM, unsigned NNODE_1D>
874 {
875 
876 };
877 
878 
879 //==================================================================
880 //Specialization for 1D DG Advection element
881 //==================================================================
882 template<unsigned NNODE_1D>
883 class DGScalarAdvectionElement<1,NNODE_1D> :
884  public QScalarAdvectionElement<1,NNODE_1D>,
885  public DGElement
886 {
887  friend class
889 
890 public:
891 
892  //There is a single required n_flux
893  unsigned required_nflux() {return 1;}
894 
895  //Calculate averages
896  void calculate_element_averages(double* &average_value)
898 
899  //Constructor
901  DGElement() {}
902 
904 
905  void build_all_faces()
906  {
907  //Make the two faces
908  Face_element_pt.resize(2);
909  //Make the face on the left
910  Face_element_pt[0] =
913  (this,-1);
914  //Make the face on the right
915  Face_element_pt[1] =
917  DGScalarAdvectionElement<1,NNODE_1D> >
918  (this,+1);
919  }
920 
921 
922  ///\short Compute the residuals for the Navier--Stokes equations;
923  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
924  void fill_in_generic_residual_contribution_flux_transport(
925  Vector<double> &residuals, DenseMatrix<double> &jacobian,
926  DenseMatrix<double> &mass_matrix, unsigned flag)
927  {
930  jacobian,
931  mass_matrix,
932  flag);
933 
934  this->add_flux_contributions_to_residuals(residuals,jacobian,flag);
935  }
936 
937 
938 };
939 
940 
941 //=======================================================================
942 /// Face geometry of the 1D DG elements
943 //=======================================================================
944 template<unsigned NNODE_1D>
946 public virtual PointElement
947 {
948  public:
950 };
951 
952 
953 
954 //==================================================================
955 ///Specialisation for 2D DG Elements
956 //==================================================================
957 template<unsigned NNODE_1D>
958 class DGScalarAdvectionElement<2,NNODE_1D> :
959  public QScalarAdvectionElement<2,NNODE_1D>,
960  public DGElement
961 {
962  friend class
964 
965 public:
966 
967  //There is a single required n_flux
968  unsigned required_nflux() {return 1;}
969 
970  //Calculate averages
971  void calculate_element_averages(double* &average_value)
973 
974  //Constructor
976  QScalarAdvectionElement<2,NNODE_1D>(),
977  DGElement() {}
978 
980 
981  void build_all_faces()
982  {
983  Face_element_pt.resize(4);
984  Face_element_pt[0] =
987  (this,2);
988  Face_element_pt[1] =
990  DGScalarAdvectionElement<2,NNODE_1D> >
991  (this,1);
992  Face_element_pt[2] =
994  DGScalarAdvectionElement<2,NNODE_1D> >
995  (this,-2);
996  Face_element_pt[3] =
998  DGScalarAdvectionElement<2,NNODE_1D> >
999  (this,-1);
1000  }
1001 
1002 
1003  ///\short Compute the residuals for the Navier--Stokes equations;
1004  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
1005  void fill_in_generic_residual_contribution_flux_transport(
1006  Vector<double> &residuals, DenseMatrix<double> &jacobian,
1007  DenseMatrix<double> &mass_matrix, unsigned flag)
1008  {
1011  jacobian,
1012  mass_matrix,
1013  flag);
1014 
1015  this->add_flux_contributions_to_residuals(residuals,jacobian,flag);
1016  }
1017 
1018 
1019 };
1020 
1021 
1022 //=======================================================================
1023 /// Face geometry of the DG elements
1024 //=======================================================================
1025 template<unsigned NNODE_1D>
1027  public virtual QElement<1,NNODE_1D>
1028 {
1029  public:
1030  FaceGeometry() : QElement<1,NNODE_1D>() {}
1031 };
1032 
1033 
1034 
1035 }
1036 
1037 #endif
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1234
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
QSpectralScalarAdvectionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
void flux(const Vector< double > &u, DenseMatrix< double > &f)
Return the flux as a function of the unknown.
ScalarAdvectionWindFctPt wind_fct_pt() const
Access function: Pointer to wind function. Const version.
QScalarAdvectionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
DGScalarAdvectionFaceElement(FiniteElement *const &element_pt, const int &face_index)
Constructor.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
ScalarAdvectionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
cstr elem_len * i
Definition: cfortran.h:607
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
double dshape_and_dtest_eulerian_flux_transport(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
A general Finite Element class.
Definition: elements.h:1271
ScalarAdvectionWindFctPt Wind_fct_pt
Function pointer to the wind function.
char t
Definition: cfortran.h:572
A Base class for DGElements.
Definition: dg_elements.h:160
General DGScalarAdvectionClass. Establish the template parameters.
void(* ScalarAdvectionWindFctPt)(const Vector< double > &x, Vector< double > &wind)
Typedef for a wind function as a possible function of position.
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2358
unsigned nflux() const
A single flux is interpolated.
FaceElement for Discontinuous Galerkin Problems.
double dshape_and_dtest_eulerian_at_knot_flux_transport(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt. Return Jacobian.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
A Rank 3 Tensor class.
Definition: matrices.h:1337
unsigned required_nflux()
Set the number of flux components.
QSpectralScalarAdvectionElement(const QSpectralScalarAdvectionElement< DIM, NNODE_1D > &dummy)
Broken copy constructor.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
static char t char * s
Definition: cfortran.h:572
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
Base class for advection equations.
QScalarAdvectionElement(const QScalarAdvectionElement< DIM, NNODE_1D > &dummy)
Broken copy constructor.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
unsigned required_nvalue(const unsigned &n) const
The number of unknowns at each node is the number of values.
double dshape_and_dtest_eulerian_at_knot_flux_transport(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt. Return Jacobian.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2215
void dflux_du(const Vector< double > &u, RankThreeTensor< double > &df_du)
Return the flux derivatives as a function of the unknowns.
Non-spectral version of the classes.
General DGScalarAdvectionClass. Establish the template parameters.
double dot(const Vector< double > &a, const Vector< double > &b)
Probably not always best/fastest because not optimised for dimension but useful...
Definition: Vector.h:289
virtual void fill_in_generic_residual_contribution_flux_transport(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
double dshape_and_dtest_eulerian_flux_transport(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
Base class for Discontinuous Galerkin Faces. These are responsible for calculating the normal fluxes ...
Definition: dg_elements.h:53
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4470
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
void calculate_element_averages(double *&average_values)
Compute the average values of the fluxes.
virtual void get_wind_scalar_adv(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Return the wind at a given position.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
void numerical_flux(const Vector< double > &n_out, const Vector< double > &u_int, const Vector< double > &u_ext, Vector< double > &flux)
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt initial_condition_pt, const double &t, Vector< double > &error, Vector< double > &norm)
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3252