unsteady_heat_elements.cc
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-2025 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// Non-inline functions for UnsteadyHeat elements
28
29
30namespace oomph
31{
32 /// 2D UnsteadyHeat elements
33
34
35 /// Default value for Alpha parameter (thermal inertia)
36 template<unsigned DIM>
38
39 /// Default value for Beta parameter (thermal conductivity)
40 template<unsigned DIM>
42
43
44 //======================================================================
45 // Set the data for the number of Variables at each node
46 //======================================================================
47 template<unsigned DIM, unsigned NNODE_1D>
49
50 //======================================================================
51 /// Compute element residual Vector and/or element Jacobian matrix
52 ///
53 /// flag=1: compute both
54 /// flag=0: compute only residual Vector
55 ///
56 /// Pure version without hanging nodes
57 //======================================================================
58 template<unsigned DIM>
62 {
63 // Find out how many nodes there are
64 unsigned n_node = nnode();
65
66 // Get continuous time from timestepper of first node
67 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
68
69 // Find the index at which the variable is stored
70 unsigned u_nodal_index = u_index_ust_heat();
71
72 // Set up memory for the shape and test functions
75
76 // Set the value of n_intpt
77 unsigned n_intpt = integral_pt()->nweight();
78
79 // Set the Vector to hold local coordinates
81
82 // Get Alpha and beta parameters number
83 double alpha_local = alpha();
84 double beta_local = beta();
85
86 // Integers to hold the local equation and unknowns
87 int local_eqn = 0, local_unknown = 0;
88
89 // Loop over the integration points
90 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
91 {
92 // Assign values of s
93 for (unsigned i = 0; i < DIM; i++) s[i] = integral_pt()->knot(ipt, i);
94
95 // Get the integral weight
96 double w = integral_pt()->weight(ipt);
97
98 // Call the derivatives of the shape and test functions
99 double J = dshape_and_dtest_eulerian_at_knot_ust_heat(
101
102 // Premultiply the weights and the Jacobian
103 double W = w * J;
104
105 // Allocate memory for local quantities and initialise to zero
106 double interpolated_u = 0.0;
107 double dudt = 0.0;
109 Vector<double> interpolated_dudx(DIM, 0.0);
111
112 // Calculate function value and derivatives:
113 // Loop over nodes
114 for (unsigned l = 0; l < n_node; l++)
115 {
116 // Calculate the value at the nodes
118 interpolated_u += u_value * psi(l);
119 dudt += du_dt_ust_heat(l) * psi(l);
120 // Loop over directions
121 for (unsigned j = 0; j < DIM; j++)
122 {
124 interpolated_dudx[j] += u_value * dpsidx(l, j);
125 }
126 }
127
128 // Mesh velocity?
129 if (!ALE_is_disabled)
130 {
131 for (unsigned l = 0; l < n_node; l++)
132 {
133 for (unsigned j = 0; j < DIM; j++)
134 {
136 }
137 }
138 }
139
140 // Get source function
141 //-------------------
142 double source = 0.0;
143 get_source_ust_heat(time, ipt, interpolated_x, source);
144
145 // Assemble residuals and Jacobian
146 //--------------------------------
147
148 // Loop over the test functions
149 for (unsigned l = 0; l < n_node; l++)
150 {
152 /*IF it's not a boundary condition*/
153 if (local_eqn >= 0)
154 {
155 // Add body force/source term and time derivative
156 residuals[local_eqn] += (alpha_local * dudt + source) * test(l) * W;
157
158 // The mesh velocity bit
159 if (!ALE_is_disabled)
160 {
161 for (unsigned k = 0; k < DIM; k++)
162 {
164 interpolated_dudx[k] * test(l) * W;
165 }
166 }
167
168 // Laplace operator
169 for (unsigned k = 0; k < DIM; k++)
170 {
172 beta_local * interpolated_dudx[k] * dtestdx(l, k) * W;
173 }
174
175
176 // Calculate the jacobian
177 //-----------------------
178 if (flag)
179 {
180 // Loop over the velocity shape functions again
181 for (unsigned l2 = 0; l2 < n_node; l2++)
182 {
184 // If at a non-zero degree of freedom add in the entry
185 if (local_unknown >= 0)
186 {
187 // Mass matrix
188 jacobian(local_eqn, local_unknown) +=
189 alpha_local * test(l) * psi(l2) *
190 node_pt(l2)->time_stepper_pt()->weight(1, 0) * W;
191
192 // Laplace operator & mesh velocity bit
193 for (unsigned i = 0; i < DIM; i++)
194 {
195 double tmp = beta_local * dtestdx(l, i);
196 if (!ALE_is_disabled)
198 jacobian(local_eqn, local_unknown) += dpsidx(l2, i) * tmp * W;
199 }
200 }
201 }
202 }
203 }
204 }
205
206
207 } // End of loop over integration points
208 }
209
210
211 //======================================================================
212 /// Compute norm of fe solution
213 //======================================================================
214 template<unsigned DIM>
216 {
217 // Initialise
218 norm = 0.0;
219
220 // Vector of local coordinates
222
223 // Vector for coordintes
225
226 // Find out how many nodes there are in the element
227 unsigned n_node = nnode();
228
230
231 // Set the value of n_intpt
232 unsigned n_intpt = integral_pt()->nweight();
233
234 // Loop over the integration points
235 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
236 {
237 // Assign values of s
238 for (unsigned i = 0; i < DIM; i++)
239 {
240 s[i] = integral_pt()->knot(ipt, i);
241 }
242
243 // Get the integral weight
244 double w = integral_pt()->weight(ipt);
245
246 // Get jacobian of mapping
247 double J = J_eulerian(s);
248
249 // Premultiply the weights and the Jacobian
250 double W = w * J;
251
252 // Get FE function value
253 double u = interpolated_u_ust_heat(s);
254
255 // Add to norm
256 norm += u * u * W;
257 }
258 }
259
260 //======================================================================
261 /// Self-test: Return 0 for OK
262 //======================================================================
263 template<unsigned DIM>
265 {
266 bool passed = true;
267
268 // Check lower-level stuff
269 if (FiniteElement::self_test() != 0)
270 {
271 passed = false;
272 }
273
274 // Return verdict
275 if (passed)
276 {
277 return 0;
278 }
279 else
280 {
281 return 1;
282 }
283 }
284
285
286 //======================================================================
287 /// Output function:
288 ///
289 /// x,y,u or x,y,z,u
290 ///
291 /// nplot points in each coordinate direction
292 //======================================================================
293 template<unsigned DIM>
295 const unsigned& nplot)
296 {
297 // Vector of local coordinates
299
300 // Tecplot header info
302
303 // Loop over plot points
305 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
306 {
307 // Get local coordinates of plot point
309
310 for (unsigned i = 0; i < DIM; i++)
311 {
312 outfile << interpolated_x(s, i) << " ";
313 }
314 outfile << interpolated_u_ust_heat(s) << std::endl;
315 }
316
317 // Write tecplot footer (e.g. FE connectivity lists)
319 }
320
321
322 //======================================================================
323 /// C-style output function:
324 ///
325 /// x,y,u or x,y,z,u
326 ///
327 /// nplot points in each coordinate direction
328 //======================================================================
329 template<unsigned DIM>
331 {
332 // Vector of local coordinates
334
335 // Tecplot header info
337
338 // Loop over plot points
340 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
341 {
342 // Get local coordinates of plot point
344
345 for (unsigned i = 0; i < DIM; i++)
346 {
347 fprintf(file_pt, "%g ", interpolated_x(s, i));
348 }
349 fprintf(file_pt, "%g \n", interpolated_u_ust_heat(s));
350 }
351
352 // Write tecplot footer (e.g. FE connectivity lists)
354 }
355
356
357 //======================================================================
358 /// Output exact solution
359 ///
360 /// Solution is provided via function pointer.
361 /// Plot at a given number of plot points.
362 ///
363 /// x,y,u_exact or x,y,z,u_exact
364 //======================================================================
365 template<unsigned DIM>
367 std::ostream& outfile,
368 const unsigned& nplot,
370 {
371 // Vector of local coordinates
373
374 // Vector for coordintes
376
377 // Tecplot header info
379
380 // Exact solution Vector (here a scalar)
382
383 // Loop over plot points
385 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
386 {
387 // Get local coordinates of plot point
389
390 // Get x position as Vector
391 interpolated_x(s, x);
392
393 // Get exact solution at this point
394 (*exact_soln_pt)(x, exact_soln);
395
396 // Output x,y,...,u_exact
397 for (unsigned i = 0; i < DIM; i++)
398 {
399 outfile << x[i] << " ";
400 }
401 outfile << exact_soln[0] << std::endl;
402 }
403
404 // Write tecplot footer (e.g. FE connectivity lists)
406 }
407
408
409 //======================================================================
410 /// Output exact solution at time t
411 ///
412 /// Solution is provided via function pointer.
413 /// Plot at a given number of plot points.
414 ///
415 /// x,y,u_exact or x,y,z,u_exact
416 //======================================================================
417 template<unsigned DIM>
419 std::ostream& outfile,
420 const unsigned& nplot,
421 const double& time,
423
424 {
425 // Vector of local coordinates
427
428 // Vector for coordintes
430
431
432 // Tecplot header info
434
435 // Exact solution Vector (here a scalar)
437
438 // Loop over plot points
440 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
441 {
442 // Get local coordinates of plot point
444
445 // Get x position as Vector
446 interpolated_x(s, x);
447
448 // Get exact solution at this point
449 (*exact_soln_pt)(time, x, exact_soln);
450
451 // Output x,y,...,u_exact
452 for (unsigned i = 0; i < DIM; i++)
453 {
454 outfile << x[i] << " ";
455 }
456 outfile << exact_soln[0] << std::endl;
457 }
458
459 // Write tecplot footer (e.g. FE connectivity lists)
461 }
462
463
464 //======================================================================
465 /// Validate against exact solution
466 ///
467 /// Solution is provided via function pointer.
468 /// Plot error at a given number of plot points.
469 ///
470 //======================================================================
471 template<unsigned DIM>
473 std::ostream& outfile,
475 double& error,
476 double& norm)
477 {
478 // Initialise
479 error = 0.0;
480 norm = 0.0;
481
482 // Vector of local coordinates
484
485 // Vector for coordintes
487
488 // Find out how many nodes there are in the element
489 unsigned n_node = nnode();
490
492
493 // Set the value of n_intpt
494 unsigned n_intpt = integral_pt()->nweight();
495
496 // Tecplot header info
497 outfile << "ZONE" << std::endl;
498
499 // Exact solution Vector (here a scalar)
501
502 // Loop over the integration points
503 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
504 {
505 // Assign values of s
506 for (unsigned i = 0; i < DIM; i++)
507 {
508 s[i] = integral_pt()->knot(ipt, i);
509 }
510
511 // Get the integral weight
512 double w = integral_pt()->weight(ipt);
513
514 // Get jacobian of mapping
515 double J = J_eulerian(s);
516
517 // Premultiply the weights and the Jacobian
518 double W = w * J;
519
520 // Get x position as Vector
521 interpolated_x(s, x);
522
523 // Get FE function value
524 double u_fe = interpolated_u_ust_heat(s);
525
526 // Get exact solution at this point
527 (*exact_soln_pt)(x, exact_soln);
528
529 // Output x,y,...,error
530 for (unsigned i = 0; i < DIM; i++)
531 {
532 outfile << x[i] << " ";
533 }
534 outfile << exact_soln[0] << " " << exact_soln[0] - u_fe << std::endl;
535
536 // Add to error and norm
537 norm += exact_soln[0] * exact_soln[0] * W;
538 error += (exact_soln[0] - u_fe) * (exact_soln[0] - u_fe) * W;
539 }
540 }
541
542
543 //======================================================================
544 /// Validate against exact solution at time t.
545 ///
546 /// Solution is provided via function pointer.
547 /// Plot error at a given number of plot points.
548 ///
549 //======================================================================
550 template<unsigned DIM>
552 std::ostream& outfile,
554 const double& time,
555 double& error,
556 double& norm)
557
558 {
559 // Initialise
560 error = 0.0;
561 norm = 0.0;
562
563 // Vector of local coordinates
565
566 // Vector for coordintes
568
569
570 // Find out how many nodes there are in the element
571 unsigned n_node = nnode();
572
574
575 // Set the value of n_intpt
576 unsigned n_intpt = integral_pt()->nweight();
577
578 // Tecplot
579 outfile << "ZONE" << std::endl;
580
581 // Exact solution Vector (here a scalar)
583
584 // Loop over the integration points
585 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
586 {
587 // Assign values of s
588 for (unsigned i = 0; i < DIM; i++)
589 {
590 s[i] = integral_pt()->knot(ipt, i);
591 }
592
593 // Get the integral weight
594 double w = integral_pt()->weight(ipt);
595
596 // Get jacobian of mapping
597 double J = J_eulerian(s);
598
599 // Premultiply the weights and the Jacobian
600 double W = w * J;
601
602 // Get x position as Vector
603 interpolated_x(s, x);
604
605 // Get FE function value
606 double u_fe = interpolated_u_ust_heat(s);
607
608 // Get exact solution at this point
609 (*exact_soln_pt)(time, x, exact_soln);
610
611 // Output x,y,...,error
612 for (unsigned i = 0; i < DIM; i++)
613 {
614 outfile << x[i] << " ";
615 }
616 outfile << exact_soln[0] << " " << exact_soln[0] - u_fe << std::endl;
617
618 // Add to error and norm
619 norm += exact_soln[0] * exact_soln[0] * W;
620 error += (exact_soln[0] - u_fe) * (exact_soln[0] - u_fe) * W;
621 }
622 }
623
624
625 //====================================================================
626 // Force build of templates
627 //====================================================================
628 template class QUnsteadyHeatElement<1, 2>;
629 template class QUnsteadyHeatElement<1, 3>;
630 template class QUnsteadyHeatElement<1, 4>;
631
632 template class QUnsteadyHeatElement<2, 2>;
633 template class QUnsteadyHeatElement<2, 3>;
634 template class QUnsteadyHeatElement<2, 4>;
635
636 template class QUnsteadyHeatElement<3, 2>;
637 template class QUnsteadyHeatElement<3, 3>;
638 template class QUnsteadyHeatElement<3, 4>;
639
640
641} // namespace oomph
static char t char * s
Definition cfortran.h:568
cstr elem_len * i
Definition cfortran.h:603
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
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
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
Definition elements.cc:4133
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition elements.h:1967
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:3165
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:3992
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:1436
unsigned nnode() const
Return the number of nodes.
Definition elements.h:2214
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition elements.h:1763
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:3152
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:3190
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
double raw_dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n. Do not use the hanging node repre...
Definition elements.h:2260
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition elements.h:2580
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation....
Definition elements.cc:1714
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:3178
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition elements.h:1769
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK.
Definition elements.cc:4470
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
static const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition shape.h:76
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
TAdvectionDiffusionReactionElement()
Constructor: Call constructors for TElement and AdvectionDiffusionReaction equations.
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Time *const & time_pt() const
Access function for the pointer to time (const version)
double & time()
Return the current value of the continuous time.
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 output(std::ostream &outfile)
Output with default number of plot points.
unsigned self_test()
Self-test: Return 0 for OK.
static double Default_beta_parameter
Static default value for the Beta parameter (thermal conductivity): One for natural scaling.
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.
static double Default_alpha_parameter
Static default value for the Alpha parameter: (thermal inertia): One for natural scaling.
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).