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