unsteady_heat_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 UnsteadyHeat elements
31 #ifndef OOMPH_UNSTEADY_HEAT_ELEMENTS_HEADER
32 #define OOMPH_UNSTEADY_HEAT_ELEMENTS_HEADER
33 
34 // Config header generated by autoconfig
35 #ifdef HAVE_CONFIG_H
36  #include <oomph-lib-config.h>
37 #endif
38 
39 
40 //OOMPH-LIB headers
41 #include "../generic/projection.h"
42 #include "../generic/nodes.h"
43 #include "../generic/Qelements.h"
44 #include "../generic/oomph_utilities.h"
45 
46 
47 namespace oomph
48 {
49 
50 /// Base class so that we don't need to know the dimension just to set the
51 /// source function!
53 {
54 public:
55  /// \short Function pointer to source function fct(t,x,f(x,t)) --
56  /// x is a Vector!
57  typedef void (*UnsteadyHeatSourceFctPt)(const double& time,
58  const Vector<double>& x,
59  double& u);
60 
61  /// Access function: Pointer to source function
63 };
64 
65 //=============================================================
66 /// A class for all isoparametric elements that solve the
67 /// UnsteadyHeat equations.
68 /// \f[
69 /// \frac{\partial^2 u}{\partial x_i^2}=\frac{\partial u}{\partial t}+f(t,x_j)
70 /// \f]
71 /// This contains the generic maths. Shape functions, geometric
72 /// mapping etc. must get implemented in derived class.
73 /// Note that this class assumes an isoparametric formulation, i.e. that
74 /// the scalar unknown is interpolated using the same shape funcitons
75 /// as the position.
76 //=============================================================
77 template <unsigned DIM>
79 {
80 
81 public:
82 
83  /// \short Function pointer to source function fct(t,x,f(x,t)) --
84  /// x is a Vector!
85  typedef void (*UnsteadyHeatSourceFctPt)(const double& time,
86  const Vector<double>& x,
87  double& u);
88 
89 
90  /// \short Constructor: Initialises the Source_fct_pt to null and
91  /// sets flag to use ALE formulation of the equations.
92  /// Also set Alpha (thermal inertia) and Beta (thermal conductivity)
93  /// parameters to defaults (both one for natural scaling)
95  {
96  //Set Alpha and Beta parameter to default (one for natural scaling of time)
99  }
100 
101 
102  /// Broken copy constructor
104  {
105  BrokenCopy::broken_copy("UnsteadyHeatEquations");
106  }
107 
108  /// Broken assignment operator
109 //Commented out broken assignment operator because this can lead to a conflict warning
110 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
111 //realise that two separate implementations of the broken function are the same and so,
112 //quite rightly, it shouts.
113  /*void operator=(const UnsteadyHeatEquations&)
114  {
115  BrokenCopy::broken_assign("UnsteadyHeatEquations");
116  }*/
117 
118  /// \short Return the index at which the unknown value
119  /// is stored. The default value, 0, is appropriate for single-physics
120  /// problems, when there is only one variable, the value that satisfies the
121  /// unsteady heat equation.
122  /// In derived multi-physics elements, this function should be overloaded
123  /// to reflect the chosen storage scheme. Note that these equations require
124  /// that the unknown is always stored at the same index at each node.
125  virtual inline unsigned u_index_ust_heat() const {return 0;}
126 
127  /// \short du/dt at local node n.
128  /// Uses suitably interpolated value for hanging nodes.
129  double du_dt_ust_heat(const unsigned &n) const
130  {
131  // Get the data's timestepper
133 
134  //Initialise dudt
135  double dudt=0.0;
136 
137  //Loop over the timesteps, if there is a non Steady timestepper
138  if (!time_stepper_pt->is_steady())
139  {
140  //Find the index at which the variable is stored
141  const unsigned u_nodal_index = u_index_ust_heat();
142 
143  // Number of timsteps (past & present)
144  const unsigned n_time = time_stepper_pt->ntstorage();
145 
146  //Add the contributions to the time derivative
147  for(unsigned t=0;t<n_time;t++)
148  {
149  dudt += time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
150  }
151  }
152  return dudt;
153  }
154 
155  /// \short Disable ALE, i.e. assert the mesh is not moving -- you do this
156  /// at your own risk!
157  void disable_ALE()
158  {
159  ALE_is_disabled=true;
160  }
161 
162 
163  /// \short (Re-)enable ALE, i.e. take possible mesh motion into account
164  /// when evaluating the time-derivative. Note: By default, ALE is
165  /// enabled, at the expense of possibly creating unnecessary work
166  /// in problems where the mesh is, in fact, stationary.
167  void enable_ALE()
168  {
169  ALE_is_disabled=false;
170  }
171 
172  /// Compute norm of fe solution
173  void compute_norm(double& norm);
174 
175  /// Output with default number of plot points
176  void output(std::ostream &outfile)
177  {
178  unsigned nplot=5;
179  output(outfile,nplot);
180  }
181 
182 
183  /// \short Output FE representation of soln: x,y,u or x,y,z,u at
184  /// n_plot^DIM plot points
185  void output(std::ostream &outfile, const unsigned &nplot);
186 
187  /// C_style output with default number of plot points
188  void output(FILE* file_pt)
189  {
190  unsigned n_plot=5;
191  output(file_pt,n_plot);
192  }
193 
194 
195  /// \short C-style output FE representation of soln: x,y,u or x,y,z,u at
196  /// n_plot^DIM plot points
197  void output(FILE* file_pt, const unsigned &n_plot);
198 
199 
200  /// Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points
201  void output_fct(std::ostream &outfile, const unsigned &nplot,
203  exact_soln_pt);
204 
205 
206  /// \short Output exact soln: x,y,u_exact or x,y,z,u_exact at
207  /// nplot^DIM plot points (time-dependent version)
208  virtual
209  void output_fct(std::ostream &outfile, const unsigned &nplot,
210  const double& time,
212  exact_soln_pt);
213 
214 
215  /// Get error against and norm of exact solution
216  void compute_error(std::ostream &outfile,
218  exact_soln_pt,
219  double& error, double& norm);
220 
221 
222  /// Get error against and norm of exact solution
223  void compute_error(std::ostream &outfile,
225  exact_soln_pt,
226  const double& time, double& error, double& norm);
227 
228 
229  /// Access function: Pointer to source function
231 
232 
233  /// Access function: Pointer to source function. Const version
235 
236 
237  /// \short Get source term at continous time t and (Eulerian) position x.
238  /// Virtual so it can be overloaded in derived multiphysics elements.
239  virtual inline void get_source_ust_heat(const double& t,
240  const unsigned& ipt,
241  const Vector<double>& x,
242  double& source) const
243  {
244  //If no source function has been set, return zero
245  if(Source_fct_pt==0) {source = 0.0;}
246  else
247  {
248  // Get source strength
249  (*Source_fct_pt)(t,x,source);
250  }
251  }
252 
253  /// Alpha parameter (thermal inertia)
254  const double &alpha() const {return *Alpha_pt;}
255 
256  /// Pointer to Alpha parameter (thermal inertia)
257  double* &alpha_pt() {return Alpha_pt;}
258 
259 
260  /// Beta parameter (thermal conductivity)
261  const double &beta() const {return *Beta_pt;}
262 
263  /// Pointer to Beta parameter (thermal conductivity)
264  double* &beta_pt() {return Beta_pt;}
265 
266  /// Get flux: flux[i] = du/dx_i
267  void get_flux(const Vector<double>& s, Vector<double>& flux) const
268  {
269  //Find out how many nodes there are in the element
270  unsigned n_node = nnode();
271 
272  //Find the index at which the variable is stored
273  unsigned u_nodal_index = u_index_ust_heat();
274 
275  //Set up memory for the shape and test functions
276  Shape psi(n_node);
277  DShape dpsidx(n_node,DIM);
278 
279  //Call the derivatives of the shape and test functions
280  dshape_eulerian(s,psi,dpsidx);
281 
282  //Initialise to zero
283  for(unsigned j=0;j<DIM;j++) {flux[j] = 0.0;}
284 
285  // Loop over nodes
286  for(unsigned l=0;l<n_node;l++)
287  {
288  //Loop over derivative directions
289  for(unsigned j=0;j<DIM;j++)
290  {
291  flux[j] += nodal_value(l,u_nodal_index)*dpsidx(l,j);
292  }
293  }
294  }
295 
296 
297  /// Compute element residual Vector (wrapper)
299  {
300  //Call the generic residuals function with flag set to 0
301  //using a dummy matrix argument
304  }
305 
306 
307  /// Compute element residual Vector and element Jacobian matrix (wrapper)
309  DenseMatrix<double> &jacobian)
310  {
311  //Call the generic routine with the flag set to 1
313  }
314 
315 
316  /// Return FE representation of function value u(s) at local coordinate s
317  inline double interpolated_u_ust_heat(const Vector<double> &s) const
318  {
319  //Find number of nodes
320  unsigned n_node = nnode();
321 
322  //Find the index at which the variable is stored
323  unsigned u_nodal_index = u_index_ust_heat();
324 
325  //Local shape function
326  Shape psi(n_node);
327 
328  //Find values of shape function
329  shape(s,psi);
330 
331  //Initialise value of u
332  double interpolated_u = 0.0;
333 
334  //Loop over the local nodes and sum
335  for(unsigned l=0;l<n_node;l++)
336  {
337  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
338  }
339 
340  return(interpolated_u);
341  }
342 
343 
344  /// \short Return FE representation of function value u(s) at local
345  /// coordinate s at previous time t (t=0: present)
346  inline double interpolated_u_ust_heat(const unsigned& t,
347  const Vector<double> &s) const
348  {
349  //Find number of nodes
350  unsigned n_node = nnode();
351 
352  //Find the index at which the variable is stored
353  unsigned u_nodal_index = u_index_ust_heat();
354 
355  //Local shape function
356  Shape psi(n_node);
357 
358  //Find values of shape function
359  shape(s,psi);
360 
361  //Initialise value of u
362  double interpolated_u = 0.0;
363 
364  //Loop over the local nodes and sum
365  for(unsigned l=0;l<n_node;l++)
366  {
367  interpolated_u += nodal_value(t,l,u_nodal_index)*psi[l];
368  }
369 
370  return(interpolated_u);
371  }
372 
373 
374  /// Return FE representation of function value du/dt(s) at local coordinate s
375  inline double interpolated_du_dt_ust_heat(const Vector<double> &s) const
376  {
377  //Find number of nodes
378  unsigned n_node = nnode();
379 
380  //Local shape function
381  Shape psi(n_node);
382 
383  //Find values of shape function
384  shape(s,psi);
385 
386  //Initialise value of du/dt
387  double interpolated_dudt = 0.0;
388 
389  //Loop over the local nodes and sum
390  for(unsigned l=0;l<n_node;l++)
391  {
392  interpolated_dudt += du_dt_ust_heat(l)*psi[l];
393  }
394 
395  return(interpolated_dudt);
396  }
397 
398 
399  /// \short Self-test: Return 0 for OK
400  unsigned self_test();
401 
402 
403 protected:
404 
405  /// \short Shape/test functions and derivs w.r.t. to global coords at
406  /// local coord. s; return Jacobian of mapping
407  virtual double dshape_and_dtest_eulerian_ust_heat(const Vector<double> &s,
408  Shape &psi,
409  DShape &dpsidx,
410  Shape &test,
411  DShape &dtestdx) const=0;
412 
413 
414  /// \short Shape/test functions and derivs w.r.t. to global coords at
415  /// integration point ipt; return Jacobian of mapping
417  const unsigned &ipt,
418  Shape &psi,
419  DShape &dpsidx,
420  Shape &test,
421  DShape &dtestdx)
422  const=0;
423 
424  /// \short Compute element residual Vector only (if flag=and/or element
425  /// Jacobian matrix
427  Vector<double> &residuals, DenseMatrix<double> &jacobian,
428  unsigned flag);
429 
430  /// Pointer to source function:
432 
433  /// \short Boolean flag to indicate if ALE formulation is disabled when
434  /// time-derivatives are computed. Only set to true if you're sure
435  /// that the mesh is stationary.
437 
438  /// Pointer to Alpha parameter (thermal inertia)
439  double *Alpha_pt;
440 
441  /// Pointer to Beta parameter (thermal conductivity)
442  double *Beta_pt;
443 
444  private:
445 
446  /// \short Static default value for the Alpha parameter:
447  /// (thermal inertia): One for natural scaling
448  static double Default_alpha_parameter;
449 
450  /// \short Static default value for the Beta parameter (thermal
451  /// conductivity): One for natural scaling
452  static double Default_beta_parameter;
453 
454 };
455 
456 
457 
458 
459 
460 
461 ///////////////////////////////////////////////////////////////////////////
462 ///////////////////////////////////////////////////////////////////////////
463 ///////////////////////////////////////////////////////////////////////////
464 
465 
466 
467 //======================================================================
468 /// QUnsteadyHeatElement elements are linear/quadrilateral/brick-shaped
469 /// UnsteadyHeat elements with isoparametric interpolation for the function.
470 //======================================================================
471 template <unsigned DIM, unsigned NNODE_1D>
472  class QUnsteadyHeatElement : public virtual QElement<DIM,NNODE_1D>,
473  public virtual UnsteadyHeatEquations<DIM>
474 {
475  private:
476 
477  /// \short Static array of ints to hold number of variables at
478  /// nodes: Initial_Nvalue[n]
479  static const unsigned Initial_Nvalue;
480 
481  public:
482 
483  ///\short Constructor: Call constructors for QElement and
484  /// UnsteadyHeat equations
485  QUnsteadyHeatElement() : QElement<DIM,NNODE_1D>(),
486  UnsteadyHeatEquations<DIM>()
487  { }
488 
489  /// Broken copy constructor
491  {
492  BrokenCopy::broken_copy("QUnsteadyHeatElement");
493  }
494 
495  /// Broken assignment operator
496  /*void operator=(const QUnsteadyHeatElement<DIM,NNODE_1D>&)
497  {
498  BrokenCopy::broken_assign("QUnsteadyHeatElement");
499  }*/
500 
501  /// \short Required # of `values' (pinned or dofs)
502  /// at node n
503  inline unsigned required_nvalue(const unsigned &n) const
504  {return Initial_Nvalue;}
505 
506  /// \short Output function:
507  /// x,y,u or x,y,z,u
508  void output(std::ostream &outfile)
510 
511 
512  /// \short Output function:
513  /// x,y,u or x,y,z,u at n_plot^DIM plot points
514  void output(std::ostream &outfile, const unsigned &n_plot)
515  {UnsteadyHeatEquations<DIM>::output(outfile,n_plot);}
516 
517 
518  /// \short C-style output function:
519  /// x,y,u or x,y,z,u
520  void output(FILE* file_pt)
522 
523 
524  /// \short C-style output function:
525  /// x,y,u or x,y,z,u at n_plot^DIM plot points
526  void output(FILE* file_pt, const unsigned &n_plot)
527  {UnsteadyHeatEquations<DIM>::output(file_pt,n_plot);}
528 
529 
530  /// \short Output function for an exact solution:
531  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
532  void output_fct(std::ostream &outfile, const unsigned &n_plot,
534  exact_soln_pt)
535  {UnsteadyHeatEquations<DIM>::output_fct(outfile,n_plot,exact_soln_pt);}
536 
537 
538 
539  /// \short Output function for a time-dependent exact solution.
540  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
541  /// (Calls the steady version)
542  void output_fct(std::ostream &outfile, const unsigned &n_plot,
543  const double& time,
545  {UnsteadyHeatEquations<DIM>::output_fct(outfile,n_plot,time,exact_soln_pt);}
546 
547 
548 
549 protected:
550 
551  /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
553  Shape &psi,
554  DShape &dpsidx,
555  Shape &test,
556  DShape &dtestdx) const;
557 
558 
559  /// \short Shape/test functions and derivs w.r.t. to global coords at
560  /// integration point ipt; return Jacobian of mapping
561  inline double dshape_and_dtest_eulerian_at_knot_ust_heat(const unsigned &ipt,
562  Shape &psi,
563  DShape &dpsidx,
564  Shape &test,
565  DShape &dtestdx)
566  const;
567 
568 };
569 
570 
571 //Inline functions:
572 
573 
574 //======================================================================
575 /// Define the shape functions and test functions and derivatives
576 /// w.r.t. global coordinates and return Jacobian of mapping.
577 ///
578 /// Galerkin: Test functions = shape functions
579 //======================================================================
580 template<unsigned DIM, unsigned NNODE_1D>
583  Shape &psi,
584  DShape &dpsidx,
585  Shape &test,
586  DShape &dtestdx) const
587 {
588  //Call the geometrical shape functions and derivatives
589  double J = this->dshape_eulerian(s,psi,dpsidx);
590 
591  //Loop over the test functions and derivatives and set them equal to the
592  //shape functions
593  for(unsigned i=0;i<NNODE_1D;i++)
594  {
595  test[i] = psi[i];
596  for(unsigned j=0;j<DIM;j++)
597  {
598  dtestdx(i,j) = dpsidx(i,j);
599  }
600  }
601 
602  //Return the jacobian
603  return J;
604 }
605 
606 
607 //======================================================================
608 /// Define the shape functions and test functions and derivatives
609 /// w.r.t. global coordinates and return Jacobian of mapping.
610 ///
611 /// Galerkin: Test functions = shape functions
612 //======================================================================
613 template<unsigned DIM,unsigned NNODE_1D>
616  const unsigned &ipt,
617  Shape &psi,
618  DShape &dpsidx,
619  Shape &test,
620  DShape &dtestdx) const
621 {
622  //Call the geometrical shape functions and derivatives
623  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
624 
625  //Set the test functions equal to the shape functions
626  //(sets internal pointers)
627  test = psi;
628  dtestdx = dpsidx;
629 
630  //Return the jacobian
631  return J;
632 }
633 
634 
635 ////////////////////////////////////////////////////////////////////////
636 ////////////////////////////////////////////////////////////////////////
637 
638 
639 
640 //=======================================================================
641 /// Face geometry for the QUnsteadyHeatElement elements: The spatial
642 /// dimension of the face elements is one lower than that of the
643 /// bulk element but they have the same number of points
644 /// along their 1D edges.
645 //=======================================================================
646 template<unsigned DIM, unsigned NNODE_1D>
647 class FaceGeometry<QUnsteadyHeatElement<DIM,NNODE_1D> >:
648  public virtual QElement<DIM-1,NNODE_1D>
649 {
650 
651  public:
652 
653  /// \short Constructor: Call the constructor for the
654  /// appropriate lower-dimensional QElement
655  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
656 
657 };
658 
659 ////////////////////////////////////////////////////////////////////////
660 ////////////////////////////////////////////////////////////////////////
661 ////////////////////////////////////////////////////////////////////////
662 
663 
664 //=======================================================================
665 /// Face geometry for the 1D QUnsteadyHeatElement elements: Point elements
666 //=======================================================================
667 template<unsigned NNODE_1D>
668 class FaceGeometry<QUnsteadyHeatElement<1,NNODE_1D> >:
669  public virtual PointElement
670 {
671 
672  public:
673 
674  /// \short Constructor: Call the constructor for the
675  /// appropriate lower-dimensional QElement
677 
678 };
679 
680 
681 ////////////////////////////////////////////////////////////////////////
682 ////////////////////////////////////////////////////////////////////////
683 ////////////////////////////////////////////////////////////////////////
684 
685 
686 
687 //==========================================================
688 /// UnsteadyHeat upgraded to become projectable
689 //==========================================================
690  template<class UNSTEADY_HEAT_ELEMENT>
692  public virtual ProjectableElement<UNSTEADY_HEAT_ELEMENT>
693  {
694 
695  public:
696 
697 
698  /// \short Constructor [this was only required explicitly
699  /// from gcc 4.5.2 onwards...]
701 
702  /// \short Specify the values associated with field fld.
703  /// The information is returned in a vector of pairs which comprise
704  /// the Data object and the value within it, that correspond to field fld.
706  {
707 
708 #ifdef PARANOID
709  if (fld!=0)
710  {
711  std::stringstream error_stream;
712  error_stream
713  << "UnsteadyHeat elements only store a single field so fld must be 0 rather"
714  << " than " << fld << std::endl;
715  throw OomphLibError(
716  error_stream.str(),
717  OOMPH_CURRENT_FUNCTION,
718  OOMPH_EXCEPTION_LOCATION);
719  }
720 #endif
721 
722  // Create the vector
723  unsigned nnod=this->nnode();
724  Vector<std::pair<Data*,unsigned> > data_values(nnod);
725 
726  // Loop over all nodes
727  for (unsigned j=0;j<nnod;j++)
728  {
729  // Add the data value associated field: The node itself
730  data_values[j]=std::make_pair(this->node_pt(j),fld);
731  }
732 
733  // Return the vector
734  return data_values;
735  }
736 
737  /// \short Number of fields to be projected: Just one
739  {
740  return 1;
741  }
742 
743  /// \short Number of history values to be stored for fld-th field.
744  /// (Note: count includes current value!)
745  unsigned nhistory_values_for_projection(const unsigned &fld)
746  {
747 #ifdef PARANOID
748  if (fld!=0)
749  {
750  std::stringstream error_stream;
751  error_stream
752  << "UnsteadyHeat elements only store a single field so fld must be 0 rather"
753  << " than " << fld << std::endl;
754  throw OomphLibError(
755  error_stream.str(),
756  OOMPH_CURRENT_FUNCTION,
757  OOMPH_EXCEPTION_LOCATION);
758  }
759 #endif
760  return this->node_pt(0)->ntstorage();
761  }
762 
763  ///\short Number of positional history values
764  /// (Note: count includes current value!)
766  {
767  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
768  }
769 
770  /// \short Return Jacobian of mapping and shape functions of field fld
771  /// at local coordinate s
772  double jacobian_and_shape_of_field(const unsigned &fld,
773  const Vector<double> &s,
774  Shape &psi)
775  {
776 #ifdef PARANOID
777  if (fld!=0)
778  {
779  std::stringstream error_stream;
780  error_stream
781  << "UnsteadyHeat elements only store a single field so fld must be 0 rather"
782  << " than " << fld << std::endl;
783  throw OomphLibError(
784  error_stream.str(),
785  OOMPH_CURRENT_FUNCTION,
786  OOMPH_EXCEPTION_LOCATION);
787  }
788 #endif
789  unsigned n_dim=this->dim();
790  unsigned n_node=this->nnode();
791  Shape test(n_node);
792  DShape dpsidx(n_node,n_dim), dtestdx(n_node,n_dim);
793  double J=this->dshape_and_dtest_eulerian_ust_heat(s,psi,dpsidx,
794  test,dtestdx);
795  return J;
796  }
797 
798 
799 
800  /// \short Return interpolated field fld at local coordinate s, at time level
801  /// t (t=0: present; t>0: history values)
802  double get_field(const unsigned &t,
803  const unsigned &fld,
804  const Vector<double>& s)
805  {
806 #ifdef PARANOID
807  if (fld!=0)
808  {
809  std::stringstream error_stream;
810  error_stream
811  << "UnsteadyHeat elements only store a single field so fld must be 0 rather"
812  << " than " << fld << std::endl;
813  throw OomphLibError(
814  error_stream.str(),
815  OOMPH_CURRENT_FUNCTION,
816  OOMPH_EXCEPTION_LOCATION);
817  }
818 #endif
819  //Find the index at which the variable is stored
820  unsigned u_nodal_index = this->u_index_ust_heat();
821 
822  //Local shape function
823  unsigned n_node=this->nnode();
824  Shape psi(n_node);
825 
826  //Find values of shape function
827  this->shape(s,psi);
828 
829  //Initialise value of u
830  double interpolated_u = 0.0;
831 
832  //Sum over the local nodes
833  for(unsigned l=0;l<n_node;l++)
834  {
835  interpolated_u += this->nodal_value(t,l,u_nodal_index)*psi[l];
836  }
837  return interpolated_u;
838  }
839 
840 
841 
842 
843  ///Return number of values in field fld: One per node
844  unsigned nvalue_of_field(const unsigned &fld)
845  {
846 #ifdef PARANOID
847  if (fld!=0)
848  {
849  std::stringstream error_stream;
850  error_stream
851  << "UnsteadyHeat elements only store a single field so fld must be 0 rather"
852  << " than " << fld << std::endl;
853  throw OomphLibError(
854  error_stream.str(),
855  OOMPH_CURRENT_FUNCTION,
856  OOMPH_EXCEPTION_LOCATION);
857  }
858 #endif
859  return this->nnode();
860  }
861 
862 
863  ///Return local equation number of value j in field fld.
864  int local_equation(const unsigned &fld,
865  const unsigned &j)
866  {
867 #ifdef PARANOID
868  if (fld!=0)
869  {
870  std::stringstream error_stream;
871  error_stream
872  << "UnsteadyHeat elements only store a single field so fld must be 0 rather"
873  << " than " << fld << std::endl;
874  throw OomphLibError(
875  error_stream.str(),
876  OOMPH_CURRENT_FUNCTION,
877  OOMPH_EXCEPTION_LOCATION);
878  }
879 #endif
880  const unsigned u_nodal_index = this->u_index_ust_heat();
881  return this->nodal_local_eqn(j,u_nodal_index);
882  }
883 
884 
885 
886 
887  /// \short Output FE representation of soln: x,y,u or x,y,z,u at
888  /// and history values at n_plot^DIM plot points
889  void output(std::ostream &outfile, const unsigned &nplot)
890  {
891  unsigned el_dim=this->dim();
892  //Vector of local coordinates
893  Vector<double> s(el_dim);
894 
895  // Tecplot header info
896  outfile << this->tecplot_zone_string(nplot);
897 
898  // Loop over plot points
899  unsigned num_plot_points=this->nplot_points(nplot);
900  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
901  {
902  // Get local coordinates of plot point
903  this->get_s_plot(iplot,nplot,s);
904 
905  for(unsigned i=0;i<el_dim;i++)
906  {
907  outfile << this->interpolated_x(s,i) << " ";
908  }
909  outfile << this->interpolated_u_ust_heat(s) << " ";
910  outfile << this->interpolated_du_dt_ust_heat(s) << " ";
911 
912 
913  // History values of coordinates
914  unsigned n_prev=this->node_pt(0)->position_time_stepper_pt()->ntstorage();
915  for (unsigned t=1;t<n_prev;t++)
916  {
917  for(unsigned i=0;i<el_dim;i++)
918  {
919  outfile << this->interpolated_x(t,s,i) << " ";
920  }
921  }
922 
923  // History values of velocities
924  n_prev=this->node_pt(0)->time_stepper_pt()->ntstorage();
925  for (unsigned t=1;t<n_prev;t++)
926  {
927  outfile << this->interpolated_u_ust_heat(t,s) << " ";
928  }
929  outfile << std::endl;
930  }
931 
932 
933 
934  // Write tecplot footer (e.g. FE connectivity lists)
935  this->write_tecplot_zone_footer(outfile,nplot);
936  }
937 
938 
939  };
940 
941 
942 //=======================================================================
943 /// Face geometry for element is the same as that for the underlying
944 /// wrapped element
945 //=======================================================================
946  template<class ELEMENT>
948  : public virtual FaceGeometry<ELEMENT>
949  {
950  public:
951  FaceGeometry() : FaceGeometry<ELEMENT>() {}
952  };
953 
954 
955 //=======================================================================
956 /// Face geometry of the Face Geometry for element is the same as
957 /// that for the underlying wrapped element
958 //=======================================================================
959  template<class ELEMENT>
961  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
962  {
963  public:
965  };
966 
967 
968 
969 
970 }
971 
972 #endif
virtual unsigned u_index_ust_heat() const
Broken assignment operator.
double du_dt_ust_heat(const unsigned &n) const
du/dt at local node n. Uses suitably interpolated value for hanging nodes.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
double *& beta_pt()
Pointer to Beta parameter (thermal conductivity)
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 broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute element residual Vector (wrapper)
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
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
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.
double dshape_and_dtest_eulerian_ust_heat(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 *& alpha_pt()
Pointer to Alpha parameter (thermal inertia)
virtual void get_source_ust_heat(const double &t, const unsigned &ipt, const Vector< double > &x, double &source) const
Get source term at continous time t and (Eulerian) position x. Virtual so it can be overloaded in der...
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 compute_norm(double &norm)
Compute norm of fe solution.
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:847
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
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
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
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 ...
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
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
const double & alpha() const
Alpha parameter (thermal inertia)
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
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
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
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
unsigned nfields_for_projection()
Number of fields to be projected: Just one.
void(* UnsteadyHeatSourceFctPt)(const double &time, const Vector< double > &x, double &u)
Function pointer to source function fct(t,x,f(x,t)) – x is a Vector!
UnsteadyHeatSourceFctPt Source_fct_pt
Pointer to source function:
void output(std::ostream &outfile, const unsigned &nplot)
Output FE representation of soln: x,y,u or x,y,z,u at and history values at n_plot^DIM plot points...
UnsteadyHeatEquations()
Constructor: Initialises the Source_fct_pt to null and sets flag to use ALE formulation of the equati...
QUnsteadyHeatElement()
Constructor: Call constructors for QElement and UnsteadyHeat equations.
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
double * Beta_pt
Pointer to Beta parameter (thermal conductivity)
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
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
UnsteadyHeatSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
unsigned self_test()
Self-test: Return 0 for OK.
static double Default_beta_parameter
Static default value for the Beta parameter (thermal conductivity): One for natural scaling...
virtual double dshape_and_dtest_eulerian_ust_heat(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...
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
double interpolated_u_ust_heat(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
static char t char * s
Definition: cfortran.h:572
void output(std::ostream &outfile)
Output with default number of plot points.
double * Alpha_pt
Pointer to Alpha parameter (thermal inertia)
const double & beta() const
Beta parameter (thermal conductivity)
QUnsteadyHeatElement(const QUnsteadyHeatElement< DIM, NNODE_1D > &dummy)
Broken copy constructor.
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!) ...
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
UnsteadyHeat upgraded to become projectable.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute element residual Vector and element Jacobian matrix (wrapper)
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
static double Default_alpha_parameter
Static default value for the Alpha parameter: (thermal inertia): One for natural scaling.
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:1383
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
UnsteadyHeatSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
virtual UnsteadyHeatSourceFctPt & source_fct_pt()=0
Access function: Pointer to source function.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
static const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
double interpolated_u_ust_heat(const unsigned &t, const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s at previous time t (t=0: presen...
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 double dshape_and_dtest_eulerian_at_knot_ust_heat(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 ...
void output(FILE *file_pt)
C_style output with default number of plot points.
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
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 void fill_in_generic_residual_contribution_ust_heat(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
void(* UnsteadyHeatSourceFctPt)(const double &time, const Vector< double > &x, double &u)
Function pointer to source function fct(t,x,f(x,t)) – x is a Vector!
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
ProjectableUnsteadyHeatElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
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...
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.
UnsteadyHeatEquations(const UnsteadyHeatEquations &dummy)
Broken copy constructor.
double interpolated_du_dt_ust_heat(const Vector< double > &s) const
Return FE representation of function value du/dt(s) at local coordinate s.
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
double dshape_and_dtest_eulerian_at_knot_ust_heat(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: flux[i] = du/dx_i.