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