space_time_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-2023 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 SpaceTimeUnsteadyHeat elements
28 
29 /// //////////////////////////////////////////////////////////////////////
30 /// //////////////////////////////////////////////////////////////////////
31 /// //////////////////////////////////////////////////////////////////////
32 
33 namespace oomph
34 {
35  //======================================================================
36  // Default parameters
37  //======================================================================
38  /// Default value for Alpha parameter (thermal inertia)
39  template<unsigned SPATIAL_DIM>
41  1.0;
42 
43  /// Default value for Beta parameter (thermal conductivity)
44  template<unsigned SPATIAL_DIM>
46  1.0;
47 
48  //======================================================================
49  // Set the data for the number of variables at each node
50  //======================================================================
51  template<unsigned SPATIAL_DIM, unsigned NNODE_1D>
52  const unsigned
54 
55  //======================================================================
56  /// Compute element residual vector and/or element Jacobian matrix
57  ///
58  /// flag=0: compute only residual vector
59  /// flag=1: compute both
60  ///
61  /// Pure version without hanging nodes
62  //======================================================================
63  template<unsigned SPATIAL_DIM>
66  Vector<double>& residuals,
67  DenseMatrix<double>& jacobian,
68  const unsigned& flag)
69  {
70  // Find out how many nodes there are
71  unsigned n_node = nnode();
72 
73  // Find the index at which the variable is stored
74  unsigned u_nodal_index = u_index_ust_heat();
75 
76  // Set up memory for the shape functions
77  Shape psi(n_node);
78 
79  // Set up memory for the test functions
80  Shape test(n_node);
81 
82  // Allocate space for the derivatives of the shape functions
83  DShape dpsidx(n_node, SPATIAL_DIM + 1);
84 
85  // Allocate space for the derivatives of the test functions
86  DShape dtestdx(n_node, SPATIAL_DIM + 1);
87 
88  // Set the value of n_intpt
89  unsigned n_intpt = integral_pt()->nweight();
90 
91  // Storage for the local coordinates
92  Vector<double> s(SPATIAL_DIM + 1, 0.0);
93 
94  // Get the Alpha parameter
95  double alpha_local = alpha();
96 
97  // Get the Beta parameter
98  double beta_local = beta();
99 
100  // Integer to hold the local equation
101  int local_eqn = 0;
102 
103  // Integer to hold the local unknowns
104  int local_unknown = 0;
105 
106  // Loop over the integration points
107  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
108  {
109  // Assign values of s
110  for (unsigned i = 0; i < SPATIAL_DIM + 1; i++)
111  {
112  // Calculate the i-th local coordinate
113  s[i] = integral_pt()->knot(ipt, i);
114  }
115 
116  // Get the integral weight
117  double w = integral_pt()->weight(ipt);
118 
119  // Call the derivatives of the shape and test functions
120  double J = dshape_and_dtest_eulerian_at_knot_ust_heat(
121  ipt, psi, dpsidx, test, dtestdx);
122 
123  // Premultiply the weights and the Jacobian
124  double W = w * J;
125 
126  // Storage for the interpolated time value
127  double interpolated_t = 0.0;
128 
129  // Storage for the interpolated solution value
130  double interpolated_u = 0.0;
131 
132  // Storage for the interpolated time-derivative of the solution
133  double interpolated_dudt = 0.0;
134 
135  // Storage for the spatial coordinates
136  Vector<double> interpolated_x(SPATIAL_DIM, 0.0);
137 
138  // Storage for the spatial derivatives of the solution
139  Vector<double> interpolated_dudx(SPATIAL_DIM, 0.0);
140 
141  // Storage for the mesh velocity
142  Vector<double> mesh_velocity(SPATIAL_DIM, 0.0);
143 
144  //-------------------------------------------------
145  // Calculate derivatives and source function value:
146  //-------------------------------------------------
147  // Loop over the nodes
148  for (unsigned l = 0; l < n_node; l++)
149  {
150  // Get the nodal value at the l-th node
151  double u_value = raw_nodal_value(l, u_nodal_index);
152 
153  // Update the interpolated time value
154  interpolated_t += raw_nodal_position(l, SPATIAL_DIM) * psi(l);
155 
156  // Loop over the coordinate directions (both spatial AND time)
157  for (unsigned j = 0; j < SPATIAL_DIM; j++)
158  {
159  // Update the interpolated x value
160  interpolated_x[j] += raw_nodal_position(l, j) * psi(l);
161 
162  // Update the interpolated du/dx_j value
163  interpolated_dudx[j] += u_value * dpsidx(l, j);
164  }
165 
166  // Update the interpolated u value
167  interpolated_u += u_value * psi(l);
168 
169  // Update the interpolated du/dt value
170  interpolated_dudt += u_value * dpsidx(l, SPATIAL_DIM);
171  } // for (unsigned l=0;l<n_node;l++)
172 
173  // Initialise the source term value
174  double source = 0.0;
175 
176  // Get the interpolated source term value
177  get_source_ust_heat(interpolated_t, ipt, interpolated_x, source);
178 
179  //---------------------------------
180  // Assemble residuals and Jacobian:
181  //---------------------------------
182  // Loop over the nodes (or equivalently the test functions)
183  for (unsigned l = 0; l < n_node; l++)
184  {
185  // Get the local equation number
186  local_eqn = nodal_local_eqn(l, u_nodal_index);
187 
188  // If it's not a boundary condition
189  if (local_eqn >= 0)
190  {
191  // Add source term and time derivative
192  residuals[local_eqn] +=
193  (source + alpha_local * interpolated_dudt) * test(l) * W;
194 
195  // Loop over the coordinate directions
196  for (unsigned k = 0; k < SPATIAL_DIM; k++)
197  {
198  // Add in the contribution from the Laplace operator
199  residuals[local_eqn] +=
200  beta_local * interpolated_dudx[k] * dtestdx(l, k) * W;
201  }
202 
203  //------------------------
204  // Calculate the Jacobian:
205  //------------------------
206  // If we also need to construct the Jacobian
207  if (flag)
208  {
209  // Loop over the velocity shape functions again
210  for (unsigned l2 = 0; l2 < n_node; l2++)
211  {
212  // Get the local equation number
213  local_unknown = nodal_local_eqn(l2, u_nodal_index);
214 
215  // If we're at a non-zero degree of freedom add in the entry
216  if (local_unknown >= 0)
217  {
218  // Add in the time derivative contribution
219  jacobian(local_eqn, local_unknown) +=
220  (alpha_local * test(l) * dpsidx(l2, SPATIAL_DIM) * W);
221 
222  // Laplace operator
223  for (unsigned i = 0; i < SPATIAL_DIM; i++)
224  {
225  // Add the test function contribution to the Jacobian
226  jacobian(local_eqn, local_unknown) +=
227  (beta_local * dpsidx(l2, i) * dtestdx(l, i) * W);
228  }
229  } // if (local_unknown>=0)
230  } // for (unsigned l2=0;l2<n_node;l2++)
231  } // if (flag)
232  } // if (local_eqn>=0)
233  } // for (unsigned l=0;l<n_node;l++)
234  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
235  } // End of fill_in_generic_residual_contribution_ust_heat
236 
237 
238  //======================================================================
239  /// Compute norm of FE solution
240  //======================================================================
241  template<unsigned SPATIAL_DIM>
243  {
244  // Initialise
245  norm = 0.0;
246 
247  // Vector of local coordinates
248  Vector<double> s(SPATIAL_DIM + 1, 0.0);
249 
250  // Vector for coordinates
251  Vector<double> x(SPATIAL_DIM + 1, 0.0);
252 
253  // Find out how many nodes there are in the element
254  unsigned n_node = nnode();
255 
256  // Allocate memory for the shape and test functions
257  Shape psi(n_node);
258 
259  // Set the value of n_intpt
260  unsigned n_intpt = integral_pt()->nweight();
261 
262  // Loop over the integration points
263  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
264  {
265  // Assign values of s
266  for (unsigned i = 0; i < SPATIAL_DIM + 1; i++)
267  {
268  // Get the i-th local coordinate at the ipt-th integration point
269  s[i] = integral_pt()->knot(ipt, i);
270  }
271 
272  // Get the integral weight
273  double w = integral_pt()->weight(ipt);
274 
275  // Get Jacobian of mapping
276  double J = J_eulerian(s);
277 
278  // Pre-multiply the weights and the Jacobian
279  double W = w * J;
280 
281  // Get FE function value
282  double u = interpolated_u_ust_heat(s);
283 
284  // Update the norm value
285  norm += u * u * W;
286  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
287  } // End of compute_norm
288 
289 
290  //======================================================================
291  /// Self-test: Return 0 for OK
292  //======================================================================
293  template<unsigned SPATIAL_DIM>
295  {
296  // Initialise the boolean variable
297  bool passed = true;
298 
299  // Check lower-level stuff
300  if (FiniteElement::self_test() != 0)
301  {
302  // If we get here then the lower-level self-tests did not pass
303  passed = false;
304  }
305 
306  // If the self-tests passed
307  if (passed)
308  {
309  // Return the value zero
310  return 0;
311  }
312  // If the self-tests didn't pass
313  else
314  {
315  // Return the value one
316  return 1;
317  }
318  } // End of self_test
319 
320 
321  //======================================================================
322  /// Output function:
323  /// x,t,u or x,y,t,u
324  /// at nplot points in each coordinate direction
325  //======================================================================
326  template<unsigned SPATIAL_DIM>
328  std::ostream& outfile, const unsigned& nplot)
329  {
330  // Vector of local coordinates
331  Vector<double> s(SPATIAL_DIM + 1, 0.0);
332 
333  // Tecplot header info
334  outfile << tecplot_zone_string(nplot);
335 
336  // Get the number of plot points
337  unsigned num_plot_points = nplot_points(nplot);
338 
339  // Loop over plot points
340  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
341  {
342  // Get local coordinates of plot point
343  get_s_plot(iplot, nplot, s);
344 
345  // Loop over the coordinate directions
346  for (unsigned i = 0; i < SPATIAL_DIM + 1; i++)
347  {
348  // Output the interpolated coordinate
349  outfile << interpolated_x(s, i) << " ";
350  }
351 
352  // Calculate the interpolated solution value
353  outfile << interpolated_u_ust_heat(s) << std::endl;
354  } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
355 
356  // Write tecplot footer (e.g. FE connectivity lists)
357  write_tecplot_zone_footer(outfile, nplot);
358  } // End of output
359 
360 
361  //======================================================================
362  /// C-style output function:
363  /// x,t,u or x,y,t,u
364  /// at nplot points in each coordinate direction
365  //======================================================================
366  template<unsigned SPATIAL_DIM>
368  FILE* file_pt, const unsigned& nplot)
369  {
370  // Vector of local coordinates
371  Vector<double> s(SPATIAL_DIM + 1, 0.0);
372 
373  // Tecplot header info
374  fprintf(file_pt, "%s", tecplot_zone_string(nplot).c_str());
375 
376  // Get the number of plot points
377  unsigned num_plot_points = nplot_points(nplot);
378 
379  // Loop over plot points
380  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
381  {
382  // Get local coordinates of plot point
383  get_s_plot(iplot, nplot, s);
384 
385  // Loop over the coordinate directions
386  for (unsigned i = 0; i < SPATIAL_DIM + 1; i++)
387  {
388  // Print the i-th coordinate value at local coordinate s
389  fprintf(file_pt, "%g ", interpolated_x(s, i));
390  }
391 
392  // Output the interpolated solution value at local coordinate s
393  fprintf(file_pt, "%g \n", interpolated_u_ust_heat(s));
394  } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
395 
396  // Write tecplot footer (e.g. FE connectivity lists)
397  write_tecplot_zone_footer(file_pt, nplot);
398  } // End of output
399 
400 
401  //======================================================================
402  /// Output exact solution at a given number of plot points:
403  /// x,t,u_exact or x,y,t,u_exact
404  /// Solution is provided via function pointer.
405  //======================================================================
406  template<unsigned SPATIAL_DIM>
408  std::ostream& outfile,
409  const unsigned& nplot,
411  {
412  // Vector of local coordinates
413  Vector<double> s(SPATIAL_DIM + 1, 0.0);
414 
415  // Vector for spatial coordinates
416  Vector<double> spatial_coordinates(SPATIAL_DIM, 0.0);
417 
418  // Tecplot header info
419  outfile << tecplot_zone_string(nplot);
420 
421  // Exact solution vector (here it's simply a scalar)
422  Vector<double> exact_soln(1, 0.0);
423 
424  // Get the number of plot points
425  unsigned num_plot_points = nplot_points(nplot);
426 
427  // Loop over plot points
428  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
429  {
430  // Get local coordinates of plot point
431  get_s_plot(iplot, nplot, s);
432 
433  // Loop over the spatial coordinates
434  for (unsigned i = 0; i < SPATIAL_DIM; i++)
435  {
436  // Assign the i-th spatial coordinate
437  spatial_coordinates[i] = interpolated_x(s, i);
438 
439  // Output the i-th coordinate at the point
440  outfile << spatial_coordinates[i] << " ";
441  }
442 
443  // Output the time value at this point
444  outfile << interpolated_x(s, SPATIAL_DIM) << " ";
445 
446  // Get the exact solution at this point
447  (*exact_soln_pt)(spatial_coordinates, exact_soln);
448 
449  // Output the exact solution at this point
450  outfile << exact_soln[0] << std::endl;
451  } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
452 
453  // Write tecplot footer (e.g. FE connectivity lists)
454  write_tecplot_zone_footer(outfile, nplot);
455  } // End of output_fct
456 
457 
458  //======================================================================
459  /// Output exact solution at a given number of plot points:
460  /// x,t,u_exact or x,y,t,u_exact
461  /// Solution is provided via function pointer.
462  //======================================================================
463  template<unsigned SPATIAL_DIM>
465  std::ostream& outfile,
466  const unsigned& nplot,
467  const double& time,
469  {
470  // Storage for the time value
471  double interpolated_t = 0.0;
472 
473  // Vector of local coordinates
474  Vector<double> s(SPATIAL_DIM + 1, 0.0);
475 
476  // Vector for spatial coordinates
477  Vector<double> spatial_coordinates(SPATIAL_DIM, 0.0);
478 
479  // Tecplot header info
480  outfile << tecplot_zone_string(nplot);
481 
482  // Exact solution vector (here it's simply a scalar)
483  Vector<double> exact_soln(1, 0.0);
484 
485  // Get the number of plot points
486  unsigned num_plot_points = nplot_points(nplot);
487 
488  // Loop over plot points
489  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
490  {
491  // Get local coordinates of plot point
492  get_s_plot(iplot, nplot, s);
493 
494  // Loop over the spatial coordinates
495  for (unsigned i = 0; i < SPATIAL_DIM; i++)
496  {
497  // Assign the i-th spatial coordinate
498  spatial_coordinates[i] = interpolated_x(s, i);
499 
500  // Output the i-th coordinate at the point
501  outfile << spatial_coordinates[i] << " ";
502  }
503 
504  // Get the time value
505  interpolated_t = interpolated_x(s, SPATIAL_DIM);
506 
507  // Output the time value at this point
508  outfile << interpolated_t << " ";
509 
510  // Get the exact solution at this point
511  (*exact_soln_pt)(interpolated_t, spatial_coordinates, exact_soln);
512 
513  // Output the exact solution at this point
514  outfile << exact_soln[0] << std::endl;
515  } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
516 
517  // Write tecplot footer (e.g. FE connectivity lists)
518  write_tecplot_zone_footer(outfile, nplot);
519  } // End of output_fct
520 
521 
522  //======================================================================
523  /// Validate against exact solution
524  ///
525  /// Solution is provided via function pointer.
526  /// Plot error at a given number of plot points.
527  //======================================================================
528  template<unsigned SPATIAL_DIM>
530  std::ostream& outfile,
532  double& error,
533  double& norm)
534  {
535  // Initialise error value
536  error = 0.0;
537 
538  // Initialise norm value
539  norm = 0.0;
540 
541  // Vector of local coordinates
542  Vector<double> s(SPATIAL_DIM + 1, 0.0);
543 
544  // Vector for spatial coordinates
545  Vector<double> spatial_coordinates(SPATIAL_DIM, 0.0);
546 
547  // Find out how many nodes there are in the element
548  unsigned n_node = nnode();
549 
550  // Initialise shape functions
551  Shape psi(n_node);
552 
553  // Set the value of n_intpt
554  unsigned n_intpt = integral_pt()->nweight();
555 
556  // Tecplot header info
557  outfile << "ZONE" << std::endl;
558 
559  // Exact solution vector (here it's simply a scalar)
560  Vector<double> exact_soln(1, 0.0);
561 
562  // Loop over the integration points
563  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
564  {
565  // Assign values of s
566  for (unsigned i = 0; i < SPATIAL_DIM + 1; i++)
567  {
568  // Get the i-th local coordinate at the ipt-th integration point
569  s[i] = integral_pt()->knot(ipt, i);
570  }
571 
572  // Get the integral weight
573  double w = integral_pt()->weight(ipt);
574 
575  // Get jacobian of mapping
576  double J = J_eulerian(s);
577 
578  // Premultiply the weights and the Jacobian
579  double W = w * J;
580 
581  // Get FE function value
582  double u_fe = interpolated_u_ust_heat(s);
583 
584  // Loop over the spatial coordinates
585  for (unsigned i = 0; i < SPATIAL_DIM; i++)
586  {
587  // Assign the i-th spatial coordinate
588  spatial_coordinates[i] = interpolated_x(s, i);
589 
590  // Output the i-th coordinate at the point
591  outfile << spatial_coordinates[i] << " ";
592  }
593 
594  // Output the i-th coordinate at this point
595  outfile << interpolated_x(s, SPATIAL_DIM) << " ";
596 
597  // Get exact solution at this point
598  (*exact_soln_pt)(spatial_coordinates, exact_soln);
599 
600  // Output the error
601  outfile << exact_soln[0] << " " << exact_soln[0] - u_fe << std::endl;
602 
603  // Add to (exact) solution norm value
604  norm += exact_soln[0] * exact_soln[0] * W;
605 
606  // Update the error norm value
607  error += (exact_soln[0] - u_fe) * (exact_soln[0] - u_fe) * W;
608  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
609  } // End of compute_error
610 
611 
612  //======================================================================
613  /// Validate against exact solution at time t.
614  ///
615  /// Solution is provided via function pointer.
616  /// Plot error at a given number of plot points.
617  //======================================================================
618  template<unsigned SPATIAL_DIM>
620  std::ostream& outfile,
622  const double& time,
623  double& error,
624  double& norm)
625  {
626  // Initialise error value
627  error = 0.0;
628 
629  // Initialise norm value
630  norm = 0.0;
631 
632  // Storage for the time value
633  double interpolated_t = 0.0;
634 
635  // Vector of local coordinates
636  Vector<double> s(SPATIAL_DIM + 1, 0.0);
637 
638  // Vector for spatial coordinates
639  Vector<double> spatial_coordinates(SPATIAL_DIM, 0.0);
640 
641  // Find out how many nodes there are in the element
642  unsigned n_node = nnode();
643 
644  // Initialise shape functions
645  Shape psi(n_node);
646 
647  // Set the value of n_intpt
648  unsigned n_intpt = integral_pt()->nweight();
649 
650  // Tecplot header info
651  outfile << "ZONE" << std::endl;
652 
653  // Exact solution vector (here it's simply a scalar)
654  Vector<double> exact_soln(1, 0.0);
655 
656  // Loop over the integration points
657  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
658  {
659  // Assign values of s
660  for (unsigned i = 0; i < SPATIAL_DIM + 1; i++)
661  {
662  s[i] = integral_pt()->knot(ipt, i);
663  }
664 
665  // Get the integral weight
666  double w = integral_pt()->weight(ipt);
667 
668  // Get jacobian of mapping
669  double J = J_eulerian(s);
670 
671  // Premultiply the weights and the Jacobian
672  double W = w * J;
673 
674  // Get FE function value
675  double u_fe = interpolated_u_ust_heat(s);
676 
677  // Loop over the spatial coordinates
678  for (unsigned i = 0; i < SPATIAL_DIM; i++)
679  {
680  // Assign the i-th spatial coordinate
681  spatial_coordinates[i] = interpolated_x(s, i);
682 
683  // Output the i-th coordinate at the point
684  outfile << spatial_coordinates[i] << " ";
685  }
686 
687  // Get the time value
688  interpolated_t = interpolated_x(s, SPATIAL_DIM);
689 
690  // Output the time value at this point
691  outfile << interpolated_t << " ";
692 
693  // Get the exact solution at this point
694  (*exact_soln_pt)(interpolated_t, spatial_coordinates, exact_soln);
695 
696  // Output the error
697  outfile << exact_soln[0] << " " << exact_soln[0] - u_fe << std::endl;
698 
699  // Add to (exact) solution norm value
700  norm += exact_soln[0] * exact_soln[0] * W;
701 
702  // Update the error norm value
703  error += (exact_soln[0] - u_fe) * (exact_soln[0] - u_fe) * W;
704  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
705  } // End of compute_error
706 
707 
708  //======================================================================
709  /// Output function:
710  /// x,t,u or x,y,t,u
711  /// at nplot points in each coordinate direction
712  //======================================================================
713  template<unsigned SPATIAL_DIM>
715  std::ofstream& file_out, const unsigned& nplot)
716  {
717  // Change the scientific format so that E is used rather than e
718  file_out.setf(std::ios_base::uppercase);
719 
720  // Make variables to hold the number of nodes and elements
721  unsigned number_of_nodes = this->nplot_points_paraview(nplot);
722 
723  // Make variables to hold the number of elements
724  unsigned total_number_of_elements = this->nsub_elements_paraview(nplot);
725 
726  //------------------
727  // File Declaration:
728  //------------------
729  // Insert the necessary lines plus header of file, and
730  // number of nodes and elements
731  file_out << "<?xml version=\"1.0\"?>\n"
732  << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
733  << "byte_order=\"LittleEndian\">\n"
734  << "<UnstructuredGrid>\n"
735  << "<Piece NumberOfPoints=\"" << number_of_nodes
736  << "\" NumberOfCells=\"" << total_number_of_elements << "\">\n";
737 
738  //------------
739  // Point Data:
740  //------------
741  // Check the number of degrees of freedom
742  unsigned ndof = this->nscalar_paraview();
743 
744  // Point data is going in here
745  file_out << "<PointData ";
746 
747  // Insert just the first scalar name, since paraview reads everything
748  // else after that as being of the same type. Get information from
749  // first element.
750  file_out << "Scalars=\"" << this->scalar_name_paraview(0) << "\">\n";
751 
752  // Loop over i scalar fields and j number of elements
753  for (unsigned i = 0; i < ndof; i++)
754  {
755  file_out << "<DataArray type=\"Float32\" "
756  << "Name=\"" << this->scalar_name_paraview(i) << "\" "
757  << "format=\"ascii\""
758  << ">\n";
759 
760  // Output the i-th scalar field with nplot plot points
761  this->scalar_value_paraview(file_out, i, nplot);
762 
763  // Close of the DataArray
764  file_out << "</DataArray>\n";
765  }
766 
767  // Close off the PointData set
768  file_out << "</PointData>\n";
769 
770  //------------------
771  // Geometric Points:
772  //------------------
773  // Always has to be 3 components for an unstructured grid
774  file_out << "<Points>\n"
775  << "<DataArray type=\"Float32\""
776  << " NumberOfComponents=\"" << 3 << "\" "
777  << "format=\"ascii\">\n";
778 
779  // Print the plot points
780  this->output_paraview(file_out, nplot);
781 
782  // Close off the geometric points set
783  file_out << "</DataArray>\n"
784  << "</Points>\n";
785 
786  //-------
787  // Cells:
788  //-------
789  file_out << "<Cells>\n"
790  << "<DataArray type=\"Int32\" Name=\""
791  << "connectivity\" format=\"ascii\">\n";
792 
793  // Make counter for keeping track of all the local elements,
794  // because Paraview requires global coordinates
795  unsigned counter = 0;
796 
797  // Write connectivity with the local elements
798  this->write_paraview_output_offset_information(file_out, nplot, counter);
799 
800  // Output header stuff
801  file_out << "</DataArray>\n"
802  << "<DataArray type=\"Int32\" "
803  << "Name=\"offsets\" format=\"ascii\">\n";
804 
805  // Make variable that holds the current offset number
806  unsigned offset_sum = 0;
807 
808  // Write the offset for the specific elements
809  this->write_paraview_offsets(file_out, nplot, offset_sum);
810 
811  // Add in header information
812  file_out << "</DataArray>\n"
813  << "<DataArray type=\"UInt8\" Name=\"types\">\n";
814 
815  // Get the type the element has
816  this->write_paraview_type(file_out, nplot);
817 
818  // Finish off the data set
819  file_out << "</DataArray>\n"
820  << "</Cells>\n";
821 
822  //--------------
823  // File Closure:
824  //--------------
825  file_out << "</Piece>\n"
826  << "</UnstructuredGrid>\n"
827  << "</VTKFile>";
828  } // End of output_element_paraview
829 
830 
831  //====================================================================
832  // Force build of templates
833  //====================================================================
837 
841 } // End of namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
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
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK.
Definition: elements.cc:4440
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
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
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.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error and norm against exact solution.
void compute_norm(double &norm)
Compute norm of FE solution.
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(std::ostream &outfile)
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.
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...