discontinuous_galerkin_space_time_unsteady_heat_mixed_order_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 SpaceTimeUnsteadyHeatMixedOrder elements
27#ifndef OOMPH_DISCONTINUOUS_GALERKIN_SPACE_TIME_UNSTEADY_HEAT_MIXED_ORDER_ELEMENTS_HEADER
28#define OOMPH_DISCONTINUOUS_GALERKIN_SPACE_TIME_UNSTEADY_HEAT_MIXED_ORDER_ELEMENTS_HEADER
29
30// Config header generated by autoconfig
31#ifdef HAVE_CONFIG_H
32#include <oomph-lib-config.h>
33#endif
34
35// Oomph-lib headers
36#include "generic/Qelements.h"
37#include "generic/shape.h"
38#include "generic/projection.h"
39
40/// //////////////////////////////////////////////////////////////////////
41/// //////////////////////////////////////////////////////////////////////
42/// //////////////////////////////////////////////////////////////////////
43
44namespace oomph
45{
46 //============================================================================
47 /// Base class so that we don't need to know the dimension just to
48 /// set the source function!
49 //============================================================================
50 class SpaceTimeUnsteadyHeatEquationsBase : public virtual FiniteElement
51 {
52 public:
53 /// Function pointer to source function fct(t,x,f(x,t)) -- x
54 /// is a Vector!
55 typedef void (*SpaceTimeUnsteadyHeatSourceFctPt)(const double& time,
56 const Vector<double>& x,
57 double& u);
58
59 /// Access function: Pointer to source function
61 };
62
63 //============================================================================
64 /// A class for all isoparametric elements that solve the
65 /// SpaceTimeUnsteadyHeatMixedOrder equations.
66 /// \f[ \frac{\partial^2 u}{\partial x_i^2}=\frac{\partial u}{\partial t}+f(t,x_j) \f]
67 /// This contains the generic maths. Shape functions, geometric
68 /// mapping etc. must get implemented in derived class.
69 /// Note that this class assumes an isoparametric formulation, i.e. that
70 /// the scalar unknown is interpolated using the same shape funcitons
71 /// as the position.
72 //============================================================================
73 template<unsigned SPATIAL_DIM>
76 {
77 public:
78 /// Constructor: Initialises the Source_fct_pt to null and sets
79 /// flag to use ALE formulation of the equations. Also, set Alpha (thermal
80 /// inertia) and Beta (thermal conductivity) parameters to defaults (both
81 /// one for natural scaling).
84 {
85 // Set Alpha parameter to default (one for natural scaling)
87
88 // Set Beta parameter to default (one for natural scaling)
90 } // End of SpaceTimeUnsteadyHeatMixedOrderEquations
91
92
93 /// Broken copy constructor
95 const SpaceTimeUnsteadyHeatMixedOrderEquations& dummy) = delete;
96
97 /// Disable ALE, i.e. assert the mesh is not moving -- you do this
98 /// at your own risk!
100 {
101 // Set the flag to true
102 ALE_is_disabled = true;
103 } // End of disable_ALE
104
105
106 /// (Re-)enable ALE, i.e. take possible mesh motion into account
107 /// when evaluating the time-derivative. Note: By default, ALE is
108 /// enabled, at the expense of possibly creating unnecessary work
109 /// in problems where the mesh is, in fact, stationary.
111 {
112 // Set the flag to false
113 ALE_is_disabled = false;
114 } // End of enable_ALE
115
116
117 /// Compute norm of FE solution
118 void compute_norm(double& norm);
119
120
121 /// Output with default number of plot points
122 void output(std::ostream& outfile)
123 {
124 // Number of plot points
125 unsigned nplot = 5;
126
127 // Output the solution
128 output(outfile, nplot);
129 } // End of output
130
131
132 /// Output FE representation of soln: x,y,u or x,y,z,u at
133 /// nplot^SPATIAL_DIM plot points
134 void output(std::ostream& outfile, const unsigned& nplot);
135
136
137 /// C_style output with default number of plot points
138 void output(FILE* file_pt)
139 {
140 // Number of plot points
141 unsigned nplot = 5;
142
143 // Output the solution
144 output(file_pt, nplot);
145 } // End of output
146
147
148 /// C-style output FE representation of soln: x,y,u or x,y,z,u at
149 /// nplot^SPATIAL_DIM plot points
150 void output(FILE* file_pt, const unsigned& nplot);
151
152
153 /// Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^SPATIAL_DIM
154 /// plot points
155 void output_fct(std::ostream& outfile,
156 const unsigned& nplot,
158
159
160 /// Output exact soln: x,y,u_exact or x,y,z,u_exact at
161 /// nplot^SPATIAL_DIM plot points (time-dependent version)
162 virtual void output_fct(
163 std::ostream& outfile,
164 const unsigned& nplot,
165 const double& time,
167
168
169 /// Get error and norm against exact solution
170 void compute_error(std::ostream& outfile,
172 double& error,
173 double& norm);
174
175
176 /// Get error and norm against exact solution
177 void compute_error(std::ostream& outfile,
179 const double& time,
180 double& error,
181 double& norm);
182
183
184 /// C-style output FE representation of soln: x,y,u or x,y,z,u at
185 /// nplot^SPATIAL_DIM plot points
186 void output_element_paraview(std::ofstream& outfile, const unsigned& nplot);
187
188
189 /// Number of scalars/fields output by this element. Reimplements
190 /// broken virtual function in base class.
191 unsigned nscalar_paraview() const
192 {
193 // Only one field to output
194 return 1;
195 } // End of nscalar_paraview
196
197
198 /// Write values of the i-th scalar field at the plot points. Needs
199 /// to be implemented for each new specific element type.
200 void scalar_value_paraview(std::ofstream& file_out,
201 const unsigned& i,
202 const unsigned& nplot) const
203 {
204#ifdef PARANOID
205 if (i != 0)
206 {
207 std::stringstream error_stream;
208 error_stream << "Space-time unsteady heat elements only store a single "
209 << "field so i must be 0 rather than " << i << std::endl;
210 throw OomphLibError(
211 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
212 }
213#endif
214
215 // Get the number of plot points
216 unsigned local_loop = this->nplot_points_paraview(nplot);
217
218 // Loop over the plot points
219 for (unsigned j = 0; j < local_loop; j++)
220 {
221 // Storage for the local coordinates
222 Vector<double> s(SPATIAL_DIM + 1);
223
224 // Get the local coordinate of the required plot point
225 this->get_s_plot(j, nplot, s);
226
227 // Output the interpolated solution value
228 file_out << this->interpolated_u_ust_heat(s) << std::endl;
229 }
230 } // End of scalar_value_paraview
231
232
233 /// Write values of the i-th scalar field at the plot points. Needs
234 /// to be implemented for each new specific element type.
236 std::ofstream& file_out,
237 const unsigned& i,
238 const unsigned& nplot,
239 FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
240 {
241#ifdef PARANOID
242 if (i != 0)
243 {
244 std::stringstream error_stream;
245 error_stream << "Space-time unsteady heat elements only store a single "
246 << "field so i must be 0 rather than " << i << std::endl;
247 throw OomphLibError(
248 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
249 }
250#endif
251
252 // Get the number of plot points
253 unsigned local_loop = this->nplot_points_paraview(nplot);
254
255 // Loop over the plot points
256 for (unsigned j = 0; j < local_loop; j++)
257 {
258 // Storage for the local coordinates
259 Vector<double> s(SPATIAL_DIM + 1);
260
261 // Storage for the global coordinates
262 Vector<double> spatial_coordinates(SPATIAL_DIM);
263
264 // Get the local coordinate of the required plot point
265 this->get_s_plot(j, nplot, s);
266
267 // Loop over the spatial coordinates
268 for (unsigned i = 0; i < SPATIAL_DIM; i++)
269 {
270 // Assign the i-th spatial coordinate
271 spatial_coordinates[i] = interpolated_x(s, i);
272 }
273
274 // Exact solution vector (here it's simply a scalar)
275 Vector<double> exact_soln(1, 0.0);
276
277 // Get the exact solution at this point
278 (*exact_soln_pt)(spatial_coordinates, exact_soln);
279
280 // Output the interpolated solution value
281 file_out << exact_soln[0] << std::endl;
282 } // for (unsigned j=0;j<local_loop;j++)
283 } // End of scalar_value_fct_paraview
284
285
286 /// Write values of the i-th scalar field at the plot points. Needs
287 /// to be implemented for each new specific element type.
289 std::ofstream& file_out,
290 const unsigned& i,
291 const unsigned& nplot,
292 const double& time,
294 {
295#ifdef PARANOID
296 if (i != 0)
297 {
298 std::stringstream error_stream;
299 error_stream << "Space-time unsteady heat elements only store a single "
300 << "field so i must be 0 rather than " << i << std::endl;
301 throw OomphLibError(
302 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
303 }
304#endif
305
306 // Get the number of plot points
307 unsigned local_loop = this->nplot_points_paraview(nplot);
308
309 // Loop over the plot points
310 for (unsigned j = 0; j < local_loop; j++)
311 {
312 // Storage for the local coordinates
313 Vector<double> s(SPATIAL_DIM + 1);
314
315 // Storage for the time value
316 double interpolated_t = 0.0;
317
318 // Storage for the global coordinates
319 Vector<double> spatial_coordinates(SPATIAL_DIM);
320
321 // Get the local coordinate of the required plot point
322 this->get_s_plot(j, nplot, s);
323
324 // Loop over the spatial coordinates
325 for (unsigned i = 0; i < SPATIAL_DIM; i++)
326 {
327 // Assign the i-th spatial coordinate
328 spatial_coordinates[i] = interpolated_x(s, i);
329 }
330
331 // Get the time value
332 interpolated_t = interpolated_x(s, SPATIAL_DIM);
333
334 // Exact solution vector (here it's simply a scalar)
335 Vector<double> exact_soln(1, 0.0);
336
337 // Get the exact solution at this point
338 (*exact_soln_pt)(interpolated_t, spatial_coordinates, exact_soln);
339
340 // Output the interpolated solution value
341 file_out << exact_soln[0] << std::endl;
342 } // for (unsigned j=0;j<local_loop;j++)
343 } // End of scalar_value_fct_paraview
344
345
346 /// Name of the i-th scalar field. Default implementation
347 /// returns V1 for the first one, V2 for the second etc.
348 std::string scalar_name_paraview(const unsigned& i) const
349 {
350 // If we're outputting the solution
351 if (i == 0)
352 {
353 // There's only one field to output
354 return "U";
355 }
356 // Never get here
357 else
358 {
359 std::stringstream error_stream;
360 error_stream << "These unsteady heat elements only store 1 field, \n"
361 << "but i is currently " << i << std::endl;
362 throw OomphLibError(
363 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
364
365 // Dummy return
366 return " ";
367 }
368 } // End of scalar_name_paraview
369
370
371 /// Access function: Pointer to source function
373 {
374 // Return the source function pointer
375 return Source_fct_pt;
376 } // End of source_fct_pt
377
378
379 /// Access function: Pointer to source function. Const version
381 {
382 // Return the source function pointer
383 return Source_fct_pt;
384 }
385
386
387 /// Get source term at continous time t and (Eulerian) position x.
388 /// Virtual so it can be overloaded in derived multi-physics elements.
389 virtual inline void get_source_ust_heat(const double& t,
390 const unsigned& ipt,
391 const Vector<double>& x,
392 double& source) const
393 {
394 // If no source function has been set, return zero
395 if (Source_fct_pt == 0)
396 {
397 // Set the source term value to zero
398 source = 0.0;
399 }
400 // Otherwise return the appropriate value
401 else
402 {
403 // Get source strength
404 (*Source_fct_pt)(t, x, source);
405 }
406 } // End of get_source_ust_heat
407
408
409 /// Alpha parameter (thermal inertia)
410 const double& alpha() const
411 {
412 // Return the value of Alpha
413 return *Alpha_pt;
414 } // End of alpha
415
416
417 /// Pointer to Alpha parameter (thermal inertia)
418 double*& alpha_pt()
419 {
420 // Return the pointer to Alpha
421 return Alpha_pt;
422 } // End of alpha_pt
423
424
425 /// Beta parameter (thermal conductivity)
426 const double& beta() const
427 {
428 // Return the pointer to Beta
429 return *Beta_pt;
430 } // End of beta
431
432
433 /// Pointer to Beta parameter (thermal conductivity)
434 double*& beta_pt()
435 {
436 // Return the pointer to Beta
437 return Beta_pt;
438 } // End of beta_pt
439
440
441 /// Get flux: flux[i]=du/dx_i
442 void get_flux(const Vector<double>& s, Vector<double>& flux) const
443 {
444 // Find out how many nodes there are in the element
445 unsigned n_node = nnode();
446
447 // Find the index at which the variable is stored
448 unsigned u_nodal_index = u_index_ust_heat();
449
450 // Set up memory for the shape and test functions
451 Shape psi(n_node);
452
453 // Set up memory for the derivatives of the shape functions
454 DShape dpsidx(n_node, SPATIAL_DIM + 1);
455
456 //------------dshape_eulerian(s,psi,dpsidx)----------------------------
457 // Find the element dimension
458 const unsigned el_dim = this->dim();
459
460 // Get the values of the shape functions and their local derivatives;
461 // temporarily stored in dpsi
462 dshape_local_ust_heat(s, psi, dpsidx);
463
464 // Allocate memory for the inverse jacobian
465 DenseMatrix<double> inverse_jacobian(el_dim);
466
467 // Now calculate the inverse jacobian
468 local_to_eulerian_mapping(dpsidx, inverse_jacobian);
469
470 // Now set the values of the derivatives to be dpsidx
471 transform_derivatives(inverse_jacobian, dpsidx);
472 //------------dshape_eulerian(s,psi,dpsidx)----------------------------
473
474 // Loop over the entries of the flux vector
475 for (unsigned j = 0; j < SPATIAL_DIM; j++)
476 {
477 // Initialise j-th flux entry to zero
478 flux[j] = 0.0;
479 }
480
481 // Loop over nodes
482 for (unsigned l = 0; l < n_node; l++)
483 {
484 // Loop over derivative directions
485 for (unsigned j = 0; j < SPATIAL_DIM; j++)
486 {
487 // Update the flux value
488 flux[j] += nodal_value(l, u_nodal_index) * dpsidx(l, j);
489 }
490 } // for (unsigned l=0;l<n_node;l++)
491 } // End of get_flux
492
493
494 /// Compute element residual Vector (wrapper)
496 {
497 // Call the generic residuals function with flag set to 0
498 // using a dummy matrix argument
501 } // End of fill_in_contribution_to_residuals
502
503
504 /// Compute element residual Vector and element Jacobian matrix (wrapper)
506 DenseMatrix<double>& jacobian)
507 {
508 // Call the generic routine with the flag set to 1
510 } // End of fill_in_contribution_to_jacobian
511
512
513 /// Return FE representation of function value u(s) at local coordinate s
514 inline double interpolated_u_ust_heat(const Vector<double>& s) const
515 {
516 // Find number of nodes
517 unsigned n_node = nnode();
518
519 // Find the index at which the variable is stored
520 unsigned u_nodal_index = u_index_ust_heat();
521
522 // Local shape function
523 Shape psi(n_node);
524
525 // Find values of the shape functions at local coordinate s
526 shape_ust_heat(s, psi);
527
528 // Initialise value of u
529 double interpolated_u = 0.0;
530
531 // Loop over the local nodes and sum
532 for (unsigned l = 0; l < n_node; l++)
533 {
534 // Update the interpolated u value
535 interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
536 }
537
538 // Return the interpolated u value
539 return interpolated_u;
540 } // End of interpolated_u_ust_heat
541
542
543 /// Return the index at which the unknown value
544 /// is stored. The default value, 0, is appropriate for single-physics
545 /// problems, when there is only one variable, the value that satisfies the
546 /// unsteady heat equation.
547 /// In derived multi-physics elements, this function should be overloaded
548 /// to reflect the chosen storage scheme. Note that these equations require
549 /// that the unknown is always stored at the same index at each node.
550 virtual inline unsigned u_index_ust_heat() const
551 {
552 // Return the default value
553 return 0;
554 } // End of u_index_ust_heat
555
556
557 /// Return FE representation of function value u(s) at local
558 /// coordinate s at previous time t (t=0: present)
559 /// DRAIG: This needs to be broken; doesn't make sense in space-time
560 /// elements!
561 inline double interpolated_u_ust_heat(const unsigned& t,
562 const Vector<double>& s) const
563 {
564 // Find number of nodes
565 unsigned n_node = nnode();
566
567 // Find the index at which the variable is stored
568 unsigned u_nodal_index = u_index_ust_heat();
569
570 // Local shape function
571 Shape psi(n_node);
572
573 // Find values of shape function
574 shape_ust_heat(s, psi);
575
576 // Initialise value of u
577 double interpolated_u = 0.0;
578
579 // Loop over the local nodes and sum
580 for (unsigned l = 0; l < n_node; l++)
581 {
582 // Update the interpolated u value
583 interpolated_u += nodal_value(t, l, u_nodal_index) * psi[l];
584 }
585
586 // Return the interpolated u value
587 return interpolated_u;
588 } // End of interpolated_u_ust_heat
589
590
591 /// Calculate du/dt at the n-th local node. Uses suitably
592 /// interpolated value for hanging nodes.
593 double du_dt_ust_heat(const unsigned& n) const
594 {
595 // Storage for the local coordinates
596 Vector<double> s(SPATIAL_DIM + 1, 0.0);
597
598 // Get the local coordinate at the n-th node
600
601 // Return the interpolated du/dt value
603 } // End of du_dt_ust_heat
604
605
606 /// Return FE representation of function value du/dt(s) at local coordinate
607 /// s
608 inline double interpolated_du_dt_ust_heat(const Vector<double>& s) const
609 {
610 // Find number of nodes
611 unsigned n_node = nnode();
612
613 // Find the index at which the variable is stored
614 unsigned u_nodal_index = u_index_ust_heat();
615
616 // Local shape function
617 Shape psi(n_node);
618
619 // Allocate space for the derivatives of the shape functions
620 DShape dpsidx(n_node, SPATIAL_DIM + 1);
621
622 //------------dshape_eulerian(s,psi,dpsidx)----------------------------
623 // Find the element dimension
624 const unsigned el_dim = this->dim();
625
626 // Get the values of the shape functions and their local derivatives;
627 // temporarily stored in dpsi
628 dshape_local_ust_heat(s, psi, dpsidx);
629
630 // Allocate memory for the inverse jacobian
631 DenseMatrix<double> inverse_jacobian(el_dim);
632
633 // Now calculate the inverse jacobian
634 local_to_eulerian_mapping(dpsidx, inverse_jacobian);
635
636 // Now set the values of the derivatives to be dpsidx
637 transform_derivatives(inverse_jacobian, dpsidx);
638 //------------dshape_eulerian(s,psi,dpsidx)----------------------------
639
640 // Initialise value of du/dt
641 double interpolated_dudt = 0.0;
642
643 // Loop over the local nodes and sum
644 for (unsigned l = 0; l < n_node; l++)
645 {
646 // Update the interpolated du/dt value
647 interpolated_dudt +=
648 nodal_value(l, u_nodal_index) * dpsidx(l, SPATIAL_DIM);
649 }
650
651 // Return the interpolated du/dt value
652 return interpolated_dudt;
653 } // End of interpolated_du_dt_ust_heat
654
655
656 /// Self-test: Return 0 for OK
657 unsigned self_test();
658
659 /// Shape/test functions and derivs w.r.t. to global coords at
660 /// local coordinate s; return Jacobian of mapping
662 const Vector<double>& s,
663 Shape& psi,
664 DShape& dpsidx,
665 Shape& test,
666 DShape& dtestdx) const = 0;
667
668
669 /// Shape/test functions and derivs w.r.t. to global coords at
670 /// integration point ipt; return Jacobian of mapping
672 const unsigned& ipt,
673 Shape& psi,
674 DShape& dpsidx,
675 Shape& test,
676 DShape& dtestdx) const = 0;
677
678 protected:
679 /// Shape functions w.r.t. to local coords
680 virtual void shape_ust_heat(const Vector<double>& s, Shape& psi) const = 0;
681
682
683 /// Shape functions & derivs. w.r.t. to local coords
685 Shape& psi,
686 DShape& dpsidx) const = 0;
687
688
689 /// Test functions & derivs. w.r.t. to local coords
691 Shape& test,
692 DShape& dtestdx) const = 0;
693
694
695 /// Compute element residual Vector only (if flag=and/or element
696 /// Jacobian matrix
698 Vector<double>& residuals,
699 DenseMatrix<double>& jacobian,
700 const unsigned& flag);
701
702 /// Pointer to source function:
704
705 /// Boolean flag to indicate if ALE formulation is disabled when
706 /// time-derivatives are computed. Only set to true if you're sure that
707 /// the mesh is stationary.
709
710 /// Pointer to Alpha parameter (thermal inertia)
711 double* Alpha_pt;
712
713 /// Pointer to Beta parameter (thermal conductivity)
714 double* Beta_pt;
715
716 private:
717 /// Static default value for the Alpha parameter (thermal inertia):
718 /// One for natural scaling
720
721 /// Static default value for the Beta parameter (thermal
722 /// conductivity): One for natural scaling
724 };
725
726
727 /// ////////////////////////////////////////////////////////////////////////
728 /// ////////////////////////////////////////////////////////////////////////
729 /// ////////////////////////////////////////////////////////////////////////
730
731
732 //=========================================================================
733 /// QUnsteadyHeatMixedOrderSpaceTimeElement elements are quadrilateral/brick-
734 /// shaped UnsteadyHeatMixedOrder elements with isoparametric interpolation
735 /// for the function.
736 //=========================================================================
737 template<unsigned SPATIAL_DIM, unsigned NNODE_1D>
739 : public virtual QElement<SPATIAL_DIM + 1, NNODE_1D>,
740 public virtual SpaceTimeUnsteadyHeatMixedOrderEquations<SPATIAL_DIM>
741 {
742 public:
743 /// Constructor: Call constructors for QElement and
744 /// SpaceTimeUnsteadyHeatMixedOrder equations
746 : QElement<SPATIAL_DIM + 1, NNODE_1D>(),
748 {
749 }
750
751 /// Broken copy constructor
754 dummy) = delete;
755
756 /// Required number of 'values' (pinned or dofs) at node n
757 inline unsigned required_nvalue(const unsigned& n) const
758 {
759 // Return the appropriate value
760 return Initial_Nvalue;
761 } // End of required_nvalue
762
763
764 /// Output function:
765 /// x,t,u or x,y,t,u
766 void output(std::ostream& outfile)
767 {
768 // Call the function in the base class
770 } // End of output
771
772
773 /// Output function:
774 /// x,t,u or x,y,t,u at n_plot^(SPATIAL_DIM+1) plot points
775 void output(std::ostream& outfile, const unsigned& n_plot)
776 {
777 // Call the function in the base class
779 n_plot);
780 } // End of output
781
782
783 /// C-style output function:
784 /// x,t,u or x,y,t,u
785 void output(FILE* file_pt)
786 {
787 // Call the function in the base class
789 } // End of output
790
791
792 /// C-style output function:
793 /// x,t,u or x,y,t,u at n_plot^(SPATIAL_DIM+1) plot points
794 void output(FILE* file_pt, const unsigned& n_plot)
795 {
796 // Call the function in the base class
798 n_plot);
799 } // End of output
800
801
802 /// Output function for an exact solution:
803 /// x,t,u_exact or x,y,t,u_exact at n_plot^(SPATIAL_DIM+1) plot points
804 void output_fct(std::ostream& outfile,
805 const unsigned& n_plot,
807 {
808 // Call the function in the base class
810 outfile, n_plot, exact_soln_pt);
811 } // End of output_fct
812
813
814 /// Output function for a time-dependent exact solution.
815 /// x,t,u_exact or x,y,t,u_exact at n_plot^(SPATIAL_DIM+1) plot points
816 /// (Calls the unsteady version)
817 void output_fct(std::ostream& outfile,
818 const unsigned& n_plot,
819 const double& time,
821 {
822 // Call the function in the base class
824 outfile, n_plot, time, exact_soln_pt);
825 } // End of output_fct
826
827
828 /// Shape/test functions & derivs. w.r.t. to global coords. Return Jacobian.
830 Shape& psi,
831 DShape& dpsidx,
832 Shape& test,
833 DShape& dtestdx) const;
834
835 /// Shape/test functions and derivs w.r.t. to global coords at
836 /// integration point ipt; return Jacobian of mapping
838 const unsigned& ipt,
839 Shape& psi,
840 DShape& dpsidx,
841 Shape& test,
842 DShape& dtestdx) const;
843
844 protected:
845 /// Shape functions w.r.t. to local coords
846 inline void shape_ust_heat(const Vector<double>& s, Shape& psi) const;
847
848
849 /// Shape functions & derivs. w.r.t. to local coords
850 inline void dshape_local_ust_heat(const Vector<double>& s,
851 Shape& psi,
852 DShape& dpsidx) const;
853
854
855 /// Test functions & derivs. w.r.t. to local coords
856 inline void dtest_local_ust_heat(const Vector<double>& s,
857 Shape& test,
858 DShape& dtestdx) const;
859
860
861 private:
862 /// Static array of ints to hold number of variables at nodes:
863 /// Initial_Nvalue[n]
864 static const unsigned Initial_Nvalue;
865 };
866
867
868 //======================================================================
869 /// Define the shape functions and test functions and derivatives
870 /// w.r.t. local coordinates
871 ///
872 /// Galerkin: Test functions=shape functions
873 //======================================================================
874 template<unsigned SPATIAL_DIM, unsigned NNODE_1D>
876 shape_ust_heat(const Vector<double>& s, Shape& psi) const
877 {
878 // Local storage (for 3D = 2D space + 1D time)
879 double psi_values[3][NNODE_1D];
880
881 // Index of the total shape function
882 unsigned index = 0;
883
884 // Call the 1D shape functions and derivatives
885 OneDimLagrange::shape<NNODE_1D>(s[0], psi_values[0]);
886 OneDimLagrange::shape<NNODE_1D>(s[1], psi_values[1]);
887
888 // Set the time discretisation
889 OneDimDiscontinuousGalerkinMixedOrderBasis::shape<NNODE_1D>(s[2],
890 psi_values[2]);
891
892 // Loop over the nodes in the third local coordinate direction
893 for (unsigned k = 0; k < NNODE_1D; k++)
894 {
895 // Loop over the nodes in the second local coordinate direction
896 for (unsigned j = 0; j < NNODE_1D; j++)
897 {
898 // Loop over the nodes in the first local coordinate direction
899 for (unsigned i = 0; i < NNODE_1D; i++)
900 {
901 // Calculate the index-th entry of psi
902 psi[index] = psi_values[0][i] * psi_values[1][j] * psi_values[2][k];
903
904 // Increment the index
905 index++;
906 }
907 } // for (unsigned j=0;j<NNODE_1D;j++)
908 } // for (unsigned k=0;k<NNODE_1D;k++)
909 } // End of shape_ust_heat
910
911
912 //======================================================================
913 /// Define the shape functions and test functions and derivatives
914 /// w.r.t. local coordinates
915 ///
916 /// Galerkin: Test functions=shape functions
917 //======================================================================
918 template<unsigned SPATIAL_DIM, unsigned NNODE_1D>
921 Shape& psi,
922 DShape& dpsidx) const
923 {
924 // Local storage (for 3D = 2D space + 1D time)
925 double psi_values[3][NNODE_1D];
926 double dpsi_values[3][NNODE_1D];
927
928 // Index of the total shape function
929 unsigned index = 0;
930
931 // Call the 1D shape functions and derivatives
932 OneDimLagrange::shape<NNODE_1D>(s[0], psi_values[0]);
933 OneDimLagrange::shape<NNODE_1D>(s[1], psi_values[1]);
934 OneDimLagrange::dshape<NNODE_1D>(s[0], dpsi_values[0]);
935 OneDimLagrange::dshape<NNODE_1D>(s[1], dpsi_values[1]);
936
937 // Set the time discretisation
938 OneDimDiscontinuousGalerkinMixedOrderBasis::shape<NNODE_1D>(s[2],
939 psi_values[2]);
940 OneDimDiscontinuousGalerkinMixedOrderBasis::dshape<NNODE_1D>(
941 s[2], dpsi_values[2]);
942
943 // Loop over the nodes in the third local coordinate direction
944 for (unsigned k = 0; k < NNODE_1D; k++)
945 {
946 // Loop over the nodes in the second local coordinate direction
947 for (unsigned j = 0; j < NNODE_1D; j++)
948 {
949 // Loop over the nodes in the first local coordinate direction
950 for (unsigned i = 0; i < NNODE_1D; i++)
951 {
952 // Calculate dpsi/ds_0
953 dpsidx(index, 0) =
954 dpsi_values[0][i] * psi_values[1][j] * psi_values[2][k];
955
956 // Calculate dpsi/ds_1
957 dpsidx(index, 1) =
958 psi_values[0][i] * dpsi_values[1][j] * psi_values[2][k];
959
960 // Calculate dpsi/ds_2
961 dpsidx(index, 2) =
962 psi_values[0][i] * psi_values[1][j] * dpsi_values[2][k];
963
964 // Calculate the index-th entry of psi
965 psi[index] = psi_values[0][i] * psi_values[1][j] * psi_values[2][k];
966
967 // Increment the index
968 index++;
969 }
970 } // for (unsigned j=0;j<NNODE_1D;j++)
971 } // for (unsigned k=0;k<NNODE_1D;k++)
972 } // End of dshape_local_ust_heat
973
974
975 //======================================================================
976 /// Define the test functions and derivatives w.r.t. local coordinates
977 //======================================================================
978 template<unsigned SPATIAL_DIM, unsigned NNODE_1D>
981 Shape& test,
982 DShape& dtestdx) const
983 {
984 // Local storage (for 3D = 2D space + 1D time)
985 double test_values[3][NNODE_1D];
986 double dtest_values[3][NNODE_1D];
987
988 // Index of the total shape function
989 unsigned index = 0;
990
991 // Call the 1D shape functions and derivatives
992 OneDimLagrange::shape<NNODE_1D>(s[0], test_values[0]);
993 OneDimLagrange::shape<NNODE_1D>(s[1], test_values[1]);
994 OneDimLagrange::dshape<NNODE_1D>(s[0], dtest_values[0]);
995 OneDimLagrange::dshape<NNODE_1D>(s[1], dtest_values[1]);
996
997 // Set the time discretisation
998 OneDimDiscontinuousGalerkinMixedOrderTest::shape<NNODE_1D>(s[2],
999 test_values[2]);
1000 OneDimDiscontinuousGalerkinMixedOrderTest::dshape<NNODE_1D>(
1001 s[2], dtest_values[2]);
1002
1003 // Loop over the nodes in the third local coordinate direction
1004 for (unsigned k = 0; k < NNODE_1D; k++)
1005 {
1006 // Loop over the nodes in the second local coordinate direction
1007 for (unsigned j = 0; j < NNODE_1D; j++)
1008 {
1009 // Loop over the nodes in the first local coordinate direction
1010 for (unsigned i = 0; i < NNODE_1D; i++)
1011 {
1012 // Calculate dtest/ds_0
1013 dtestdx(index, 0) =
1014 dtest_values[0][i] * test_values[1][j] * test_values[2][k];
1015
1016 // Calculate dtest/ds_1
1017 dtestdx(index, 1) =
1018 test_values[0][i] * dtest_values[1][j] * test_values[2][k];
1019
1020 // Calculate dtest/ds_2
1021 dtestdx(index, 2) =
1022 test_values[0][i] * test_values[1][j] * dtest_values[2][k];
1023
1024 // Calculate the index-th entry of test
1025 test[index] =
1026 test_values[0][i] * test_values[1][j] * test_values[2][k];
1027
1028 // Increment the index
1029 index++;
1030 }
1031 } // for (unsigned j=0;j<NNODE_1D;j++)
1032 } // for (unsigned k=0;k<NNODE_1D;k++)
1033 } // End of dtest_local_ust_heat
1034
1035
1036 //======================================================================
1037 /// Define the shape functions and test functions and derivatives
1038 /// w.r.t. global coordinates and return Jacobian of mapping.
1039 ///
1040 /// Galerkin: Test functions=shape functions
1041 //======================================================================
1042 template<unsigned SPATIAL_DIM, unsigned NNODE_1D>
1045 Shape& psi,
1046 DShape& dpsidx,
1047 Shape& test,
1048 DShape& dtestdx) const
1049 {
1050 //--------------------------
1051 // Call the shape functions:
1052 //--------------------------
1053 // Find the element dimension
1054 const unsigned el_dim = this->dim();
1055
1056 // Make sure we're using 3D space-time elements
1057 if (el_dim != 3)
1058 {
1059 // Create an output stream
1060 std::ostringstream error_message_stream;
1061
1062 // Create an error message
1063 error_message_stream << "Need 3D space-time elements for this to work!"
1064 << std::endl;
1065
1066 // Throw the error message
1067 throw OomphLibError(error_message_stream.str(),
1068 OOMPH_CURRENT_FUNCTION,
1069 OOMPH_EXCEPTION_LOCATION);
1070 }
1071
1072 // Compute the geometric shape functions and also first derivatives
1073 // w.r.t. local coordinates at local coordinate s
1074 dshape_local_ust_heat(s, psi, dpsidx);
1075
1076 // Allocate memory for the inverse jacobian
1077 DenseMatrix<double> inverse_jacobian(el_dim);
1078
1079 // Now calculate the inverse jacobian
1080 const double det =
1081 this->local_to_eulerian_mapping(dpsidx, inverse_jacobian);
1082
1083 // Now set the values of the derivatives to be dpsidx
1084 this->transform_derivatives(inverse_jacobian, dpsidx);
1085
1086 //-------------------------
1087 // Call the test functions:
1088 //-------------------------
1089 // Compute the geometric test functions and also first derivatives
1090 // w.r.t. local coordinates at local coordinate s
1091 dtest_local_ust_heat(s, test, dtestdx);
1092
1093 // Transform derivatives from dtest/ds to dtest/dx
1094 this->transform_derivatives(inverse_jacobian, dtestdx);
1095
1096 // Return the determinant value
1097 return det;
1098 } // End of dshape_and_dtest_eulerian_ust_heat
1099
1100
1101 //======================================================================
1102 /// Define the shape functions and test functions and derivatives
1103 /// w.r.t. global coordinates and return Jacobian of mapping.
1104 ///
1105 /// Galerkin: Test functions=shape functions
1106 //======================================================================
1107 template<unsigned SPATIAL_DIM, unsigned NNODE_1D>
1110 Shape& psi,
1111 DShape& dpsidx,
1112 Shape& test,
1113 DShape& dtestdx) const
1114 {
1115 // Find the element dimension
1116 const unsigned el_dim = SPATIAL_DIM + 1;
1117
1118 // Storage for the local coordinates of the integration point
1119 Vector<double> s(el_dim, 0.0);
1120
1121 // Set the local coordinate
1122 for (unsigned i = 0; i < el_dim; i++)
1123 {
1124 // Calculate the i-th local coordinate at the ipt-th knot point
1125 s[i] = this->integral_pt()->knot(ipt, i);
1126 }
1127
1128 // Return the Jacobian of the geometrical shape functions and derivatives
1129 return dshape_and_dtest_eulerian_ust_heat(s, psi, dpsidx, test, dtestdx);
1130 } // End of dshape_and_dtest_eulerian_at_knot_ust_heat
1131
1132
1133 /// /////////////////////////////////////////////////////////////////////
1134 /// /////////////////////////////////////////////////////////////////////
1135 /// /////////////////////////////////////////////////////////////////////
1136
1137
1138 //=======================================================================
1139 /// Face geometry for the QUnsteadyHeatMixedOrderSpaceTimeElement elements:
1140 /// The spatial dimension of the face elements is one lower than that of the
1141 /// bulk element but they have the same number of points along their 1D edges.
1142 //=======================================================================
1143 template<unsigned SPATIAL_DIM, unsigned NNODE_1D>
1145 QUnsteadyHeatMixedOrderSpaceTimeElement<SPATIAL_DIM, NNODE_1D>>
1146 : public virtual QElement<SPATIAL_DIM, NNODE_1D>
1147 {
1148 public:
1149 /// Constructor: Call the constructor for the appropriate
1150 /// lower-dimensional QElement
1151 FaceGeometry() : QElement<SPATIAL_DIM, NNODE_1D>() {}
1152 };
1153
1154 /// /////////////////////////////////////////////////////////////////////
1155 /// /////////////////////////////////////////////////////////////////////
1156 /// /////////////////////////////////////////////////////////////////////
1157
1158 //=======================================================================
1159 /// Face geometry for the 1D QUnsteadyHeatMixedOrderSpaceTimeElement
1160 /// elements: Point elements
1161 //=======================================================================
1162 template<unsigned NNODE_1D>
1164 : public virtual PointElement
1165 {
1166 public:
1167 /// Constructor: Call the constructor for the appropriate
1168 /// lower-dimensional QElement
1170 };
1171
1172
1173 /// /////////////////////////////////////////////////////////////////////
1174 /// /////////////////////////////////////////////////////////////////////
1175 /// /////////////////////////////////////////////////////////////////////
1176
1177
1178 //==========================================================
1179 /// SpaceTimeUnsteadyHeatMixedOrder upgraded to become projectable
1180 //==========================================================
1181 template<class UNSTEADY_HEAT_ELEMENT>
1183 : public virtual ProjectableElement<UNSTEADY_HEAT_ELEMENT>
1184 {
1185 public:
1186 /// Constructor [this was only required explicitly
1187 /// from gcc 4.5.2 onwards...]
1189
1190
1191 /// Specify the values associated with field fld. The information
1192 /// is returned in a vector of pairs which comprise the Data object and
1193 /// the value within it, that correspond to field fld.
1195 {
1196#ifdef PARANOID
1197 // If we're not dealing with the first field
1198 if (fld != 0)
1199 {
1200 // Create a stringstream object to create an error message
1201 std::stringstream error_stream;
1202
1203 // Create the error string
1204 error_stream
1205 << "SpaceTimeUnsteadyHeatMixedOrder elements only store a single "
1206 << "field so fld must be 0 rather than " << fld << std::endl;
1207
1208 // Throw an error
1209 throw OomphLibError(
1210 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1211 }
1212#endif
1213
1214 // The number of nodes in this element
1215 unsigned nnod = this->nnode();
1216
1217 // Storage for the pairs
1218 Vector<std::pair<Data*, unsigned>> data_values(nnod);
1219
1220 // Loop over all nodes
1221 for (unsigned j = 0; j < nnod; j++)
1222 {
1223 // Add the data value and associated field: The node itself
1224 data_values[j] = std::make_pair(this->node_pt(j), fld);
1225 }
1226
1227 // Return the vector
1228 return data_values;
1229 } // End of data_values_of_field
1230
1231
1232 /// Number of fields to be projected: Just one
1234 {
1235 // Return the appropriate value
1236 return 1;
1237 } // End of nfields_for_projection
1238
1239
1240 /// Number of history values to be stored for fld-th field.
1241 /// (Note: count includes current value!)
1242 unsigned nhistory_values_for_projection(const unsigned& fld)
1243 {
1244#ifdef PARANOID
1245 // If we're not dealing with the first field
1246 if (fld != 0)
1247 {
1248 // Create a stringstream object to create an error message
1249 std::stringstream error_stream;
1250
1251 // Create the error string
1252 error_stream
1253 << "SpaceTimeUnsteadyHeatMixedOrder elements only store a single "
1254 << "field so fld must be 0 rather than " << fld << std::endl;
1255
1256 // Throw an error
1257 throw OomphLibError(
1258 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1259 }
1260#endif
1261
1262 // Return the number of stored values
1263 return this->node_pt(0)->ntstorage();
1264 } // End of nhistory_values_for_projection
1265
1266
1267 /// Number of positional history values (Note: count includes
1268 /// current value!)
1270 {
1271 // Return the number of history values stored by the position timestepper
1272 return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
1273 } // End of nhistory_values_for_coordinate_projection
1274
1275
1276 /// Return Jacobian of mapping and shape functions of field fld
1277 /// at local coordinate s
1278 double jacobian_and_shape_of_field(const unsigned& fld,
1279 const Vector<double>& s,
1280 Shape& psi)
1281 {
1282#ifdef PARANOID
1283 // If we're not dealing with the first field
1284 if (fld != 0)
1285 {
1286 // Create a stringstream object to create an error message
1287 std::stringstream error_stream;
1288
1289 // Create the error string
1290 error_stream
1291 << "SpaceTimeUnsteadyHeatMixedOrder elements only store a single "
1292 << "field so fld must be 0 rather than " << fld << std::endl;
1293
1294 // Throw an error
1295 throw OomphLibError(
1296 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1297 }
1298#endif
1299
1300 // Get the number of dimensions in the element
1301 unsigned n_dim = this->dim();
1302
1303 // Get the number of nodes in the element
1304 unsigned n_node = this->nnode();
1305
1306 // Allocate space for the test functions
1307 Shape test(n_node);
1308
1309 // Allocate space for the derivatives of the shape functions
1310 DShape dpsidx(n_node, n_dim);
1311
1312 // Allocate space for the derivatives of the test functions
1313 DShape dtestdx(n_node, n_dim);
1314
1315 // Calculate the shape functions and their derivatives at the local
1316 // coordinate s (and the same for the test functions). On top of this
1317 // calculate the determinant of the Jacobian
1318 double J =
1319 this->dshape_and_dtest_eulerian_ust_heat(s, psi, dpsidx, test, dtestdx);
1320
1321 // Return the determinant of the Jacobian
1322 return J;
1323 } // End of jacobian_and_shape_of_field
1324
1325
1326 /// Return interpolated field fld at local coordinate s, at time
1327 /// level t (t=0: present; t>0: history values)
1328 double get_field(const unsigned& t,
1329 const unsigned& fld,
1330 const Vector<double>& s)
1331 {
1332#ifdef PARANOID
1333 // If we're not dealing with the first field
1334 if (fld != 0)
1335 {
1336 // Create a stringstream object to create an error message
1337 std::stringstream error_stream;
1338
1339 // Create the error string
1340 error_stream
1341 << "SpaceTimeUnsteadyHeatMixedOrder elements only store a single "
1342 << "field so fld must be 0 rather than " << fld << std::endl;
1343
1344 // Throw an error
1345 throw OomphLibError(
1346 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1347 }
1348#endif
1349
1350 // Find the index at which the variable is stored
1351 unsigned u_nodal_index = this->u_index_ust_heat();
1352
1353 // Get the number of nodes in the element
1354 unsigned n_node = this->nnode();
1355
1356 // Local shape function
1357 Shape psi(n_node);
1358
1359 // Find values of shape function
1360 this->shape(s, psi);
1361
1362 // Initialise value of u
1363 double interpolated_u = 0.0;
1364
1365 // Loop over the local nodes
1366 for (unsigned l = 0; l < n_node; l++)
1367 {
1368 // Update the interpolated solution value
1369 interpolated_u += this->nodal_value(t, l, u_nodal_index) * psi[l];
1370 }
1371
1372 // Return the interpolated solution value
1373 return interpolated_u;
1374 } // End of get_field
1375
1376
1377 /// Return number of values in field fld: One per node
1378 unsigned nvalue_of_field(const unsigned& fld)
1379 {
1380#ifdef PARANOID
1381 // If we're not dealing with the first field
1382 if (fld != 0)
1383 {
1384 // Create a stringstream object to create an error message
1385 std::stringstream error_stream;
1386
1387 // Create the error string
1388 error_stream
1389 << "SpaceTimeUnsteadyHeatMixedOrder elements only store a single "
1390 << "field so fld must be 0 rather than " << fld << std::endl;
1391
1392 // Throw an error
1393 throw OomphLibError(
1394 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1395 }
1396#endif
1397
1398 // Return the number of nodes in the element
1399 return this->nnode();
1400 } // End of nvalue_of_field
1401
1402
1403 /// Return local equation number of value j in field fld.
1404 int local_equation(const unsigned& fld, const unsigned& j)
1405 {
1406#ifdef PARANOID
1407 // If we're not dealing with the first field
1408 if (fld != 0)
1409 {
1410 // Create a stringstream object to create an error message
1411 std::stringstream error_stream;
1412
1413 // Create the error string
1414 error_stream << "SpaceTimeUnsteadyHeatMixedOrder elements only store a "
1415 << "single field so fld must be 0 rather than " << fld
1416 << std::endl;
1417
1418 // Throw an error
1419 throw OomphLibError(
1420 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1421 }
1422#endif
1423
1424 // Get the nodal index of the unknown
1425 const unsigned u_nodal_index = this->u_index_ust_heat();
1426
1427 // Output the local equation number
1428 return this->nodal_local_eqn(j, u_nodal_index);
1429 } // End of local_equation
1430
1431
1432 /// Output FE representation of soln: x,t,u or x,y,t,u
1433 /// at n_plot^(SPATIAL_DIM+1) plot points
1434 void output(std::ostream& outfile, const unsigned& nplot)
1435 {
1436 // Get the dimension of the element
1437 unsigned el_dim = this->dim();
1438
1439 // Vector of local coordinates
1440 Vector<double> s(el_dim, 0.0);
1441
1442 // Tecplot header info
1443 outfile << this->tecplot_zone_string(nplot);
1444
1445 // Get the number of plot points
1446 unsigned num_plot_points = this->nplot_points(nplot);
1447
1448 // Loop over plot points
1449 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1450 {
1451 // Get local coordinates of plot point
1452 this->get_s_plot(iplot, nplot, s);
1453
1454 // Loop over the coordinate directions
1455 for (unsigned i = 0; i < el_dim; i++)
1456 {
1457 // Output the interpolated coordinates
1458 outfile << this->interpolated_x(s, i) << " ";
1459 }
1460
1461 // Output the interpolated value of u(s)
1462 outfile << this->interpolated_u_ust_heat(s) << " ";
1463
1464 // Output the interpolated value of du/dt(s)
1465 outfile << this->interpolated_du_dt_ust_heat(s) << " ";
1466
1467 // History values of coordinates
1468 unsigned n_prev =
1470
1471 // Loop over the previous timesteps
1472 for (unsigned t = 1; t < n_prev; t++)
1473 {
1474 // Loop over the coordinate directions
1475 for (unsigned i = 0; i < el_dim; i++)
1476 {
1477 // Output the coordinates
1478 outfile << this->interpolated_x(t, s, i) << " ";
1479 }
1480 } // for (unsigned t=1;t<n_prev;t++)
1481
1482 // Number of history values of velocities
1483 n_prev = this->node_pt(0)->time_stepper_pt()->ntstorage();
1484
1485 // Loop over the previous timesteps
1486 for (unsigned t = 1; t < n_prev; t++)
1487 {
1488 // Output the solution
1489 outfile << this->interpolated_u_ust_heat(t, s) << " ";
1490 }
1491
1492 // Finish the line
1493 outfile << std::endl;
1494 } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
1495
1496 // Write tecplot footer (e.g. FE connectivity lists)
1497 this->write_tecplot_zone_footer(outfile, nplot);
1498 } // End of output
1499 };
1500
1501
1502 //=======================================================================
1503 /// Face geometry for element is the same as that for the underlying
1504 /// wrapped element
1505 //=======================================================================
1506 template<class ELEMENT>
1508 : public virtual FaceGeometry<ELEMENT>
1509 {
1510 public:
1511 FaceGeometry() : FaceGeometry<ELEMENT>() {}
1512 };
1513
1514
1515 //=======================================================================
1516 /// Face geometry of the Face Geometry for element is the same as
1517 /// that for the underlying wrapped element
1518 //=======================================================================
1519 template<class ELEMENT>
1522 : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
1523 {
1524 public:
1526 };
1527} // End of namespace oomph
1528
1529#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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
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:2862
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual)
Definition: elements.h:1842
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 void transform_derivatives(const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
Convert derivative w.r.t.local coordinates to derivatives w.r.t the coordinates used to assemble the ...
Definition: elements.cc:2833
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 double local_to_eulerian_mapping(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates, given the derivatives of the shape function...
Definition: elements.h:1508
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
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 *& 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
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
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 nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
ProjectableUnsteadyHeatMixedOrderSpaceTimeElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
void output(std::ostream &outfile, const unsigned &nplot)
Output FE representation of soln: x,t,u or x,y,t,u at n_plot^(SPATIAL_DIM+1) plot points.
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!...
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 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 ...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
QUnsteadyHeatMixedOrderSpaceTimeElement()
Constructor: Call constructors for QElement and SpaceTimeUnsteadyHeatMixedOrder equations.
void dtest_local_ust_heat(const Vector< double > &s, Shape &test, DShape &dtestdx) const
Test functions & derivs. w.r.t. to local coords.
unsigned required_nvalue(const unsigned &n) const
Required number of 'values' (pinned or dofs) at node n.
void dshape_local_ust_heat(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Shape functions & derivs. w.r.t. to local coords.
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 ...
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,t,u or x,y,t,u at n_plot^(SPATIAL_DIM+1) plot points.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: x,t,u or x,y,t,u at n_plot^(SPATIAL_DIM+1) plot points.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: x,t,u_exact or x,y,t,u_exact at n_plot^(SPATIAL_DIM+1) plot po...
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,t,u_exact or x,y,t,u_exact at n_plot^(SPATIAL_...
QUnsteadyHeatMixedOrderSpaceTimeElement(const QUnsteadyHeatMixedOrderSpaceTimeElement< SPATIAL_DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
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.
void shape_ust_heat(const Vector< double > &s, Shape &psi) const
Shape functions w.r.t. to local coords.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
Base class so that we don't need to know the dimension just to set the source function!
void(* SpaceTimeUnsteadyHeatSourceFctPt)(const double &time, const Vector< double > &x, double &u)
Function pointer to source function fct(t,x,f(x,t)) – x is a Vector!
virtual SpaceTimeUnsteadyHeatSourceFctPt & source_fct_pt()=0
Access function: Pointer to source function.
A class for all isoparametric elements that solve the SpaceTimeUnsteadyHeatMixedOrder equations.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
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)
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...
SpaceTimeUnsteadyHeatSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
double interpolated_u_ust_heat(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class.
void scalar_value_fct_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
SpaceTimeUnsteadyHeatMixedOrderEquations(const SpaceTimeUnsteadyHeatMixedOrderEquations &dummy)=delete
Broken copy constructor.
void output_element_paraview(std::ofstream &outfile, const unsigned &nplot)
C-style output FE representation of soln: x,y,u or x,y,z,u at nplot^SPATIAL_DIM plot points.
static double Default_alpha_parameter
Static default value for the Alpha parameter (thermal inertia): One for natural scaling.
double interpolated_du_dt_ust_heat(const Vector< double > &s) const
Return FE representation of function value du/dt(s) at local coordinate s.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
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...
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 ...
static double Default_beta_parameter
Static default value for the Beta parameter (thermal conductivity): One for natural scaling.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute element residual Vector (wrapper)
SpaceTimeUnsteadyHeatMixedOrderEquations()
Constructor: Initialises the Source_fct_pt to null and sets flag to use ALE formulation of the equati...
virtual void dshape_local_ust_heat(const Vector< double > &s, Shape &psi, DShape &dpsidx) const =0
Shape functions & derivs. w.r.t. to local coords.
void output(std::ostream &outfile)
Output with default number of plot points.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one,...
virtual void shape_ust_heat(const Vector< double > &s, Shape &psi) const =0
Shape functions w.r.t. to local coords.
void scalar_value_fct_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
virtual void dtest_local_ust_heat(const Vector< double > &s, Shape &test, DShape &dtestdx) const =0
Test functions & derivs. w.r.t. to local coords.
double du_dt_ust_heat(const unsigned &n) const
Calculate du/dt at the n-th local node. Uses suitably interpolated value for hanging nodes.
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 coordinate s; return Jacobian of map...
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative....
void output(FILE *file_pt)
C_style output with default number of plot points.
virtual void fill_in_generic_residual_contribution_ust_heat(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
virtual unsigned u_index_ust_heat() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...
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...
SpaceTimeUnsteadyHeatSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error and norm against exact solution.
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^SPATIAL_DIM plot points.
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...