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-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Non-inline functions for NS elements
28 
29 /// //////////////////////////////////////////////////////////////////////
30 /// //////////////////////////////////////////////////////////////////////
31 /// //////////////////////////////////////////////////////////////////////
32 
33 namespace 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>
195  void SpaceTimeNavierStokesEquations<DIM>::compute_norm(Vector<double>& norm)
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>
928  void SpaceTimeNavierStokesEquations<DIM>::output_veloc(std::ostream& outfile,
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>
982  void SpaceTimeNavierStokesEquations<DIM>::output(std::ostream& outfile,
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>
1091  void SpaceTimeNavierStokesEquations<DIM>::full_output(std::ostream& outfile,
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>
1994  fill_in_generic_residual_contribution_nst(Vector<double>& residuals,
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 Strouhal value 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  // Update the i-th mesh velocity component
2163  mesh_velocity[i] += raw_nodal_position(l, i) * dpsifdx(l, DIM);
2164  }
2165  } // for (unsigned l=0;l<n_node;l++)
2166  } // if (!ALE_is_disabled)
2167 
2168  // Allocate space for the body force
2169  Vector<double> body_force(DIM, 0.0);
2170 
2171  // Get the user-defined body force term
2172  get_body_force_nst(interpolated_t, ipt, s, interpolated_x, body_force);
2173 
2174  // Get the user-defined source function
2175  double source = get_source_nst(interpolated_t, ipt, interpolated_x);
2176 
2177  //---------------------------------
2178  // Assemble residuals and Jacobian:
2179  //---------------------------------
2180  //--------------------
2181  // Momentum equations:
2182  //--------------------
2183  // Loop over the test functions
2184  for (unsigned l = 0; l < n_node; l++)
2185  {
2186  // Loop over the velocity components
2187  for (unsigned i = 0; i < DIM; i++)
2188  {
2189  // Get the local equation number associated with this unknown and node
2190  local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
2191 
2192  // If it's not a boundary condition
2193  if (local_eqn >= 0)
2194  {
2195  // Add the user-defined body force terms
2196  residuals[local_eqn] += body_force[i] * testf[l] * W;
2197 
2198  // Add the gravitational body force term
2199  residuals[local_eqn] += scaled_re_inv_fr * testf[l] * G[i] * W;
2200 
2201  // Add the pressure gradient term
2202  residuals[local_eqn] += interpolated_p * dtestfdx(l, i) * W;
2203 
2204  // Add in the contribution from the time derivative
2205  residuals[local_eqn] -=
2206  scaled_re_st * interpolated_dudt[i] * testf[l] * W;
2207 
2208  // If ALE is enabled
2209  if (!ALE_is_disabled)
2210  {
2211  // Loop over the coordinate directions
2212  for (unsigned k = 0; k < DIM; k++)
2213  {
2214  // Add in the mesh velocity contribution
2215  residuals[local_eqn] +=
2216  (scaled_re_st * mesh_velocity[k] * interpolated_dudx(i, k) *
2217  testf[l] * W);
2218  }
2219  } // if (!ALE_is_disabled)
2220 
2221  // Loop over the coordinate directions
2222  for (unsigned k = 0; k < DIM; k++)
2223  {
2224  // Add in the convective term contribution
2225  residuals[local_eqn] -= (scaled_re * interpolated_u[k] *
2226  interpolated_dudx(i, k) * testf[l] * W);
2227  }
2228 
2229  // Loop over the coordinate directions
2230  for (unsigned k = 0; k < DIM; k++)
2231  {
2232  // Add in the stress tensor terms:
2233  // NOTE: The viscosity ratio needs to go in here to ensure
2234  // continuity of normal stress is satisfied even in flows
2235  // with zero pressure gradient!
2236  residuals[local_eqn] -= ((interpolated_dudx(i, k) +
2237  Gamma[i] * interpolated_dudx(k, i)) *
2238  visc_ratio * dtestfdx(l, k) * W);
2239  }
2240 
2241  //------------------------
2242  // Calculate the Jacobian:
2243  //------------------------
2244  // If we also need to construct the Jacobian
2245  if (flag)
2246  {
2247  // Loop over the velocity shape functions again
2248  for (unsigned l2 = 0; l2 < n_node; l2++)
2249  {
2250  // Loop over the velocity components again
2251  for (unsigned i2 = 0; i2 < DIM; i2++)
2252  {
2253  // Get the local equation number associated with this unknown
2254  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
2255 
2256  // If we're at a proper degree of freedom
2257  if (local_unknown >= 0)
2258  {
2259  // Add the contribution to the elemental matrix
2260  jacobian(local_eqn, local_unknown) -=
2261  (visc_ratio * Gamma[i] * dpsifdx(l2, i) *
2262  dtestfdx(l, i2) * W);
2263 
2264  // Now add in the inertial terms
2265  jacobian(local_eqn, local_unknown) -=
2266  (scaled_re * psif[l2] * interpolated_dudx(i, i2) *
2267  testf[l] * W);
2268 
2269  // Extra component if i2=i
2270  if (i2 == i)
2271  {
2272  // If we also need to construct the mass matrix (only
2273  // diagonal entries)
2274  if (flag == 2)
2275  {
2276  // NOTE: This is positive because the mass matrix is
2277  // taken to the other side of the equation when
2278  // formulating the generalised eigenproblem.
2279  mass_matrix(local_eqn, local_unknown) +=
2280  (scaled_re_st * psif[l2] * testf[l] * W);
2281  }
2282 
2283  // Add in the time-derivative contribution
2284  jacobian(local_eqn, local_unknown) -=
2285  (scaled_re_st * dpsifdx(l2, DIM) * testf[l] * W);
2286 
2287  // If ALE is enabled
2288  if (!ALE_is_disabled)
2289  {
2290  // Loop over the velocity components
2291  for (unsigned k = 0; k < DIM; k++)
2292  {
2293  // Add in the mesh velocity contribution
2294  jacobian(local_eqn, local_unknown) +=
2295  (scaled_re_st * mesh_velocity[k] * dpsifdx(l2, k) *
2296  testf[l] * W);
2297  }
2298  } // if (!ALE_is_disabled)
2299 
2300  // Loop over the velocity components
2301  for (unsigned k = 0; k < DIM; k++)
2302  {
2303  // Add in the convective term contribution
2304  jacobian(local_eqn, local_unknown) -=
2305  (scaled_re * interpolated_u[k] * dpsifdx(l2, k) *
2306  testf[l] * W);
2307  }
2308 
2309  // Loop over the velocity components
2310  for (unsigned k = 0; k < DIM; k++)
2311  {
2312  // Add in the velocity gradient terms
2313  jacobian(local_eqn, local_unknown) -=
2314  (visc_ratio * dpsifdx(l2, k) * dtestfdx(l, k) * W);
2315  }
2316  } // if (i2==i)
2317  } // if (local_unknown>=0)
2318  } // for (unsigned i2=0;i2<DIM;i2++)
2319  } // for (unsigned l2=0;l2<n_node;l2++)
2320 
2321  // Loop over pressure shape functions
2322  for (unsigned l2 = 0; l2 < n_pres; l2++)
2323  {
2324  // Get the local equation number associated with this degree of
2325  // freedom
2326  local_unknown = p_local_eqn(l2);
2327 
2328  // If we are at a proper degree of freedom
2329  if (local_unknown >= 0)
2330  {
2331  // Add in the pressure gradient contribution
2332  jacobian(local_eqn, local_unknown) +=
2333  psip[l2] * dtestfdx(l, i) * W;
2334  }
2335  } // for (unsigned l2=0;l2<n_pres;l2++)
2336 
2337  //------------------------------------
2338  // Calculate external data information
2339  //------------------------------------
2340  // If we're storing the Strouhal number as external data then we
2341  // need to calculate dr/d(St) in the elemental Jacobian
2342  if (Strouhal_is_stored_as_external_data)
2343  {
2344  // The index of the external data (there's only one!)
2345  unsigned data_index = 0;
2346 
2347  // The index of the unknown value stored in the external data
2348  unsigned value_index = 0;
2349 
2350  // Get the local equation number associated with the extra
2351  // unknown
2352  local_unknown =
2353  this->external_local_eqn(data_index, value_index);
2354 
2355  // If we're at a non-zero degree of freedom add in the entry
2356  if (local_unknown >= 0)
2357  {
2358  // Add in the contribution from the time derivative
2359  jacobian(local_eqn, local_unknown) -=
2360  (scaled_dre_st_dst * interpolated_dudt[i] * testf[l] * W);
2361 
2362  // If ALE is enabled
2363  if (!ALE_is_disabled)
2364  {
2365  // Loop over the coordinate directions
2366  for (unsigned k = 0; k < DIM; k++)
2367  {
2368  // Add in the mesh velocity contribution
2369  jacobian(local_eqn, local_unknown) +=
2370  (scaled_dre_st_dst * mesh_velocity[k] *
2371  interpolated_dudx(i, k) * testf[l] * W);
2372  }
2373  } // if (!ALE_is_disabled)
2374  } // if (local_unknown>=0)
2375  } // if (Strouhal_is_stored_as_external_data)
2376  } // if (flag)
2377  } // if (local_eqn>=0)
2378  } // for (unsigned i=0;i<DIM;i++)
2379  } // for (unsigned l=0;l<n_node;l++)
2380 
2381  //---------------------
2382  // Continuity equation:
2383  //---------------------
2384  // Loop over the shape functions
2385  for (unsigned l = 0; l < n_pres; l++)
2386  {
2387  // Get the local equation number associated with the pressure dof
2388  local_eqn = p_local_eqn(l);
2389 
2390  // If it's not a boundary condition
2391  if (local_eqn >= 0)
2392  {
2393  // Add in the source term contribution
2394  residuals[local_eqn] -= source * testp[l] * W;
2395 
2396  // Loop over velocity components
2397  for (unsigned k = 0; k < DIM; k++)
2398  {
2399  // Add in the velocity gradient terms
2400  residuals[local_eqn] += interpolated_dudx(k, k) * testp[l] * W;
2401  }
2402 
2403  //------------------------
2404  // Calculate the Jacobian:
2405  //------------------------
2406  // If we also need to construct the Jacobian
2407  if (flag)
2408  {
2409  // Loop over the velocity shape functions
2410  for (unsigned l2 = 0; l2 < n_node; l2++)
2411  {
2412  // Loop over velocity components
2413  for (unsigned i2 = 0; i2 < DIM; i2++)
2414  {
2415  // Get the local equation number associated with this node
2416  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
2417 
2418  // If we're at a proper degree of freedom
2419  if (local_unknown >= 0)
2420  {
2421  // Add in the velocity gradient contribution
2422  jacobian(local_eqn, local_unknown) +=
2423  dpsifdx(l2, i2) * testp[l] * W;
2424  }
2425  } // for (unsigned i2=0;i2<DIM;i2++)
2426  } // for (unsigned l2=0;l2<n_node;l2++)
2427  } // if (flag)
2428  } // if (local_eqn>=0)
2429  } // for (unsigned l=0;l<n_pres;l++)
2430  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
2431  } // End of fill_in_generic_residual_contribution_nst
2432 
2433  //================================================================
2434  /// Compute the derivatives of the residuals for the Navier-Stokes
2435  /// equations with respect to a parameter;
2436  /// flag=2 or 1 or 0: do (or don't) compute the Jacobian and
2437  /// mass matrix as well
2438  //================================================================
2439  template<unsigned DIM>
2442  double* const& parameter_pt,
2443  Vector<double>& dres_dparam,
2444  DenseMatrix<double>& djac_dparam,
2445  DenseMatrix<double>& dmass_matrix_dparam,
2446  const unsigned& flag)
2447  {
2448  // Throw an error
2449  throw OomphLibError("Not yet implemented\n",
2450  OOMPH_CURRENT_FUNCTION,
2451  OOMPH_EXCEPTION_LOCATION);
2452  } // End of fill_in_generic_dresidual_contribution_nst
2453 
2454  //==================================================================
2455  /// Compute the hessian tensor vector products required to
2456  /// perform continuation of bifurcations analytically
2457  //==================================================================
2458  template<unsigned DIM>
2461  Vector<double> const& Y,
2462  DenseMatrix<double> const& C,
2463  DenseMatrix<double>& product)
2464  {
2465  // Throw an error
2466  throw OomphLibError("Not yet implemented\n",
2467  OOMPH_CURRENT_FUNCTION,
2468  OOMPH_EXCEPTION_LOCATION);
2469  } // End of fill_in_contribution_to_hessian_vector_products
2470 
2471  //======================================================================
2472  /// Compute derivatives of elemental residual vector with respect
2473  /// to nodal coordinates.
2474  /// dresidual_dnodal_coordinates(l,i,j)=d res(l) / dX_{ij}
2475  /// Overloads the FD-based version in the FE base class.
2476  /// DRAIG: This needs doing carefully if the ALE nodes aren't fixed!!!
2477  //======================================================================
2478  template<unsigned DIM>
2480  RankThreeTensor<double>& dresidual_dnodal_coordinates)
2481  {
2482  // Throw a warning
2483  throw OomphLibError("Space-time update needs to be checked!",
2484  OOMPH_CURRENT_FUNCTION,
2485  OOMPH_EXCEPTION_LOCATION);
2486 
2487  // Return immediately if there are no dofs
2488  if (ndof() == 0)
2489  {
2490  return;
2491  }
2492 
2493  // Determine number of nodes in element
2494  const unsigned n_node = nnode();
2495 
2496  // Determine number of pressure dofs in element
2497  const unsigned n_pres = npres_nst();
2498 
2499  // Find the indices at which the local velocities are stored
2500  unsigned u_nodal_index[DIM];
2501  for (unsigned i = 0; i < DIM; i++)
2502  {
2503  u_nodal_index[i] = u_index_nst(i);
2504  }
2505 
2506  // Set up memory for the shape and test functions
2507  Shape psif(n_node);
2508  Shape testf(n_node);
2509  DShape dpsifdx(n_node, DIM + 1);
2510  DShape dtestfdx(n_node, DIM + 1);
2511 
2512  // Set up memory for pressure shape and test functions
2513  Shape psip(n_pres), testp(n_pres);
2514 
2515  // Deriatives of shape fct derivatives w.r.t. nodal coords
2516  RankFourTensor<double> d_dpsifdx_dX(DIM, n_node, n_node, DIM);
2517  RankFourTensor<double> d_dtestfdx_dX(DIM, n_node, n_node, DIM);
2518 
2519  // Derivative of Jacobian of mapping w.r.t. to nodal coords
2520  DenseMatrix<double> dJ_dX(DIM, n_node);
2521 
2522  // Derivatives of derivative of u w.r.t. nodal coords
2523  RankFourTensor<double> d_dudx_dX(DIM, n_node, DIM, DIM);
2524 
2525  // Derivatives of nodal velocities w.r.t. nodal coords:
2526  // Assumption: Interaction only local via no-slip so
2527  // X_ij only affects U_ij.
2528  DenseMatrix<double> d_U_dX(DIM, n_node, 0.0);
2529 
2530  // Determine the number of integration points
2531  const unsigned n_intpt = integral_pt()->nweight();
2532 
2533  // Vector to hold local coordinates
2534  Vector<double> s(DIM + 1, 0.0);
2535 
2536  // Get physical variables from element
2537  // (Reynolds number must be multiplied by the density ratio)
2538  double scaled_re = re() * density_ratio();
2539  double scaled_re_st = re_st() * density_ratio();
2540  double scaled_re_inv_fr = re_invfr() * density_ratio();
2541  double visc_ratio = viscosity_ratio();
2542  Vector<double> G = g();
2543 
2544  // FD step
2546 
2547  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
2548  // Assumption: Interaction only local via no-slip so
2549  // X_ij only affects U_ij.
2550  bool element_has_node_with_aux_node_update_fct = false;
2551  for (unsigned q = 0; q < n_node; q++)
2552  {
2553  Node* nod_pt = node_pt(q);
2554 
2555  // Only compute if there's a node-update fct involved
2556  if (nod_pt->has_auxiliary_node_update_fct_pt())
2557  {
2558  element_has_node_with_aux_node_update_fct = true;
2559 
2560  // Current nodal velocity
2561  Vector<double> u_ref(DIM);
2562  for (unsigned i = 0; i < DIM; i++)
2563  {
2564  u_ref[i] = *(nod_pt->value_pt(u_nodal_index[i]));
2565  }
2566 
2567  // FD
2568  for (unsigned p = 0; p < DIM; p++)
2569  {
2570  // Make backup
2571  double backup = nod_pt->x(p);
2572 
2573  // Do FD step. No node update required as we're
2574  // attacking the coordinate directly...
2575  nod_pt->x(p) += eps_fd;
2576 
2577  // Do auxiliary node update (to apply no-slip)
2578  nod_pt->perform_auxiliary_node_update_fct();
2579 
2580  // Evaluate
2581  d_U_dX(p, q) =
2582  (*(nod_pt->value_pt(u_nodal_index[p])) - u_ref[p]) / eps_fd;
2583 
2584  // Reset
2585  nod_pt->x(p) = backup;
2586 
2587  // Do auxiliary node update (to apply no slip)
2588  nod_pt->perform_auxiliary_node_update_fct();
2589  }
2590  }
2591  }
2592 
2593  // Integer to store the local equation number
2594  int local_eqn = 0;
2595 
2596  // Loop over the integration points
2597  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2598  {
2599  // Assign values of s
2600  for (unsigned i = 0; i < DIM + 1; i++)
2601  {
2602  s[i] = integral_pt()->knot(ipt, i);
2603  }
2604 
2605  // Get the integral weight
2606  const double w = integral_pt()->weight(ipt);
2607 
2608  // Call the derivatives of the shape and test functions
2609  const double J = dshape_and_dtest_eulerian_at_knot_nst(ipt,
2610  psif,
2611  dpsifdx,
2612  d_dpsifdx_dX,
2613  testf,
2614  dtestfdx,
2615  d_dtestfdx_dX,
2616  dJ_dX);
2617 
2618  // Call the pressure shape and test functions
2619  pshape_nst(s, psip, testp);
2620 
2621  // Storage for the interpolated time value
2622  double interpolated_t = 0.0;
2623 
2624  // Storage for the interpolated pressure value
2625  double interpolated_p = 0.0;
2626 
2627  // Storage for the spatial coordinates
2628  Vector<double> interpolated_x(DIM, 0.0);
2629 
2630  // Storage for the interpolated velocity components
2631  Vector<double> interpolated_u(DIM, 0.0);
2632 
2633  // Storage for the interpolated time-derivative of the velocities
2634  Vector<double> interpolated_dudt(DIM, 0.0);
2635 
2636  // Storage for the mesh velocity
2637  Vector<double> mesh_velocity(DIM, 0.0);
2638 
2639  // Storage for the spatial derivatives of the velocity components
2640  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
2641 
2642  // Calculate pressure
2643  for (unsigned l = 0; l < n_pres; l++)
2644  {
2645  // Update the current approximation to the interpolated pressure
2646  interpolated_p += p_nst(l) * psip[l];
2647  }
2648 
2649  //-------------------------------------------------------------------
2650  // Calculate velocities, derivatives, source function and body force:
2651  //-------------------------------------------------------------------
2652  // Loop over nodes
2653  for (unsigned l = 0; l < n_node; l++)
2654  {
2655  // Update the interpolated time value
2656  interpolated_t += raw_nodal_position(l, DIM) * psif(l);
2657 
2658  // Loop over coordinate directions
2659  for (unsigned i = 0; i < DIM; i++)
2660  {
2661  // Get the nodal value
2662  double u_value = raw_nodal_value(l, u_nodal_index[i]);
2663 
2664  // Update the i-th interpolated velocity component
2665  interpolated_u[i] += u_value * psif[l];
2666 
2667  // Update the i-th interpolated coordinate value
2668  interpolated_x[i] += raw_nodal_position(l, i) * psif[l];
2669 
2670  // Update the interpolated du_i/dt value
2671  interpolated_dudt[i] += u_value * dpsifdx(l, DIM);
2672 
2673  // Loop over derivative directions
2674  for (unsigned j = 0; j < DIM; j++)
2675  {
2676  // Update the interpolated du_i/dx_j value
2677  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
2678  }
2679  } // for (unsigned i=0;i<DIM;i++)
2680  } // for (unsigned l=0;l<n_node;l++)
2681 
2682  // If ALE is enabled
2683  if (!ALE_is_disabled)
2684  {
2685  // Loop over nodes
2686  for (unsigned l = 0; l < n_node; l++)
2687  {
2688  // Loop over directions
2689  for (unsigned i = 0; i < DIM; i++)
2690  {
2691  // Update the i-th mesh velocity component
2692  mesh_velocity[i] += raw_nodal_position(l, i) * dpsifdx(l, DIM);
2693  }
2694  } // for (unsigned l=0;l<n_node;l++)
2695  } // if (!ALE_is_disabled)
2696 
2697  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
2698  // DRAIG: CHECK
2699  for (unsigned q = 0; q < n_node; q++)
2700  {
2701  // Loop over coordinate directions
2702  for (unsigned p = 0; p < DIM; p++)
2703  {
2704  for (unsigned i = 0; i < DIM; i++)
2705  {
2706  for (unsigned k = 0; k < DIM; k++)
2707  {
2708  double aux = 0.0;
2709 
2710  for (unsigned j = 0; j < n_node; j++)
2711  {
2712  aux += raw_nodal_value(j, u_nodal_index[i]) *
2713  d_dpsifdx_dX(p, q, j, k);
2714  }
2715 
2716  //
2717  d_dudx_dX(p, q, i, k) = aux;
2718  }
2719  }
2720  }
2721  }
2722 
2723  // Allocate space for the body force
2724  Vector<double> body_force(DIM);
2725 
2726  // Get the user-defined body force term
2727  get_body_force_nst(interpolated_t, ipt, s, interpolated_x, body_force);
2728 
2729  // Get the user-defined source function
2730  const double source = get_source_nst(interpolated_t, ipt, interpolated_x);
2731 
2732  // Note: Can use raw values and avoid bypassing hanging information
2733  // because this is the non-refineable version!
2734 
2735  // Allocate space for the gradient of the body force function
2736  DenseMatrix<double> d_body_force_dx(DIM, DIM, 0.0);
2737 
2738  // Get gradient of body force function
2739  get_body_force_gradient_nst(
2740  interpolated_t, ipt, s, interpolated_x, d_body_force_dx);
2741 
2742  // Allocate space for the gradient of the source function
2743  Vector<double> source_gradient(DIM, 0.0);
2744 
2745  // Get gradient of source function
2746  get_source_gradient_nst(
2747  interpolated_t, ipt, interpolated_x, source_gradient);
2748 
2749  // Get weight of actual nodal position in computation of mesh velocity
2750  // from positional time stepper
2751  const double pos_time_weight =
2752  (node_pt(0)->position_time_stepper_pt()->weight(1, 0));
2753 
2754  // Get weight of actual nodal value in computation of mesh velocity
2755  // from value time stepper
2756  const double val_time_weight =
2757  node_pt(0)->time_stepper_pt()->weight(1, 0);
2758 
2759  //---------------------------------
2760  // Assemble residuals and Jacobian:
2761  //---------------------------------
2762  //--------------------
2763  // Momentum equations:
2764  //--------------------
2765  // Loop over the test functions
2766  for (unsigned l = 0; l < n_node; l++)
2767  {
2768  // Loop over coordinate directions
2769  for (unsigned i = 0; i < DIM; i++)
2770  {
2771  // Get the local equation number associated with this degree of
2772  // freedom
2773  local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
2774 
2775  // If it's not a boundary condition
2776  if (local_eqn >= 0)
2777  {
2778  // Loop over coordinate directions
2779  for (unsigned p = 0; p < DIM; p++)
2780  {
2781  // Loop over nodes
2782  for (unsigned q = 0; q < n_node; q++)
2783  {
2784  //-----------------------------------
2785  // Residual x derivative of Jacobian:
2786  //-----------------------------------
2787  // Add the user-defined body force terms
2788  double sum = body_force[i] * testf[l];
2789 
2790  // Add the gravitational body force term
2791  sum += scaled_re_inv_fr * testf[l] * G[i];
2792 
2793  // Add the pressure gradient term
2794  sum += interpolated_p * dtestfdx(l, i);
2795 
2796  // Loop over the coordinate directions
2797  for (unsigned k = 0; k < DIM; k++)
2798  {
2799  // Add in the stress tensor contribution
2800  // NOTE: The viscosity ratio needs to go in here to ensure
2801  // continuity of normal stress is satisfied even in flows
2802  // with zero pressure gradient!
2803  sum -= (visc_ratio *
2804  (interpolated_dudx(i, k) +
2805  Gamma[i] * interpolated_dudx(k, i)) *
2806  dtestfdx(l, k));
2807  }
2808 
2809  // Add in the contribution from the time derivative
2810  sum -= scaled_re_st * interpolated_dudt[i] * testf[l];
2811 
2812  // Loop over the coordinate directions
2813  for (unsigned k = 0; k < DIM; k++)
2814  {
2815  // Add in convective term contribution
2816  sum -= scaled_re * interpolated_u[k] *
2817  interpolated_dudx(i, k) * testf[l];
2818  }
2819 
2820  // If ALE is enabled
2821  if (!ALE_is_disabled)
2822  {
2823  // Loop over the coordinate directions
2824  for (unsigned k = 0; k < DIM; k++)
2825  {
2826  // Add in the mesh velocity contribution
2827  sum += scaled_re_st * mesh_velocity[k] *
2828  interpolated_dudx(i, k) * testf[l];
2829  }
2830  }
2831 
2832  // Multiply through by derivative of Jacobian and integration
2833  // weight
2834  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2835  sum * dJ_dX(p, q) * w;
2836 
2837  //----------------------------------
2838  // Derivative of residual x Jacobian
2839  //----------------------------------
2840  // Body force
2841  sum = d_body_force_dx(i, p) * psif(q) * testf(l);
2842 
2843  // Pressure gradient term
2844  sum += interpolated_p * d_dtestfdx_dX(p, q, l, i);
2845 
2846  // Viscous term
2847  for (unsigned k = 0; k < DIM; k++)
2848  {
2849  sum -= (visc_ratio * ((interpolated_dudx(i, k) +
2850  Gamma[i] * interpolated_dudx(k, i)) *
2851  d_dtestfdx_dX(p, q, l, k) +
2852  (d_dudx_dX(p, q, i, k) +
2853  Gamma[i] * d_dudx_dX(p, q, k, i)) *
2854  dtestfdx(l, k)));
2855  }
2856 
2857  // Convective terms, including mesh velocity
2858  for (unsigned k = 0; k < DIM; k++)
2859  {
2860  double tmp = scaled_re * interpolated_u[k];
2861  if (!ALE_is_disabled)
2862  {
2863  tmp -= scaled_re_st * mesh_velocity[k];
2864  }
2865  sum -= tmp * d_dudx_dX(p, q, i, k) * testf(l);
2866  }
2867 
2868  if (!ALE_is_disabled)
2869  {
2870  sum += scaled_re_st * pos_time_weight * psif(q) *
2871  interpolated_dudx(i, p) * testf(l);
2872  }
2873 
2874  // Multiply through by Jacobian and integration weight
2875  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2876 
2877  // Derivs w.r.t. to nodal velocities
2878  // ---------------------------------
2879  if (element_has_node_with_aux_node_update_fct)
2880  {
2881  sum =
2882  -visc_ratio * Gamma[i] * dpsifdx(q, i) * dtestfdx(l, p) -
2883  scaled_re * psif(q) * interpolated_dudx(i, p) * testf(l);
2884  if (i == p)
2885  {
2886  sum -= scaled_re_st * val_time_weight * psif(q) * testf(l);
2887  for (unsigned k = 0; k < DIM; k++)
2888  {
2889  sum -= visc_ratio * dpsifdx(q, k) * dtestfdx(l, k);
2890  double tmp = scaled_re * interpolated_u[k];
2891  if (!ALE_is_disabled)
2892  tmp -= scaled_re_st * mesh_velocity[k];
2893  sum -= tmp * dpsifdx(q, k) * testf(l);
2894  }
2895  }
2896  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2897  sum * d_U_dX(p, q) * J * w;
2898  }
2899  }
2900  }
2901  }
2902  }
2903  } // End of loop over test functions
2904 
2905 
2906  // CONTINUITY EQUATION
2907  // -------------------
2908 
2909  // Loop over the shape functions
2910  for (unsigned l = 0; l < n_pres; l++)
2911  {
2912  local_eqn = p_local_eqn(l);
2913 
2914  // If not a boundary conditions
2915  if (local_eqn >= 0)
2916  {
2917  // Loop over coordinate directions
2918  for (unsigned p = 0; p < DIM; p++)
2919  {
2920  // Loop over nodes
2921  for (unsigned q = 0; q < n_node; q++)
2922  {
2923  // Residual x deriv of Jacobian
2924  //-----------------------------
2925 
2926  // Source term
2927  double aux = -source;
2928 
2929  // Loop over velocity components
2930  for (unsigned k = 0; k < DIM; k++)
2931  {
2932  aux += interpolated_dudx(k, k);
2933  }
2934 
2935  // Multiply through by deriv of Jacobian and integration weight
2936  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2937  aux * dJ_dX(p, q) * testp[l] * w;
2938 
2939  // Derivative of residual x Jacobian
2940  //----------------------------------
2941 
2942  // Loop over velocity components
2943  aux = -source_gradient[p] * psif(q);
2944  for (unsigned k = 0; k < DIM; k++)
2945  {
2946  aux += d_dudx_dX(p, q, k, k);
2947  }
2948  // Multiply through by Jacobian and integration weight
2949  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2950  aux * testp[l] * J * w;
2951 
2952  // Derivs w.r.t. to nodal velocities
2953  //---------------------------------
2954  if (element_has_node_with_aux_node_update_fct)
2955  {
2956  aux = dpsifdx(q, p) * testp(l);
2957  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2958  aux * d_U_dX(p, q) * J * w;
2959  }
2960  }
2961  }
2962  }
2963  }
2964  } // End of loop over integration points
2965  } // End of get_dresidual_dnodal_coordinates
2966 
2967  /// ////////////////////////////////////////////////////////////////////////
2968  /// ////////////////////////////////////////////////////////////////////////
2969  /// ////////////////////////////////////////////////////////////////////////
2970 
2971  //---------------------------------
2972  // Space-time Taylor-Hood elements:
2973  //---------------------------------
2974  // Set the data for the number of Variables at each node
2975  template<>
2977  3, 2, 3, 2, 2, 2, 3, 2, 3, 2, 2, 2, 2, 2,
2978  2, 2, 2, 2, 3, 2, 3, 2, 2, 2, 3, 2, 3};
2979 
2980  // Set the data for the pressure conversion array
2981  template<>
2983  0, 2, 6, 8, 18, 20, 24, 26};
2984 
2985  //=========================================================================
2986  /// Add to the set \c paired_load_data pairs containing
2987  /// - the pointer to a Data object
2988  /// and
2989  /// - the index of the value in that Data object
2990  /// .
2991  /// for all values (pressures, velocities) that affect the
2992  /// load computed in the \c get_load(...) function.
2993  //=========================================================================
2994  template<unsigned DIM>
2996  std::set<std::pair<Data*, unsigned>>& paired_load_data)
2997  {
2998  // Allocate storage for the indices of the velocity components
2999  unsigned u_index[DIM];
3000 
3001  // Loop over the velocity components
3002  for (unsigned i = 0; i < DIM; i++)
3003  {
3004  // Get the index at which the i-th velocity component is stored
3005  u_index[i] = this->u_index_nst(i);
3006  }
3007 
3008  // Get the number of nodes in this element
3009  unsigned n_node = this->nnode();
3010 
3011  // Loop over the nodes
3012  for (unsigned n = 0; n < n_node; n++)
3013  {
3014  // Loop over the velocity components and add pointer to their data
3015  // and indices to the vectors
3016  for (unsigned i = 0; i < DIM; i++)
3017  {
3018  // Add in the node and equation number pair
3019  paired_load_data.insert(std::make_pair(this->node_pt(n), u_index[i]));
3020  }
3021  } // for (unsigned n=0;n<n_node;n++)
3022 
3023  // Identify the pressure data
3024  this->identify_pressure_data(paired_load_data);
3025  } // End of identify_load_data
3026 
3027  //=========================================================================
3028  /// Add to the set \c paired_pressure_data pairs containing
3029  /// - the pointer to a Data object
3030  /// and
3031  /// - the index of the value in that Data object
3032  /// .
3033  /// for pressure values that affect the
3034  /// load computed in the \c get_load(...) function.,
3035  //=========================================================================
3036  template<unsigned DIM>
3038  std::set<std::pair<Data*, unsigned>>& paired_pressure_data)
3039  {
3040  // Find the index at which the pressure is stored
3041  unsigned p_index = static_cast<unsigned>(this->p_nodal_index_nst());
3042 
3043  // Get the number of pressure degrees of freedom
3044  unsigned n_pres = npres_nst();
3045 
3046  // Loop over the pressure data
3047  for (unsigned l = 0; l < n_pres; l++)
3048  {
3049  // The DIM-th entry in each nodal data is the pressure, which affects
3050  // the traction
3051  paired_pressure_data.insert(
3052  std::make_pair(this->node_pt(Pconv[l]), p_index));
3053  }
3054  } // End of identify_pressure_data
3055 
3056  /*
3057  //============================================================================
3058  /// Create a list of pairs for all unknowns in this element,
3059  /// so the first entry in each pair contains the global equation
3060  /// number of the unknown, while the second one contains the number
3061  /// of the "DOF type" that this unknown is associated with.
3062  /// (Function can obviously only be called if the equation numbering
3063  /// scheme has been set up.)
3064  //============================================================================
3065  template<unsigned DIM>
3066  void QTaylorHoodSpaceTimeElement<DIM>::get_dof_numbers_for_unknowns(
3067  std::list<std::pair<unsigned long, unsigned> >& dof_lookup_list) const
3068  {
3069  // Get the number of nodes in this element
3070  unsigned n_node=this->nnode();
3071 
3072  // Temporary pair (used to store dof lookup prior to being added to list)
3073  std::pair<unsigned,unsigned> dof_lookup;
3074 
3075  // Loop over the nodes
3076  for (unsigned n=0;n<n_node;n++)
3077  {
3078  // Find the number of values at this node
3079  unsigned n_value=this->required_nvalue(n);
3080 
3081  // Loop over these values
3082  for (unsigned i_value=0;i_value<n_value;i_value++)
3083  {
3084  // Determine the local eqn number associated with this value
3085  int local_eqn_number=this->nodal_local_eqn(n,i_value);
3086 
3087  // Ignore pinned values - far away degrees of freedom resulting
3088  // from hanging nodes can be ignored since these are be dealt
3089  // with by the element containing their master nodes
3090  if (local_eqn_number>=0)
3091  {
3092  // Store dof lookup in temporary pair; the global equation number
3093  // is the first entry in pair
3094  dof_lookup.first=this->eqn_number(local_eqn_number);
3095 
3096  // Set dof numbers; dof number is the second entry in pair
3097  // DRAIG: Uncomment whichever one you want to use. Setting
3098  // all dof numbers to 0 means you don't distinguish between
3099  // velocity and pressure component; just aggregate everything.
3100  // In contrast, setting the dof number to i_value means the
3101  // dofs associated with the first velocity component, second
3102  // velocity component and the pressure are all separated
3103  dof_lookup.second=i_value;
3104 
3105  // Add to list
3106  dof_lookup_list.push_front(dof_lookup);
3107  }
3108  } // for (unsigned v=0;v<nv;v++)
3109  } // for (unsigned n=0;n<n_node;n++)
3110  } // End of get_dof_numbers_for_unknowns
3111  */
3112 
3113  //====================================================================
3114  /// Force build of templates
3115  //====================================================================
3116  template class SpaceTimeNavierStokesEquations<2>;
3117  template class QTaylorHoodSpaceTimeElement<2>;
3118 } // 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
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
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void identify_load_data(std::set< std::pair< Data *, unsigned >> &paired_load_data)
Add to the set paired_load_data pairs containing.
void identify_pressure_data(std::set< std::pair< Data *, unsigned >> &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...