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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Non-inline functions for UnsteadyHeat elements
31 #include "unsteady_heat_elements.h"
32 
33 
34 namespace oomph
35 {
36 
37 ///2D UnsteadyHeat elements
38 
39 
40 /// Default value for Alpha parameter (thermal inertia)
41 template<unsigned DIM>
42 double UnsteadyHeatEquations<DIM>::Default_alpha_parameter=1.0;
43 
44 /// Default value for Beta parameter (thermal conductivity)
45 template<unsigned DIM>
46 double UnsteadyHeatEquations<DIM>::Default_beta_parameter=1.0;
47 
48 
49 //======================================================================
50 // Set the data for the number of Variables at each node
51 //======================================================================
52 template<unsigned DIM, unsigned NNODE_1D>
53 const unsigned QUnsteadyHeatElement<DIM,NNODE_1D>::Initial_Nvalue = 1;
54 
55 //======================================================================
56 /// Compute element residual Vector and/or element Jacobian matrix
57 ///
58 /// flag=1: compute both
59 /// flag=0: compute only residual Vector
60 ///
61 /// Pure version without hanging nodes
62 //======================================================================
63 template <unsigned DIM>
66  DenseMatrix<double> &jacobian,
67  unsigned flag)
68 {
69  //Find out how many nodes there are
70  unsigned n_node = nnode();
71 
72  // Get continuous time from timestepper of first node
73  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
74 
75  //Find the index at which the variable is stored
76  unsigned u_nodal_index = u_index_ust_heat();
77 
78  //Set up memory for the shape and test functions
79  Shape psi(n_node), test(n_node);
80  DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
81 
82  //Set the value of n_intpt
83  unsigned n_intpt = integral_pt()->nweight();
84 
85  //Set the Vector to hold local coordinates
86  Vector<double> s(DIM);
87 
88  //Get Alpha and beta parameters number
89  double alpha_local = alpha();
90  double beta_local = beta();
91 
92  //Integers to hold the local equation and unknowns
93  int local_eqn=0, local_unknown=0;
94 
95  //Loop over the integration points
96  for(unsigned ipt=0;ipt<n_intpt;ipt++)
97  {
98 
99  //Assign values of s
100  for(unsigned i=0;i<DIM;i++) s[i] = integral_pt()->knot(ipt,i);
101 
102  //Get the integral weight
103  double w = integral_pt()->weight(ipt);
104 
105  //Call the derivatives of the shape and test functions
106  double J =
107  dshape_and_dtest_eulerian_at_knot_ust_heat(ipt,psi,dpsidx,test,dtestdx);
108 
109  //Premultiply the weights and the Jacobian
110  double W = w*J;
111 
112  //Allocate memory for local quantities and initialise to zero
113  double interpolated_u=0.0;
114  double dudt=0.0;
115  Vector<double> interpolated_x(DIM,0.0);
116  Vector<double> interpolated_dudx(DIM,0.0);
117  Vector<double> mesh_velocity(DIM,0.0);
118 
119  //Calculate function value and derivatives:
120  // Loop over nodes
121  for(unsigned l=0;l<n_node;l++)
122  {
123  //Calculate the value at the nodes
124  double u_value = raw_nodal_value(l,u_nodal_index);
125  interpolated_u += u_value*psi(l);
126  dudt += du_dt_ust_heat(l)*psi(l);
127  // Loop over directions
128  for(unsigned j=0;j<DIM;j++)
129  {
130  interpolated_x[j] += raw_nodal_position(l,j)*psi(l);
131  interpolated_dudx[j] += u_value*dpsidx(l,j);
132  }
133  }
134 
135  // Mesh velocity?
136  if (!ALE_is_disabled)
137  {
138  for(unsigned l=0;l<n_node;l++)
139  {
140  for(unsigned j=0;j<DIM;j++)
141  {
142  mesh_velocity[j] += raw_dnodal_position_dt(l,j)*psi(l);
143  }
144  }
145  }
146 
147  //Get source function
148  //-------------------
149  double source=0.0;
150  get_source_ust_heat(time,ipt,interpolated_x,source);
151 
152  // Assemble residuals and Jacobian
153  //--------------------------------
154 
155  // Loop over the test functions
156  for(unsigned l=0;l<n_node;l++)
157  {
158  local_eqn = nodal_local_eqn(l,u_nodal_index);
159  /*IF it's not a boundary condition*/
160  if(local_eqn >= 0)
161  {
162 
163  // Add body force/source term and time derivative
164  residuals[local_eqn] += (alpha_local*dudt+source)*test(l)*W;
165 
166  // The mesh velocity bit
167  if (!ALE_is_disabled)
168  {
169  for(unsigned k=0;k<DIM;k++)
170  {
171  residuals[local_eqn] -= alpha_local*
172  mesh_velocity[k]*interpolated_dudx[k]*test(l)*W;
173  }
174  }
175 
176  // Laplace operator
177  for(unsigned k=0;k<DIM;k++)
178  {
179  residuals[local_eqn] += beta_local*
180  interpolated_dudx[k]*dtestdx(l,k)*W;
181  }
182 
183 
184  // Calculate the jacobian
185  //-----------------------
186  if(flag)
187  {
188  //Loop over the velocity shape functions again
189  for(unsigned l2=0;l2<n_node;l2++)
190  {
191  local_unknown = nodal_local_eqn(l2,u_nodal_index);
192  //If at a non-zero degree of freedom add in the entry
193  if(local_unknown >= 0)
194  {
195  // Mass matrix
196  jacobian(local_eqn,local_unknown)
197  += alpha_local*test(l)*psi(l2)*
198  node_pt(l2)->time_stepper_pt()->weight(1,0)*W;
199 
200  // Laplace operator & mesh velocity bit
201  for(unsigned i=0;i<DIM;i++)
202  {
203  double tmp=beta_local*dtestdx(l,i);
204  if (!ALE_is_disabled) tmp-=alpha_local*mesh_velocity[i]*test(l);
205  jacobian(local_eqn,local_unknown) +=
206  dpsidx(l2,i)*tmp*W;
207  }
208  }
209  }
210  }
211  }
212  }
213 
214 
215  } // End of loop over integration points
216 }
217 
218 
219 
220 //======================================================================
221  /// Compute norm of fe solution
222 //======================================================================
223 template <unsigned DIM>
225 {
226 
227  // Initialise
228  norm=0.0;
229 
230  //Vector of local coordinates
231  Vector<double> s(DIM);
232 
233  // Vector for coordintes
234  Vector<double> x(DIM);
235 
236  //Find out how many nodes there are in the element
237  unsigned n_node = nnode();
238 
239  Shape psi(n_node);
240 
241  //Set the value of n_intpt
242  unsigned n_intpt = integral_pt()->nweight();
243 
244  //Loop over the integration points
245  for(unsigned ipt=0;ipt<n_intpt;ipt++)
246  {
247 
248  //Assign values of s
249  for(unsigned i=0;i<DIM;i++)
250  {
251  s[i] = integral_pt()->knot(ipt,i);
252  }
253 
254  //Get the integral weight
255  double w = integral_pt()->weight(ipt);
256 
257  // Get jacobian of mapping
258  double J=J_eulerian(s);
259 
260  //Premultiply the weights and the Jacobian
261  double W = w*J;
262 
263  // Get FE function value
264  double u=interpolated_u_ust_heat(s);
265 
266  // Add to norm
267  norm+=u*u*W;
268 
269  }
270  }
271 
272 //======================================================================
273 /// Self-test: Return 0 for OK
274 //======================================================================
275 template <unsigned DIM>
277 {
278 
279  bool passed=true;
280 
281  // Check lower-level stuff
282  if (FiniteElement::self_test()!=0)
283  {
284  passed=false;
285  }
286 
287  // Return verdict
288  if (passed)
289  {
290  return 0;
291  }
292  else
293  {
294  return 1;
295  }
296 
297 }
298 
299 
300 
301 //======================================================================
302 /// Output function:
303 ///
304 /// x,y,u or x,y,z,u
305 ///
306 /// nplot points in each coordinate direction
307 //======================================================================
308 template <unsigned DIM>
309 void UnsteadyHeatEquations<DIM>::output(std::ostream &outfile,
310  const unsigned &nplot)
311 {
312  //Vector of local coordinates
313  Vector<double> s(DIM);
314 
315  // Tecplot header info
316  outfile << tecplot_zone_string(nplot);
317 
318  // Loop over plot points
319  unsigned num_plot_points=nplot_points(nplot);
320  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
321  {
322  // Get local coordinates of plot point
323  get_s_plot(iplot,nplot,s);
324 
325  for(unsigned i=0;i<DIM;i++)
326  {
327  outfile << interpolated_x(s,i) << " ";
328  }
329  outfile << interpolated_u_ust_heat(s) << std::endl;
330  }
331 
332  // Write tecplot footer (e.g. FE connectivity lists)
333  write_tecplot_zone_footer(outfile,nplot);
334 
335 }
336 
337 
338 
339 //======================================================================
340 /// C-style output function:
341 ///
342 /// x,y,u or x,y,z,u
343 ///
344 /// nplot points in each coordinate direction
345 //======================================================================
346 template <unsigned DIM>
348  const unsigned &nplot)
349 {
350  //Vector of local coordinates
351  Vector<double> s(DIM);
352 
353  // Tecplot header info
354  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
355 
356  // Loop over plot points
357  unsigned num_plot_points=nplot_points(nplot);
358  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
359  {
360 
361  // Get local coordinates of plot point
362  get_s_plot(iplot,nplot,s);
363 
364  for(unsigned i=0;i<DIM;i++)
365  {
366  fprintf(file_pt,"%g ",interpolated_x(s,i));
367 
368  }
369  fprintf(file_pt,"%g \n",interpolated_u_ust_heat(s));
370  }
371 
372  // Write tecplot footer (e.g. FE connectivity lists)
373  write_tecplot_zone_footer(file_pt,nplot);
374 
375 }
376 
377 
378 
379 //======================================================================
380  /// Output exact solution
381  ///
382  /// Solution is provided via function pointer.
383  /// Plot at a given number of plot points.
384  ///
385  /// x,y,u_exact or x,y,z,u_exact
386 //======================================================================
387 template <unsigned DIM>
388 void UnsteadyHeatEquations<DIM>::output_fct(std::ostream &outfile,
389  const unsigned &nplot,
391 {
392 
393  //Vector of local coordinates
394  Vector<double> s(DIM);
395 
396  // Vector for coordintes
397  Vector<double> x(DIM);
398 
399  // Tecplot header info
400  outfile << tecplot_zone_string(nplot);
401 
402  // Exact solution Vector (here a scalar)
403  Vector<double> exact_soln(1);
404 
405  // Loop over plot points
406  unsigned num_plot_points=nplot_points(nplot);
407  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
408  {
409 
410  // Get local coordinates of plot point
411  get_s_plot(iplot,nplot,s);
412 
413  // Get x position as Vector
414  interpolated_x(s,x);
415 
416  // Get exact solution at this point
417  (*exact_soln_pt)(x,exact_soln);
418 
419  //Output x,y,...,u_exact
420  for(unsigned i=0;i<DIM;i++)
421  {
422  outfile << x[i] << " ";
423  }
424  outfile << exact_soln[0] << std::endl;
425 
426  }
427 
428  // Write tecplot footer (e.g. FE connectivity lists)
429  write_tecplot_zone_footer(outfile,nplot);
430 
431 }
432 
433 
434 
435 //======================================================================
436 /// Output exact solution at time t
437 ///
438 /// Solution is provided via function pointer.
439 /// Plot at a given number of plot points.
440 ///
441 /// x,y,u_exact or x,y,z,u_exact
442 //======================================================================
443 template <unsigned DIM>
444 void UnsteadyHeatEquations<DIM>::output_fct(std::ostream &outfile,
445  const unsigned &nplot,
446  const double& time,
448 
449 {
450 
451  //Vector of local coordinates
452  Vector<double> s(DIM);
453 
454  // Vector for coordintes
455  Vector<double> x(DIM);
456 
457 
458  // Tecplot header info
459  outfile << tecplot_zone_string(nplot);
460 
461  // Exact solution Vector (here a scalar)
462  Vector<double> exact_soln(1);
463 
464  // Loop over plot points
465  unsigned num_plot_points=nplot_points(nplot);
466  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
467  {
468 
469  // Get local coordinates of plot point
470  get_s_plot(iplot,nplot,s);
471 
472  // Get x position as Vector
473  interpolated_x(s,x);
474 
475  // Get exact solution at this point
476  (*exact_soln_pt)(time,x,exact_soln);
477 
478  //Output x,y,...,u_exact
479  for(unsigned i=0;i<DIM;i++)
480  {
481  outfile << x[i] << " ";
482  }
483  outfile << exact_soln[0] << std::endl;
484 
485  }
486 
487  // Write tecplot footer (e.g. FE connectivity lists)
488  write_tecplot_zone_footer(outfile,nplot);
489 
490 }
491 
492 
493 
494 //======================================================================
495  /// Validate against exact solution
496  ///
497  /// Solution is provided via function pointer.
498  /// Plot error at a given number of plot points.
499  ///
500 //======================================================================
501 template <unsigned DIM>
502 void UnsteadyHeatEquations<DIM>::compute_error(std::ostream &outfile,
504  double& error, double& norm)
505 {
506 
507  // Initialise
508  error=0.0;
509  norm=0.0;
510 
511  //Vector of local coordinates
512  Vector<double> s(DIM);
513 
514  // Vector for coordintes
515  Vector<double> x(DIM);
516 
517  //Find out how many nodes there are in the element
518  unsigned n_node = nnode();
519 
520  Shape psi(n_node);
521 
522  //Set the value of n_intpt
523  unsigned n_intpt = integral_pt()->nweight();
524 
525  // Tecplot header info
526  outfile << "ZONE" << std::endl;
527 
528  // Exact solution Vector (here a scalar)
529  Vector<double> exact_soln(1);
530 
531  //Loop over the integration points
532  for(unsigned ipt=0;ipt<n_intpt;ipt++)
533  {
534 
535  //Assign values of s
536  for(unsigned i=0;i<DIM;i++)
537  {
538  s[i] = integral_pt()->knot(ipt,i);
539  }
540 
541  //Get the integral weight
542  double w = integral_pt()->weight(ipt);
543 
544  // Get jacobian of mapping
545  double J=J_eulerian(s);
546 
547  //Premultiply the weights and the Jacobian
548  double W = w*J;
549 
550  // Get x position as Vector
551  interpolated_x(s,x);
552 
553  // Get FE function value
554  double u_fe = interpolated_u_ust_heat(s);
555 
556  // Get exact solution at this point
557  (*exact_soln_pt)(x,exact_soln);
558 
559  //Output x,y,...,error
560  for(unsigned i=0;i<DIM;i++)
561  {
562  outfile << x[i] << " ";
563  }
564  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
565 
566  // Add to error and norm
567  norm+=exact_soln[0]*exact_soln[0]*W;
568  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
569  }
570 }
571 
572 
573 
574 
575 //======================================================================
576 /// Validate against exact solution at time t.
577 ///
578 /// Solution is provided via function pointer.
579 /// Plot error at a given number of plot points.
580 ///
581 //======================================================================
582 template<unsigned DIM>
583 void UnsteadyHeatEquations<DIM>::compute_error(std::ostream &outfile,
585  const double& time, double& error, double& norm)
586 
587 {
588 
589  // Initialise
590  error=0.0;
591  norm=0.0;
592 
593  //Vector of local coordinates
594  Vector<double> s(DIM);
595 
596  // Vector for coordintes
597  Vector<double> x(DIM);
598 
599 
600  //Find out how many nodes there are in the element
601  unsigned n_node = nnode();
602 
603  Shape psi(n_node);
604 
605  //Set the value of n_intpt
606  unsigned n_intpt = integral_pt()->nweight();
607 
608  // Tecplot
609  outfile << "ZONE" << std::endl;
610 
611  // Exact solution Vector (here a scalar)
612  Vector<double> exact_soln(1);
613 
614  //Loop over the integration points
615  for(unsigned ipt=0;ipt<n_intpt;ipt++)
616  {
617 
618  //Assign values of s
619  for(unsigned i=0;i<DIM;i++)
620  {
621  s[i] = integral_pt()->knot(ipt,i);
622  }
623 
624  //Get the integral weight
625  double w = integral_pt()->weight(ipt);
626 
627  // Get jacobian of mapping
628  double J=J_eulerian(s);
629 
630  //Premultiply the weights and the Jacobian
631  double W = w*J;
632 
633  // Get x position as Vector
634  interpolated_x(s,x);
635 
636  // Get FE function value
637  double u_fe = interpolated_u_ust_heat(s);
638 
639  // Get exact solution at this point
640  (*exact_soln_pt)(time,x,exact_soln);
641 
642  //Output x,y,...,error
643  for(unsigned i=0;i<DIM;i++)
644  {
645  outfile << x[i] << " ";
646  }
647  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
648 
649  // Add to error and norm
650  norm+=exact_soln[0]*exact_soln[0]*W;
651  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
652 
653  }
654 }
655 
656 
657 //====================================================================
658 // Force build of templates
659 //====================================================================
660 template class QUnsteadyHeatElement<1,2>;
661 template class QUnsteadyHeatElement<1,3>;
662 template class QUnsteadyHeatElement<1,4>;
663 
664 template class QUnsteadyHeatElement<2,2>;
665 template class QUnsteadyHeatElement<2,3>;
666 template class QUnsteadyHeatElement<2,4>;
667 
668 template class QUnsteadyHeatElement<3,2>;
669 template class QUnsteadyHeatElement<3,3>;
670 template class QUnsteadyHeatElement<3,4>;
671 
672 
673 }
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points.
void compute_norm(double &norm)
Compute norm of fe solution.
cstr elem_len * i
Definition: cfortran.h:607
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
Definition: elements.cc:4247
unsigned self_test()
Self-test: Return 0 for OK.
static char t char * s
Definition: cfortran.h:572
void output(std::ostream &outfile)
Output with default number of plot points.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
virtual 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.