discontinuous_galerkin_equal_order_pressure_spacetime_navier_stokes_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 NS elements
28
29/// //////////////////////////////////////////////////////////////////////
30/// //////////////////////////////////////////////////////////////////////
31/// //////////////////////////////////////////////////////////////////////
32
33namespace oomph
34{
35 /// Navier-Stokes equations static data
36 template<unsigned DIM>
37 Vector<double> SpaceTimeNavierStokesEquations<DIM>::Gamma(DIM, 1.0);
38
39 //===================================================================
40 /// "Magic" negative number that indicates that the pressure is
41 /// not stored at a node. This cannot be -1 because that represents
42 /// the positional hanging scheme in the hanging_pt object of nodes.
43 //===================================================================
44 template<unsigned DIM>
46
47 /// Navier-Stokes equations static data
48 template<unsigned DIM>
50 0.0;
51
52 /// Navier-Stokes equations static data
53 template<unsigned DIM>
55 1.0;
56
57 /// Navier-Stokes equations default gravity vector
58 template<unsigned DIM>
60 DIM, 0.0);
61
62 //===================================================================
63 /// Compute the diagonal of the velocity/pressure mass matrices.
64 /// If which one=0, both are computed, otherwise only the pressure
65 /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
66 /// LSC version of the preconditioner only needs that one).
67 //===================================================================
68 template<unsigned DIM>
71 Vector<double>& press_mass_diag,
72 Vector<double>& veloc_mass_diag,
73 const unsigned& which_one)
74 {
75 // Resize and initialise
76 unsigned n_dof = ndof();
77
78 // If space needs to be assigned for the pressure mass matrix
79 if ((which_one == 0) || (which_one == 1))
80 {
81 // Assign the appropriate amount of space
82 press_mass_diag.assign(n_dof, 0.0);
83 }
84
85 // If space needs to be assigned for the velocity mass matrix
86 if ((which_one == 0) || (which_one == 2))
87 {
88 // Assign the appropriate amount of space
89 veloc_mass_diag.assign(n_dof, 0.0);
90 }
91
92 // Find out how many nodes there are
93 unsigned n_node = nnode();
94
95 // Storage for the local coordinates
96 Vector<double> s(DIM + 1, 0.0);
97
98 // Storage for the local velocity indices
99 Vector<unsigned> u_nodal_index(DIM, 0.0);
100
101 // Find the indices at which the local velocities are stored
102 for (unsigned i = 0; i < DIM; i++)
103 {
104 // Calculate the i-th local velocity component
105 u_nodal_index[i] = this->u_index_nst(i);
106 }
107
108 // Set up memory for velocity shape functions
109 Shape psi(n_node);
110
111 // Find number of pressure dofs
112 unsigned n_pres = npres_nst();
113
114 // Pressure shape function
115 Shape psi_p(n_pres);
116
117 // Number of integration points
118 unsigned n_intpt = integral_pt()->nweight();
119
120 // Integer to store the local equations no
121 int local_eqn = 0;
122
123 // Loop over the integration points
124 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
125 {
126 // Get the integral weight
127 double w = integral_pt()->weight(ipt);
128
129 // Get determinant of Jacobian of the mapping
130 double J = J_eulerian_at_knot(ipt);
131
132 // Assign values of s
133 for (unsigned i = 0; i < DIM + 1; i++)
134 {
135 // Calculate the i-th local coordinate at the ipt-th integration point
136 s[i] = integral_pt()->knot(ipt, i);
137 }
138
139 // Premultiply weights and Jacobian
140 double W = w * J;
141
142 // Do we want the velocity one?
143 if ((which_one == 0) || (which_one == 2))
144 {
145 // Get the velocity shape functions
146 shape_at_knot(ipt, psi);
147
148 // Loop over the velocity shape functions
149 for (unsigned l = 0; l < n_node; l++)
150 {
151 // Loop over the velocity components
152 for (unsigned i = 0; i < DIM; i++)
153 {
154 // Get the local equation number
155 local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
156
157 // If not a boundary condition
158 if (local_eqn >= 0)
159 {
160 // Add the contribution
161 veloc_mass_diag[local_eqn] += pow(psi[l], 2) * W;
162 }
163 } // for (unsigned i=0;i<n_dim;i++)
164 } // for (unsigned l=0;l<n_node;l++)
165 } // if ((which_one==0)||(which_one==2))
166
167 // Do we want the pressure one?
168 if ((which_one == 0) || (which_one == 1))
169 {
170 // Get the pressure shape functions
171 pshape_nst(s, psi_p);
172
173 // Loop over the veclocity shape functions
174 for (unsigned l = 0; l < n_pres; l++)
175 {
176 // Get equation number
177 local_eqn = p_local_eqn(l);
178
179 // If not a boundary condition
180 if (local_eqn >= 0)
181 {
182 // Add the contribution
183 press_mass_diag[local_eqn] += pow(psi_p[l], 2) * W;
184 }
185 } // for (unsigned l=0;l<n_pres;l++)
186 } // if ((which_one==0)||(which_one==1))
187 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
188 } // End of get_pressure_and_velocity_mass_matrix_diagonal
189
190
191 //======================================================================
192 /// Compute the vector norm of FEM solution
193 //======================================================================
194 template<unsigned DIM>
196 {
197 // Resize the solution norm vector
198 norm.resize(DIM + 1, 0.0);
199
200 // Vector of local coordinates
201 Vector<double> s(DIM + 1, 0.0);
202
203 // Set the value of n_intpt
204 unsigned n_intpt = integral_pt()->nweight();
205
206 // Loop over the integration points
207 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
208 {
209 // Assign values of s
210 for (unsigned i = 0; i < DIM + 1; i++)
211 {
212 s[i] = integral_pt()->knot(ipt, i);
213 }
214
215 // Get the integral weight
216 double w = integral_pt()->weight(ipt);
217
218 // Get Jacobian of mapping
219 double J = J_eulerian(s);
220
221 // Premultiply the weights and the Jacobian
222 double W = w * J;
223
224 // Compute the velocity norm
225 for (unsigned i = 0; i < DIM; i++)
226 {
227 // Update the solution norm value
228 norm[i] += pow(interpolated_u_nst(s, i), 2) * W;
229 }
230
231 // Update the pressure norm value
232 norm[DIM] += pow(interpolated_p_nst(s), 2) * W;
233 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
234 } // End of compute_norm
235
236
237 //======================================================================
238 /// Validate against exact velocity solution at given time. Solution is
239 /// provided via function pointer. Plot at a given number of plot points
240 /// and compute L2 error and L2 norm of velocity solution over element.
241 //=======================================================================
242 template<unsigned DIM>
244 std::ostream& outfile,
246 const double& time,
247 double& error,
248 double& norm)
249 {
250 // Initialise the error norm value
251 error = 0.0;
252
253 // Initialise the solution norm value
254 norm = 0.0;
255
256 // Storage for the time value
257 double interpolated_t = 0.0;
258
259 // Vector of local coordinates
260 Vector<double> s(DIM + 1, 0.0);
261
262 // Vector for the spatial coordinates
263 Vector<double> x(DIM, 0.0);
264
265 // Set the value of n_intpt
266 unsigned n_intpt = integral_pt()->nweight();
267
268 // Exact solution Vector (u,v,[w],p)
269 Vector<double> exact_soln(DIM + 1, 0.0);
270
271 // Loop over the integration points
272 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
273 {
274 // Assign values of s
275 for (unsigned i = 0; i < DIM + 1; i++)
276 {
277 s[i] = integral_pt()->knot(ipt, i);
278 }
279
280 // Assign values of x
281 for (unsigned i = 0; i < DIM; i++)
282 {
283 // Get x position as Vector
284 x[i] = interpolated_x(s, i);
285 }
286
287 // Get the time value
288 interpolated_t = interpolated_x(s, DIM);
289
290 // Get exact solution at this point
291 (*exact_soln_pt)(interpolated_t, x, exact_soln);
292
293 // Get the integral weight
294 double w = integral_pt()->weight(ipt);
295
296 // Get Jacobian of mapping
297 double J = J_eulerian(s);
298
299 // Premultiply the weights and the Jacobian
300 double W = w * J;
301
302 // Velocity error
303 for (unsigned i = 0; i < DIM; i++)
304 {
305 // Update the solution norm value
306 norm += exact_soln[i] * exact_soln[i] * W;
307
308 // Update the error norm value
309 error += pow(exact_soln[i] - interpolated_u_nst(s, i), 2) * W;
310 }
311
312 // ------ DRAIG: REMOVE ----------------------------------------
313 // Update the solution norm value
314 norm += pow(exact_soln[DIM], 2) * W;
315
316 // Update the error norm value
317 error += pow(exact_soln[DIM] - interpolated_p_nst(s), 2) * W;
318 // ------ DRAIG: REMOVE ----------------------------------------
319 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
320
321 //------------------------------------------------
322 // Output the error at the appropriate plot points
323 //------------------------------------------------
324 // Just output at the default number of plot points
325 unsigned n_plot = 5;
326
327 // Tecplot header info
328 outfile << tecplot_zone_string(n_plot);
329
330 // How many plot points are there in total?
331 unsigned num_plot_points = nplot_points(n_plot);
332
333 // Loop over plot points
334 for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
335 {
336 // Get local coordinates of plot point
337 get_s_plot(i_plot, n_plot, s);
338
339 // Loop over the spatial coordinates
340 for (unsigned i = 0; i < DIM; i++)
341 {
342 // Assign the i-th spatial coordinate
343 x[i] = interpolated_x(s, i);
344 }
345
346 // Get the time value
347 interpolated_t = interpolated_x(s, DIM);
348
349 // Get exact solution at this point
350 (*exact_soln_pt)(interpolated_t, x, exact_soln);
351
352 // Output x,y,...
353 for (unsigned i = 0; i < DIM; i++)
354 {
355 // Output the i-th spatial coordinate
356 outfile << x[i] << " ";
357 }
358
359 // Output the time value at this point
360 outfile << interpolated_t << " ";
361
362 // Output u_error,v_error
363 for (unsigned i = 0; i < DIM; i++)
364 {
365 // Output the error in the i-th velocity component at this point
366 outfile << exact_soln[i] - interpolated_u_nst(s, i) << " ";
367 }
368
369 // Output the error in the pressure field at this point
370 outfile << exact_soln[DIM] - interpolated_p_nst(s) << " ";
371
372 // End the line
373 outfile << std::endl;
374 } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
375
376 // Write tecplot footer (e.g. FE connectivity lists)
377 write_tecplot_zone_footer(outfile, n_plot);
378 } // End of compute_error
379
380
381 //======================================================================
382 /// Validate against exact velocity solution at given time. Solution is
383 /// provided via function pointer. Plot at a given number of plot points
384 /// and compute L2 error and L2 norm of velocity solution over element.
385 //=======================================================================
386 template<unsigned DIM>
388 std::ostream& outfile,
390 const double& time,
391 Vector<double>& error,
392 Vector<double>& norm)
393 {
394 // Resize the error norm vector
395 error.resize(DIM + 1, 0.0);
396
397 // Resize the solution norm vector
398 norm.resize(DIM + 1, 0.0);
399
400 // Storage for the time value
401 double interpolated_t = 0.0;
402
403 // Vector of local coordinates
404 Vector<double> s(DIM + 1, 0.0);
405
406 // Vector for the spatial coordinates
407 Vector<double> x(DIM, 0.0);
408
409 // Set the value of n_intpt
410 unsigned n_intpt = integral_pt()->nweight();
411
412 // Exact solution Vector (u,v,[w],p)
413 Vector<double> exact_soln(DIM + 1, 0.0);
414
415 // Loop over the integration points
416 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
417 {
418 // Assign values of s
419 for (unsigned i = 0; i < DIM + 1; i++)
420 {
421 s[i] = integral_pt()->knot(ipt, i);
422 }
423
424 // Assign values of x
425 for (unsigned i = 0; i < DIM; i++)
426 {
427 // Get x position as Vector
428 x[i] = interpolated_x(s, i);
429 }
430
431 // Get the time value
432 interpolated_t = interpolated_x(s, DIM);
433
434 // Get exact solution at this point
435 (*exact_soln_pt)(interpolated_t, x, exact_soln);
436
437 // Get the integral weight
438 double w = integral_pt()->weight(ipt);
439
440 // Get Jacobian of mapping
441 double J = J_eulerian(s);
442
443 // Premultiply the weights and the Jacobian
444 double W = w * J;
445
446 // Velocity error
447 for (unsigned i = 0; i < DIM; i++)
448 {
449 // Update the solution norm value
450 norm[i] += exact_soln[i] * exact_soln[i] * W;
451
452 // Update the error norm value
453 error[i] += pow(exact_soln[i] - interpolated_u_nst(s, i), 2) * W;
454 }
455
456 // ------ DRAIG: REMOVE ----------------------------------------
457 // Update the solution norm value
458 norm[DIM] += pow(exact_soln[DIM], 2) * W;
459
460 // Update the error norm value
461 error[DIM] += pow(exact_soln[DIM] - interpolated_p_nst(s), 2) * W;
462 // ------ DRAIG: REMOVE ----------------------------------------
463 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
464
465 //------------------------------------------------
466 // Output the error at the appropriate plot points
467 //------------------------------------------------
468 // Just output at the default number of plot points
469 unsigned n_plot = 5;
470
471 // Tecplot header info
472 outfile << tecplot_zone_string(n_plot);
473
474 // How many plot points are there in total?
475 unsigned num_plot_points = nplot_points(n_plot);
476
477 // Loop over plot points
478 for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
479 {
480 // Get local coordinates of plot point
481 get_s_plot(i_plot, n_plot, s);
482
483 // Loop over the spatial coordinates
484 for (unsigned i = 0; i < DIM; i++)
485 {
486 // Assign the i-th spatial coordinate
487 x[i] = interpolated_x(s, i);
488 }
489
490 // Get the time value
491 interpolated_t = interpolated_x(s, DIM);
492
493 // Get exact solution at this point
494 (*exact_soln_pt)(interpolated_t, x, exact_soln);
495
496 // Output x,y,...
497 for (unsigned i = 0; i < DIM; i++)
498 {
499 // Output the i-th spatial coordinate
500 outfile << x[i] << " ";
501 }
502
503 // Output the time value at this point
504 outfile << interpolated_t << " ";
505
506 // Output u_error,v_error
507 for (unsigned i = 0; i < DIM; i++)
508 {
509 // Output the error in the i-th velocity component at this point
510 outfile << exact_soln[i] - interpolated_u_nst(s, i) << " ";
511 }
512
513 // Output the error in the pressure field at this point
514 outfile << exact_soln[DIM] - interpolated_p_nst(s) << " ";
515
516 // End the line
517 outfile << std::endl;
518 } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
519
520 // Write tecplot footer (e.g. FE connectivity lists)
521 write_tecplot_zone_footer(outfile, n_plot);
522 } // End of compute_error
523
524
525 //======================================================================
526 /// Validate against exact velocity solution at given time.
527 /// Solution is provided via function pointer.
528 /// Compute L2 error and L2 norm of velocity solution over element.
529 //=======================================================================
530 template<unsigned DIM>
533 const double& time,
534 double& error,
535 double& norm)
536 {
537 // Initialise the error norm value
538 error = 0.0;
539
540 // Initialise the solution norm value
541 norm = 0.0;
542
543 // Storage for the time value
544 double interpolated_t = 0.0;
545
546 // Vector of local coordinates
547 Vector<double> s(DIM + 1, 0.0);
548
549 // Vector for the spatial coordinates
550 Vector<double> x(DIM, 0.0);
551
552 // Set the value of n_intpt
553 unsigned n_intpt = integral_pt()->nweight();
554
555 // Exact solution Vector (u,v,[w],p)
556 Vector<double> exact_soln(DIM + 1, 0.0);
557
558 // Loop over the integration points
559 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
560 {
561 // Assign values of s
562 for (unsigned i = 0; i < DIM + 1; i++)
563 {
564 s[i] = integral_pt()->knot(ipt, i);
565 }
566
567 // Assign values of x
568 for (unsigned i = 0; i < DIM; i++)
569 {
570 // Get x position as Vector
571 x[i] = interpolated_x(s, i);
572 }
573
574 // Get the time value
575 interpolated_t = interpolated_x(s, DIM);
576
577 // Get exact solution at this point
578 (*exact_soln_pt)(interpolated_t, x, exact_soln);
579
580 // Get the integral weight
581 double w = integral_pt()->weight(ipt);
582
583 // Get Jacobian of mapping
584 double J = J_eulerian(s);
585
586 // Premultiply the weights and the Jacobian
587 double W = w * J;
588
589 // Velocity error
590 for (unsigned i = 0; i < DIM; i++)
591 {
592 // Update the solution norm value
593 norm += exact_soln[i] * exact_soln[i] * W;
594
595 // Update the error norm value
596 error += pow(exact_soln[i] - interpolated_u_nst(s, i), 2) * W;
597 }
598 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
599 } // End of compute_error
600
601 //======================================================================
602 /// Validate against exact velocity solution Solution is provided via a
603 /// function pointer. Compute L2 error and L2 norm of velocity solution
604 /// over element.
605 //=======================================================================
606 template<unsigned DIM>
609 double& error,
610 double& norm)
611 {
612 // Initialise the error norm value
613 error = 0.0;
614
615 // Initialise the solution norm value
616 norm = 0.0;
617
618 // Vector of local coordinates
619 Vector<double> s(DIM + 1, 0.0);
620
621 // Vector for the spatial coordinates
622 Vector<double> x(DIM, 0.0);
623
624 // Set the value of n_intpt
625 unsigned n_intpt = integral_pt()->nweight();
626
627 // Exact solution Vector (u,v,[w],p)
628 Vector<double> exact_soln(DIM + 1, 0.0);
629
630 // Loop over the integration points
631 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
632 {
633 // Assign values of s
634 for (unsigned i = 0; i < DIM + 1; i++)
635 {
636 s[i] = integral_pt()->knot(ipt, i);
637 }
638
639 // Assign values of x
640 for (unsigned i = 0; i < DIM; i++)
641 {
642 // Get x position as Vector
643 x[i] = interpolated_x(s, i);
644 }
645
646 // Get exact solution at this point
647 (*exact_soln_pt)(x, exact_soln);
648
649 // Get the integral weight
650 double w = integral_pt()->weight(ipt);
651
652 // Get Jacobian of mapping
653 double J = J_eulerian(s);
654
655 // Premultiply the weights and the Jacobian
656 double W = w * J;
657
658 // Velocity error
659 for (unsigned i = 0; i < DIM; i++)
660 {
661 // Update the solution norm value
662 norm += exact_soln[i] * exact_soln[i] * W;
663
664 // Update the error norm value
665 error += pow(exact_soln[i] - interpolated_u_nst(s, i), 2) * W;
666 }
667 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
668 } // End of compute_error
669
670 //======================================================================
671 /// Validate against exact velocity solution
672 /// Solution is provided via function pointer.
673 /// Plot at a given number of plot points and compute L2 error
674 /// and L2 norm of velocity solution over element.
675 //=======================================================================
676 template<unsigned DIM>
678 std::ostream& outfile,
680 double& error,
681 double& norm)
682 {
683 // Initialise the error norm value
684 error = 0.0;
685
686 // Initialise the solution norm value
687 norm = 0.0;
688
689 // Storage for the time value
690 double interpolated_t = 0.0;
691
692 // Vector of local coordinates
693 Vector<double> s(DIM + 1, 0.0);
694
695 // Vector for the spatial coordinates
696 Vector<double> x(DIM, 0.0);
697
698 // Set the value of n_intpt
699 unsigned n_intpt = integral_pt()->nweight();
700
701 // Output the tecplot header
702 outfile << "ZONE" << std::endl;
703
704 // Exact solution Vector (u,v,[w],p)
705 Vector<double> exact_soln(DIM + 1, 0.0);
706
707 // Loop over the integration points
708 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
709 {
710 // Assign values of s
711 for (unsigned i = 0; i < DIM + 1; i++)
712 {
713 s[i] = integral_pt()->knot(ipt, i);
714 }
715
716 // Assign values of x
717 for (unsigned i = 0; i < DIM; i++)
718 {
719 // Get x position as Vector
720 x[i] = interpolated_x(s, i);
721 }
722
723 // Get the time value
724 interpolated_t = interpolated_x(s, DIM);
725
726 // Get exact solution at this point
727 (*exact_soln_pt)(x, exact_soln);
728
729 // Get the integral weight
730 double w = integral_pt()->weight(ipt);
731
732 // Get Jacobian of mapping
733 double J = J_eulerian(s);
734
735 // Premultiply the weights and the Jacobian
736 double W = w * J;
737
738 // Velocity error
739 for (unsigned i = 0; i < DIM; i++)
740 {
741 // Update the solution norm value
742 norm += exact_soln[i] * exact_soln[i] * W;
743
744 // Update the error norm value
745 error += pow(exact_soln[i] - interpolated_u_nst(s, i), 2) * W;
746 }
747
748 // Output x,y,...,u_exact
749 for (unsigned i = 0; i < DIM; i++)
750 {
751 // Output the i-th coordinate value
752 outfile << x[i] << " ";
753 }
754
755 // Output the time value at this point
756 outfile << interpolated_t << " ";
757
758 // Output x,y,u_error,v_error
759 for (unsigned i = 0; i < DIM; i++)
760 {
761 // Output the error in the i-th velocity component at this point
762 outfile << exact_soln[i] - interpolated_u_nst(s, i) << " ";
763 }
764
765 // Finish the line off
766 outfile << std::endl;
767 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
768 } // End of compute_error
769
770 //======================================================================
771 /// Output "exact" solution
772 /// Solution is provided via function pointer.
773 /// Plot at a given number of plot points.
774 /// Function prints as many components as are returned in solution Vector.
775 //=======================================================================
776 template<unsigned DIM>
778 std::ostream& outfile,
779 const unsigned& n_plot,
781 {
782 // Storage for the time value
783 double interpolated_t = 0.0;
784
785 // Vector of local coordinates
786 Vector<double> s(DIM + 1, 0.0);
787
788 // Vector for coordinates
789 Vector<double> x(DIM, 0.0);
790
791 // Tecplot header info
792 outfile << tecplot_zone_string(n_plot);
793
794 // Exact solution Vector
795 Vector<double> exact_soln;
796
797 // How many plot points are there in total?
798 unsigned num_plot_points = nplot_points(n_plot);
799
800 // Loop over plot points
801 for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
802 {
803 // Get local coordinates of plot point
804 get_s_plot(i_plot, n_plot, s);
805
806 // Loop over the spatial coordinates
807 for (unsigned i = 0; i < DIM; i++)
808 {
809 // Assign the i-th spatial coordinate
810 x[i] = interpolated_x(s, i);
811 }
812
813 // Get the time value
814 interpolated_t = interpolated_x(s, DIM);
815
816 // Get exact solution at this point
817 (*exact_soln_pt)(x, exact_soln);
818
819 // Output x,y,...
820 for (unsigned i = 0; i < DIM; i++)
821 {
822 // Output the i-th spatial coordinate
823 outfile << x[i] << " ";
824 }
825
826 // Output the time value at this point
827 outfile << interpolated_t << " ";
828
829 // Output "exact solution"
830 for (unsigned i = 0; i < exact_soln.size(); i++)
831 {
832 // Output the i-th (exact) velocity component value
833 outfile << exact_soln[i] << " ";
834 }
835
836 // End the line
837 outfile << std::endl;
838 } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
839
840 // Write tecplot footer (e.g. FE connectivity lists)
841 write_tecplot_zone_footer(outfile, n_plot);
842 } // End of output_fct
843
844 //======================================================================
845 /// Output "exact" solution at a given time
846 /// Solution is provided via function pointer.
847 /// Plot at a given number of plot points.
848 /// Function prints as many components as are returned in solution Vector.
849 //=======================================================================
850 template<unsigned DIM>
852 std::ostream& outfile,
853 const unsigned& n_plot,
854 const double& time,
856 {
857 // Storage for the time value
858 double interpolated_t = 0.0;
859
860 // Vector of local coordinates
861 Vector<double> s(DIM + 1, 0.0);
862
863 // Vector for coordinates
864 Vector<double> x(DIM, 0.0);
865
866 // Tecplot header info
867 outfile << tecplot_zone_string(n_plot);
868
869 // Exact solution Vector
870 Vector<double> exact_soln;
871
872 // How many plot points are there in total?
873 unsigned num_plot_points = nplot_points(n_plot);
874
875 // Loop over plot points
876 for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
877 {
878 // Get local coordinates of plot point
879 get_s_plot(i_plot, n_plot, s);
880
881 // Loop over the spatial coordinates
882 for (unsigned i = 0; i < DIM; i++)
883 {
884 // Assign the i-th spatial coordinate
885 x[i] = interpolated_x(s, i);
886 }
887
888 // Get the time value
889 interpolated_t = interpolated_x(s, DIM);
890
891 // Get exact solution at this point
892 (*exact_soln_pt)(interpolated_t, x, exact_soln);
893
894 // Output x,y,...
895 for (unsigned i = 0; i < DIM; i++)
896 {
897 // Output the i-th spatial coordinate value
898 outfile << x[i] << " ";
899 }
900
901 // Output the time value at this point
902 outfile << interpolated_t << " ";
903
904 // Output "exact solution"
905 for (unsigned i = 0; i < exact_soln.size(); i++)
906 {
907 // Output the i-th (exact) velocity component value
908 outfile << exact_soln[i] << " ";
909 }
910
911 // End the line
912 outfile << std::endl;
913 } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
914
915 // Write tecplot footer (e.g. FE connectivity lists)
916 write_tecplot_zone_footer(outfile, n_plot);
917 } // End of output_fct
918
919 //==============================================================
920 /// Output function: Velocities only
921 /// x,y,[z],u,v,[w]
922 /// in tecplot format at specified previous timestep (t=0: present;
923 /// t>0: previous timestep). Specified number of plot points in each
924 /// coordinate direction.
925 /// DRAIG: Should be broken!
926 //==============================================================
927 template<unsigned DIM>
929 const unsigned& n_plot,
930 const unsigned& t)
931 {
932 // Find number of nodes
933 unsigned n_node = nnode();
934
935 // Local shape function
936 Shape psi(n_node);
937
938 // Vectors of local coordinates
939 Vector<double> s(DIM + 1, 0.0);
940
941 // Tecplot header info
942 outfile << tecplot_zone_string(n_plot);
943
944 // How many plot points are there?
945 unsigned num_plot_points = nplot_points(n_plot);
946
947 // Loop over plot points
948 for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
949 {
950 // Get the local coordinates of plot point
951 get_s_plot(i_plot, n_plot, s);
952
953 // Coordinates
954 for (unsigned i = 0; i < DIM + 1; i++)
955 {
956 // Output the i-th spatial coordinate value
957 outfile << interpolated_x(s, i) << " ";
958 }
959
960 // Velocities
961 for (unsigned i = 0; i < DIM; i++)
962 {
963 // Output the i-th velocity component
964 outfile << interpolated_u_nst(s, i) << " ";
965 }
966
967 // End the line
968 outfile << std::endl;
969 } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
970
971 // Write tecplot footer (e.g. FE connectivity lists)
972 write_tecplot_zone_footer(outfile, n_plot);
973 } // End of output_veloc
974
975 //==============================================================
976 /// Output function:
977 /// x,y,[z],u,v,[w],p
978 /// in tecplot format. Specified number of plot points in each
979 /// coordinate direction.
980 //==============================================================
981 template<unsigned DIM>
983 const unsigned& n_plot)
984 {
985 // Find number of nodes
986 unsigned n_node = nnode();
987
988 // Local shape function
989 Shape psi(n_node);
990
991 // Vectors of local coordinates
992 Vector<double> s(DIM + 1, 0.0);
993
994 // Tecplot header info
995 outfile << tecplot_zone_string(n_plot);
996
997 // How many plot points are there?
998 unsigned num_plot_points = nplot_points(n_plot);
999
1000 // Loop over plot points
1001 for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1002 {
1003 // Get the local coordinates of plot point
1004 get_s_plot(i_plot, n_plot, s);
1005
1006 // Coordinates
1007 for (unsigned i = 0; i < DIM + 1; i++)
1008 {
1009 // Output the i-th spatial coordinate value
1010 outfile << interpolated_x(s, i) << " ";
1011 }
1012
1013 // Velocities
1014 for (unsigned i = 0; i < DIM; i++)
1015 {
1016 // Output the i-th velocity component
1017 outfile << interpolated_u_nst(s, i) << " ";
1018 }
1019
1020 // Pressure
1021 outfile << interpolated_p_nst(s) << " ";
1022
1023 // End the line
1024 outfile << std::endl;
1025 } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
1026
1027 // Add an extra line
1028 outfile << std::endl;
1029
1030 // Write tecplot footer (e.g. FE connectivity lists)
1031 write_tecplot_zone_footer(outfile, n_plot);
1032 } // End of output
1033
1034 //==============================================================
1035 /// C-style output function:
1036 /// x,y,[z],u,v,[w],p
1037 /// in tecplot format. Specified number of plot points in each
1038 /// coordinate direction.
1039 //==============================================================
1040 template<unsigned DIM>
1042 const unsigned& n_plot)
1043 {
1044 // Vector of local coordinates
1045 Vector<double> s(DIM + 1, 0.0);
1046
1047 // Tecplot header info
1048 fprintf(file_pt, "%s", tecplot_zone_string(n_plot).c_str());
1049
1050 // How many plot points in total?
1051 unsigned num_plot_points = nplot_points(n_plot);
1052
1053 // Loop over plot points
1054 for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1055 {
1056 // Get the local coordinates associated with this plot point
1057 get_s_plot(i_plot, n_plot, s);
1058
1059 // Loop over the coordinate directions
1060 for (unsigned i = 0; i < DIM + 1; i++)
1061 {
1062 // Output the i-th coordinate value to file
1063 fprintf(file_pt, "%g ", interpolated_x(s, i));
1064 }
1065
1066 // Loop over the velocity components
1067 for (unsigned i = 0; i < DIM; i++)
1068 {
1069 // Output the i-th velocity component to file
1070 fprintf(file_pt, "%g ", interpolated_u_nst(s, i));
1071 }
1072
1073 // Pressure
1074 fprintf(file_pt, "%g \n", interpolated_p_nst(s));
1075 }
1076
1077 // Finish the line off
1078 fprintf(file_pt, "\n");
1079
1080 // Write tecplot footer (e.g. FE connectivity lists)
1081 write_tecplot_zone_footer(file_pt, n_plot);
1082 } // End of output
1083
1084 //==============================================================
1085 /// Full output function:
1086 /// x,y,t,u,v,p,du/dt,dv/dt,dissipation
1087 /// in tecplot format. Specified number of plot points in each
1088 /// coordinate direction
1089 //==============================================================
1090 template<unsigned DIM>
1092 const unsigned& n_plot)
1093 {
1094 // Vector of local coordinates
1095 Vector<double> s(DIM + 1, 0.0);
1096
1097 // Tecplot header info
1098 outfile << tecplot_zone_string(n_plot);
1099
1100 // Find out how many nodes there are
1101 unsigned n_node = nnode();
1102
1103 // Set up memory for the shape functions
1104 Shape psif(n_node);
1105
1106 // Set up memory for the shape function derivatives
1107 DShape dpsifdx(n_node, DIM);
1108
1109 // How many plot points in total?
1110 unsigned num_plot_points = nplot_points(n_plot);
1111
1112 // Loop over plot points
1113 for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1114 {
1115 // Get the local coordinates of plot point
1116 get_s_plot(i_plot, n_plot, s);
1117
1118 // Call the derivatives of the shape and test functions
1119 dshape_eulerian(s, psif, dpsifdx);
1120
1121 // Allocate storage for the mesh velocity
1122 Vector<double> mesh_velocity(DIM, 0.0);
1123
1124 // Allocate storage for the acceleration
1125 Vector<double> du_dt(DIM, 0.0);
1126
1127 // Allocate storage for the ALE acceleration
1128 Vector<double> du_dt_ALE(DIM, 0.0);
1129
1130 // Allocate storage for the velocity derivatives
1131 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1132
1133 //--------------------------------------
1134 // Calculate velocities and derivatives:
1135 //--------------------------------------
1136 // Loop over directions
1137 for (unsigned i = 0; i < DIM; i++)
1138 {
1139 // Get the index at which velocity i is stored
1140 unsigned u_nodal_index = u_index_nst(i);
1141
1142 // Loop over nodes
1143 for (unsigned l = 0; l < n_node; l++)
1144 {
1145 // Get the nodal value
1146 double u_value = nodal_value(l, u_nodal_index);
1147
1148 // Update the i-th acceleration component
1149 du_dt[i] += u_value * dpsifdx(l, DIM);
1150
1151 // Update the i-th mesh velocity component
1152 mesh_velocity[i] += nodal_position(l, i) * dpsifdx(l, DIM);
1153
1154 // Loop over derivative directions for velocity gradients
1155 for (unsigned j = 0; j < DIM; j++)
1156 {
1157 // Update the value of du_i/dx_j
1158 interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1159 }
1160 } // for (unsigned l=0;l<n_node;l++)
1161 } // for (unsigned i=0;i<DIM;i++)
1162
1163 //---------------------------------------------
1164 // Get du/dt in ALE form (incl. mesh velocity):
1165 //---------------------------------------------
1166 // Loop over the coordinate directions
1167 for (unsigned i = 0; i < DIM; i++)
1168 {
1169 // Store the i-th acceleration component
1170 du_dt_ALE[i] = du_dt[i];
1171
1172 // Loop over the coordinate directions
1173 for (unsigned k = 0; k < DIM; k++)
1174 {
1175 // Take the mesh velocity into account
1176 du_dt_ALE[i] -= mesh_velocity[k] * interpolated_dudx(i, k);
1177 }
1178 } // for (unsigned i=0;i<DIM;i++)
1179
1180 // Output the coordinates
1181 for (unsigned i = 0; i < DIM + 1; i++)
1182 {
1183 // Output the i-th coordinate value
1184 outfile << interpolated_x(s, i) << " ";
1185 }
1186
1187 // Output the velocity components
1188 for (unsigned i = 0; i < DIM; i++)
1189 {
1190 // Output the i-th velocity component
1191 outfile << interpolated_u_nst(s, i) << " ";
1192 }
1193
1194 // Output the pressure
1195 outfile << interpolated_p_nst(s) << " ";
1196
1197 // Output the acceleration
1198 for (unsigned i = 0; i < DIM; i++)
1199 {
1200 // Output the ALE acceleration term
1201 outfile << du_dt_ALE[i] << " ";
1202 }
1203
1204 // Dissipation
1205 outfile << dissipation(s) << " ";
1206
1207 // End the line
1208 outfile << std::endl;
1209 } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
1210
1211 // Write tecplot footer (e.g. FE connectivity lists)
1212 write_tecplot_zone_footer(outfile, n_plot);
1213 } // End of full_output
1214
1215
1216 //==============================================================
1217 /// Output function for vorticity.
1218 /// x,y,[z],[omega_x,omega_y,[and/or omega_z]]
1219 /// in tecplot format. Specified number of plot points in each
1220 /// coordinate direction.
1221 //==============================================================
1222 template<unsigned DIM>
1224 std::ostream& outfile, const unsigned& n_plot)
1225 {
1226 // Vector of local coordinates
1227 Vector<double> s(DIM + 1, 0.0);
1228
1229 // Create vorticity vector of the required size
1230 Vector<double> vorticity;
1231
1232 // If we're in 2D the vorticity field is a scalar field
1233 if (DIM == 2)
1234 {
1235 // Resize the vorticity vector
1236 vorticity.resize(1);
1237 }
1238 // If we're in 3D the vorticity field is a vector field
1239 else if (DIM == 3)
1240 {
1241 // Resize the vorticity vector
1242 vorticity.resize(3);
1243 }
1244 // If we're in 1D
1245 else
1246 {
1247 std::string error_message =
1248 "Can't output vorticity in 1D - in fact, what ";
1249 error_message += "do you mean\nby the 1D Navier-Stokes equations?\n";
1250 throw OomphLibError(
1251 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1252 }
1253
1254 // Tecplot header info
1255 outfile << tecplot_zone_string(n_plot);
1256
1257 // How many plot points are there in total?
1258 unsigned num_plot_points = nplot_points(n_plot);
1259
1260 // Loop over plot points
1261 for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1262 {
1263 // Get local coordinates of plot point
1264 get_s_plot(i_plot, n_plot, s);
1265
1266 // Coordinates
1267 for (unsigned i = 0; i < DIM + 1; i++)
1268 {
1269 // Output the i-th coordinate value
1270 outfile << interpolated_x(s, i) << " ";
1271 }
1272
1273 // Get vorticity
1274 get_vorticity(s, vorticity);
1275
1276 // If we're in 2D
1277 if (DIM == 2)
1278 {
1279 // Output the vorticity field value
1280 outfile << vorticity[0];
1281 }
1282 // If we're in 3D
1283 else
1284 {
1285 // Output the vorticity field
1286 outfile << vorticity[0] << " " << vorticity[1] << " " << vorticity[2]
1287 << " ";
1288 }
1289
1290 // Finish the line off
1291 outfile << std::endl;
1292 }
1293
1294 // Add in an extra line
1295 outfile << std::endl;
1296
1297 // Write tecplot footer (e.g. FE connectivity lists)
1298 write_tecplot_zone_footer(outfile, n_plot);
1299 } // End of output_vorticity
1300
1301 //==============================================================
1302 /// Return integral of dissipation over element
1303 //==============================================================
1304 template<unsigned DIM>
1306 {
1307 // Initialise the elemental dissipation value
1308 double diss = 0.0;
1309
1310 // Set the value of n_intpt
1311 unsigned n_intpt = integral_pt()->nweight();
1312
1313 // Set the Vector to hold local coordinates
1314 Vector<double> s(DIM + 1, 0.0);
1315
1316 // Loop over the integration points
1317 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1318 {
1319 // Assign values of s
1320 for (unsigned i = 0; i < DIM + 1; i++)
1321 {
1322 // Calculate the i-th local coordinate value
1323 s[i] = integral_pt()->knot(ipt, i);
1324 }
1325
1326 // Get the integral weight
1327 double w = integral_pt()->weight(ipt);
1328
1329 // Get Jacobian of mapping
1330 double J = J_eulerian(s);
1331
1332 // Storage for the strain rate matrix
1333 DenseMatrix<double> strainrate(DIM, DIM);
1334
1335 // Get strain rate matrix
1336 strain_rate(s, strainrate);
1337
1338 // Initialise the local dissipation
1339 double local_diss = 0.0;
1340
1341 // Loop over the coordinate directions
1342 for (unsigned i = 0; i < DIM; i++)
1343 {
1344 // Loop over the coordinate directions
1345 for (unsigned j = 0; j < DIM; j++)
1346 {
1347 // Update the local dissipation value
1348 local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
1349 }
1350 } // for (unsigned i=0;i<DIM;i++)
1351
1352 // Update the elemental dissipation value
1353 diss += local_diss * w * J;
1354 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
1355
1356 // Return the elemental dissipation value
1357 return diss;
1358 } // End of dissipation
1359
1360 //==============================================================
1361 /// Compute traction (on the viscous scale) exerted onto
1362 /// the fluid at local coordinate s. N has to be outer unit normal
1363 /// to the fluid.
1364 //==============================================================
1365 template<unsigned DIM>
1367 const Vector<double>& s, const Vector<double>& N, Vector<double>& traction)
1368 {
1369 // Allocate space for the strain rate matrix
1370 DenseMatrix<double> strainrate(DIM, DIM);
1371
1372 // Get velocity gradients
1373 strain_rate(s, strainrate);
1374
1375 // Get pressure
1376 double press = interpolated_p_nst(s);
1377
1378 // Loop over traction components
1379 for (unsigned i = 0; i < DIM; i++)
1380 {
1381 // Add in the pressure contribution
1382 traction[i] = -press * N[i];
1383
1384 // Loop over the strain rate entries
1385 for (unsigned j = 0; j < DIM; j++)
1386 {
1387 // Add in the strain rate contribution
1388 traction[i] += 2.0 * strainrate(i, j) * N[j];
1389 }
1390 } // for (unsigned i=0;i<DIM;i++)
1391 } // End of get_traction
1392
1393 //==============================================================
1394 /// Compute traction (on the viscous scale) exerted onto
1395 /// the fluid at local coordinate s, decomposed into pressure and
1396 /// normal and tangential viscous components.
1397 /// N has to be outer unit normal to the fluid.
1398 //==============================================================
1399 template<unsigned DIM>
1401 const Vector<double>& s,
1402 const Vector<double>& N,
1403 Vector<double>& traction_p,
1404 Vector<double>& traction_visc_n,
1405 Vector<double>& traction_visc_t)
1406 {
1407 // Allocate space for the traction components
1408 Vector<double> traction_visc(DIM);
1409
1410 // Allocate space for the velocity gradients
1411 DenseMatrix<double> strainrate(DIM, DIM);
1412
1413 // Calculate the strain rate at the local coordinate s
1414 strain_rate(s, strainrate);
1415
1416 // Get pressure
1417 double press = interpolated_p_nst(s);
1418
1419 // Loop over traction components
1420 for (unsigned i = 0; i < DIM; i++)
1421 {
1422 // Get the pressure contribution
1423 traction_p[i] = -press * N[i];
1424
1425 // Loop over the coordinate directions
1426 for (unsigned j = 0; j < DIM; j++)
1427 {
1428 // Add in the viscous stress contribution
1429 traction_visc[i] += 2.0 * strainrate(i, j) * N[j];
1430 }
1431
1432 // Get the normal component of the viscous stress
1433 traction_visc_n[i] = traction_visc[i] * N[i];
1434
1435 // Get the tangential component of the viscous stress
1436 traction_visc_t[i] = traction_visc[i] - traction_visc_n[i];
1437 } // for (unsigned i=0;i<DIM;i++)
1438 } // End of get_traction
1439
1440 //==============================================================
1441 /// Return dissipation at local coordinate s
1442 //==============================================================
1443 template<unsigned DIM>
1445 const Vector<double>& s) const
1446 {
1447 // Get strain rate matrix
1448 DenseMatrix<double> strainrate(DIM, DIM);
1449 strain_rate(s, strainrate);
1450
1451 // Initialise
1452 double local_diss = 0.0;
1453 for (unsigned i = 0; i < DIM; i++)
1454 {
1455 for (unsigned j = 0; j < DIM; j++)
1456 {
1457 local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
1458 }
1459 }
1460
1461 return local_diss;
1462 }
1463
1464 //==============================================================
1465 /// Get strain-rate tensor: (1/2)*(du_i/dx_j+du_j/dx_i)
1466 //==============================================================
1467 template<unsigned DIM>
1469 const Vector<double>& s, DenseMatrix<double>& strainrate) const
1470 {
1471#ifdef PARANOID
1472 if ((strainrate.ncol() != DIM) || (strainrate.nrow() != DIM))
1473 {
1474 std::ostringstream error_message;
1475 error_message << "The strain rate has incorrect dimensions "
1476 << strainrate.ncol() << " x " << strainrate.nrow()
1477 << " not " << DIM << std::endl;
1478 throw OomphLibError(
1479 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1480 }
1481#endif
1482
1483 // Velocity gradient matrix
1484 DenseMatrix<double> dudx(DIM);
1485
1486 // Find out how many nodes there are in the element
1487 unsigned n_node = nnode();
1488
1489 // Set up memory for the shape and test functions
1490 Shape psi(n_node);
1491
1492 // Set up memory for the shape and test function derivatives
1493 DShape dpsidx(n_node, DIM + 1);
1494
1495 // Evaluate the shape functions and shape function derivatives at this point
1496 dshape_eulerian(s, psi, dpsidx);
1497
1498 // Initialise to zero
1499 dudx.initialise(0.0);
1500
1501 // Loop over velocity components
1502 for (unsigned i = 0; i < DIM; i++)
1503 {
1504 // Get the index at which the i-th velocity is stored
1505 unsigned u_nodal_index = u_index_nst(i);
1506
1507 // Loop over derivative directions
1508 for (unsigned j = 0; j < DIM; j++)
1509 {
1510 // Loop over nodes
1511 for (unsigned l = 0; l < n_node; l++)
1512 {
1513 // Update the value of du_i/dx_j
1514 dudx(i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
1515 }
1516 } // for (unsigned j=0;j<DIM;j++)
1517 } // for (unsigned i=0;i<DIM;i++)
1518
1519 // Loop over velocity components
1520 for (unsigned i = 0; i < DIM; i++)
1521 {
1522 // Loop over derivative directions
1523 for (unsigned j = 0; j < DIM; j++)
1524 {
1525 // Calculate the (i,j)-th strain rate entry
1526 strainrate(i, j) = 0.5 * (dudx(i, j) + dudx(j, i));
1527 }
1528 } // for (unsigned i=0;i<DIM;i++)
1529 } // End of strain_rate
1530
1531 //==============================================================
1532 /// Compute 2D vorticity vector at local coordinate s (return in
1533 /// one and only component of vorticity vector
1534 //==============================================================
1535 template<>
1537 const Vector<double>& s, Vector<double>& vorticity) const
1538 {
1539#ifdef PARANOID
1540 if (vorticity.size() != 1)
1541 {
1542 std::ostringstream error_message;
1543 error_message << "Computation of vorticity in 2D requires a 1D vector\n"
1544 << "which contains the only non-zero component of the\n"
1545 << "vorticity vector. You've passed a vector of size "
1546 << vorticity.size() << std::endl;
1547
1548 throw OomphLibError(
1549 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1550 }
1551#endif
1552
1553 // Specify spatial dimension
1554 unsigned DIM = 2;
1555
1556 // Velocity gradient matrix
1557 DenseMatrix<double> dudx(DIM, DIM, 0.0);
1558
1559 // Find out how many nodes there are in the element
1560 unsigned n_node = nnode();
1561
1562 // Set up memory for the shape and test functions
1563 Shape psi(n_node);
1564
1565 // Set up memory for the shape and test function derivatives
1566 DShape dpsidx(n_node, DIM + 1);
1567
1568 // Evaluate the shape functions and shape function derivatives at this point
1569 dshape_eulerian(s, psi, dpsidx);
1570
1571 // Initialise to zero
1572 dudx.initialise(0.0);
1573
1574 // Loop over velocity components
1575 for (unsigned i = 0; i < DIM; i++)
1576 {
1577 // Get the index at which the i-th velocity is stored
1578 unsigned u_nodal_index = u_index_nst(i);
1579
1580 // Loop over derivative directions
1581 for (unsigned j = 0; j < DIM; j++)
1582 {
1583 // Loop over nodes
1584 for (unsigned l = 0; l < n_node; l++)
1585 {
1586 // Update the value of du_i/dx_j
1587 dudx(i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
1588 }
1589 } // for (unsigned j=0;j<DIM;j++)
1590 } // for (unsigned i=0;i<DIM;i++)
1591
1592 // Z-component of vorticity
1593 vorticity[0] = dudx(1, 0) - dudx(0, 1);
1594 } // End of get_vorticity
1595
1596 //==============================================================
1597 /// Compute 2D vorticity vector at local coordinate s (return in
1598 /// one and only component of vorticity vector as a double
1599 //==============================================================
1600 template<>
1602 double& vorticity) const
1603 {
1604 // Create a vector to store the vorticity
1605 Vector<double> vort(1, 0.0);
1606
1607 // Calculate the vorticity
1608 get_vorticity(s, vort);
1609
1610 // Assign the vorticity
1611 vorticity = vort[0];
1612 } // End of get_vorticity
1613
1614 //==============================================================
1615 /// Get integral of kinetic energy over element:
1616 /// Note that this is the "raw" kinetic energy in the sense
1617 /// that the density ratio has not been included. In problems
1618 /// with two or more fluids the user will have to remember
1619 /// to premultiply certain elements by the appropriate density
1620 /// ratio.
1621 //==============================================================
1622 template<unsigned DIM>
1624 {
1625 // Initialise the elemental kinetic energy value
1626 double kin_en = 0.0;
1627
1628 // Set the value of n_intpt
1629 unsigned n_intpt = integral_pt()->nweight();
1630
1631 // Set the Vector to hold local coordinates
1632 Vector<double> s(DIM + 1, 0.0);
1633
1634 // Loop over the integration points
1635 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1636 {
1637 // Assign values of s
1638 for (unsigned i = 0; i < DIM + 1; i++)
1639 {
1640 s[i] = integral_pt()->knot(ipt, i);
1641 }
1642
1643 // Get the integral weight
1644 double w = integral_pt()->weight(ipt);
1645
1646 // Get Jacobian of mapping
1647 double J = J_eulerian(s);
1648
1649 // Loop over directions
1650 double veloc_squared = 0.0;
1651 for (unsigned i = 0; i < DIM; i++)
1652 {
1653 veloc_squared += interpolated_u_nst(s, i) * interpolated_u_nst(s, i);
1654 }
1655
1656 kin_en += 0.5 * veloc_squared * w * J;
1657 }
1658
1659 return kin_en;
1660
1661 } // End of kin_energy
1662
1663
1664 //==========================================================================
1665 /// Get integral of time derivative of kinetic energy over element:
1666 //==========================================================================
1667 template<unsigned DIM>
1669 {
1670 // Initialise
1671 double d_kin_en_dt = 0.0;
1672
1673 // Set the value of n_intpt
1674 const unsigned n_intpt = integral_pt()->nweight();
1675
1676 // Get the number of nodes
1677 const unsigned n_node = this->nnode();
1678
1679 // Storage for the shape function
1680 Shape psi(n_node);
1681 DShape dpsidx(n_node, DIM + 1);
1682
1683 // Get the value at which the velocities are stored
1684 unsigned u_index[DIM];
1685 for (unsigned i = 0; i < DIM; i++)
1686 {
1687 u_index[i] = this->u_index_nst(i);
1688 }
1689
1690 // Loop over the integration points
1691 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1692 {
1693 // Get the jacobian of the mapping
1694 double J = dshape_eulerian_at_knot(ipt, psi, dpsidx);
1695
1696 // Get the integral weight
1697 double w = integral_pt()->weight(ipt);
1698
1699 // Now work out the velocity and the time derivative
1700 Vector<double> interpolated_u(DIM, 0.0);
1701 Vector<double> interpolated_dudt(DIM, 0.0);
1702
1703 // Loop over the shape functions
1704 for (unsigned l = 0; l < n_node; l++)
1705 {
1706 // Loop over the dimensions
1707 for (unsigned i = 0; i < DIM; i++)
1708 {
1709 interpolated_u[i] += nodal_value(l, u_index[i]) * psi(l);
1710 interpolated_dudt[i] += du_dt_nst(l, u_index[i]) * psi(l);
1711 }
1712 }
1713
1714 // Get mesh velocity bit
1715 if (!ALE_is_disabled)
1716 {
1717 Vector<double> mesh_velocity(DIM, 0.0);
1718 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1719
1720 // Loop over nodes again
1721 for (unsigned l = 0; l < n_node; l++)
1722 {
1723 for (unsigned i = 0; i < DIM; i++)
1724 {
1725 mesh_velocity[i] += this->dnodal_position_dt(l, i) * psi(l);
1726
1727 for (unsigned j = 0; j < DIM; j++)
1728 {
1729 interpolated_dudx(i, j) +=
1730 this->nodal_value(l, u_index[i]) * dpsidx(l, j);
1731 }
1732 }
1733 }
1734
1735 // Subtract mesh velocity from du_dt
1736 for (unsigned i = 0; i < DIM; i++)
1737 {
1738 for (unsigned k = 0; k < DIM; k++)
1739 {
1740 interpolated_dudt[i] -= mesh_velocity[k] * interpolated_dudx(i, k);
1741 }
1742 }
1743 }
1744
1745
1746 // Loop over directions and add up u du/dt terms
1747 double sum = 0.0;
1748 for (unsigned i = 0; i < DIM; i++)
1749 {
1750 sum += interpolated_u[i] * interpolated_dudt[i];
1751 }
1752
1753 d_kin_en_dt += sum * w * J;
1754 }
1755
1756 return d_kin_en_dt;
1757
1758 } // End of d_kin_energy_dt
1759
1760
1761 //==============================================================
1762 /// Return pressure integrated over the element
1763 //==============================================================
1764 template<unsigned DIM>
1766 {
1767 // Initialise
1768 double press_int = 0;
1769
1770 // Set the value of n_intpt
1771 unsigned n_intpt = integral_pt()->nweight();
1772
1773 // Set the Vector to hold local coordinates
1774 Vector<double> s(DIM + 1, 0.0);
1775
1776 // Loop over the integration points
1777 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1778 {
1779 // Assign values of s
1780 for (unsigned i = 0; i < DIM + 1; i++)
1781 {
1782 s[i] = integral_pt()->knot(ipt, i);
1783 }
1784
1785 // Get the integral weight
1786 double w = integral_pt()->weight(ipt);
1787
1788 // Get Jacobian of mapping
1789 double J = J_eulerian(s);
1790
1791 // Premultiply the weights and the Jacobian
1792 double W = w * J;
1793
1794 // Get pressure
1795 double press = interpolated_p_nst(s);
1796
1797 // Add
1798 press_int += press * W;
1799 }
1800
1801 return press_int;
1802 }
1803
1804
1805 //==============================================================
1806 /// Compute the residuals for the associated pressure advection
1807 /// diffusion problem. Used by the Fp preconditioner.
1808 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
1809 //==============================================================
1810 template<unsigned DIM>
1813 Vector<double>& residuals,
1814 DenseMatrix<double>& jacobian,
1815 const unsigned& flag)
1816 {
1817 OomphLibWarning("I'm not sure this is correct yet...",
1818 OOMPH_CURRENT_FUNCTION,
1819 OOMPH_EXCEPTION_LOCATION);
1820
1821 // Return immediately if there are no dofs
1822 if (ndof() == 0) return;
1823
1824 // Find out how many nodes there are
1825 unsigned n_node = nnode();
1826
1827 // Find out how many pressure dofs there are
1828 unsigned n_pres = npres_nst();
1829
1830 // Find the indices at which the local velocities are stored
1831 unsigned u_nodal_index[DIM];
1832 for (unsigned i = 0; i < DIM; i++)
1833 {
1834 u_nodal_index[i] = u_index_nst(i);
1835 }
1836
1837 // Set up memory for the velocity shape fcts
1838 Shape psif(n_node);
1839 DShape dpsidx(n_node, DIM + 1);
1840
1841 // Set up memory for pressure shape and test functions
1842 Shape psip(n_pres);
1843 Shape testp(n_pres);
1844 DShape dpsip(n_pres, DIM);
1845 DShape dtestp(n_pres, DIM);
1846
1847 // Number of integration points
1848 unsigned n_intpt = integral_pt()->nweight();
1849
1850 // Set the Vector to hold local coordinates
1851 Vector<double> s(DIM + 1, 0.0);
1852
1853 // Get Physical Variables from Element
1854 // Reynolds number must be multiplied by the density ratio
1855 double scaled_re = re() * density_ratio();
1856
1857 // Integers to store the local equations and unknowns
1858 int local_eqn = 0;
1859 int local_unknown = 0;
1860
1861 // Loop over the integration points
1862 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1863 {
1864 // Assign values of s
1865 for (unsigned i = 0; i < DIM + 1; i++)
1866 {
1867 s[i] = integral_pt()->knot(ipt, i);
1868 }
1869
1870 // Get the integral weight
1871 double w = integral_pt()->weight(ipt);
1872
1873 // Call the derivatives of the velocity shape functions
1874 // (Derivs not needed but they are free)
1875 double J = this->dshape_eulerian_at_knot(ipt, psif, dpsidx);
1876
1877 // Call the pressure shape and test functions
1878 this->dpshape_and_dptest_eulerian_nst(s, psip, dpsip, testp, dtestp);
1879
1880 // Premultiply the weights and the Jacobian
1881 double W = w * J;
1882
1883 // Storage for the (Eulerian) coordinates
1884 Vector<double> x(DIM + 1, 0.0);
1885
1886 // Calculate local values of the pressure and velocity components
1887 Vector<double> interpolated_u(DIM, 0.0);
1888 Vector<double> interpolated_dpdx(DIM, 0.0);
1889
1890 // Calculate pressure gradient
1891 for (unsigned l = 0; l < n_pres; l++)
1892 {
1893 for (unsigned i = 0; i < DIM; i++)
1894 {
1895 interpolated_dpdx[i] += p_nst(l) * dpsip(l, i);
1896 }
1897 }
1898
1899 // Loop over the velocity components
1900 for (unsigned i = 0; i < DIM; i++)
1901 {
1902 interpolated_u[i] = interpolated_u_nst(s, i);
1903 }
1904
1905 // Loop over coordinate directions
1906 for (unsigned i = 0; i < DIM + 1; i++)
1907 {
1908 x[i] = interpolated_x(s, i);
1909 }
1910
1911 // Source function (for validaton only)
1912 double source = 0.0;
1913 if (Press_adv_diff_source_fct_pt != 0)
1914 {
1915 source = Press_adv_diff_source_fct_pt(x);
1916 }
1917
1918 // Loop over the shape functions
1919 for (unsigned l = 0; l < n_pres; l++)
1920 {
1921 local_eqn = p_local_eqn(l);
1922
1923 // If not a boundary conditions
1924 if (local_eqn >= 0)
1925 {
1926 residuals[local_eqn] -= source * testp[l] * W;
1927
1928 for (unsigned k = 0; k < DIM; k++)
1929 {
1930 residuals[local_eqn] +=
1931 interpolated_dpdx[k] *
1932 (scaled_re * interpolated_u[k] * testp[l] + dtestp(l, k)) * W;
1933 }
1934
1935 // Jacobian too?
1936 if (flag)
1937 {
1938 // Loop over the shape functions
1939 for (unsigned l2 = 0; l2 < n_pres; l2++)
1940 {
1941 local_unknown = p_local_eqn(l2);
1942
1943 // If not a boundary conditions
1944 if (local_unknown >= 0)
1945 {
1946 if ((int(eqn_number(local_eqn)) != Pinned_fp_pressure_eqn) &&
1947 (int(eqn_number(local_unknown)) != Pinned_fp_pressure_eqn))
1948 {
1949 for (unsigned k = 0; k < DIM; k++)
1950 {
1951 jacobian(local_eqn, local_unknown) +=
1952 dtestp(l2, k) *
1953 (scaled_re * interpolated_u[k] * testp[l] +
1954 dtestp(l, k)) *
1955 W;
1956 }
1957 }
1958 else
1959 {
1960 if ((int(eqn_number(local_eqn)) == Pinned_fp_pressure_eqn) &&
1961 (int(eqn_number(local_unknown)) ==
1962 Pinned_fp_pressure_eqn))
1963 {
1964 //
1965 jacobian(local_eqn, local_unknown) = 1.0;
1966 }
1967 } // if
1968 // ((int(eqn_number(local_eqn))!=Pinned_fp_pressure_eqn)&&...
1969 } // if (local_unknown>=0)
1970 } // for (unsigned l2=0;l2<n_pres;l2++)
1971 } // if (flag)
1972 } // if (local_eqn>=0)
1973 } // for (unsigned l=0;l<n_pres;l++)
1974 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
1975
1976 // Now add boundary contributions from Robin BCs
1977 unsigned nrobin = Pressure_advection_diffusion_robin_element_pt.size();
1978 for (unsigned e = 0; e < nrobin; e++)
1979 {
1980 Pressure_advection_diffusion_robin_element_pt[e]
1981 ->fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(
1982 residuals, jacobian, flag);
1983 }
1984 } // End of fill_in_generic_pressure_advection_diffusion_contribution_nst
1985
1986
1987 //==============================================================
1988 /// Compute the residuals for the Navier-Stokes
1989 /// equations; flag=1 (or 0): do (or don't) compute the
1990 /// Jacobian as well.
1991 //==============================================================
1992 template<unsigned DIM>
1995 DenseMatrix<double>& jacobian,
1996 DenseMatrix<double>& mass_matrix,
1997 const unsigned& flag)
1998 {
1999 // Return immediately if there are no dofs
2000 if (ndof() == 0) return;
2001
2002 // Find out how many nodes there are
2003 unsigned n_node = nnode();
2004
2005 // Find out how many pressure dofs there are
2006 unsigned n_pres = npres_nst();
2007
2008 // Allocate storage for the indices of the velocity components
2009 unsigned u_nodal_index[DIM];
2010
2011 // Loop over the velocity components
2012 for (unsigned i = 0; i < DIM; i++)
2013 {
2014 // Find the index at which the i-th local velocity is stored
2015 u_nodal_index[i] = u_index_nst(i);
2016 }
2017
2018 // Set up memory for the shape functions
2019 Shape psif(n_node);
2020
2021 // Set up memory for the test functions
2022 Shape testf(n_node);
2023
2024 // Allocate space for the derivatives of the shape functions
2025 DShape dpsifdx(n_node, DIM + 1);
2026
2027 // Allocate space for the derivatives of the test functions
2028 DShape dtestfdx(n_node, DIM + 1);
2029
2030 // Set up memory for pressure shape functions
2031 Shape psip(n_pres);
2032
2033 // Set up memory for pressure test functions
2034 Shape testp(n_pres);
2035
2036 // Number of integration points
2037 unsigned n_intpt = integral_pt()->nweight();
2038
2039 // Set the Vector to hold local coordinates
2040 Vector<double> s(DIM + 1, 0.0);
2041
2042 //-------------------------------------
2043 // Get physical variables from element:
2044 //-------------------------------------
2045 // Reynolds number must be multiplied by the density ratio
2046 double scaled_re = re() * density_ratio();
2047
2048 // Get the scaled Strouhal value
2049 double scaled_re_st = re() * st() * density_ratio();
2050
2051 // Get the scaled Womersley number differentiated w.r.t. the Strouhal number
2052 double scaled_dre_st_dst = re() * density_ratio();
2053
2054 // Get the scaled Reynolds / Froude number
2055 double scaled_re_inv_fr = re_invfr() * density_ratio();
2056
2057 // Get the viscosity ratio
2058 double visc_ratio = viscosity_ratio();
2059
2060 // Get the gravity vector
2061 Vector<double> G = g();
2062
2063 // Integer to store the local equation number
2064 int local_eqn = 0;
2065
2066 // Integer to store the local unknown number
2067 int local_unknown = 0;
2068
2069 // Loop over the integration points
2070 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2071 {
2072 // Assign values of s
2073 for (unsigned i = 0; i < DIM + 1; i++)
2074 {
2075 // Calculate the i-th local coordinate
2076 s[i] = integral_pt()->knot(ipt, i);
2077 }
2078
2079 // Get the integral weight
2080 double w = integral_pt()->weight(ipt);
2081
2082 // Call the derivatives of the shape and test functions
2083 double J = dshape_and_dtest_eulerian_at_knot_nst(
2084 ipt, psif, dpsifdx, testf, dtestfdx);
2085
2086 // Call the pressure shape and test functions
2087 pshape_nst(s, psip, testp);
2088
2089 // Pre-multiply the weights and the Jacobian
2090 double W = w * J;
2091
2092 // Storage for the interpolated time value
2093 double interpolated_t = 0.0;
2094
2095 // Storage for the interpolated pressure value
2096 double interpolated_p = 0.0;
2097
2098 // Storage for the spatial coordinates
2099 Vector<double> interpolated_x(DIM, 0.0);
2100
2101 // Storage for the interpolated velocity components
2102 Vector<double> interpolated_u(DIM, 0.0);
2103
2104 // Storage for the interpolated time-derivative of the velocities
2105 Vector<double> interpolated_dudt(DIM, 0.0);
2106
2107 // Storage for the mesh velocity
2108 Vector<double> mesh_velocity(DIM, 0.0);
2109
2110 // Storage for the spatial derivatives of the velocity components
2111 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
2112
2113 // Calculate pressure
2114 for (unsigned l = 0; l < n_pres; l++)
2115 {
2116 // Update the current approximation to the interpolated pressure
2117 interpolated_p += p_nst(l) * psip[l];
2118 }
2119
2120 //-------------------------------------------------------------------
2121 // Calculate velocities, derivatives, source function and body force:
2122 //-------------------------------------------------------------------
2123 // Loop over nodes
2124 for (unsigned l = 0; l < n_node; l++)
2125 {
2126 // Update the interpolated time value
2127 interpolated_t += raw_nodal_position(l, DIM) * psif(l);
2128
2129 // Loop over coordinate directions
2130 for (unsigned i = 0; i < DIM; i++)
2131 {
2132 // Get the nodal value
2133 double u_value = raw_nodal_value(l, u_nodal_index[i]);
2134
2135 // Update the i-th interpolated velocity component
2136 interpolated_u[i] += u_value * psif[l];
2137
2138 // Update the i-th interpolated coordinate value
2139 interpolated_x[i] += raw_nodal_position(l, i) * psif[l];
2140
2141 // Update the interpolated du_i/dt value
2142 interpolated_dudt[i] += u_value * dpsifdx(l, DIM);
2143
2144 // Loop over derivative directions
2145 for (unsigned j = 0; j < DIM; j++)
2146 {
2147 // Update the interpolated du_i/dx_j value
2148 interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
2149 }
2150 } // for (unsigned i=0;i<DIM;i++)
2151 } // for (unsigned l=0;l<n_node;l++)
2152
2153 // If ALE is enabled
2154 if (!ALE_is_disabled)
2155 {
2156 // Loop over nodes
2157 for (unsigned l = 0; l < n_node; l++)
2158 {
2159 // Loop over directions
2160 for (unsigned i = 0; i < DIM; i++)
2161 {
2162 // DRAIG: Check to see if the velocity is actually consistent
2163 // with the imposed mesh velocity...
2164 // Update the i-th mesh velocity component
2165 mesh_velocity[i] += raw_nodal_position(l, i) * dpsifdx(l, DIM);
2166 }
2167 } // for (unsigned l=0;l<n_node;l++)
2168 } // if (!ALE_is_disabled)
2169
2170 // Allocate space for the body force
2171 Vector<double> body_force(DIM, 0.0);
2172
2173 // Get the user-defined body force term
2174 get_body_force_nst(interpolated_t, ipt, s, interpolated_x, body_force);
2175
2176 // Get the user-defined source function
2177 double source = get_source_nst(interpolated_t, ipt, interpolated_x);
2178
2179 //---------------------------------
2180 // Assemble residuals and Jacobian:
2181 //---------------------------------
2182 //--------------------
2183 // Momentum equations:
2184 //--------------------
2185 // Loop over the test functions
2186 for (unsigned l = 0; l < n_node; l++)
2187 {
2188 // Loop over the velocity components
2189 for (unsigned i = 0; i < DIM; i++)
2190 {
2191 // Get the local equation number associated with this unknown and node
2192 local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
2193
2194 // If it's not a boundary condition
2195 if (local_eqn >= 0)
2196 {
2197 // Add the user-defined body force terms
2198 residuals[local_eqn] += body_force[i] * testf[l] * W;
2199
2200 // Add the gravitational body force term
2201 residuals[local_eqn] += scaled_re_inv_fr * testf[l] * G[i] * W;
2202
2203 // Add the pressure gradient term
2204 residuals[local_eqn] += interpolated_p * dtestfdx(l, i) * W;
2205
2206 // Add in the contribution from the time derivative
2207 residuals[local_eqn] -=
2208 scaled_re_st * interpolated_dudt[i] * testf[l] * W;
2209
2210 // If ALE is enabled
2211 if (!ALE_is_disabled)
2212 {
2213 // Loop over the coordinate directions
2214 for (unsigned k = 0; k < DIM; k++)
2215 {
2216 // Add in the mesh velocity contribution
2217 residuals[local_eqn] +=
2218 (scaled_re_st * mesh_velocity[k] * interpolated_dudx(i, k) *
2219 testf[l] * W);
2220 }
2221 } // if (!ALE_is_disabled)
2222
2223 // Loop over the coordinate directions
2224 for (unsigned k = 0; k < DIM; k++)
2225 {
2226 // Add in the convective term contribution
2227 residuals[local_eqn] -= (scaled_re * interpolated_u[k] *
2228 interpolated_dudx(i, k) * testf[l] * W);
2229 }
2230
2231 // Loop over the coordinate directions
2232 for (unsigned k = 0; k < DIM; k++)
2233 {
2234 // Add in the stress tensor terms:
2235 // NOTE: The viscosity ratio needs to go in here to ensure
2236 // continuity of normal stress is satisfied even in flows
2237 // with zero pressure gradient!
2238 residuals[local_eqn] -= ((interpolated_dudx(i, k) +
2239 Gamma[i] * interpolated_dudx(k, i)) *
2240 visc_ratio * dtestfdx(l, k) * W);
2241 }
2242
2243 //------------------------
2244 // Calculate the Jacobian:
2245 //------------------------
2246 // If we also need to construct the Jacobian
2247 if (flag)
2248 {
2249 // Loop over the velocity shape functions again
2250 for (unsigned l2 = 0; l2 < n_node; l2++)
2251 {
2252 // Loop over the velocity components again
2253 for (unsigned i2 = 0; i2 < DIM; i2++)
2254 {
2255 // Get the local equation number associated with this unknown
2256 local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
2257
2258 // If we're at a proper degree of freedom
2259 if (local_unknown >= 0)
2260 {
2261 // Add the contribution to the elemental matrix
2262 jacobian(local_eqn, local_unknown) -=
2263 (visc_ratio * Gamma[i] * dpsifdx(l2, i) *
2264 dtestfdx(l, i2) * W);
2265
2266 // Now add in the inertial terms
2267 jacobian(local_eqn, local_unknown) -=
2268 (scaled_re * psif[l2] * interpolated_dudx(i, i2) *
2269 testf[l] * W);
2270
2271 // Extra component if i2=i
2272 if (i2 == i)
2273 {
2274 // If we also need to construct the mass matrix (only
2275 // diagonal entries)
2276 if (flag == 2)
2277 {
2278 // NOTE: This is positive because the mass matrix is
2279 // taken to the other side of the equation when
2280 // formulating the generalised eigenproblem.
2281 mass_matrix(local_eqn, local_unknown) +=
2282 (scaled_re_st * psif[l2] * testf[l] * W);
2283 }
2284
2285 // Add in the time-derivative contribution
2286 jacobian(local_eqn, local_unknown) -=
2287 (scaled_re_st * dpsifdx(l2, DIM) * testf[l] * W);
2288
2289 // If ALE is enabled
2290 if (!ALE_is_disabled)
2291 {
2292 // Loop over the velocity components
2293 for (unsigned k = 0; k < DIM; k++)
2295 // Add in the mesh velocity contribution
2296 jacobian(local_eqn, local_unknown) +=
2297 (scaled_re_st * mesh_velocity[k] * dpsifdx(l2, k) *
2298 testf[l] * W);
2299 }
2300 } // if (!ALE_is_disabled)
2301
2302 // Loop over the velocity components
2303 for (unsigned k = 0; k < DIM; k++)
2304 {
2305 // Add in the convective term contribution
2306 jacobian(local_eqn, local_unknown) -=
2307 (scaled_re * interpolated_u[k] * dpsifdx(l2, k) *
2308 testf[l] * W);
2309 }
2310
2311 // Loop over the velocity components
2312 for (unsigned k = 0; k < DIM; k++)
2313 {
2314 // Add in the velocity gradient terms
2315 jacobian(local_eqn, local_unknown) -=
2316 (visc_ratio * dpsifdx(l2, k) * dtestfdx(l, k) * W);
2317 }
2318 } // if (i2==i)
2319 } // if (local_unknown>=0)
2320 } // for (unsigned i2=0;i2<DIM;i2++)
2321 } // for (unsigned l2=0;l2<n_node;l2++)
2322
2323 // Loop over pressure shape functions
2324 for (unsigned l2 = 0; l2 < n_pres; l2++)
2325 {
2326 // Get the local equation number associated with this degree of
2327 // freedom
2328 local_unknown = p_local_eqn(l2);
2329
2330 // If we are at a proper degree of freedom
2331 if (local_unknown >= 0)
2332 {
2333 // Add in the pressure gradient contribution
2334 jacobian(local_eqn, local_unknown) +=
2335 psip[l2] * dtestfdx(l, i) * W;
2336 }
2337 } // for (unsigned l2=0;l2<n_pres;l2++)
2338
2339 //------------------------------------
2340 // Calculate external data information
2341 //------------------------------------
2342 // If we're storing the Strouhal number as external data then we
2343 // need to calculate dr/d(St) in the elemental Jacobian
2344 if (Strouhal_is_stored_as_external_data)
2345 {
2346 // The index of the external data (there's only one!)
2347 unsigned data_index = 0;
2348
2349 // The index of the unknown value stored in the external data
2350 unsigned value_index = 0;
2351
2352 // Get the local equation number associated with the extra
2353 // unknown
2354 local_unknown =
2355 this->external_local_eqn(data_index, value_index);
2356
2357 // If we're at a non-zero degree of freedom add in the entry
2358 if (local_unknown >= 0)
2359 {
2360 // Add in the contribution from the time derivative
2361 jacobian(local_eqn, local_unknown) -=
2362 (scaled_dre_st_dst * interpolated_dudt[i] * testf[l] * W);
2363
2364 // If ALE is enabled
2365 if (!ALE_is_disabled)
2366 {
2367 // Loop over the coordinate directions
2368 for (unsigned k = 0; k < DIM; k++)
2369 {
2370 // Add in the mesh velocity contribution
2371 jacobian(local_eqn, local_unknown) +=
2372 (scaled_dre_st_dst * mesh_velocity[k] *
2373 interpolated_dudx(i, k) * testf[l] * W);
2374 }
2375 } // if (!ALE_is_disabled)
2376 } // if (local_unknown>=0)
2377 } // if (Strouhal_is_stored_as_external_data)
2378 } // if (flag)
2379 } // if (local_eqn>=0)
2380 } // for (unsigned i=0;i<DIM;i++)
2381 } // for (unsigned l=0;l<n_node;l++)
2382
2383 //---------------------
2384 // Continuity equation:
2385 //---------------------
2386 // Loop over the shape functions
2387 for (unsigned l = 0; l < n_pres; l++)
2388 {
2389 // Get the local equation number associated with the pressure dof
2390 local_eqn = p_local_eqn(l);
2391
2392 // If it's not a boundary condition
2393 if (local_eqn >= 0)
2394 {
2395 // Add in the source term contribution
2396 residuals[local_eqn] -= source * testp[l] * W;
2397
2398 // Loop over velocity components
2399 for (unsigned k = 0; k < DIM; k++)
2400 {
2401 // Add in the velocity gradient terms
2402 residuals[local_eqn] += interpolated_dudx(k, k) * testp[l] * W;
2403 }
2404
2405 //------------------------
2406 // Calculate the Jacobian:
2407 //------------------------
2408 // If we also need to construct the Jacobian
2409 if (flag)
2410 {
2411 // Loop over the velocity shape functions
2412 for (unsigned l2 = 0; l2 < n_node; l2++)
2413 {
2414 // Loop over velocity components
2415 for (unsigned i2 = 0; i2 < DIM; i2++)
2416 {
2417 // Get the local equation number associated with this node
2418 local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
2419
2420 // If we're at a proper degree of freedom
2421 if (local_unknown >= 0)
2422 {
2423 // Add in the velocity gradient contribution
2424 jacobian(local_eqn, local_unknown) +=
2425 dpsifdx(l2, i2) * testp[l] * W;
2426 }
2427 } // for (unsigned i2=0;i2<DIM;i2++)
2428 } // for (unsigned l2=0;l2<n_node;l2++)
2429 } // if (flag)
2430 } // if (local_eqn>=0)
2431 } // for (unsigned l=0;l<n_pres;l++)
2432 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
2433 } // End of fill_in_generic_residual_contribution_nst
2434
2435 //================================================================
2436 /// Compute the derivatives of the residuals for the Navier-Stokes
2437 /// equations with respect to a parameter;
2438 /// flag=2 or 1 or 0: do (or don't) compute the Jacobian and
2439 /// mass matrix as well
2440 //================================================================
2441 template<unsigned DIM>
2444 double* const& parameter_pt,
2445 Vector<double>& dres_dparam,
2446 DenseMatrix<double>& djac_dparam,
2447 DenseMatrix<double>& dmass_matrix_dparam,
2448 const unsigned& flag)
2449 {
2450 // Throw an error
2451 throw OomphLibError("Not yet implemented\n",
2452 OOMPH_CURRENT_FUNCTION,
2453 OOMPH_EXCEPTION_LOCATION);
2454 } // End of fill_in_generic_dresidual_contribution_nst
2455
2456 //==================================================================
2457 /// Compute the hessian tensor vector products required to
2458 /// perform continuation of bifurcations analytically
2459 //==================================================================
2460 template<unsigned DIM>
2463 Vector<double> const& Y,
2464 DenseMatrix<double> const& C,
2465 DenseMatrix<double>& product)
2466 {
2467 // Throw an error
2468 throw OomphLibError("Not yet implemented\n",
2469 OOMPH_CURRENT_FUNCTION,
2470 OOMPH_EXCEPTION_LOCATION);
2471 } // End of fill_in_contribution_to_hessian_vector_products
2472
2473 //======================================================================
2474 /// Compute derivatives of elemental residual vector with respect
2475 /// to nodal coordinates.
2476 /// dresidual_dnodal_coordinates(l,i,j)=d res(l) / dX_{ij}
2477 /// Overloads the FD-based version in the FE base class.
2478 /// DRAIG: This needs doing carefully if the ALE nodes aren't fixed!!!
2479 //======================================================================
2480 template<unsigned DIM>
2482 RankThreeTensor<double>& dresidual_dnodal_coordinates)
2483 {
2484 // Throw a warning
2485 throw OomphLibError("Space-time update needs to be checked!",
2486 OOMPH_CURRENT_FUNCTION,
2487 OOMPH_EXCEPTION_LOCATION);
2488
2489 // Return immediately if there are no dofs
2490 if (ndof() == 0)
2491 {
2492 return;
2493 }
2494
2495 // Determine number of nodes in element
2496 const unsigned n_node = nnode();
2497
2498 // Determine number of pressure dofs in element
2499 const unsigned n_pres = npres_nst();
2500
2501 // Find the indices at which the local velocities are stored
2502 unsigned u_nodal_index[DIM];
2503 for (unsigned i = 0; i < DIM; i++)
2504 {
2505 u_nodal_index[i] = u_index_nst(i);
2506 }
2507
2508 // Set up memory for the shape and test functions
2509 Shape psif(n_node);
2510 Shape testf(n_node);
2511 DShape dpsifdx(n_node, DIM + 1);
2512 DShape dtestfdx(n_node, DIM + 1);
2513
2514 // Set up memory for pressure shape and test functions
2515 Shape psip(n_pres), testp(n_pres);
2516
2517 // Deriatives of shape fct derivatives w.r.t. nodal coords
2518 RankFourTensor<double> d_dpsifdx_dX(DIM, n_node, n_node, DIM);
2519 RankFourTensor<double> d_dtestfdx_dX(DIM, n_node, n_node, DIM);
2520
2521 // Derivative of Jacobian of mapping w.r.t. to nodal coords
2522 DenseMatrix<double> dJ_dX(DIM, n_node);
2523
2524 // Derivatives of derivative of u w.r.t. nodal coords
2525 RankFourTensor<double> d_dudx_dX(DIM, n_node, DIM, DIM);
2526
2527 // Derivatives of nodal velocities w.r.t. nodal coords:
2528 // Assumption: Interaction only local via no-slip so
2529 // X_ij only affects U_ij.
2530 DenseMatrix<double> d_U_dX(DIM, n_node, 0.0);
2531
2532 // Determine the number of integration points
2533 const unsigned n_intpt = integral_pt()->nweight();
2534
2535 // Vector to hold local coordinates
2536 Vector<double> s(DIM + 1, 0.0);
2537
2538 // Get physical variables from element
2539 // (Reynolds number must be multiplied by the density ratio)
2540 double scaled_re = re() * density_ratio();
2541 double scaled_re_st = re_st() * density_ratio();
2542 double scaled_re_inv_fr = re_invfr() * density_ratio();
2543 double visc_ratio = viscosity_ratio();
2544 Vector<double> G = g();
2545
2546 // FD step
2548
2549 // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
2550 // Assumption: Interaction only local via no-slip so
2551 // X_ij only affects U_ij.
2552 bool element_has_node_with_aux_node_update_fct = false;
2553 for (unsigned q = 0; q < n_node; q++)
2554 {
2555 Node* nod_pt = node_pt(q);
2556
2557 // Only compute if there's a node-update fct involved
2559 {
2560 element_has_node_with_aux_node_update_fct = true;
2561
2562 // Current nodal velocity
2563 Vector<double> u_ref(DIM);
2564 for (unsigned i = 0; i < DIM; i++)
2565 {
2566 u_ref[i] = *(nod_pt->value_pt(u_nodal_index[i]));
2567 }
2568
2569 // FD
2570 for (unsigned p = 0; p < DIM; p++)
2571 {
2572 // Make backup
2573 double backup = nod_pt->x(p);
2574
2575 // Do FD step. No node update required as we're
2576 // attacking the coordinate directly...
2577 nod_pt->x(p) += eps_fd;
2578
2579 // Do auxiliary node update (to apply no-slip)
2581
2582 // Evaluate
2583 d_U_dX(p, q) =
2584 (*(nod_pt->value_pt(u_nodal_index[p])) - u_ref[p]) / eps_fd;
2585
2586 // Reset
2587 nod_pt->x(p) = backup;
2588
2589 // Do auxiliary node update (to apply no slip)
2591 }
2592 }
2593 }
2594
2595 // Integer to store the local equation number
2596 int local_eqn = 0;
2597
2598 // Loop over the integration points
2599 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2600 {
2601 // Assign values of s
2602 for (unsigned i = 0; i < DIM + 1; i++)
2603 {
2604 s[i] = integral_pt()->knot(ipt, i);
2605 }
2606
2607 // Get the integral weight
2608 const double w = integral_pt()->weight(ipt);
2609
2610 // Call the derivatives of the shape and test functions
2611 const double J = dshape_and_dtest_eulerian_at_knot_nst(ipt,
2612 psif,
2613 dpsifdx,
2614 d_dpsifdx_dX,
2615 testf,
2616 dtestfdx,
2617 d_dtestfdx_dX,
2618 dJ_dX);
2619
2620 // Call the pressure shape and test functions
2621 pshape_nst(s, psip, testp);
2622
2623 // Storage for the interpolated time value
2624 double interpolated_t = 0.0;
2625
2626 // Storage for the interpolated pressure value
2627 double interpolated_p = 0.0;
2628
2629 // Storage for the spatial coordinates
2630 Vector<double> interpolated_x(DIM, 0.0);
2631
2632 // Storage for the interpolated velocity components
2633 Vector<double> interpolated_u(DIM, 0.0);
2634
2635 // Storage for the interpolated time-derivative of the velocities
2636 Vector<double> interpolated_dudt(DIM, 0.0);
2637
2638 // Storage for the mesh velocity
2639 Vector<double> mesh_velocity(DIM, 0.0);
2640
2641 // Storage for the spatial derivatives of the velocity components
2642 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
2643
2644 // Calculate pressure
2645 for (unsigned l = 0; l < n_pres; l++)
2646 {
2647 // Update the current approximation to the interpolated pressure
2648 interpolated_p += p_nst(l) * psip[l];
2649 }
2650
2651 //-------------------------------------------------------------------
2652 // Calculate velocities, derivatives, source function and body force:
2653 //-------------------------------------------------------------------
2654 // Loop over nodes
2655 for (unsigned l = 0; l < n_node; l++)
2656 {
2657 // Update the interpolated time value
2658 interpolated_t += raw_nodal_position(l, DIM) * psif(l);
2659
2660 // Loop over coordinate directions
2661 for (unsigned i = 0; i < DIM; i++)
2662 {
2663 // Get the nodal value
2664 double u_value = raw_nodal_value(l, u_nodal_index[i]);
2665
2666 // Update the i-th interpolated velocity component
2667 interpolated_u[i] += u_value * psif[l];
2668
2669 // Update the i-th interpolated coordinate value
2670 interpolated_x[i] += raw_nodal_position(l, i) * psif[l];
2671
2672 // Update the interpolated du_i/dt value
2673 interpolated_dudt[i] += u_value * dpsifdx(l, DIM);
2674
2675 // Loop over derivative directions
2676 for (unsigned j = 0; j < DIM; j++)
2677 {
2678 // Update the interpolated du_i/dx_j value
2679 interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
2680 }
2681 } // for (unsigned i=0;i<DIM;i++)
2682 } // for (unsigned l=0;l<n_node;l++)
2683
2684 // If ALE is enabled
2685 if (!ALE_is_disabled)
2686 {
2687 // Loop over nodes
2688 for (unsigned l = 0; l < n_node; l++)
2689 {
2690 // Loop over directions
2691 for (unsigned i = 0; i < DIM; i++)
2692 {
2693 // Update the i-th mesh velocity component
2694 mesh_velocity[i] += raw_nodal_position(l, i) * dpsifdx(l, DIM);
2695 }
2696 } // for (unsigned l=0;l<n_node;l++)
2697 } // if (!ALE_is_disabled)
2698
2699 // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
2700 // DRAIG: CHECK
2701 for (unsigned q = 0; q < n_node; q++)
2702 {
2703 // Loop over coordinate directions
2704 for (unsigned p = 0; p < DIM; p++)
2705 {
2706 for (unsigned i = 0; i < DIM; i++)
2707 {
2708 for (unsigned k = 0; k < DIM; k++)
2709 {
2710 double aux = 0.0;
2711
2712 for (unsigned j = 0; j < n_node; j++)
2713 {
2714 aux += raw_nodal_value(j, u_nodal_index[i]) *
2715 d_dpsifdx_dX(p, q, j, k);
2716 }
2717
2718 //
2719 d_dudx_dX(p, q, i, k) = aux;
2720 }
2721 }
2722 }
2723 }
2724
2725 // Allocate space for the body force
2726 Vector<double> body_force(DIM);
2727
2728 // Get the user-defined body force term
2729 get_body_force_nst(interpolated_t, ipt, s, interpolated_x, body_force);
2730
2731 // Get the user-defined source function
2732 const double source = get_source_nst(interpolated_t, ipt, interpolated_x);
2733
2734 // Note: Can use raw values and avoid bypassing hanging information
2735 // because this is the non-refineable version!
2736
2737 // Allocate space for the gradient of the body force function
2738 DenseMatrix<double> d_body_force_dx(DIM, DIM, 0.0);
2739
2740 // Get gradient of body force function
2741 get_body_force_gradient_nst(
2742 interpolated_t, ipt, s, interpolated_x, d_body_force_dx);
2743
2744 // Allocate space for the gradient of the source function
2745 Vector<double> source_gradient(DIM, 0.0);
2746
2747 // Get gradient of source function
2748 get_source_gradient_nst(
2749 interpolated_t, ipt, interpolated_x, source_gradient);
2750
2751 // Get weight of actual nodal position in computation of mesh velocity
2752 // from positional time stepper
2753 const double pos_time_weight =
2754 (node_pt(0)->position_time_stepper_pt()->weight(1, 0));
2755
2756 // Get weight of actual nodal value in computation of mesh velocity
2757 // from value time stepper
2758 const double val_time_weight =
2759 node_pt(0)->time_stepper_pt()->weight(1, 0);
2760
2761 //---------------------------------
2762 // Assemble residuals and Jacobian:
2763 //---------------------------------
2764 //--------------------
2765 // Momentum equations:
2766 //--------------------
2767 // Loop over the test functions
2768 for (unsigned l = 0; l < n_node; l++)
2769 {
2770 // Loop over coordinate directions
2771 for (unsigned i = 0; i < DIM; i++)
2772 {
2773 // Get the local equation number associated with this degree of
2774 // freedom
2775 local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
2776
2777 // If it's not a boundary condition
2778 if (local_eqn >= 0)
2779 {
2780 // Loop over coordinate directions
2781 for (unsigned p = 0; p < DIM; p++)
2782 {
2783 // Loop over nodes
2784 for (unsigned q = 0; q < n_node; q++)
2785 {
2786 //-----------------------------------
2787 // Residual x derivative of Jacobian:
2788 //-----------------------------------
2789 // Add the user-defined body force terms
2790 double sum = body_force[i] * testf[l];
2791
2792 // Add the gravitational body force term
2793 sum += scaled_re_inv_fr * testf[l] * G[i];
2794
2795 // Add the pressure gradient term
2796 sum += interpolated_p * dtestfdx(l, i);
2797
2798 // Loop over the coordinate directions
2799 for (unsigned k = 0; k < DIM; k++)
2800 {
2801 // Add in the stress tensor contribution
2802 // NOTE: The viscosity ratio needs to go in here to ensure
2803 // continuity of normal stress is satisfied even in flows
2804 // with zero pressure gradient!
2805 sum -= (visc_ratio *
2806 (interpolated_dudx(i, k) +
2807 Gamma[i] * interpolated_dudx(k, i)) *
2808 dtestfdx(l, k));
2809 }
2810
2811 // Add in the contribution from the time derivative
2812 sum -= scaled_re_st * interpolated_dudt[i] * testf[l];
2813
2814 // Loop over the coordinate directions
2815 for (unsigned k = 0; k < DIM; k++)
2816 {
2817 // Add in convective term contribution
2818 sum -= scaled_re * interpolated_u[k] *
2819 interpolated_dudx(i, k) * testf[l];
2820 }
2821
2822 // If ALE is enabled
2823 if (!ALE_is_disabled)
2824 {
2825 // Loop over the coordinate directions
2826 for (unsigned k = 0; k < DIM; k++)
2827 {
2828 // Add in the mesh velocity contribution
2829 sum += scaled_re_st * mesh_velocity[k] *
2830 interpolated_dudx(i, k) * testf[l];
2831 }
2832 }
2833
2834 // Multiply through by derivative of Jacobian and integration
2835 // weight
2836 dresidual_dnodal_coordinates(local_eqn, p, q) +=
2837 sum * dJ_dX(p, q) * w;
2838
2839 //----------------------------------
2840 // Derivative of residual x Jacobian
2841 //----------------------------------
2842 // Body force
2843 sum = d_body_force_dx(i, p) * psif(q) * testf(l);
2844
2845 // Pressure gradient term
2846 sum += interpolated_p * d_dtestfdx_dX(p, q, l, i);
2847
2848 // Viscous term
2849 for (unsigned k = 0; k < DIM; k++)
2850 {
2851 sum -= (visc_ratio * ((interpolated_dudx(i, k) +
2852 Gamma[i] * interpolated_dudx(k, i)) *
2853 d_dtestfdx_dX(p, q, l, k) +
2854 (d_dudx_dX(p, q, i, k) +
2855 Gamma[i] * d_dudx_dX(p, q, k, i)) *
2856 dtestfdx(l, k)));
2857 }
2858
2859 // Convective terms, including mesh velocity
2860 for (unsigned k = 0; k < DIM; k++)
2861 {
2862 double tmp = scaled_re * interpolated_u[k];
2863 if (!ALE_is_disabled)
2864 {
2865 tmp -= scaled_re_st * mesh_velocity[k];
2866 }
2867 sum -= tmp * d_dudx_dX(p, q, i, k) * testf(l);
2868 }
2869
2870 if (!ALE_is_disabled)
2871 {
2872 sum += scaled_re_st * pos_time_weight * psif(q) *
2873 interpolated_dudx(i, p) * testf(l);
2874 }
2875
2876 // Multiply through by Jacobian and integration weight
2877 dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2878
2879 // Derivs w.r.t. to nodal velocities
2880 // ---------------------------------
2881 if (element_has_node_with_aux_node_update_fct)
2882 {
2883 sum =
2884 -visc_ratio * Gamma[i] * dpsifdx(q, i) * dtestfdx(l, p) -
2885 scaled_re * psif(q) * interpolated_dudx(i, p) * testf(l);
2886 if (i == p)
2887 {
2888 sum -= scaled_re_st * val_time_weight * psif(q) * testf(l);
2889 for (unsigned k = 0; k < DIM; k++)
2890 {
2891 sum -= visc_ratio * dpsifdx(q, k) * dtestfdx(l, k);
2892 double tmp = scaled_re * interpolated_u[k];
2893 if (!ALE_is_disabled)
2894 tmp -= scaled_re_st * mesh_velocity[k];
2895 sum -= tmp * dpsifdx(q, k) * testf(l);
2896 }
2897 }
2898 dresidual_dnodal_coordinates(local_eqn, p, q) +=
2899 sum * d_U_dX(p, q) * J * w;
2900 }
2901 }
2902 }
2903 }
2904 }
2905 } // End of loop over test functions
2906
2907
2908 // CONTINUITY EQUATION
2909 // -------------------
2910
2911 // Loop over the shape functions
2912 for (unsigned l = 0; l < n_pres; l++)
2913 {
2914 local_eqn = p_local_eqn(l);
2915
2916 // If not a boundary conditions
2917 if (local_eqn >= 0)
2918 {
2919 // Loop over coordinate directions
2920 for (unsigned p = 0; p < DIM; p++)
2921 {
2922 // Loop over nodes
2923 for (unsigned q = 0; q < n_node; q++)
2924 {
2925 // Residual x deriv of Jacobian
2926 //-----------------------------
2927
2928 // Source term
2929 double aux = -source;
2930
2931 // Loop over velocity components
2932 for (unsigned k = 0; k < DIM; k++)
2933 {
2934 aux += interpolated_dudx(k, k);
2935 }
2936
2937 // Multiply through by deriv of Jacobian and integration weight
2938 dresidual_dnodal_coordinates(local_eqn, p, q) +=
2939 aux * dJ_dX(p, q) * testp[l] * w;
2940
2941 // Derivative of residual x Jacobian
2942 //----------------------------------
2943
2944 // Loop over velocity components
2945 aux = -source_gradient[p] * psif(q);
2946 for (unsigned k = 0; k < DIM; k++)
2947 {
2948 aux += d_dudx_dX(p, q, k, k);
2949 }
2950 // Multiply through by Jacobian and integration weight
2951 dresidual_dnodal_coordinates(local_eqn, p, q) +=
2952 aux * testp[l] * J * w;
2953
2954 // Derivs w.r.t. to nodal velocities
2955 //---------------------------------
2956 if (element_has_node_with_aux_node_update_fct)
2957 {
2958 aux = dpsifdx(q, p) * testp(l);
2959 dresidual_dnodal_coordinates(local_eqn, p, q) +=
2960 aux * d_U_dX(p, q) * J * w;
2961 }
2962 }
2963 }
2964 }
2965 }
2966 } // End of loop over integration points
2967 } // End of get_dresidual_dnodal_coordinates
2968
2969
2970 /// ////////////////////////////////////////////////////////////////////////
2971 /// ////////////////////////////////////////////////////////////////////////
2972 /// ////////////////////////////////////////////////////////////////////////
2973
2974 //---------------------------------
2975 // Space-time Taylor-Hood elements:
2976 //---------------------------------
2977 // Set the data for the number of Variables at each node
2978 template<>
2980 3, 2, 3, 2, 2, 2, 3, 2, 3, 3, 2, 3, 2, 2,
2981 2, 3, 2, 3, 3, 2, 3, 2, 2, 2, 3, 2, 3};
2982
2983 // Set the data for the pressure conversion array
2984 template<>
2986 0, 2, 6, 8, 9, 11, 15, 17, 18, 20, 24, 26};
2987
2988 //=========================================================================
2989 /// Add to the set \c paired_load_data pairs containing
2990 /// - the pointer to a Data object
2991 /// and
2992 /// - the index of the value in that Data object
2993 /// .
2994 /// for all values (pressures, velocities) that affect the
2995 /// load computed in the \c get_load(...) function.
2996 //=========================================================================
2997 template<unsigned DIM>
2999 std::set<std::pair<Data*, unsigned>>& paired_load_data)
3000 {
3001 // Allocate storage for the indices of the velocity components
3002 unsigned u_index[DIM];
3003
3004 // Loop over the velocity components
3005 for (unsigned i = 0; i < DIM; i++)
3006 {
3007 // Get the index at which the i-th velocity component is stored
3008 u_index[i] = this->u_index_nst(i);
3009 }
3010
3011 // Get the number of nodes in this element
3012 unsigned n_node = this->nnode();
3013
3014 // Loop over the nodes
3015 for (unsigned n = 0; n < n_node; n++)
3016 {
3017 // Loop over the velocity components and add pointer to their data
3018 // and indices to the vectors
3019 for (unsigned i = 0; i < DIM; i++)
3020 {
3021 // Add in the node and equation number pair
3022 paired_load_data.insert(std::make_pair(this->node_pt(n), u_index[i]));
3023 }
3024 } // for (unsigned n=0;n<n_node;n++)
3025
3026 // Identify the pressure data
3027 this->identify_pressure_data(paired_load_data);
3028 } // End of identify_load_data
3029
3030 //=========================================================================
3031 /// Add to the set \c paired_pressure_data pairs containing
3032 /// - the pointer to a Data object
3033 /// and
3034 /// - the index of the value in that Data object
3035 /// .
3036 /// for pressure values that affect the
3037 /// load computed in the \c get_load(...) function.,
3038 //=========================================================================
3039 template<unsigned DIM>
3041 std::set<std::pair<Data*, unsigned>>& paired_pressure_data)
3042 {
3043 // Find the index at which the pressure is stored
3044 unsigned p_index = static_cast<unsigned>(this->p_nodal_index_nst());
3045
3046 // Get the number of pressure degrees of freedom
3047 unsigned n_pres = npres_nst();
3048
3049 // Loop over the pressure data
3050 for (unsigned l = 0; l < n_pres; l++)
3051 {
3052 // The DIM-th entry in each nodal data is the pressure, which affects
3053 // the traction
3054 paired_pressure_data.insert(
3055 std::make_pair(this->node_pt(Pconv[l]), p_index));
3056 }
3057 } // End of identify_pressure_data
3058
3059 //============================================================================
3060 /// Create a list of pairs for all unknowns in this element,
3061 /// so the first entry in each pair contains the global equation
3062 /// number of the unknown, while the second one contains the number
3063 /// of the "DOF type" that this unknown is associated with.
3064 /// (Function can obviously only be called if the equation numbering
3065 /// scheme has been set up.)
3066 //============================================================================
3067 template<unsigned DIM>
3069 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
3070 {
3071 // Get the number of nodes in this element
3072 unsigned n_node = this->nnode();
3073
3074 // Temporary pair (used to store dof lookup prior to being added to list)
3075 std::pair<unsigned, unsigned> dof_lookup;
3076
3077 // Loop over the nodes
3078 for (unsigned n = 0; n < n_node; n++)
3079 {
3080 // Find the number of values at this node
3081 unsigned n_value = this->required_nvalue(n);
3082
3083 // Loop over these values
3084 for (unsigned i_value = 0; i_value < n_value; i_value++)
3085 {
3086 // Determine the local eqn number associated with this value
3087 int local_eqn_number = this->nodal_local_eqn(n, i_value);
3088
3089 // Ignore pinned values - far away degrees of freedom resulting
3090 // from hanging nodes can be ignored since these are be dealt
3091 // with by the element containing their master nodes
3092 if (local_eqn_number >= 0)
3093 {
3094 // Store dof lookup in temporary pair; the global equation number
3095 // is the first entry in pair
3096 dof_lookup.first = this->eqn_number(local_eqn_number);
3097
3098 // Set dof numbers; dof number is the second entry in pair
3099 // DRAIG: Uncomment whichever one you want to use. Setting
3100 // all dof numbers to 0 means you don't distinguish between
3101 // velocity and pressure component; just aggregate everything.
3102 // In contrast, setting the dof number to i_value means the
3103 // dofs associated with the first velocity component, second
3104 // velocity component and the pressure are all separated
3105 // dof_lookup.second=i_value;
3106 dof_lookup.second = 0;
3107
3108 // Add to list
3109 dof_lookup_list.push_front(dof_lookup);
3110 }
3111 } // for (unsigned v=0;v<nv;v++)
3112 } // for (unsigned n=0;n<n_node;n++)
3113 } // End of get_dof_numbers_for_unknowns
3114
3115
3116 //====================================================================
3117 /// Force build of templates
3118 //====================================================================
3119 template class SpaceTimeNavierStokesEquations<2>;
3120 template class QTaylorHoodSpaceTimeElement<2>;
3121} // namespace oomph
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
Definition: nodes.h:324
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:485
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:491
void initialise(const T &val)
Initialize all values in the matrix to val.
Definition: matrices.h:514
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
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1198
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
void perform_auxiliary_node_update_fct()
Execute auxiliary update function (if any) – this can be used to update any nodal values following th...
Definition: nodes.h:1615
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
bool has_auxiliary_node_update_fct_pt()
Boolean to indicate if node has a pointer to and auxiliary update function.
Definition: nodes.h:1606
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1701
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed,...
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not N.B. This needs to be public so th...
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)
Compute traction (on the viscous scale) exerted onto the fluid at local coordinate s....
void full_output(std::ostream &outfile)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format....
void compute_norm(Vector< double > &norm)
Compute the vector norm of the FEM solution.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution specified via function pointer at a given number of plot points....
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
void fill_in_contribution_to_hessian_vector_products(Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Compute the hessian tensor vector products required to perform continuation of bifurcations analytica...
virtual void fill_in_generic_pressure_advection_diffusion_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
virtual void fill_in_generic_dresidual_contribution_nst(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, const unsigned &flag)
Compute the derivatives of the residuals for the Navier-Stokes equations with respect to a parameter ...
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
void output(std::ostream &outfile)
Output function: x,y,t,u,v,p in tecplot format. The default number of plot points is five.
double d_kin_energy_dt() const
Get integral of time derivative of kinetic energy over element.
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
virtual void fill_in_generic_residual_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, const unsigned &flag)
Compute the residuals for the Navier-Stokes equations. Flag=1 (or 0): do (or don't) compute the Jacob...
void output_veloc(std::ostream &outfile, const unsigned &nplot, const unsigned &t)
Output function: x,y,t,u,v in tecplot format. Use n_plot points in each coordinate direction at times...
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Validate against exact solution at given time Solution is provided via function pointer....
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
void get_vorticity(const Vector< double > &s, Vector< double > &vorticity) const
Compute the vorticity vector at local coordinate s.
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: 1/2 (du_i/dx_j+du_j/dx_i)
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
void output_vorticity(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,t,omega in tecplot format. nplot points in each coordinate direction.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...