mixed_order_petrov_galerkin_space_time_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>
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>
45  int
47  -100;
48 
49  /// Navier-Stokes equations static data
50  template<unsigned DIM>
51  double SpaceTimeNavierStokesMixedOrderEquations<
53 
54  /// Navier-Stokes equations static data
55  template<unsigned DIM>
56  double SpaceTimeNavierStokesMixedOrderEquations<
57  DIM>::Default_Physical_Ratio_Value = 1.0;
58 
59  /// Navier-Stokes equations default gravity vector
60  template<unsigned DIM>
61  Vector<double> SpaceTimeNavierStokesMixedOrderEquations<
62  DIM>::Default_Gravity_vector(DIM, 0.0);
63 
64  //===================================================================
65  /// Compute the diagonal of the velocity/pressure mass matrices.
66  /// If which one=0, both are computed, otherwise only the pressure
67  /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
68  /// LSC version of the preconditioner only needs that one).
69  //===================================================================
70  template<unsigned DIM>
73  Vector<double>& press_mass_diag,
74  Vector<double>& veloc_mass_diag,
75  const unsigned& which_one)
76  {
77  // Resize and initialise
78  unsigned n_dof = ndof();
79 
80  // If space needs to be assigned for the pressure mass matrix
81  if ((which_one == 0) || (which_one == 1))
82  {
83  // Assign the appropriate amount of space
84  press_mass_diag.assign(n_dof, 0.0);
85  }
86 
87  // If space needs to be assigned for the velocity mass matrix
88  if ((which_one == 0) || (which_one == 2))
89  {
90  // Assign the appropriate amount of space
91  veloc_mass_diag.assign(n_dof, 0.0);
92  }
93 
94  // Find out how many nodes there are
95  unsigned n_node = nnode();
96 
97  // Storage for the local coordinates
98  Vector<double> s(DIM + 1, 0.0);
99 
100  // Storage for the local velocity indices
101  Vector<unsigned> u_nodal_index(DIM, 0.0);
102 
103  // Find the indices at which the local velocities are stored
104  for (unsigned i = 0; i < DIM; i++)
105  {
106  // Calculate the i-th local velocity component
107  u_nodal_index[i] = this->u_index_nst(i);
108  }
109 
110  // Set up memory for velocity shape functions
111  Shape psi(n_node);
112 
113  // Find number of pressure dofs
114  unsigned n_pres = npres_nst();
115 
116  // Pressure shape function
117  Shape psi_p(n_pres);
118 
119  // Number of integration points
120  unsigned n_intpt = integral_pt()->nweight();
121 
122  // Integer to store the local equations no
123  int local_eqn = 0;
124 
125  // Loop over the integration points
126  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
127  {
128  // Get the integral weight
129  double w = integral_pt()->weight(ipt);
130 
131  // Get determinant of Jacobian of the mapping
132  double J = J_eulerian_at_knot(ipt);
133 
134  // Assign values of s
135  for (unsigned i = 0; i < DIM + 1; i++)
136  {
137  // Calculate the i-th local coordinate at the ipt-th integration point
138  s[i] = integral_pt()->knot(ipt, i);
139  }
140 
141  // Premultiply weights and Jacobian
142  double W = w * J;
143 
144  // Do we want the velocity one?
145  if ((which_one == 0) || (which_one == 2))
146  {
147  // Get the velocity shape functions
148  shape_at_knot(ipt, psi);
149 
150  // Loop over the velocity shape functions
151  for (unsigned l = 0; l < n_node; l++)
152  {
153  // Loop over the velocity components
154  for (unsigned i = 0; i < DIM; i++)
155  {
156  // Get the local equation number
157  local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
158 
159  // If not a boundary condition
160  if (local_eqn >= 0)
161  {
162  // Add the contribution
163  veloc_mass_diag[local_eqn] += pow(psi[l], 2) * W;
164  }
165  } // for (unsigned i=0;i<n_dim;i++)
166  } // for (unsigned l=0;l<n_node;l++)
167  } // if ((which_one==0)||(which_one==2))
168 
169  // Do we want the pressure one?
170  if ((which_one == 0) || (which_one == 1))
171  {
172  // Get the pressure shape functions
173  pshape_nst(s, psi_p);
174 
175  // Loop over the veclocity shape functions
176  for (unsigned l = 0; l < n_pres; l++)
177  {
178  // Get equation number
179  local_eqn = p_local_eqn(l);
180 
181  // If not a boundary condition
182  if (local_eqn >= 0)
183  {
184  // Add the contribution
185  press_mass_diag[local_eqn] += pow(psi_p[l], 2) * W;
186  }
187  } // for (unsigned l=0;l<n_pres;l++)
188  } // if ((which_one==0)||(which_one==1))
189  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
190  } // End of get_pressure_and_velocity_mass_matrix_diagonal
191 
192 
193  //======================================================================
194  /// Compute the vector norm of FEM solution
195  //======================================================================
196  template<unsigned DIM>
198  Vector<double>& norm)
199  {
200  // Resize the solution norm vector
201  norm.resize(DIM + 1, 0.0);
202 
203  // Vector of local coordinates
204  Vector<double> s(DIM + 1, 0.0);
205 
206  // Set the value of n_intpt
207  unsigned n_intpt = integral_pt()->nweight();
208 
209  // Loop over the integration points
210  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
211  {
212  // Assign values of s
213  for (unsigned i = 0; i < DIM + 1; i++)
214  {
215  s[i] = integral_pt()->knot(ipt, i);
216  }
217 
218  // Get the integral weight
219  double w = integral_pt()->weight(ipt);
220 
221  // Get Jacobian of mapping
222  double J = J_eulerian(s);
223 
224  // Premultiply the weights and the Jacobian
225  double W = w * J;
226 
227  // Compute the velocity norm
228  for (unsigned i = 0; i < DIM; i++)
229  {
230  // Update the solution norm value
231  norm[i] += pow(interpolated_u_nst(s, i), 2) * W;
232  }
233 
234  // Update the pressure norm value
235  norm[DIM] += pow(interpolated_p_nst(s), 2) * W;
236  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
237  } // End of compute_norm
238 
239 
240  //======================================================================
241  /// Validate against exact velocity solution at given time. Solution is
242  /// provided via function pointer. Plot at a given number of plot points
243  /// and compute L2 error and L2 norm of velocity solution over element.
244  //=======================================================================
245  template<unsigned DIM>
247  std::ostream& outfile,
249  const double& time,
250  double& error,
251  double& norm)
252  {
253  // Initialise the error norm value
254  error = 0.0;
255 
256  // Initialise the solution norm value
257  norm = 0.0;
258 
259  // Storage for the time value
260  double interpolated_t = 0.0;
261 
262  // Vector of local coordinates
263  Vector<double> s(DIM + 1, 0.0);
264 
265  // Vector for the spatial coordinates
266  Vector<double> x(DIM, 0.0);
267 
268  // Set the value of n_intpt
269  unsigned n_intpt = integral_pt()->nweight();
270 
271  // Exact solution Vector (u,v,[w],p)
272  Vector<double> exact_soln(DIM + 1, 0.0);
273 
274  // Loop over the integration points
275  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
276  {
277  // Assign values of s
278  for (unsigned i = 0; i < DIM + 1; i++)
279  {
280  s[i] = integral_pt()->knot(ipt, i);
281  }
282 
283  // Assign values of x
284  for (unsigned i = 0; i < DIM; i++)
285  {
286  // Get x position as Vector
287  x[i] = interpolated_x(s, i);
288  }
289 
290  // Get the time value
291  interpolated_t = interpolated_x(s, DIM);
292 
293  // Get exact solution at this point
294  (*exact_soln_pt)(interpolated_t, x, exact_soln);
295 
296  // Get the integral weight
297  double w = integral_pt()->weight(ipt);
298 
299  // Get Jacobian of mapping
300  double J = J_eulerian(s);
301 
302  // Premultiply the weights and the Jacobian
303  double W = w * J;
304 
305  // Velocity error
306  for (unsigned i = 0; i < DIM; i++)
307  {
308  // Update the solution norm value
309  norm += exact_soln[i] * exact_soln[i] * W;
310 
311  // Update the error norm value
312  error += pow(exact_soln[i] - interpolated_u_nst(s, i), 2) * W;
313  }
314 
315  // ------ DRAIG: REMOVE ----------------------------------------
316  // Update the solution norm value
317  norm += pow(exact_soln[DIM], 2) * W;
318 
319  // Update the error norm value
320  error += pow(exact_soln[DIM] - interpolated_p_nst(s), 2) * W;
321  // ------ DRAIG: REMOVE ----------------------------------------
322  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
323 
324  //------------------------------------------------
325  // Output the error at the appropriate plot points
326  //------------------------------------------------
327  // Just output at the default number of plot points
328  unsigned n_plot = 5;
329 
330  // Tecplot header info
331  outfile << tecplot_zone_string(n_plot);
332 
333  // How many plot points are there in total?
334  unsigned num_plot_points = nplot_points(n_plot);
335 
336  // Loop over plot points
337  for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
338  {
339  // Get local coordinates of plot point
340  get_s_plot(i_plot, n_plot, s);
341 
342  // Loop over the spatial coordinates
343  for (unsigned i = 0; i < DIM; i++)
344  {
345  // Assign the i-th spatial coordinate
346  x[i] = interpolated_x(s, i);
347  }
348 
349  // Get the time value
350  interpolated_t = interpolated_x(s, DIM);
351 
352  // Get exact solution at this point
353  (*exact_soln_pt)(interpolated_t, x, exact_soln);
354 
355  // Output x,y,...
356  for (unsigned i = 0; i < DIM; i++)
357  {
358  // Output the i-th spatial coordinate
359  outfile << x[i] << " ";
360  }
361 
362  // Output the time value at this point
363  outfile << interpolated_t << " ";
364 
365  // Output u_error,v_error
366  for (unsigned i = 0; i < DIM; i++)
367  {
368  // Output the error in the i-th velocity component at this point
369  outfile << exact_soln[i] - interpolated_u_nst(s, i) << " ";
370  }
371 
372  // Output the error in the pressure field at this point
373  outfile << exact_soln[DIM] - interpolated_p_nst(s) << " ";
374 
375  // End the line
376  outfile << std::endl;
377  } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
378 
379  // Write tecplot footer (e.g. FE connectivity lists)
380  write_tecplot_zone_footer(outfile, n_plot);
381  } // End of compute_error
382 
383 
384  //======================================================================
385  /// Validate against exact velocity solution at given time. Solution is
386  /// provided via function pointer. Plot at a given number of plot points
387  /// and compute L2 error and L2 norm of velocity solution over element.
388  //=======================================================================
389  template<unsigned DIM>
391  std::ostream& outfile,
393  const double& time,
394  Vector<double>& error,
395  Vector<double>& norm)
396  {
397  // Resize the error norm vector
398  error.resize(DIM + 1, 0.0);
399 
400  // Resize the solution norm vector
401  norm.resize(DIM + 1, 0.0);
402 
403  // Storage for the time value
404  double interpolated_t = 0.0;
405 
406  // Vector of local coordinates
407  Vector<double> s(DIM + 1, 0.0);
408 
409  // Vector for the spatial coordinates
410  Vector<double> x(DIM, 0.0);
411 
412  // Set the value of n_intpt
413  unsigned n_intpt = integral_pt()->nweight();
414 
415  // Exact solution Vector (u,v,[w],p)
416  Vector<double> exact_soln(DIM + 1, 0.0);
417 
418  // Loop over the integration points
419  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
420  {
421  // Assign values of s
422  for (unsigned i = 0; i < DIM + 1; i++)
423  {
424  s[i] = integral_pt()->knot(ipt, i);
425  }
426 
427  // Assign values of x
428  for (unsigned i = 0; i < DIM; i++)
429  {
430  // Get x position as Vector
431  x[i] = interpolated_x(s, i);
432  }
433 
434  // Get the time value
435  interpolated_t = interpolated_x(s, DIM);
436 
437  // Get exact solution at this point
438  (*exact_soln_pt)(interpolated_t, x, exact_soln);
439 
440  // Get the integral weight
441  double w = integral_pt()->weight(ipt);
442 
443  // Get Jacobian of mapping
444  double J = J_eulerian(s);
445 
446  // Premultiply the weights and the Jacobian
447  double W = w * J;
448 
449  // Velocity error
450  for (unsigned i = 0; i < DIM; i++)
451  {
452  // Update the solution norm value
453  norm[i] += exact_soln[i] * exact_soln[i] * W;
454 
455  // Update the error norm value
456  error[i] += pow(exact_soln[i] - interpolated_u_nst(s, i), 2) * W;
457  }
458 
459  // ------ DRAIG: REMOVE ----------------------------------------
460  // Update the solution norm value
461  norm[DIM] += pow(exact_soln[DIM], 2) * W;
462 
463  // Update the error norm value
464  error[DIM] += pow(exact_soln[DIM] - interpolated_p_nst(s), 2) * W;
465  // ------ DRAIG: REMOVE ----------------------------------------
466  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
467 
468  //------------------------------------------------
469  // Output the error at the appropriate plot points
470  //------------------------------------------------
471  // Just output at the default number of plot points
472  unsigned n_plot = 5;
473 
474  // Tecplot header info
475  outfile << tecplot_zone_string(n_plot);
476 
477  // How many plot points are there in total?
478  unsigned num_plot_points = nplot_points(n_plot);
479 
480  // Loop over plot points
481  for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
482  {
483  // Get local coordinates of plot point
484  get_s_plot(i_plot, n_plot, s);
485 
486  // Loop over the spatial coordinates
487  for (unsigned i = 0; i < DIM; i++)
488  {
489  // Assign the i-th spatial coordinate
490  x[i] = interpolated_x(s, i);
491  }
492 
493  // Get the time value
494  interpolated_t = interpolated_x(s, DIM);
495 
496  // Get exact solution at this point
497  (*exact_soln_pt)(interpolated_t, x, exact_soln);
498 
499  // Output x,y,...
500  for (unsigned i = 0; i < DIM; i++)
501  {
502  // Output the i-th spatial coordinate
503  outfile << x[i] << " ";
504  }
505 
506  // Output the time value at this point
507  outfile << interpolated_t << " ";
508 
509  // Output u_error,v_error
510  for (unsigned i = 0; i < DIM; i++)
511  {
512  // Output the error in the i-th velocity component at this point
513  outfile << exact_soln[i] - interpolated_u_nst(s, i) << " ";
514  }
515 
516  // Output the error in the pressure field at this point
517  outfile << exact_soln[DIM] - interpolated_p_nst(s) << " ";
518 
519  // End the line
520  outfile << std::endl;
521  } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
522 
523  // Write tecplot footer (e.g. FE connectivity lists)
524  write_tecplot_zone_footer(outfile, n_plot);
525  } // End of compute_error
526 
527 
528  //======================================================================
529  /// Validate against exact velocity solution at given time.
530  /// Solution is provided via function pointer.
531  /// Compute L2 error and L2 norm of velocity solution over element.
532  //=======================================================================
533  template<unsigned DIM>
536  const double& time,
537  double& error,
538  double& norm)
539  {
540  // Initialise the error norm value
541  error = 0.0;
542 
543  // Initialise the solution norm value
544  norm = 0.0;
545 
546  // Storage for the time value
547  double interpolated_t = 0.0;
548 
549  // Vector of local coordinates
550  Vector<double> s(DIM + 1, 0.0);
551 
552  // Vector for the spatial coordinates
553  Vector<double> x(DIM, 0.0);
554 
555  // Set the value of n_intpt
556  unsigned n_intpt = integral_pt()->nweight();
557 
558  // Exact solution Vector (u,v,[w],p)
559  Vector<double> exact_soln(DIM + 1, 0.0);
560 
561  // Loop over the integration points
562  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
563  {
564  // Assign values of s
565  for (unsigned i = 0; i < DIM + 1; i++)
566  {
567  s[i] = integral_pt()->knot(ipt, i);
568  }
569 
570  // Assign values of x
571  for (unsigned i = 0; i < DIM; i++)
572  {
573  // Get x position as Vector
574  x[i] = interpolated_x(s, i);
575  }
576 
577  // Get the time value
578  interpolated_t = interpolated_x(s, DIM);
579 
580  // Get exact solution at this point
581  (*exact_soln_pt)(interpolated_t, x, exact_soln);
582 
583  // Get the integral weight
584  double w = integral_pt()->weight(ipt);
585 
586  // Get Jacobian of mapping
587  double J = J_eulerian(s);
588 
589  // Premultiply the weights and the Jacobian
590  double W = w * J;
591 
592  // Velocity error
593  for (unsigned i = 0; i < DIM; i++)
594  {
595  // Update the solution norm value
596  norm += exact_soln[i] * exact_soln[i] * W;
597 
598  // Update the error norm value
599  error += pow(exact_soln[i] - interpolated_u_nst(s, i), 2) * W;
600  }
601  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
602  } // End of compute_error
603 
604  //======================================================================
605  /// Validate against exact velocity solution Solution is provided via a
606  /// function pointer. Compute L2 error and L2 norm of velocity solution
607  /// over element.
608  //=======================================================================
609  template<unsigned DIM>
612  double& error,
613  double& norm)
614  {
615  // Initialise the error norm value
616  error = 0.0;
617 
618  // Initialise the solution norm value
619  norm = 0.0;
620 
621  // Vector of local coordinates
622  Vector<double> s(DIM + 1, 0.0);
623 
624  // Vector for the spatial coordinates
625  Vector<double> x(DIM, 0.0);
626 
627  // Set the value of n_intpt
628  unsigned n_intpt = integral_pt()->nweight();
629 
630  // Exact solution Vector (u,v,[w],p)
631  Vector<double> exact_soln(DIM + 1, 0.0);
632 
633  // Loop over the integration points
634  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
635  {
636  // Assign values of s
637  for (unsigned i = 0; i < DIM + 1; i++)
638  {
639  s[i] = integral_pt()->knot(ipt, i);
640  }
641 
642  // Assign values of x
643  for (unsigned i = 0; i < DIM; i++)
644  {
645  // Get x position as Vector
646  x[i] = interpolated_x(s, i);
647  }
648 
649  // Get exact solution at this point
650  (*exact_soln_pt)(x, exact_soln);
651 
652  // Get the integral weight
653  double w = integral_pt()->weight(ipt);
654 
655  // Get Jacobian of mapping
656  double J = J_eulerian(s);
657 
658  // Premultiply the weights and the Jacobian
659  double W = w * J;
660 
661  // Velocity error
662  for (unsigned i = 0; i < DIM; i++)
663  {
664  // Update the solution norm value
665  norm += exact_soln[i] * exact_soln[i] * W;
666 
667  // Update the error norm value
668  error += pow(exact_soln[i] - interpolated_u_nst(s, i), 2) * W;
669  }
670  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
671  } // End of compute_error
672 
673  //======================================================================
674  /// Validate against exact velocity solution
675  /// Solution is provided via function pointer.
676  /// Plot at a given number of plot points and compute L2 error
677  /// and L2 norm of velocity solution over element.
678  //=======================================================================
679  template<unsigned DIM>
681  std::ostream& outfile,
683  double& error,
684  double& norm)
685  {
686  // Initialise the error norm value
687  error = 0.0;
688 
689  // Initialise the solution norm value
690  norm = 0.0;
691 
692  // Storage for the time value
693  double interpolated_t = 0.0;
694 
695  // Vector of local coordinates
696  Vector<double> s(DIM + 1, 0.0);
697 
698  // Vector for the spatial coordinates
699  Vector<double> x(DIM, 0.0);
700 
701  // Set the value of n_intpt
702  unsigned n_intpt = integral_pt()->nweight();
703 
704  // Output the tecplot header
705  outfile << "ZONE" << std::endl;
706 
707  // Exact solution Vector (u,v,[w],p)
708  Vector<double> exact_soln(DIM + 1, 0.0);
709 
710  // Loop over the integration points
711  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
712  {
713  // Assign values of s
714  for (unsigned i = 0; i < DIM + 1; i++)
715  {
716  s[i] = integral_pt()->knot(ipt, i);
717  }
718 
719  // Assign values of x
720  for (unsigned i = 0; i < DIM; i++)
721  {
722  // Get x position as Vector
723  x[i] = interpolated_x(s, i);
724  }
725 
726  // Get the time value
727  interpolated_t = interpolated_x(s, DIM);
728 
729  // Get exact solution at this point
730  (*exact_soln_pt)(x, exact_soln);
731 
732  // Get the integral weight
733  double w = integral_pt()->weight(ipt);
734 
735  // Get Jacobian of mapping
736  double J = J_eulerian(s);
737 
738  // Premultiply the weights and the Jacobian
739  double W = w * J;
740 
741  // Velocity error
742  for (unsigned i = 0; i < DIM; i++)
743  {
744  // Update the solution norm value
745  norm += exact_soln[i] * exact_soln[i] * W;
746 
747  // Update the error norm value
748  error += pow(exact_soln[i] - interpolated_u_nst(s, i), 2) * W;
749  }
750 
751  // Output x,y,...,u_exact
752  for (unsigned i = 0; i < DIM; i++)
753  {
754  // Output the i-th coordinate value
755  outfile << x[i] << " ";
756  }
757 
758  // Output the time value at this point
759  outfile << interpolated_t << " ";
760 
761  // Output x,y,u_error,v_error
762  for (unsigned i = 0; i < DIM; i++)
763  {
764  // Output the error in the i-th velocity component at this point
765  outfile << exact_soln[i] - interpolated_u_nst(s, i) << " ";
766  }
767 
768  // Finish the line off
769  outfile << std::endl;
770  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
771  } // End of compute_error
772 
773  //======================================================================
774  /// Output "exact" solution
775  /// Solution is provided via function pointer.
776  /// Plot at a given number of plot points.
777  /// Function prints as many components as are returned in solution Vector.
778  //=======================================================================
779  template<unsigned DIM>
781  std::ostream& outfile,
782  const unsigned& n_plot,
784  {
785  // Storage for the time value
786  double interpolated_t = 0.0;
787 
788  // Vector of local coordinates
789  Vector<double> s(DIM + 1, 0.0);
790 
791  // Vector for coordinates
792  Vector<double> x(DIM, 0.0);
793 
794  // Tecplot header info
795  outfile << tecplot_zone_string(n_plot);
796 
797  // Exact solution Vector
798  Vector<double> exact_soln;
799 
800  // How many plot points are there in total?
801  unsigned num_plot_points = nplot_points(n_plot);
802 
803  // Loop over plot points
804  for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
805  {
806  // Get local coordinates of plot point
807  get_s_plot(i_plot, n_plot, s);
808 
809  // Loop over the spatial coordinates
810  for (unsigned i = 0; i < DIM; i++)
811  {
812  // Assign the i-th spatial coordinate
813  x[i] = interpolated_x(s, i);
814  }
815 
816  // Get the time value
817  interpolated_t = interpolated_x(s, DIM);
818 
819  // Get exact solution at this point
820  (*exact_soln_pt)(x, exact_soln);
821 
822  // Output x,y,...
823  for (unsigned i = 0; i < DIM; i++)
824  {
825  // Output the i-th spatial coordinate
826  outfile << x[i] << " ";
827  }
828 
829  // Output the time value at this point
830  outfile << interpolated_t << " ";
831 
832  // Output "exact solution"
833  for (unsigned i = 0; i < exact_soln.size(); i++)
834  {
835  // Output the i-th (exact) velocity component value
836  outfile << exact_soln[i] << " ";
837  }
838 
839  // End the line
840  outfile << std::endl;
841  } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
842 
843  // Write tecplot footer (e.g. FE connectivity lists)
844  write_tecplot_zone_footer(outfile, n_plot);
845  } // End of output_fct
846 
847  //======================================================================
848  /// Output "exact" solution at a given time
849  /// Solution is provided via function pointer.
850  /// Plot at a given number of plot points.
851  /// Function prints as many components as are returned in solution Vector.
852  //=======================================================================
853  template<unsigned DIM>
855  std::ostream& outfile,
856  const unsigned& n_plot,
857  const double& time,
859  {
860  // Storage for the time value
861  double interpolated_t = 0.0;
862 
863  // Vector of local coordinates
864  Vector<double> s(DIM + 1, 0.0);
865 
866  // Vector for coordinates
867  Vector<double> x(DIM, 0.0);
868 
869  // Tecplot header info
870  outfile << tecplot_zone_string(n_plot);
871 
872  // Exact solution Vector
873  Vector<double> exact_soln;
874 
875  // How many plot points are there in total?
876  unsigned num_plot_points = nplot_points(n_plot);
877 
878  // Loop over plot points
879  for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
880  {
881  // Get local coordinates of plot point
882  get_s_plot(i_plot, n_plot, s);
883 
884  // Loop over the spatial coordinates
885  for (unsigned i = 0; i < DIM; i++)
886  {
887  // Assign the i-th spatial coordinate
888  x[i] = interpolated_x(s, i);
889  }
890 
891  // Get the time value
892  interpolated_t = interpolated_x(s, DIM);
893 
894  // Get exact solution at this point
895  (*exact_soln_pt)(interpolated_t, x, exact_soln);
896 
897  // Output x,y,...
898  for (unsigned i = 0; i < DIM; i++)
899  {
900  // Output the i-th spatial coordinate value
901  outfile << x[i] << " ";
902  }
903 
904  // Output the time value at this point
905  outfile << interpolated_t << " ";
906 
907  // Output "exact solution"
908  for (unsigned i = 0; i < exact_soln.size(); i++)
909  {
910  // Output the i-th (exact) velocity component value
911  outfile << exact_soln[i] << " ";
912  }
913 
914  // End the line
915  outfile << std::endl;
916  } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
917 
918  // Write tecplot footer (e.g. FE connectivity lists)
919  write_tecplot_zone_footer(outfile, n_plot);
920  } // End of output_fct
921 
922  //==============================================================
923  /// Output function: Velocities only
924  /// x,y,[z],u,v,[w]
925  /// in tecplot format at specified previous timestep (t=0: present;
926  /// t>0: previous timestep). Specified number of plot points in each
927  /// coordinate direction.
928  /// DRAIG: Should be broken!
929  //==============================================================
930  template<unsigned DIM>
932  std::ostream& outfile, const unsigned& n_plot, const unsigned& t)
933  {
934  // Find number of nodes
935  unsigned n_node = nnode();
936 
937  // Local shape function
938  Shape psi(n_node);
939 
940  // Vectors of local coordinates
941  Vector<double> s(DIM + 1, 0.0);
942 
943  // Tecplot header info
944  outfile << tecplot_zone_string(n_plot);
945 
946  // How many plot points are there?
947  unsigned num_plot_points = nplot_points(n_plot);
948 
949  // Loop over plot points
950  for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
951  {
952  // Get the local coordinates of plot point
953  get_s_plot(i_plot, n_plot, s);
954 
955  // Coordinates
956  for (unsigned i = 0; i < DIM + 1; i++)
957  {
958  // Output the i-th spatial coordinate value
959  outfile << interpolated_x(s, i) << " ";
960  }
961 
962  // Velocities
963  for (unsigned i = 0; i < DIM; i++)
964  {
965  // Output the i-th velocity component
966  outfile << interpolated_u_nst(s, i) << " ";
967  }
968 
969  // End the line
970  outfile << std::endl;
971  } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
972 
973  // Write tecplot footer (e.g. FE connectivity lists)
974  write_tecplot_zone_footer(outfile, n_plot);
975  } // End of output_veloc
976 
977  //==============================================================
978  /// Output function:
979  /// x,y,[z],u,v,[w],p
980  /// in tecplot format. Specified number of plot points in each
981  /// coordinate direction.
982  //==============================================================
983  template<unsigned DIM>
985  std::ostream& outfile, const unsigned& n_plot)
986  {
987  // Find number of nodes
988  unsigned n_node = nnode();
989 
990  // Local shape function
991  Shape psi(n_node);
992 
993  // Vectors of local coordinates
994  Vector<double> s(DIM + 1, 0.0);
995 
996  // Tecplot header info
997  outfile << tecplot_zone_string(n_plot);
998 
999  // How many plot points are there?
1000  unsigned num_plot_points = nplot_points(n_plot);
1001 
1002  // Loop over plot points
1003  for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1004  {
1005  // Get the local coordinates of plot point
1006  get_s_plot(i_plot, n_plot, s);
1007 
1008  // Coordinates
1009  for (unsigned i = 0; i < DIM + 1; i++)
1010  {
1011  // Output the i-th spatial coordinate value
1012  outfile << interpolated_x(s, i) << " ";
1013  }
1014 
1015  // Velocities
1016  for (unsigned i = 0; i < DIM; i++)
1017  {
1018  // Output the i-th velocity component
1019  outfile << interpolated_u_nst(s, i) << " ";
1020  }
1021 
1022  // Pressure
1023  outfile << interpolated_p_nst(s) << " ";
1024 
1025  // End the line
1026  outfile << std::endl;
1027  } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
1028 
1029  // Add an extra line
1030  outfile << std::endl;
1031 
1032  // Write tecplot footer (e.g. FE connectivity lists)
1033  write_tecplot_zone_footer(outfile, n_plot);
1034  } // End of output
1035 
1036  //==============================================================
1037  /// C-style output function:
1038  /// x,y,[z],u,v,[w],p
1039  /// in tecplot format. Specified number of plot points in each
1040  /// coordinate direction.
1041  //==============================================================
1042  template<unsigned DIM>
1044  FILE* file_pt, const unsigned& n_plot)
1045  {
1046  // Vector of local coordinates
1047  Vector<double> s(DIM + 1, 0.0);
1048 
1049  // Tecplot header info
1050  fprintf(file_pt, "%s", tecplot_zone_string(n_plot).c_str());
1051 
1052  // How many plot points in total?
1053  unsigned num_plot_points = nplot_points(n_plot);
1054 
1055  // Loop over plot points
1056  for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1057  {
1058  // Get the local coordinates associated with this plot point
1059  get_s_plot(i_plot, n_plot, s);
1060 
1061  // Loop over the coordinate directions
1062  for (unsigned i = 0; i < DIM + 1; i++)
1063  {
1064  // Output the i-th coordinate value to file
1065  fprintf(file_pt, "%g ", interpolated_x(s, i));
1066  }
1067 
1068  // Loop over the velocity components
1069  for (unsigned i = 0; i < DIM; i++)
1070  {
1071  // Output the i-th velocity component to file
1072  fprintf(file_pt, "%g ", interpolated_u_nst(s, i));
1073  }
1074 
1075  // Pressure
1076  fprintf(file_pt, "%g \n", interpolated_p_nst(s));
1077  }
1078 
1079  // Finish the line off
1080  fprintf(file_pt, "\n");
1081 
1082  // Write tecplot footer (e.g. FE connectivity lists)
1083  write_tecplot_zone_footer(file_pt, n_plot);
1084  } // End of output
1085 
1086  //==============================================================
1087  /// Full output function:
1088  /// x,y,t,u,v,p,du/dt,dv/dt,dissipation
1089  /// in tecplot format. Specified number of plot points in each
1090  /// coordinate direction
1091  //==============================================================
1092  template<unsigned DIM>
1094  std::ostream& outfile, const unsigned& n_plot)
1095  {
1096  // Vector of local coordinates
1097  Vector<double> s(DIM + 1, 0.0);
1098 
1099  // Tecplot header info
1100  outfile << tecplot_zone_string(n_plot);
1101 
1102  // Find out how many nodes there are
1103  unsigned n_node = nnode();
1104 
1105  // Set up memory for the shape functions
1106  Shape psif(n_node);
1107 
1108  // Set up memory for the shape function derivatives
1109  DShape dpsifdx(n_node, DIM);
1110 
1111  // How many plot points in total?
1112  unsigned num_plot_points = nplot_points(n_plot);
1113 
1114  // Loop over plot points
1115  for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1116  {
1117  // Get the local coordinates of plot point
1118  get_s_plot(i_plot, n_plot, s);
1119 
1120  // DRAIG: Fix
1121  throw OomphLibError("This needs to be fixed...",
1122  OOMPH_CURRENT_FUNCTION,
1123  OOMPH_EXCEPTION_LOCATION);
1124 
1125  // Call the derivatives of the shape and test functions
1126  dshape_eulerian(s, psif, dpsifdx);
1127 
1128  // Allocate storage for the mesh velocity
1129  Vector<double> mesh_velocity(DIM, 0.0);
1130 
1131  // Allocate storage for the acceleration
1132  Vector<double> du_dt(DIM, 0.0);
1133 
1134  // Allocate storage for the ALE acceleration
1135  Vector<double> du_dt_ALE(DIM, 0.0);
1136 
1137  // Allocate storage for the velocity derivatives
1138  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1139 
1140  //--------------------------------------
1141  // Calculate velocities and derivatives:
1142  //--------------------------------------
1143  // Loop over directions
1144  for (unsigned i = 0; i < DIM; i++)
1145  {
1146  // Get the index at which velocity i is stored
1147  unsigned u_nodal_index = u_index_nst(i);
1148 
1149  // Loop over nodes
1150  for (unsigned l = 0; l < n_node; l++)
1151  {
1152  // Get the nodal value
1153  double u_value = nodal_value(l, u_nodal_index);
1154 
1155  // Update the i-th acceleration component
1156  du_dt[i] += u_value * dpsifdx(l, DIM);
1157 
1158  // Update the i-th mesh velocity component
1159  mesh_velocity[i] += nodal_position(l, i) * dpsifdx(l, DIM);
1160 
1161  // Loop over derivative directions for velocity gradients
1162  for (unsigned j = 0; j < DIM; j++)
1163  {
1164  // Update the value of du_i/dx_j
1165  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1166  }
1167  } // for (unsigned l=0;l<n_node;l++)
1168  } // for (unsigned i=0;i<DIM;i++)
1169 
1170  //---------------------------------------------
1171  // Get du/dt in ALE form (incl. mesh velocity):
1172  //---------------------------------------------
1173  // Loop over the coordinate directions
1174  for (unsigned i = 0; i < DIM; i++)
1175  {
1176  // Store the i-th acceleration component
1177  du_dt_ALE[i] = du_dt[i];
1178 
1179  // Loop over the coordinate directions
1180  for (unsigned k = 0; k < DIM; k++)
1181  {
1182  // Take the mesh velocity into account
1183  du_dt_ALE[i] -= mesh_velocity[k] * interpolated_dudx(i, k);
1184  }
1185  } // for (unsigned i=0;i<DIM;i++)
1186 
1187  // Output the coordinates
1188  for (unsigned i = 0; i < DIM + 1; i++)
1189  {
1190  // Output the i-th coordinate value
1191  outfile << interpolated_x(s, i) << " ";
1192  }
1193 
1194  // Output the velocity components
1195  for (unsigned i = 0; i < DIM; i++)
1196  {
1197  // Output the i-th velocity component
1198  outfile << interpolated_u_nst(s, i) << " ";
1199  }
1200 
1201  // Output the pressure
1202  outfile << interpolated_p_nst(s) << " ";
1203 
1204  // Output the acceleration
1205  for (unsigned i = 0; i < DIM; i++)
1206  {
1207  // Output the ALE acceleration term
1208  outfile << du_dt_ALE[i] << " ";
1209  }
1210 
1211  // Dissipation
1212  outfile << dissipation(s) << " ";
1213 
1214  // End the line
1215  outfile << std::endl;
1216  } // for (unsigned i_plot=0;i_plot<num_plot_points;i_plot++)
1217 
1218  // Write tecplot footer (e.g. FE connectivity lists)
1219  write_tecplot_zone_footer(outfile, n_plot);
1220  } // End of full_output
1221 
1222 
1223  //==============================================================
1224  /// Output function for vorticity.
1225  /// x,y,[z],[omega_x,omega_y,[and/or omega_z]]
1226  /// in tecplot format. Specified number of plot points in each
1227  /// coordinate direction.
1228  //==============================================================
1229  template<unsigned DIM>
1231  std::ostream& outfile, const unsigned& n_plot)
1232  {
1233  // Vector of local coordinates
1234  Vector<double> s(DIM + 1, 0.0);
1235 
1236  // Create vorticity vector of the required size
1237  Vector<double> vorticity;
1238 
1239  // If we're in 2D the vorticity field is a scalar field
1240  if (DIM == 2)
1241  {
1242  // Resize the vorticity vector
1243  vorticity.resize(1);
1244  }
1245  // If we're in 3D the vorticity field is a vector field
1246  else if (DIM == 3)
1247  {
1248  // Resize the vorticity vector
1249  vorticity.resize(3);
1250  }
1251  // If we're in 1D
1252  else
1253  {
1254  std::string error_message =
1255  "Can't output vorticity in 1D - in fact, what ";
1256  error_message += "do you mean\nby the 1D Navier-Stokes equations?\n";
1257  throw OomphLibError(
1258  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1259  }
1260 
1261  // Tecplot header info
1262  outfile << tecplot_zone_string(n_plot);
1263 
1264  // How many plot points are there in total?
1265  unsigned num_plot_points = nplot_points(n_plot);
1266 
1267  // Loop over plot points
1268  for (unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1269  {
1270  // Get local coordinates of plot point
1271  get_s_plot(i_plot, n_plot, s);
1272 
1273  // Coordinates
1274  for (unsigned i = 0; i < DIM + 1; i++)
1275  {
1276  // Output the i-th coordinate value
1277  outfile << interpolated_x(s, i) << " ";
1278  }
1279 
1280  // Get vorticity
1281  get_vorticity(s, vorticity);
1282 
1283  // If we're in 2D
1284  if (DIM == 2)
1285  {
1286  // Output the vorticity field value
1287  outfile << vorticity[0];
1288  }
1289  // If we're in 3D
1290  else
1291  {
1292  // Output the vorticity field
1293  outfile << vorticity[0] << " " << vorticity[1] << " " << vorticity[2]
1294  << " ";
1295  }
1296 
1297  // Finish the line off
1298  outfile << std::endl;
1299  }
1300 
1301  // Add in an extra line
1302  outfile << std::endl;
1303 
1304  // Write tecplot footer (e.g. FE connectivity lists)
1305  write_tecplot_zone_footer(outfile, n_plot);
1306  } // End of output_vorticity
1307 
1308  //==============================================================
1309  /// Return integral of dissipation over element
1310  //==============================================================
1311  template<unsigned DIM>
1313  {
1314  // Initialise the elemental dissipation value
1315  double diss = 0.0;
1316 
1317  // Set the value of n_intpt
1318  unsigned n_intpt = integral_pt()->nweight();
1319 
1320  // Set the Vector to hold local coordinates
1321  Vector<double> s(DIM + 1, 0.0);
1322 
1323  // Loop over the integration points
1324  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1325  {
1326  // Assign values of s
1327  for (unsigned i = 0; i < DIM + 1; i++)
1328  {
1329  // Calculate the i-th local coordinate value
1330  s[i] = integral_pt()->knot(ipt, i);
1331  }
1332 
1333  // Get the integral weight
1334  double w = integral_pt()->weight(ipt);
1335 
1336  // Get Jacobian of mapping
1337  double J = J_eulerian(s);
1338 
1339  // Storage for the strain rate matrix
1340  DenseMatrix<double> strainrate(DIM, DIM);
1341 
1342  // Get strain rate matrix
1343  strain_rate(s, strainrate);
1344 
1345  // Initialise the local dissipation
1346  double local_diss = 0.0;
1347 
1348  // Loop over the coordinate directions
1349  for (unsigned i = 0; i < DIM; i++)
1350  {
1351  // Loop over the coordinate directions
1352  for (unsigned j = 0; j < DIM; j++)
1353  {
1354  // Update the local dissipation value
1355  local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
1356  }
1357  } // for (unsigned i=0;i<DIM;i++)
1358 
1359  // Update the elemental dissipation value
1360  diss += local_diss * w * J;
1361  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
1362 
1363  // Return the elemental dissipation value
1364  return diss;
1365  } // End of dissipation
1366 
1367  //==============================================================
1368  /// Compute traction (on the viscous scale) exerted onto
1369  /// the fluid at local coordinate s. N has to be outer unit normal
1370  /// to the fluid.
1371  //==============================================================
1372  template<unsigned DIM>
1374  const Vector<double>& s, const Vector<double>& N, Vector<double>& traction)
1375  {
1376  // Allocate space for the strain rate matrix
1377  DenseMatrix<double> strainrate(DIM, DIM);
1378 
1379  // Get velocity gradients
1380  strain_rate(s, strainrate);
1381 
1382  // Get pressure
1383  double press = interpolated_p_nst(s);
1384 
1385  // Loop over traction components
1386  for (unsigned i = 0; i < DIM; i++)
1387  {
1388  // Add in the pressure contribution
1389  traction[i] = -press * N[i];
1390 
1391  // Loop over the strain rate entries
1392  for (unsigned j = 0; j < DIM; j++)
1393  {
1394  // Add in the strain rate contribution
1395  traction[i] += 2.0 * strainrate(i, j) * N[j];
1396  }
1397  } // for (unsigned i=0;i<DIM;i++)
1398  } // End of get_traction
1399 
1400  //==============================================================
1401  /// Compute traction (on the viscous scale) exerted onto
1402  /// the fluid at local coordinate s, decomposed into pressure and
1403  /// normal and tangential viscous components.
1404  /// N has to be outer unit normal to the fluid.
1405  //==============================================================
1406  template<unsigned DIM>
1408  const Vector<double>& s,
1409  const Vector<double>& N,
1410  Vector<double>& traction_p,
1411  Vector<double>& traction_visc_n,
1412  Vector<double>& traction_visc_t)
1413  {
1414  // Allocate space for the traction components
1415  Vector<double> traction_visc(DIM);
1416 
1417  // Allocate space for the velocity gradients
1418  DenseMatrix<double> strainrate(DIM, DIM);
1419 
1420  // Calculate the strain rate at the local coordinate s
1421  strain_rate(s, strainrate);
1422 
1423  // Get pressure
1424  double press = interpolated_p_nst(s);
1425 
1426  // Loop over traction components
1427  for (unsigned i = 0; i < DIM; i++)
1428  {
1429  // Get the pressure contribution
1430  traction_p[i] = -press * N[i];
1431 
1432  // Loop over the coordinate directions
1433  for (unsigned j = 0; j < DIM; j++)
1434  {
1435  // Add in the viscous stress contribution
1436  traction_visc[i] += 2.0 * strainrate(i, j) * N[j];
1437  }
1438 
1439  // Get the normal component of the viscous stress
1440  traction_visc_n[i] = traction_visc[i] * N[i];
1441 
1442  // Get the tangential component of the viscous stress
1443  traction_visc_t[i] = traction_visc[i] - traction_visc_n[i];
1444  } // for (unsigned i=0;i<DIM;i++)
1445  } // End of get_traction
1446 
1447  //==============================================================
1448  /// Return dissipation at local coordinate s
1449  //==============================================================
1450  template<unsigned DIM>
1452  const Vector<double>& s) const
1453  {
1454  // Get strain rate matrix
1455  DenseMatrix<double> strainrate(DIM, DIM);
1456  strain_rate(s, strainrate);
1457 
1458  // Initialise
1459  double local_diss = 0.0;
1460  for (unsigned i = 0; i < DIM; i++)
1461  {
1462  for (unsigned j = 0; j < DIM; j++)
1463  {
1464  local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
1465  }
1466  }
1467 
1468  return local_diss;
1469  }
1470 
1471  //==============================================================
1472  /// Get strain-rate tensor: (1/2)*(du_i/dx_j+du_j/dx_i)
1473  //==============================================================
1474  template<unsigned DIM>
1476  const Vector<double>& s, DenseMatrix<double>& strainrate) const
1477  {
1478 #ifdef PARANOID
1479  if ((strainrate.ncol() != DIM) || (strainrate.nrow() != DIM))
1480  {
1481  std::ostringstream error_message;
1482  error_message << "The strain rate has incorrect dimensions "
1483  << strainrate.ncol() << " x " << strainrate.nrow()
1484  << " not " << DIM << std::endl;
1485  throw OomphLibError(
1486  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1487  }
1488 #endif
1489 
1490  // Velocity gradient matrix
1491  DenseMatrix<double> dudx(DIM);
1492 
1493  // Find out how many nodes there are in the element
1494  unsigned n_node = nnode();
1495 
1496  // Set up memory for the shape and test functions
1497  Shape psi(n_node);
1498 
1499  // Set up memory for the shape and test function derivatives
1500  DShape dpsidx(n_node, DIM + 1);
1501 
1502  // Evaluate the shape functions and shape function derivatives at this point
1503  dshape_eulerian(s, psi, dpsidx);
1504 
1505  // Initialise to zero
1506  dudx.initialise(0.0);
1507 
1508  // Loop over velocity components
1509  for (unsigned i = 0; i < DIM; i++)
1510  {
1511  // Get the index at which the i-th velocity is stored
1512  unsigned u_nodal_index = u_index_nst(i);
1513 
1514  // Loop over derivative directions
1515  for (unsigned j = 0; j < DIM; j++)
1516  {
1517  // Loop over nodes
1518  for (unsigned l = 0; l < n_node; l++)
1519  {
1520  // Update the value of du_i/dx_j
1521  dudx(i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
1522  }
1523  } // for (unsigned j=0;j<DIM;j++)
1524  } // for (unsigned i=0;i<DIM;i++)
1525 
1526  // Loop over velocity components
1527  for (unsigned i = 0; i < DIM; i++)
1528  {
1529  // Loop over derivative directions
1530  for (unsigned j = 0; j < DIM; j++)
1531  {
1532  // Calculate the (i,j)-th strain rate entry
1533  strainrate(i, j) = 0.5 * (dudx(i, j) + dudx(j, i));
1534  }
1535  } // for (unsigned i=0;i<DIM;i++)
1536  } // End of strain_rate
1537 
1538  //==============================================================
1539  /// Compute 2D vorticity vector at local coordinate s (return in
1540  /// one and only component of vorticity vector
1541  //==============================================================
1542  template<>
1544  const Vector<double>& s, Vector<double>& vorticity) const
1545  {
1546 #ifdef PARANOID
1547  if (vorticity.size() != 1)
1548  {
1549  // Create an output stream
1550  std::ostringstream error_message;
1551 
1552  // Create an error message
1553  error_message << "Computation of vorticity in 2D requires a 1D vector\n"
1554  << "which contains the only non-zero component of the\n"
1555  << "vorticity vector. You've passed a vector of size "
1556  << vorticity.size() << std::endl;
1557 
1558  // Throw an error
1559  throw OomphLibError(
1560  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1561  }
1562 #endif
1563 
1564  // Specify spatial dimension
1565  unsigned dim = 2;
1566 
1567  // Velocity gradient matrix
1568  DenseMatrix<double> dudx(dim, dim, 0.0);
1569 
1570  // Find out how many nodes there are in the element
1571  unsigned n_node = nnode();
1572 
1573  // Set up memory for the shape and test functions
1574  Shape psi(n_node);
1575 
1576  // Set up memory for the shape and test function derivatives
1577  DShape dpsidx(n_node, dim + 1);
1578 
1579  // Evaluate the shape functions and shape function derivatives at this point
1580  dshape_eulerian(s, psi, dpsidx);
1581 
1582  // Initialise to zero
1583  dudx.initialise(0.0);
1584 
1585  // Loop over velocity components
1586  for (unsigned i = 0; i < dim; i++)
1587  {
1588  // Get the index at which the i-th velocity is stored
1589  unsigned u_nodal_index = u_index_nst(i);
1590 
1591  // Loop over derivative directions
1592  for (unsigned j = 0; j < dim; j++)
1593  {
1594  // Loop over nodes
1595  for (unsigned l = 0; l < n_node; l++)
1596  {
1597  // Update the value of du_i/dx_j
1598  dudx(i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
1599  }
1600  } // for (unsigned j=0;j<DIM;j++)
1601  } // for (unsigned i=0;i<DIM;i++)
1602 
1603  // Z-component of vorticity
1604  vorticity[0] = dudx(1, 0) - dudx(0, 1);
1605  } // End of get_vorticity
1606 
1607  //==============================================================
1608  /// Compute 2D vorticity vector at local coordinate s (return in
1609  /// one and only component of vorticity vector as a double
1610  //==============================================================
1611  template<>
1613  const Vector<double>& s, double& vorticity) const
1614  {
1615  // Create a vector to store the vorticity
1616  Vector<double> vort(1, 0.0);
1617 
1618  // Calculate the vorticity
1619  get_vorticity(s, vort);
1620 
1621  // Assign the vorticity
1622  vorticity = vort[0];
1623  } // End of get_vorticity
1624 
1625  //==============================================================
1626  /// Get integral of kinetic energy over element:
1627  /// Note that this is the "raw" kinetic energy in the sense
1628  /// that the density ratio has not been included. In problems
1629  /// with two or more fluids the user will have to remember
1630  /// to premultiply certain elements by the appropriate density
1631  /// ratio.
1632  //==============================================================
1633  template<unsigned DIM>
1635  {
1636  // Initialise the elemental kinetic energy value
1637  double kin_en = 0.0;
1638 
1639  // Set the value of n_intpt
1640  unsigned n_intpt = integral_pt()->nweight();
1641 
1642  // Set the Vector to hold local coordinates
1643  Vector<double> s(DIM + 1, 0.0);
1644 
1645  // Loop over the integration points
1646  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1647  {
1648  // Assign values of s
1649  for (unsigned i = 0; i < DIM + 1; i++)
1650  {
1651  s[i] = integral_pt()->knot(ipt, i);
1652  }
1653 
1654  // Get the integral weight
1655  double w = integral_pt()->weight(ipt);
1656 
1657  // Get Jacobian of mapping
1658  double J = J_eulerian(s);
1659 
1660  // Loop over directions
1661  double veloc_squared = 0.0;
1662  for (unsigned i = 0; i < DIM; i++)
1663  {
1664  veloc_squared += interpolated_u_nst(s, i) * interpolated_u_nst(s, i);
1665  }
1666 
1667  kin_en += 0.5 * veloc_squared * w * J;
1668  }
1669 
1670  return kin_en;
1671 
1672  } // End of kin_energy
1673 
1674 
1675  //==========================================================================
1676  /// Get integral of time derivative of kinetic energy over element:
1677  //==========================================================================
1678  template<unsigned DIM>
1680  {
1681  // Initialise
1682  double d_kin_en_dt = 0.0;
1683 
1684  // Set the value of n_intpt
1685  const unsigned n_intpt = integral_pt()->nweight();
1686 
1687  // Get the number of nodes
1688  const unsigned n_node = this->nnode();
1689 
1690  // Storage for the shape function
1691  Shape psi(n_node);
1692  DShape dpsidx(n_node, DIM + 1);
1693 
1694  // Get the value at which the velocities are stored
1695  unsigned u_index[DIM];
1696  for (unsigned i = 0; i < DIM; i++)
1697  {
1698  u_index[i] = this->u_index_nst(i);
1699  }
1700 
1701  // Loop over the integration points
1702  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1703  {
1704  // Get the jacobian of the mapping
1705  double J = dshape_eulerian_at_knot(ipt, psi, dpsidx);
1706 
1707  // Get the integral weight
1708  double w = integral_pt()->weight(ipt);
1709 
1710  // Now work out the velocity and the time derivative
1711  Vector<double> interpolated_u(DIM, 0.0);
1712  Vector<double> interpolated_dudt(DIM, 0.0);
1713 
1714  // Loop over the shape functions
1715  for (unsigned l = 0; l < n_node; l++)
1716  {
1717  // Loop over the dimensions
1718  for (unsigned i = 0; i < DIM; i++)
1719  {
1720  interpolated_u[i] += nodal_value(l, u_index[i]) * psi(l);
1721  interpolated_dudt[i] += du_dt_nst(l, u_index[i]) * psi(l);
1722  }
1723  }
1724 
1725  // Get mesh velocity bit
1726  if (!ALE_is_disabled)
1727  {
1728  Vector<double> mesh_velocity(DIM, 0.0);
1729  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1730 
1731  // Loop over nodes again
1732  for (unsigned l = 0; l < n_node; l++)
1733  {
1734  for (unsigned i = 0; i < DIM; i++)
1735  {
1736  mesh_velocity[i] += this->dnodal_position_dt(l, i) * psi(l);
1737 
1738  for (unsigned j = 0; j < DIM; j++)
1739  {
1740  interpolated_dudx(i, j) +=
1741  this->nodal_value(l, u_index[i]) * dpsidx(l, j);
1742  }
1743  }
1744  }
1745 
1746  // Subtract mesh velocity from du_dt
1747  for (unsigned i = 0; i < DIM; i++)
1748  {
1749  for (unsigned k = 0; k < DIM; k++)
1750  {
1751  interpolated_dudt[i] -= mesh_velocity[k] * interpolated_dudx(i, k);
1752  }
1753  }
1754  }
1755 
1756 
1757  // Loop over directions and add up u du/dt terms
1758  double sum = 0.0;
1759  for (unsigned i = 0; i < DIM; i++)
1760  {
1761  sum += interpolated_u[i] * interpolated_dudt[i];
1762  }
1763 
1764  d_kin_en_dt += sum * w * J;
1765  }
1766 
1767  return d_kin_en_dt;
1768 
1769  } // End of d_kin_energy_dt
1770 
1771 
1772  //==============================================================
1773  /// Return pressure integrated over the element
1774  //==============================================================
1775  template<unsigned DIM>
1777  const
1778  {
1779  // Initialise
1780  double press_int = 0;
1781 
1782  // Set the value of n_intpt
1783  unsigned n_intpt = integral_pt()->nweight();
1784 
1785  // Set the Vector to hold local coordinates
1786  Vector<double> s(DIM + 1, 0.0);
1787 
1788  // Loop over the integration points
1789  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1790  {
1791  // Assign values of s
1792  for (unsigned i = 0; i < DIM + 1; i++)
1793  {
1794  s[i] = integral_pt()->knot(ipt, i);
1795  }
1796 
1797  // Get the integral weight
1798  double w = integral_pt()->weight(ipt);
1799 
1800  // Get Jacobian of mapping
1801  double J = J_eulerian(s);
1802 
1803  // Premultiply the weights and the Jacobian
1804  double W = w * J;
1805 
1806  // Get pressure
1807  double press = interpolated_p_nst(s);
1808 
1809  // Add
1810  press_int += press * W;
1811  }
1812 
1813  return press_int;
1814  }
1815 
1816 
1817  //==============================================================
1818  /// Compute the residuals for the associated pressure advection
1819  /// diffusion problem. Used by the Fp preconditioner.
1820  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
1821  //==============================================================
1822  template<unsigned DIM>
1825  Vector<double>& residuals,
1826  DenseMatrix<double>& jacobian,
1827  const unsigned& flag)
1828  {
1829  // Throw an error
1830  throw OomphLibError("Not yet implemented\n",
1831  OOMPH_CURRENT_FUNCTION,
1832  OOMPH_EXCEPTION_LOCATION);
1833  } // End of fill_in_generic_pressure_advection_diffusion_contribution_nst
1834 
1835 
1836  //==============================================================
1837  /// Compute the residuals for the Navier-Stokes
1838  /// equations; flag=1 (or 0): do (or don't) compute the
1839  /// Jacobian as well.
1840  //==============================================================
1841  template<unsigned DIM>
1844  DenseMatrix<double>& jacobian,
1845  DenseMatrix<double>& mass_matrix,
1846  const unsigned& flag)
1847  {
1848  // Return immediately if there are no dofs
1849  if (ndof() == 0) return;
1850 
1851  // Find out how many nodes there are
1852  unsigned n_node = nnode();
1853 
1854  // Find out how many pressure dofs there are
1855  unsigned n_pres = npres_nst();
1856 
1857  // Allocate storage for the indices of the velocity components
1858  unsigned u_nodal_index[DIM];
1859 
1860  // Loop over the velocity components
1861  for (unsigned i = 0; i < DIM; i++)
1862  {
1863  // Find the index at which the i-th local velocity is stored
1864  u_nodal_index[i] = u_index_nst(i);
1865  }
1866 
1867  // Set up memory for the shape functions
1868  Shape psif(n_node);
1869 
1870  // Set up memory for the test functions
1871  Shape testf(n_node);
1872 
1873  // Allocate space for the derivatives of the shape functions
1874  DShape dpsifdx(n_node, DIM + 1);
1875 
1876  // Allocate space for the derivatives of the test functions
1877  DShape dtestfdx(n_node, DIM + 1);
1878 
1879  // Set up memory for pressure shape functions
1880  Shape psip(n_pres);
1881 
1882  // Set up memory for pressure test functions
1883  Shape testp(n_pres);
1884 
1885  // Number of integration points
1886  unsigned n_intpt = integral_pt()->nweight();
1887 
1888  // Set the Vector to hold local coordinates
1889  Vector<double> s(DIM + 1, 0.0);
1890 
1891  //-------------------------------------
1892  // Get physical variables from element:
1893  //-------------------------------------
1894  // Reynolds number must be multiplied by the density ratio
1895  double scaled_re = re() * density_ratio();
1896 
1897  // Get the scaled Strouhal value
1898  double scaled_re_st = re_st() * density_ratio();
1899 
1900  // Get the scaled Womersley number differentiated w.r.t. the Womersley
1901  // number
1902  double scaled_dre_st_dre_st = density_ratio();
1903 
1904  // Get the scaled Reynolds / Froude number
1905  double scaled_re_inv_fr = re_invfr() * density_ratio();
1906 
1907  // Get the viscosity ratio
1908  double visc_ratio = viscosity_ratio();
1909 
1910  // Get the gravity vector
1911  Vector<double> G = g();
1912 
1913  // Integer to store the local equation number
1914  int local_eqn = 0;
1915 
1916  // Integer to store the local unknown number
1917  int local_unknown = 0;
1918 
1919  // Loop over the integration points
1920  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1921  {
1922  // Assign values of s
1923  for (unsigned i = 0; i < DIM + 1; i++)
1924  {
1925  // Calculate the i-th local coordinate
1926  s[i] = integral_pt()->knot(ipt, i);
1927  }
1928 
1929  // Get the integral weight
1930  double w = integral_pt()->weight(ipt);
1931 
1932  // Call the derivatives of the shape and test functions
1933  double J = dshape_and_dtest_eulerian_at_knot_nst(
1934  ipt, psif, dpsifdx, testf, dtestfdx);
1935 
1936  // Call the pressure shape and test functions
1937  pshape_nst(s, psip, testp);
1938 
1939  // Pre-multiply the weights and the Jacobian
1940  double W = w * J;
1941 
1942  // Storage for the interpolated time value
1943  double interpolated_t = 0.0;
1944 
1945  // Storage for the interpolated pressure value
1946  double interpolated_p = 0.0;
1947 
1948  // Storage for the spatial coordinates
1949  Vector<double> interpolated_x(DIM, 0.0);
1950 
1951  // Storage for the interpolated velocity components
1952  Vector<double> interpolated_u(DIM, 0.0);
1953 
1954  // Storage for the interpolated time-derivative of the velocities
1955  Vector<double> interpolated_dudt(DIM, 0.0);
1956 
1957  // Storage for the mesh velocity
1958  Vector<double> mesh_velocity(DIM, 0.0);
1959 
1960  // Storage for the spatial derivatives of the velocity components
1961  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1962 
1963  // Calculate pressure
1964  for (unsigned l = 0; l < n_pres; l++)
1965  {
1966  // Update the current approximation to the interpolated pressure
1967  interpolated_p += p_nst(l) * psip[l];
1968  }
1969 
1970  //-------------------------------------------------------------------
1971  // Calculate velocities, derivatives, source function and body force:
1972  //-------------------------------------------------------------------
1973  // Loop over nodes
1974  for (unsigned l = 0; l < n_node; l++)
1975  {
1976  // Update the interpolated time value
1977  interpolated_t += raw_nodal_position(l, DIM) * psif(l);
1978 
1979  // Loop over coordinate directions
1980  for (unsigned i = 0; i < DIM; i++)
1981  {
1982  // Get the nodal value
1983  double u_value = raw_nodal_value(l, u_nodal_index[i]);
1984 
1985  // Update the i-th interpolated velocity component
1986  interpolated_u[i] += u_value * psif[l];
1987 
1988  // Update the i-th interpolated coordinate value
1989  interpolated_x[i] += raw_nodal_position(l, i) * psif[l];
1990 
1991  // Update the interpolated du_i/dt value
1992  interpolated_dudt[i] += u_value * dpsifdx(l, DIM);
1993 
1994  // Loop over derivative directions
1995  for (unsigned j = 0; j < DIM; j++)
1996  {
1997  // Update the interpolated du_i/dx_j value
1998  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1999  }
2000  } // for (unsigned i=0;i<DIM;i++)
2001  } // for (unsigned l=0;l<n_node;l++)
2002 
2003  // If ALE is enabled
2004  if (!ALE_is_disabled)
2005  {
2006  // Loop over nodes
2007  for (unsigned l = 0; l < n_node; l++)
2008  {
2009  // Loop over directions
2010  for (unsigned i = 0; i < DIM; i++)
2011  {
2012  // Update the i-th mesh velocity component
2013  mesh_velocity[i] += raw_nodal_position(l, i) * dpsifdx(l, DIM);
2014  }
2015  } // for (unsigned l=0;l<n_node;l++)
2016  } // if (!ALE_is_disabled)
2017 
2018  // Allocate space for the body force
2019  Vector<double> body_force(DIM, 0.0);
2020 
2021  // Get the user-defined body force term
2022  get_body_force_nst(interpolated_t, ipt, s, interpolated_x, body_force);
2023 
2024  // Get the user-defined source function
2025  double source = get_source_nst(interpolated_t, ipt, interpolated_x);
2026 
2027  //---------------------------------
2028  // Assemble residuals and Jacobian:
2029  //---------------------------------
2030  //--------------------
2031  // Momentum equations:
2032  //--------------------
2033  // Loop over the test functions
2034  for (unsigned l = 0; l < n_node; l++)
2035  {
2036  // Loop over the velocity components
2037  for (unsigned i = 0; i < DIM; i++)
2038  {
2039  // Get the local equation number associated with this unknown and node
2040  local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
2041 
2042  // If it's not a boundary condition
2043  if (local_eqn >= 0)
2044  {
2045  // Add the user-defined body force terms
2046  residuals[local_eqn] += body_force[i] * testf[l] * W;
2047 
2048  // Add the gravitational body force term
2049  residuals[local_eqn] += scaled_re_inv_fr * testf[l] * G[i] * W;
2050 
2051  // Add the pressure gradient term
2052  residuals[local_eqn] += interpolated_p * dtestfdx(l, i) * W;
2053 
2054  // Add in the contribution from the time derivative
2055  residuals[local_eqn] -=
2056  scaled_re_st * interpolated_dudt[i] * testf[l] * W;
2057 
2058  // If ALE is enabled
2059  if (!ALE_is_disabled)
2060  {
2061  // Loop over the coordinate directions
2062  for (unsigned k = 0; k < DIM; k++)
2063  {
2064  // Add in the mesh velocity contribution
2065  residuals[local_eqn] +=
2066  (scaled_re_st * mesh_velocity[k] * interpolated_dudx(i, k) *
2067  testf[l] * W);
2068  }
2069  } // if (!ALE_is_disabled)
2070 
2071  // Loop over the coordinate directions
2072  for (unsigned k = 0; k < DIM; k++)
2073  {
2074  // Add in the convective term contribution
2075  residuals[local_eqn] -= (scaled_re * interpolated_u[k] *
2076  interpolated_dudx(i, k) * testf[l] * W);
2077  }
2078 
2079  // Loop over the coordinate directions
2080  for (unsigned k = 0; k < DIM; k++)
2081  {
2082  // Add in the stress tensor terms:
2083  // NOTE: The viscosity ratio needs to go in here to ensure
2084  // continuity of normal stress is satisfied even in flows
2085  // with zero pressure gradient!
2086  residuals[local_eqn] -= ((interpolated_dudx(i, k) +
2087  Gamma[i] * interpolated_dudx(k, i)) *
2088  visc_ratio * dtestfdx(l, k) * W);
2089  }
2090 
2091  //------------------------
2092  // Calculate the Jacobian:
2093  //------------------------
2094  // If we also need to construct the Jacobian
2095  if (flag)
2096  {
2097  // Loop over the velocity shape functions again
2098  for (unsigned l2 = 0; l2 < n_node; l2++)
2099  {
2100  // Loop over the velocity components again
2101  for (unsigned i2 = 0; i2 < DIM; i2++)
2102  {
2103  // Get the local equation number associated with this unknown
2104  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
2105 
2106  // If we're at a proper degree of freedom
2107  if (local_unknown >= 0)
2108  {
2109  // Add the contribution to the elemental matrix
2110  jacobian(local_eqn, local_unknown) -=
2111  (visc_ratio * Gamma[i] * dpsifdx(l2, i) *
2112  dtestfdx(l, i2) * W);
2113 
2114  // Now add in the inertial terms
2115  jacobian(local_eqn, local_unknown) -=
2116  (scaled_re * psif[l2] * interpolated_dudx(i, i2) *
2117  testf[l] * W);
2118 
2119  // Extra component if i2=i
2120  if (i2 == i)
2121  {
2122  // If we also need to construct the mass matrix (only
2123  // diagonal entries)
2124  if (flag == 2)
2125  {
2126  // NOTE: This is positive because the mass matrix is
2127  // taken to the other side of the equation when
2128  // formulating the generalised eigenproblem.
2129  mass_matrix(local_eqn, local_unknown) +=
2130  (scaled_re_st * psif[l2] * testf[l] * W);
2131  }
2132 
2133  // Add in the time-derivative contribution
2134  jacobian(local_eqn, local_unknown) -=
2135  (scaled_re_st * dpsifdx(l2, DIM) * testf[l] * W);
2136 
2137  // If ALE is enabled
2138  if (!ALE_is_disabled)
2139  {
2140  // Loop over the velocity components
2141  for (unsigned k = 0; k < DIM; k++)
2142  {
2143  // Add in the mesh velocity contribution
2144  jacobian(local_eqn, local_unknown) +=
2145  (scaled_re_st * mesh_velocity[k] * dpsifdx(l2, k) *
2146  testf[l] * W);
2147  }
2148  } // if (!ALE_is_disabled)
2149 
2150  // Loop over the velocity components
2151  for (unsigned k = 0; k < DIM; k++)
2152  {
2153  // Add in the convective term contribution
2154  jacobian(local_eqn, local_unknown) -=
2155  (scaled_re * interpolated_u[k] * dpsifdx(l2, k) *
2156  testf[l] * W);
2157  }
2158 
2159  // Loop over the velocity components
2160  for (unsigned k = 0; k < DIM; k++)
2161  {
2162  // Add in the velocity gradient terms
2163  jacobian(local_eqn, local_unknown) -=
2164  (visc_ratio * dpsifdx(l2, k) * dtestfdx(l, k) * W);
2165  }
2166  } // if (i2==i)
2167  } // if (local_unknown>=0)
2168  } // for (unsigned i2=0;i2<DIM;i2++)
2169  } // for (unsigned l2=0;l2<n_node;l2++)
2170 
2171  // Loop over pressure shape functions
2172  for (unsigned l2 = 0; l2 < n_pres; l2++)
2173  {
2174  // Get the local equation number associated with this degree of
2175  // freedom
2176  local_unknown = p_local_eqn(l2);
2177 
2178  // If we are at a proper degree of freedom
2179  if (local_unknown >= 0)
2180  {
2181  // Add in the pressure gradient contribution
2182  jacobian(local_eqn, local_unknown) +=
2183  psip[l2] * dtestfdx(l, i) * W;
2184  }
2185  } // for (unsigned l2=0;l2<n_pres;l2++)
2186 
2187  //------------------------------------
2188  // Calculate external data information
2189  //------------------------------------
2190  // If we're storing the Strouhal number as external data then we
2191  // need to calculate dr/d(St) in the elemental Jacobian
2192  if (ReynoldsStrouhal_is_stored_as_external_data)
2193  {
2194  // The index of the external data (there's only one!)
2195  unsigned data_index = 0;
2196 
2197  // The index of the unknown value stored in the external data
2198  unsigned value_index = 0;
2199 
2200  // Get the local equation number associated with the extra
2201  // unknown
2202  local_unknown =
2203  this->external_local_eqn(data_index, value_index);
2204 
2205  // If we're at a non-zero degree of freedom add in the entry
2206  if (local_unknown >= 0)
2207  {
2208  // Add in the contribution from the time derivative
2209  jacobian(local_eqn, local_unknown) -=
2210  (scaled_dre_st_dre_st * interpolated_dudt[i] * testf[l] *
2211  W);
2212 
2213  // If ALE is enabled
2214  if (!ALE_is_disabled)
2215  {
2216  // Loop over the coordinate directions
2217  for (unsigned k = 0; k < DIM; k++)
2218  {
2219  // Add in the mesh velocity contribution
2220  jacobian(local_eqn, local_unknown) +=
2221  (scaled_dre_st_dre_st * mesh_velocity[k] *
2222  interpolated_dudx(i, k) * testf[l] * W);
2223  }
2224  } // if (!ALE_is_disabled)
2225  } // if (local_unknown>=0)
2226  } // if (Strouhal_is_stored_as_external_data)
2227  } // if (flag)
2228  } // if (local_eqn>=0)
2229  } // for (unsigned i=0;i<DIM;i++)
2230  } // for (unsigned l=0;l<n_node;l++)
2231 
2232  //---------------------
2233  // Continuity equation:
2234  //---------------------
2235  // Loop over the shape functions
2236  for (unsigned l = 0; l < n_pres; l++)
2237  {
2238  // Get the local equation number associated with the pressure dof
2239  local_eqn = p_local_eqn(l);
2240 
2241  // If it's not a boundary condition
2242  if (local_eqn >= 0)
2243  {
2244  // Add in the source term contribution
2245  residuals[local_eqn] -= source * testp[l] * W;
2246 
2247  // Loop over velocity components
2248  for (unsigned k = 0; k < DIM; k++)
2249  {
2250  // Add in the velocity gradient terms
2251  residuals[local_eqn] += interpolated_dudx(k, k) * testp[l] * W;
2252  }
2253 
2254  //------------------------
2255  // Calculate the Jacobian:
2256  //------------------------
2257  // If we also need to construct the Jacobian
2258  if (flag)
2259  {
2260  // Loop over the velocity shape functions
2261  for (unsigned l2 = 0; l2 < n_node; l2++)
2262  {
2263  // Loop over velocity components
2264  for (unsigned i2 = 0; i2 < DIM; i2++)
2265  {
2266  // Get the local equation number associated with this node
2267  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
2268 
2269  // If we're at a proper degree of freedom
2270  if (local_unknown >= 0)
2271  {
2272  // Add in the velocity gradient contribution
2273  jacobian(local_eqn, local_unknown) +=
2274  dpsifdx(l2, i2) * testp[l] * W;
2275  }
2276  } // for (unsigned i2=0;i2<DIM;i2++)
2277  } // for (unsigned l2=0;l2<n_node;l2++)
2278  } // if (flag)
2279  } // if (local_eqn>=0)
2280  } // for (unsigned l=0;l<n_pres;l++)
2281  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
2282  } // End of fill_in_generic_residual_contribution_nst
2283 
2284  //================================================================
2285  /// Compute the derivatives of the residuals for the Navier-Stokes
2286  /// equations with respect to a parameter;
2287  /// flag=2 or 1 or 0: do (or don't) compute the Jacobian and
2288  /// mass matrix as well
2289  //================================================================
2290  template<unsigned DIM>
2293  double* const& parameter_pt,
2294  Vector<double>& dres_dparam,
2295  DenseMatrix<double>& djac_dparam,
2296  DenseMatrix<double>& dmass_matrix_dparam,
2297  const unsigned& flag)
2298  {
2299  // Throw an error
2300  throw OomphLibError("Not yet implemented\n",
2301  OOMPH_CURRENT_FUNCTION,
2302  OOMPH_EXCEPTION_LOCATION);
2303  } // End of fill_in_generic_dresidual_contribution_nst
2304 
2305  //==================================================================
2306  /// Compute the hessian tensor vector products required to
2307  /// perform continuation of bifurcations analytically
2308  //==================================================================
2309  template<unsigned DIM>
2312  Vector<double> const& Y,
2313  DenseMatrix<double> const& C,
2314  DenseMatrix<double>& product)
2315  {
2316  // Throw an error
2317  throw OomphLibError("Not yet implemented\n",
2318  OOMPH_CURRENT_FUNCTION,
2319  OOMPH_EXCEPTION_LOCATION);
2320  } // End of fill_in_contribution_to_hessian_vector_products
2321 
2322  //======================================================================
2323  /// Compute derivatives of elemental residual vector with respect
2324  /// to nodal coordinates.
2325  /// dresidual_dnodal_coordinates(l,i,j)=d res(l) / dX_{ij}
2326  /// Overloads the FD-based version in the FE base class.
2327  /// DRAIG: This needs doing carefully if the ALE nodes aren't fixed!!!
2328  //======================================================================
2329  template<unsigned DIM>
2331  DIM>::get_dresidual_dnodal_coordinates(RankThreeTensor<double>&
2332  dresidual_dnodal_coordinates)
2333  {
2334  // Throw a warning
2335  throw OomphLibError("Space-time update needs to be checked!",
2336  OOMPH_CURRENT_FUNCTION,
2337  OOMPH_EXCEPTION_LOCATION);
2338  } // End of get_dresidual_dnodal_coordinates
2339 
2340 
2341  /// ////////////////////////////////////////////////////////////////////////
2342  /// ////////////////////////////////////////////////////////////////////////
2343  /// ////////////////////////////////////////////////////////////////////////
2344 
2345  //---------------------------------
2346  // Space-time Taylor-Hood elements:
2347  //---------------------------------
2348  // Set the data for the number of Variables at each node
2349  template<>
2351  {3, 2, 3, 2, 2, 2, 3, 2, 3, 2, 2, 2, 2, 2,
2352  2, 2, 2, 2, 3, 2, 3, 2, 2, 2, 3, 2, 3};
2353 
2354  // Set the data for the pressure conversion array
2355  template<>
2357  0, 2, 6, 8, 18, 20, 24, 26};
2358 
2359  //=========================================================================
2360  /// Add to the set \c paired_load_data pairs containing
2361  /// - the pointer to a Data object
2362  /// and
2363  /// - the index of the value in that Data object
2364  /// .
2365  /// for all values (pressures, velocities) that affect the
2366  /// load computed in the \c get_load(...) function.
2367  //=========================================================================
2368  template<unsigned DIM>
2370  std::set<std::pair<Data*, unsigned>>& paired_load_data)
2371  {
2372  // Allocate storage for the indices of the velocity components
2373  unsigned u_index[DIM];
2374 
2375  // Loop over the velocity components
2376  for (unsigned i = 0; i < DIM; i++)
2377  {
2378  // Get the index at which the i-th velocity component is stored
2379  u_index[i] = this->u_index_nst(i);
2380  }
2381 
2382  // Get the number of nodes in this element
2383  unsigned n_node = this->nnode();
2384 
2385  // Loop over the nodes
2386  for (unsigned n = 0; n < n_node; n++)
2387  {
2388  // Loop over the velocity components and add pointer to their data
2389  // and indices to the vectors
2390  for (unsigned i = 0; i < DIM; i++)
2391  {
2392  // Add in the node and equation number pair
2393  paired_load_data.insert(std::make_pair(this->node_pt(n), u_index[i]));
2394  }
2395  } // for (unsigned n=0;n<n_node;n++)
2396 
2397  // Identify the pressure data
2398  this->identify_pressure_data(paired_load_data);
2399  } // End of identify_load_data
2400 
2401  //=========================================================================
2402  /// Add to the set \c paired_pressure_data pairs containing
2403  /// - the pointer to a Data object
2404  /// and
2405  /// - the index of the value in that Data object
2406  /// .
2407  /// for pressure values that affect the
2408  /// load computed in the \c get_load(...) function.,
2409  //=========================================================================
2410  template<unsigned DIM>
2412  std::set<std::pair<Data*, unsigned>>& paired_pressure_data)
2413  {
2414  // Find the index at which the pressure is stored
2415  unsigned p_index = static_cast<unsigned>(this->p_nodal_index_nst());
2416 
2417  // Get the number of pressure degrees of freedom
2418  unsigned n_pres = npres_nst();
2419 
2420  // Loop over the pressure data
2421  for (unsigned l = 0; l < n_pres; l++)
2422  {
2423  // The DIM-th entry in each nodal data is the pressure, which affects
2424  // the traction
2425  paired_pressure_data.insert(
2426  std::make_pair(this->node_pt(Pconv[l]), p_index));
2427  }
2428  } // End of identify_pressure_data
2429 
2430  //====================================================================
2431  /// Force build of templates
2432  //====================================================================
2435 } // namespace oomph
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
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
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void identify_pressure_data(std::set< std::pair< Data *, unsigned >> &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
void identify_load_data(std::set< std::pair< Data *, unsigned >> &paired_load_data)
Add to the set paired_load_data pairs containing.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h: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 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_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 ...
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....
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
double d_kin_energy_dt() const
Get integral of time derivative of kinetic energy over element.
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 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....
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...
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,...
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 compute_norm(Vector< double > &norm)
Compute the vector norm of the FEM solution.
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 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 output_vorticity(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,t,omega in tecplot format. nplot points in each coordinate direction.
void get_vorticity(const Vector< double > &s, Vector< double > &vorticity) const
Compute the vorticity vector at local coordinate s.
void output(std::ostream &outfile)
Output function: x,y,t,u,v,p in tecplot format. The default number of plot points is five.
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 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....
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
double Default_Physical_Constant_Value
Default value for physical constants.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...