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-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// Non-inline functions for SpaceTimeUnsteadyHeat elements
28
29/// //////////////////////////////////////////////////////////////////////
30/// //////////////////////////////////////////////////////////////////////
31/// //////////////////////////////////////////////////////////////////////
32
33namespace 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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...