advection_diffusion_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 file for Advection Diffusion elements
31 #ifndef OOMPH_ADV_DIFF_ELEMENTS_HEADER
32 #define OOMPH_ADV_DIFF_ELEMENTS_HEADER
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37  #include <oomph-lib-config.h>
38 #endif
39 
40 //OOMPH-LIB headers
41 #include "../generic/nodes.h"
42 #include "../generic/Qelements.h"
43 #include "../generic/oomph_utilities.h"
44 
45 namespace oomph
46 {
47 
48 //=============================================================
49 /// \short A class for all elements that solve the Advection
50 /// Diffusion equations using isoparametric elements.
51 /// \f[
52 /// \frac{\partial^2 u}{\partial x_i^2} =
53 /// Pe w_i(x_k) \frac{\partial u}{\partial x_i} + f(x_j)
54 /// \f]
55 /// This contains the generic maths. Shape functions, geometric
56 /// mapping etc. must get implemented in derived class.
57 //=============================================================
58 template <unsigned DIM>
60 {
61 
62 public:
63 
64  /// \short Function pointer to source function fct(x,f(x)) --
65  /// x is a Vector!
67  double& f);
68 
69 
70  /// \short Function pointer to wind function fct(x,w(x)) --
71  /// x is a Vector!
72  typedef void (*AdvectionDiffusionWindFctPt)(const Vector<double>& x,
73  Vector<double>& wind);
74 
75 
76  /// \short Constructor: Initialise the Source_fct_pt and Wind_fct_pt
77  /// to null and set (pointer to) Peclet number to default
79  ALE_is_disabled(false)
80  {
81  //Set Peclet number to default
83  //Set Peclet Strouhal number to default
85  }
86 
87  /// Broken copy constructor
89  {
90  BrokenCopy::broken_copy("AdvectionDiffusionEquations");
91  }
92 
93  /// Broken assignment operator
95  {
96  BrokenCopy::broken_assign("AdvectionDiffusionEquations");
97  }
98 
99  /// \short Return the index at which the unknown value
100  /// is stored. The default value, 0, is appropriate for single-physics
101  /// problems, when there is only one variable, the value that satisfies
102  /// the advection-diffusion equation.
103  /// In derived multi-physics elements, this function should be overloaded
104  /// to reflect the chosen storage scheme. Note that these equations require
105  /// that the unknown is always stored at the same index at each node.
106  virtual inline unsigned u_index_adv_diff() const {return 0;}
107 
108 /// \short du/dt at local node n.
109  /// Uses suitably interpolated value for hanging nodes.
110  double du_dt_adv_diff(const unsigned &n) const
111  {
112  // Get the data's timestepper
114 
115  //Initialise dudt
116  double dudt=0.0;
117  //Loop over the timesteps, if there is a non Steady timestepper
118  if (!time_stepper_pt->is_steady())
119  {
120  //Find the index at which the variable is stored
121  const unsigned u_nodal_index = u_index_adv_diff();
122 
123  // Number of timsteps (past & present)
124  const unsigned n_time = time_stepper_pt->ntstorage();
125 
126  for(unsigned t=0;t<n_time;t++)
127  {
128  dudt += time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
129  }
130  }
131  return dudt;
132  }
133 
134  /// \short Disable ALE, i.e. assert the mesh is not moving -- you do this
135  /// at your own risk!
136  void disable_ALE()
137  {
138  ALE_is_disabled=true;
139  }
140 
141 
142  /// \short (Re-)enable ALE, i.e. take possible mesh motion into account
143  /// when evaluating the time-derivative. Note: By default, ALE is
144  /// enabled, at the expense of possibly creating unnecessary work
145  /// in problems where the mesh is, in fact, stationary.
146  void enable_ALE()
147  {
148  ALE_is_disabled=false;
149  }
150 
151 /// \short Number of scalars/fields output by this element. Reimplements
152  /// broken virtual function in base class.
153  unsigned nscalar_paraview() const
154  {
155  return DIM+1;
156  }
157 
158  /// \short Write values of the i-th scalar field at the plot points. Needs
159  /// to be implemented for each new specific element type.
160  void scalar_value_paraview(std::ofstream& file_out,
161  const unsigned& i,
162  const unsigned& nplot) const
163  {
164  //Vector of local coordinates
165  Vector<double> s(DIM);
166 
167  // Loop over plot points
168  unsigned num_plot_points=nplot_points_paraview(nplot);
169  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
170  {
171 
172  // Get local coordinates of plot point
173  get_s_plot(iplot,nplot,s);
174 
175  // Get Eulerian coordinate of plot point
176  Vector<double> x(DIM);
177  interpolated_x(s,x);
178 
179  if(i<DIM)
180  {
181  // Get the wind
182  Vector<double> wind(DIM);
183 
184  // Dummy ipt argument needed... ?
185  unsigned ipt=0;
186  get_wind_adv_diff(ipt,s,x,wind);
187 
188  file_out << wind[i] << std::endl;
189  }
190  // Advection Diffusion
191  else if(i==DIM)
192  {
193  file_out << interpolated_u_adv_diff(s) << std::endl;
194  }
195  // Never get here
196  else
197  {
198  std::stringstream error_stream;
199  error_stream
200  << "Advection Diffusion Elements only store " << DIM+1 << " fields "
201  << std::endl;
202  throw OomphLibError(
203  error_stream.str(),
204  OOMPH_CURRENT_FUNCTION,
205  OOMPH_EXCEPTION_LOCATION);
206  }
207  }
208  }
209 
210  /// \short Name of the i-th scalar field. Default implementation
211  /// returns V1 for the first one, V2 for the second etc. Can (should!) be
212  /// overloaded with more meaningful names in specific elements.
213  std::string scalar_name_paraview(const unsigned& i) const
214  {
215  // Winds
216  if(i<DIM)
217  {
218  return "Wind "+StringConversion::to_string(i);
219  }
220  // Advection Diffusion field
221  else if(i==DIM)
222  {
223  return "Advection Diffusion";
224  }
225  // Never get here
226  else
227  {
228  std::stringstream error_stream;
229  error_stream
230  << "Advection Diffusion Elements only store " << DIM+1 << " fields"
231  << std::endl;
232  throw OomphLibError(
233  error_stream.str(),
234  OOMPH_CURRENT_FUNCTION,
235  OOMPH_EXCEPTION_LOCATION);
236  // Dummy return
237  return " ";
238  }
239  }
240 
241  /// Output with default number of plot points
242  void output(std::ostream &outfile)
243  {
244  unsigned nplot=5;
245  output(outfile,nplot);
246  }
247 
248  /// \short Output FE representation of soln: x,y,u or x,y,z,u at
249  /// nplot^DIM plot points
250  void output(std::ostream &outfile, const unsigned &nplot);
251 
252 
253  /// C_style output with default number of plot points
254  void output(FILE* file_pt)
255  {
256  unsigned n_plot=5;
257  output(file_pt,n_plot);
258  }
259 
260  /// \short C-style output FE representation of soln: x,y,u or x,y,z,u at
261  /// n_plot^DIM plot points
262  void output(FILE* file_pt, const unsigned &n_plot);
263 
264 
265  /// Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points
266  void output_fct(std::ostream &outfile, const unsigned &nplot,
268  exact_soln_pt);
269 
270  /// \short Output exact soln: x,y,u_exact or x,y,z,u_exact at
271  /// nplot^DIM plot points (dummy time-dependent version to
272  /// keep intel compiler happy)
273  virtual void output_fct(std::ostream &outfile, const unsigned &nplot,
274  const double& time,
276  {
277  throw OomphLibError(
278  "There is no time-dependent output_fct() for Advection Diffusion elements",
279  OOMPH_CURRENT_FUNCTION,
280  OOMPH_EXCEPTION_LOCATION);
281  }
282 
283 
284  /// Get error against and norm of exact solution
285  void compute_error(std::ostream &outfile,
287  exact_soln_pt, double& error, double& norm);
288 
289 
290  /// Dummy, time dependent error checker
291  void compute_error(std::ostream &outfile,
293  exact_soln_pt,
294  const double& time, double& error, double& norm)
295  {
296  throw OomphLibError(
297  "No time-dependent compute_error() for Advection Diffusion elements",
298  OOMPH_CURRENT_FUNCTION,
299  OOMPH_EXCEPTION_LOCATION);
300  }
301 
302 
303  /// Access function: Pointer to source function
305 
306 
307  /// Access function: Pointer to source function. Const version
309 
310 
311  /// Access function: Pointer to wind function
313 
314 
315  /// Access function: Pointer to wind function. Const version
317 
318  /// Peclet number
319  const double &pe() const {return *Pe_pt;}
320 
321  /// Pointer to Peclet number
322  double* &pe_pt() {return Pe_pt;}
323 
324  /// Peclet number multiplied by Strouhal number
325  const double &pe_st() const {return *PeSt_pt;}
326 
327  /// Pointer to Peclet number multipled by Strouha number
328  double* &pe_st_pt() {return PeSt_pt;}
329 
330  /// \short Get source term at (Eulerian) position x. This function is
331  /// virtual to allow overloading in multi-physics problems where
332  /// the strength of the source function might be determined by
333  /// another system of equations
334  inline virtual void get_source_adv_diff(const unsigned& ipt,
335  const Vector<double>& x,
336  double& source) const
337  {
338  //If no source function has been set, return zero
339  if(Source_fct_pt==0) {source = 0.0;}
340  else
341  {
342  // Get source strength
343  (*Source_fct_pt)(x,source);
344  }
345  }
346 
347  /// \short Get wind at (Eulerian) position x and/or local coordinate s.
348  /// This function is
349  /// virtual to allow overloading in multi-physics problems where
350  /// the wind function might be determined by
351  /// another system of equations
352  inline virtual void get_wind_adv_diff(const unsigned& ipt,
353  const Vector<double> &s,
354  const Vector<double>& x,
355  Vector<double>& wind) const
356  {
357  //If no wind function has been set, return zero
358  if(Wind_fct_pt==0)
359  {
360  for(unsigned i=0;i<DIM;i++) {wind[i]= 0.0;}
361  }
362  else
363  {
364  // Get wind
365  (*Wind_fct_pt)(x,wind);
366  }
367  }
368 
369  /// Get flux: \f$\mbox{flux}[i] = \mbox{d}u / \mbox{d}x_i \f$
370  void get_flux(const Vector<double>& s, Vector<double>& flux) const
371  {
372  //Find out how many nodes there are in the element
373  unsigned n_node = nnode();
374 
375  //Get the nodal index at which the unknown is stored
376  unsigned u_nodal_index = u_index_adv_diff();
377 
378  //Set up memory for the shape and test functions
379  Shape psi(n_node);
380  DShape dpsidx(n_node,DIM);
381 
382  //Call the derivatives of the shape and test functions
383  dshape_eulerian(s,psi,dpsidx);
384 
385  //Initialise to zero
386  for(unsigned j=0;j<DIM;j++) {flux[j] = 0.0;}
387 
388  // Loop over nodes
389  for(unsigned l=0;l<n_node;l++)
390  {
391  //Loop over derivative directions
392  for(unsigned j=0;j<DIM;j++)
393  {
394  flux[j] += nodal_value(l,u_nodal_index)*dpsidx(l,j);
395  }
396  }
397  }
398 
399 
400  /// Add the element's contribution to its residual vector (wrapper)
402  {
403  //Call the generic residuals function with flag set to 0 and using
404  //a dummy matrix
408  }
409 
410 
411  /// \short Add the element's contribution to its residual vector and
412  /// the element Jacobian matrix (wrapper)
414  DenseMatrix<double> &jacobian)
415  {
416  //Call the generic routine with the flag set to 1
418  residuals,jacobian,GeneralisedElement::Dummy_matrix,1);
419  }
420 
421 
422  /// Add the element's contribution to its residuals vector,
423  /// jacobian matrix and mass matrix
425  Vector<double> &residuals, DenseMatrix<double> &jacobian,
426  DenseMatrix<double> &mass_matrix)
427  {
428  //Call the generic routine with the flag set to 2
430  jacobian,mass_matrix,2);
431  }
432 
433 
434  /// Return FE representation of function value u(s) at local coordinate s
435  inline double interpolated_u_adv_diff(const Vector<double> &s) const
436  {
437  //Find number of nodes
438  unsigned n_node = nnode();
439 
440  //Get the nodal index at which the unknown is stored
441  unsigned u_nodal_index = u_index_adv_diff();
442 
443  //Local shape function
444  Shape psi(n_node);
445 
446  //Find values of shape function
447  shape(s,psi);
448 
449  //Initialise value of u
450  double interpolated_u = 0.0;
451 
452  //Loop over the local nodes and sum
453  for(unsigned l=0;l<n_node;l++)
454  {
455  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
456  }
457 
458  return(interpolated_u);
459  }
460 
461 
462  ///\short Return derivative of u at point s with respect to all data
463  ///that can affect its value.
464  ///In addition, return the global equation numbers corresponding to the
465  ///data. This is virtual so that it can be overloaded in the
466  ///refineable version
468  const Vector<double> &s, Vector<double> &du_ddata,
469  Vector<unsigned> &global_eqn_number)
470  {
471  //Find number of nodes
472  const unsigned n_node = nnode();
473 
474  //Get the nodal index at which the unknown is stored
475  const unsigned u_nodal_index = u_index_adv_diff();
476 
477  //Local shape function
478  Shape psi(n_node);
479 
480  //Find values of shape function
481  shape(s,psi);
482 
483  //Find the number of dofs associated with interpolated u
484  unsigned n_u_dof=0;
485  for(unsigned l=0;l<n_node;l++)
486  {
487  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
488  //If it's positive add to the count
489  if (global_eqn >=0) {++n_u_dof;}
490  }
491 
492  //Now resize the storage schemes
493  du_ddata.resize(n_u_dof,0.0);
494  global_eqn_number.resize(n_u_dof,0);
495 
496  //Loop over the nodes again and set the derivatives
497  unsigned count=0;
498  for(unsigned l=0;l<n_node;l++)
499  {
500  //Get the global equation number
501  int global_eqn=this->node_pt(l)->eqn_number(u_nodal_index);
502  //If it's positive
503  if (global_eqn >= 0)
504  {
505  //Set the global equation number
506  global_eqn_number[count] = global_eqn;
507  //Set the derivative with respect to the unknown
508  du_ddata[count] = psi[l];
509  //Increase the counter
510  ++count;
511  }
512  }
513  }
514 
515 
516  /// \short Self-test: Return 0 for OK
517  unsigned self_test();
518 
519 protected:
520 
521 
522  /// \short Shape/test functions and derivs w.r.t. to global coords at
523  /// local coord. s; return Jacobian of mapping
524  virtual double dshape_and_dtest_eulerian_adv_diff(const Vector<double> &s,
525  Shape &psi,
526  DShape &dpsidx,
527  Shape &test,
528  DShape &dtestdx) const=0;
529 
530  /// \short Shape/test functions and derivs w.r.t. to global coords at
531  /// integration point ipt; return Jacobian of mapping
533  const unsigned &ipt,
534  Shape &psi,
535  DShape &dpsidx,
536  Shape &test,
537  DShape &dtestdx)
538  const=0;
539 
540  /// \short Add the element's contribution to its residual vector only
541  /// (if flag=and/or element Jacobian matrix
543  Vector<double> &residuals, DenseMatrix<double> &jacobian,
544  DenseMatrix<double> &mass_matrix, unsigned flag);
545 
546  /// Pointer to global Peclet number
547  double *Pe_pt;
548 
549  /// Pointer to global Peclet number multiplied by Strouhal number
550  double *PeSt_pt;
551 
552  /// Pointer to source function:
554 
555  /// Pointer to wind function:
557 
558  /// \short Boolean flag to indicate if ALE formulation is disabled when
559  /// time-derivatives are computed. Only set to true if you're sure
560  /// that the mesh is stationary.
562 
563  private:
564 
565  /// Static default value for the Peclet number
566  static double Default_peclet_number;
567 
568 
569 };
570 
571 
572 ///////////////////////////////////////////////////////////////////////////
573 ///////////////////////////////////////////////////////////////////////////
574 ///////////////////////////////////////////////////////////////////////////
575 
576 
577 
578 //======================================================================
579 /// \short QAdvectionDiffusionElement elements are
580 /// linear/quadrilateral/brick-shaped Advection Diffusion elements with
581 /// isoparametric interpolation for the function.
582 //======================================================================
583 template <unsigned DIM, unsigned NNODE_1D>
584  class QAdvectionDiffusionElement : public virtual QElement<DIM,NNODE_1D>,
585  public virtual AdvectionDiffusionEquations<DIM>
586 {
587 
588 private:
589 
590  /// \short Static array of ints to hold number of variables at
591  /// nodes: Initial_Nvalue[n]
592  static const unsigned Initial_Nvalue;
593 
594  public:
595 
596 
597  ///\short Constructor: Call constructors for QElement and
598  /// Advection Diffusion equations
599  QAdvectionDiffusionElement() : QElement<DIM,NNODE_1D>(),
601  { }
602 
603  /// Broken copy constructor
605  dummy)
606  {
607  BrokenCopy::broken_copy("QAdvectionDiffusionElement");
608  }
609 
610  /// Broken assignment operator
612  {
613  BrokenCopy::broken_assign("QAdvectionDiffusionElement");
614  }
615 
616  /// \short Required # of `values' (pinned or dofs)
617  /// at node n
618  inline unsigned required_nvalue(const unsigned &n) const
619  {return Initial_Nvalue;}
620 
621  /// \short Output function:
622  /// x,y,u or x,y,z,u
623  void output(std::ostream &outfile)
625 
626  /// \short Output function:
627  /// x,y,u or x,y,z,u at n_plot^DIM plot points
628  void output(std::ostream &outfile, const unsigned &n_plot)
630 
631 
632  /// \short C-style output function:
633  /// x,y,u or x,y,z,u
634  void output(FILE* file_pt)
635  {
637  }
638 
639  /// \short C-style output function:
640  /// x,y,u or x,y,z,u at n_plot^DIM plot points
641  void output(FILE* file_pt, const unsigned &n_plot)
642  {
644  }
645 
646  /// \short Output function for an exact solution:
647  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
648  void output_fct(std::ostream &outfile, const unsigned &n_plot,
650  exact_soln_pt)
651  {AdvectionDiffusionEquations<DIM>::output_fct(outfile,n_plot,exact_soln_pt);}
652 
653 
654  /// \short Output function for a time-dependent exact solution.
655  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
656  /// (Calls the steady version)
657  void output_fct(std::ostream &outfile, const unsigned &n_plot,
658  const double& time,
660  exact_soln_pt)
661  {
663  output_fct(outfile,n_plot,time,exact_soln_pt);
664  }
665 
666 
667 protected:
668 
669  /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
671  Shape &psi,
672  DShape &dpsidx,
673  Shape &test,
674  DShape &dtestdx) const;
675 
676  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
677  /// integration point ipt. Return Jacobian.
678  inline double dshape_and_dtest_eulerian_at_knot_adv_diff(const unsigned& ipt,
679  Shape &psi,
680  DShape &dpsidx,
681  Shape &test,
682  DShape &dtestdx)
683  const;
684 
685 };
686 
687 //Inline functions:
688 
689 
690 //======================================================================
691 /// \short Define the shape functions and test functions and derivatives
692 /// w.r.t. global coordinates and return Jacobian of mapping.
693 ///
694 /// Galerkin: Test functions = shape functions
695 //======================================================================
696 template<unsigned DIM, unsigned NNODE_1D>
699  Shape &psi,
700  DShape &dpsidx,
701  Shape &test,
702  DShape &dtestdx) const
703 {
704  //Call the geometrical shape functions and derivatives
705  double J = this->dshape_eulerian(s,psi,dpsidx);
706 
707  //Loop over the test functions and derivatives and set them equal to the
708  //shape functions
709  for(unsigned i=0;i<NNODE_1D;i++)
710  {
711  test[i] = psi[i];
712  for(unsigned j=0;j<DIM;j++)
713  {
714  dtestdx(i,j) = dpsidx(i,j);
715  }
716  }
717 
718  //Return the jacobian
719  return J;
720 }
721 
722 
723 
724 //======================================================================
725 /// Define the shape functions and test functions and derivatives
726 /// w.r.t. global coordinates and return Jacobian of mapping.
727 ///
728 /// Galerkin: Test functions = shape functions
729 //======================================================================
730 template<unsigned DIM, unsigned NNODE_1D>
733  const unsigned &ipt,
734  Shape &psi,
735  DShape &dpsidx,
736  Shape &test,
737  DShape &dtestdx) const
738 {
739  //Call the geometrical shape functions and derivatives
740  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
741 
742  //Set the test functions equal to the shape functions (pointer copy)
743  test = psi;
744  dtestdx = dpsidx;
745 
746  //Return the jacobian
747  return J;
748 }
749 
750 
751 ////////////////////////////////////////////////////////////////////////
752 ////////////////////////////////////////////////////////////////////////
753 ////////////////////////////////////////////////////////////////////////
754 
755 
756 
757 //=======================================================================
758 /// \short Face geometry for the QAdvectionDiffusionElement elements:
759 /// The spatial dimension of the face elements is one lower than that
760 /// of the bulk element but they have the same number of points along
761 /// their 1D edges.
762 //=======================================================================
763 template<unsigned DIM, unsigned NNODE_1D>
765  public virtual QElement<DIM-1,NNODE_1D>
766 {
767 
768  public:
769 
770  /// \short Constructor: Call the constructor for the
771  /// appropriate lower-dimensional QElement
772  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
773 
774 };
775 
776 
777 
778 ////////////////////////////////////////////////////////////////////////
779 ////////////////////////////////////////////////////////////////////////
780 ////////////////////////////////////////////////////////////////////////
781 
782 
783 //=======================================================================
784 /// Face geometry for the 1D QAdvectionDiffusion elements: Point elements
785 //=======================================================================
786 template<unsigned NNODE_1D>
788  public virtual PointElement
789 {
790 
791  public:
792 
793  /// \short Constructor: Call the constructor for the
794  /// appropriate lower-dimensional QElement
796 
797 };
798 
799 }
800 
801 #endif
virtual void fill_in_generic_residual_contribution_adv_diff(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add the element's contribution to its residual vector only (if flag=and/or element Jacobian matrix...
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
double dshape_and_dtest_eulerian_at_knot_adv_diff(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 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 nplot points in each "coordinate direc...
Definition: elements.h:2938
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
virtual void get_wind_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Get wind at (Eulerian) position x and/or local coordinate s. This function is virtual to allow overlo...
void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot ...
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
AdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
AdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
void(* AdvectionDiffusionWindFctPt)(const Vector< double > &x, Vector< double > &wind)
Function pointer to wind function fct(x,w(x)) – x is a Vector!
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
cstr elem_len * i
Definition: cfortran.h:607
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points...
void(* AdvectionDiffusionSourceFctPt)(const Vector< double > &x, double &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
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
QAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Advection Diffusion element...
A general Finite Element class.
Definition: elements.h:1271
const double & pe() const
Peclet number.
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative. Note: By default, ALE is enabled, at the expense of possibly creating unnecessary work in problems where the mesh is, in fact, stationary.
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:557
char t
Definition: cfortran.h:572
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
AdvectionDiffusionSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
virtual void output_fct(std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points (dummy time-dependent versio...
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:375
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
double * Pe_pt
Pointer to global Peclet number.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points.
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2700
double du_dt_adv_diff(const unsigned &n) const
du/dt at local node n.
A class for all elements that solve the Advection Diffusion equations using isoparametric elements...
virtual unsigned u_index_adv_diff() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...
virtual void dinterpolated_u_adv_diff_ddata(const Vector< double > &s, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
Return derivative of u at point s with respect to all data that can affect its value. In addition, return the global equation numbers corresponding to the data. This is virtual so that it can be overloaded in the refineable version.
AdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
double *& pe_pt()
Pointer to Peclet number.
double dshape_and_dtest_eulerian_adv_diff(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.
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
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
QAdvectionDiffusionElement(const QAdvectionDiffusionElement< DIM, NNODE_1D > &dummy)
Broken copy constructor.
virtual double dshape_and_dtest_eulerian_adv_diff(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at local coord. s; return Jacobian of mapping...
static double Default_peclet_number
Static default value for the Peclet number.
static char t char * s
Definition: cfortran.h:572
AdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and the element Jacobian matrix (wrapper) ...
virtual void get_source_adv_diff(const unsigned &ipt, const Vector< double > &x, double &source) const
Get source term at (Eulerian) position x. This function is virtual to allow overloading in multi-phys...
unsigned required_nvalue(const unsigned &n) const
Required # of `values' (pinned or dofs) at node n.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
AdvectionDiffusionEquations(const AdvectionDiffusionEquations &dummy)
Broken copy constructor.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
double interpolated_u_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
void operator=(const AdvectionDiffusionEquations &)
Broken assignment operator.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void output(std::ostream &outfile)
Output with default number of plot points.
void output(FILE *file_pt)
C_style output with default number of plot points.
unsigned self_test()
Self-test: Return 0 for OK.
virtual double dshape_and_dtest_eulerian_at_knot_adv_diff(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
AdvectionDiffusionWindFctPt wind_fct_pt() const
Access function: Pointer to wind function. Const version.
AdvectionDiffusionEquations()
Constructor: Initialise the Source_fct_pt and Wind_fct_pt to null and set (pointer to) Peclet number ...
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
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
QAdvectionDiffusionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
static const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
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.
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.
const double & pe_st() const
Peclet number multiplied by Strouhal number.
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
void operator=(const QAdvectionDiffusionElement< DIM, NNODE_1D > &)
Broken assignment operator.
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types)...