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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Non-inline functions for NS elements
31 
32 #include "navier_stokes_elements.h"
33 
34 
35 namespace oomph
36 {
37 
38 
39 
40 
41 
42 
43 
44 /////////////////////////////////////////////////////////////////////////
45 /////////////////////////////////////////////////////////////////////////
46 /////////////////////////////////////////////////////////////////////////
47 
48 
49 /// Navier--Stokes equations static data
50 template<unsigned DIM>
51 Vector<double> NavierStokesEquations<DIM>::Gamma(DIM,1.0);
52 
53 //=================================================================
54 /// "Magic" negative number that indicates that the pressure is
55 /// not stored at a node. This cannot be -1 because that represents
56 /// the positional hanging scheme in the hanging_pt object of nodes
57 //=================================================================
58 template<unsigned DIM>
59 int NavierStokesEquations<DIM>::Pressure_not_stored_at_node = -100;
60 
61 /// Navier--Stokes equations static data
62 template<unsigned DIM>
64 
65 /// Navier--Stokes equations static data
66 template<unsigned DIM>
67 double NavierStokesEquations<DIM>::Default_Physical_Ratio_Value = 1.0;
68 
69 /// Navier-Stokes equations default gravity vector
70 template<unsigned DIM>
72 
73 
74 
75 
76 //===================================================================
77  /// Compute the diagonal of the velocity/pressure mass matrices.
78  /// If which one=0, both are computed, otherwise only the pressure
79  /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
80  /// LSC version of the preconditioner only needs that one)
81 //===================================================================
82  template<unsigned DIM>
85  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
86  const unsigned& which_one)
87  {
88  // Resize and initialise
89  unsigned n_dof=ndof();
90 
91  if ((which_one==0)||(which_one==1)) press_mass_diag.assign(n_dof,0.0);
92  if ((which_one==0)||(which_one==2)) veloc_mass_diag.assign(n_dof,0.0);
93 
94  // find out how many nodes there are
95  unsigned n_node = nnode();
96 
97  // find number of spatial dimensions
98  unsigned n_dim = this->dim();
99 
100  // Local coordinates
101  Vector<double> s(n_dim);
102 
103  // find the indices at which the local velocities are stored
104  Vector<unsigned> u_nodal_index(n_dim);
105  for(unsigned i=0; i<n_dim; i++)
106  {
107  u_nodal_index[i] = this->u_index_nst(i);
108  }
109 
110  //Set up memory for veloc 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 
129  //Get the integral weight
130  double w = integral_pt()->weight(ipt);
131 
132  //Get determinant of Jacobian of the mapping
133  double J = J_eulerian_at_knot(ipt);
134 
135  //Assign values of s
136  for(unsigned i=0;i<n_dim;i++)
137  {
138  s[i] = integral_pt()->knot(ipt,i);
139  }
140 
141  //Premultiply weights and Jacobian
142  double W = w*J;
143 
144 
145 
146  // Do we want the velocity one?
147  if ((which_one==0)||(which_one==2))
148  {
149 
150  //Get the velocity shape functions
151  shape_at_knot(ipt,psi);
152 
153  //Loop over the veclocity shape functions
154  for(unsigned l=0; l<n_node; l++)
155  {
156  //Loop over the velocity components
157  for(unsigned i=0; i<n_dim; i++)
158  {
159  local_eqn = nodal_local_eqn(l,u_nodal_index[i]);
160 
161  //If not a boundary condition
162  if(local_eqn >= 0)
163  {
164  //Add the contribution
165  veloc_mass_diag[local_eqn] += pow(psi[l],2) * W;
166  }
167  }
168  }
169  }
170 
171  // Do we want the pressure one?
172  if ((which_one==0)||(which_one==1))
173  {
174  //Get the pressure shape functions
175  pshape_nst(s,psi_p);
176 
177  //Loop over the veclocity shape functions
178  for(unsigned l=0; l<n_pres; l++)
179  {
180  // Get equation number
181  local_eqn = p_local_eqn(l);
182 
183  //If not a boundary condition
184  if(local_eqn >= 0)
185  {
186  //Add the contribution
187  press_mass_diag[local_eqn] += pow(psi_p[l],2) * W;
188  }
189  }
190  }
191 
192  }
193  }
194 
195 
196 //======================================================================
197 /// Validate against exact velocity solution at given time.
198 /// Solution is provided via function pointer.
199 /// Plot at a given number of plot points and compute L2 error
200 /// and L2 norm of velocity solution over element.
201 //=======================================================================
202 template<unsigned DIM>
204 compute_error(std::ostream &outfile,
206  const double& time,
207  double& error, double& norm)
208 {
209  error=0.0;
210  norm=0.0;
211 
212  //Vector of local coordinates
213  Vector<double> s(DIM);
214 
215  // Vector for coordintes
216  Vector<double> x(DIM);
217 
218  //Set the value of n_intpt
219  unsigned n_intpt = integral_pt()->nweight();
220 
221  outfile << "ZONE" << std::endl;
222 
223  // Exact solution Vector (u,v,[w],p)
224  Vector<double> exact_soln(DIM+1);
225 
226  //Loop over the integration points
227  for(unsigned ipt=0;ipt<n_intpt;ipt++)
228  {
229 
230  //Assign values of s
231  for(unsigned i=0;i<DIM;i++)
232  {
233  s[i] = integral_pt()->knot(ipt,i);
234  }
235 
236  //Get the integral weight
237  double w = integral_pt()->weight(ipt);
238 
239  // Get jacobian of mapping
240  double J=J_eulerian(s);
241 
242  //Premultiply the weights and the Jacobian
243  double W = w*J;
244 
245  // Get x position as Vector
246  interpolated_x(s,x);
247 
248  // Get exact solution at this point
249  (*exact_soln_pt)(time,x,exact_soln);
250 
251  // Velocity error
252  for(unsigned i=0;i<DIM;i++)
253  {
254  norm+=exact_soln[i]*exact_soln[i]*W;
255  error+=(exact_soln[i]-interpolated_u_nst(s,i))*
256  (exact_soln[i]-interpolated_u_nst(s,i))*W;
257  }
258 
259  //Output x,y,...,u_exact
260  for(unsigned i=0;i<DIM;i++)
261  {
262  outfile << x[i] << " ";
263  }
264 
265  //Output x,y,[z],u_error,v_error,[w_error]
266  for(unsigned i=0;i<DIM;i++)
267  {
268  outfile << exact_soln[i]-interpolated_u_nst(s,i) << " ";
269  }
270 
271  outfile << std::endl;
272 
273  }
274 }
275 
276 //======================================================================
277 /// Validate against exact velocity solution
278 /// Solution is provided via function pointer.
279 /// Plot at a given number of plot points and compute L2 error
280 /// and L2 norm of velocity solution over element.
281 //=======================================================================
282 template<unsigned DIM>
284  std::ostream &outfile,
286  double& error, double& norm)
287 {
288 
289  error=0.0;
290  norm=0.0;
291 
292  //Vector of local coordinates
293  Vector<double> s(DIM);
294 
295  // Vector for coordintes
296  Vector<double> x(DIM);
297 
298  //Set the value of n_intpt
299  unsigned n_intpt = integral_pt()->nweight();
300 
301 
302  outfile << "ZONE" << std::endl;
303 
304  // Exact solution Vector (u,v,[w],p)
305  Vector<double> exact_soln(DIM+1);
306 
307  //Loop over the integration points
308  for(unsigned ipt=0;ipt<n_intpt;ipt++)
309  {
310 
311  //Assign values of s
312  for(unsigned i=0;i<DIM;i++)
313  {
314  s[i] = integral_pt()->knot(ipt,i);
315  }
316 
317  //Get the integral weight
318  double w = integral_pt()->weight(ipt);
319 
320  // Get jacobian of mapping
321  double J=J_eulerian(s);
322 
323  //Premultiply the weights and the Jacobian
324  double W = w*J;
325 
326  // Get x position as Vector
327  interpolated_x(s,x);
328 
329  // Get exact solution at this point
330  (*exact_soln_pt)(x,exact_soln);
331 
332  // Velocity error
333  for(unsigned i=0;i<DIM;i++)
334  {
335  norm+=exact_soln[i]*exact_soln[i]*W;
336  error+=(exact_soln[i]-interpolated_u_nst(s,i))*
337  (exact_soln[i]-interpolated_u_nst(s,i))*W;
338  }
339 
340  //Output x,y,...,u_exact
341  for(unsigned i=0;i<DIM;i++)
342  {
343  outfile << x[i] << " ";
344  }
345 
346  //Output x,y,[z],u_error,v_error,[w_error]
347  for(unsigned i=0;i<DIM;i++)
348  {
349  outfile << exact_soln[i]-interpolated_u_nst(s,i) << " ";
350  }
351  outfile << std::endl;
352  }
353 }
354 
355 //======================================================================
356 /// Output "exact" solution
357 /// Solution is provided via function pointer.
358 /// Plot at a given number of plot points.
359 /// Function prints as many components as are returned in solution Vector.
360 //=======================================================================
361 template<unsigned DIM>
362 void NavierStokesEquations<DIM>::output_fct(std::ostream &outfile,
363  const unsigned &nplot,
365 {
366 
367  //Vector of local coordinates
368  Vector<double> s(DIM);
369 
370  // Vector for coordintes
371  Vector<double> x(DIM);
372 
373  // Tecplot header info
374  outfile << tecplot_zone_string(nplot);
375 
376  // Exact solution Vector
377  Vector<double> exact_soln;
378 
379  // Loop over plot points
380  unsigned num_plot_points=nplot_points(nplot);
381  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
382  {
383 
384  // Get local coordinates of plot point
385  get_s_plot(iplot,nplot,s);
386 
387  // Get x position as Vector
388  interpolated_x(s,x);
389 
390  // Get exact solution at this point
391  (*exact_soln_pt)(x,exact_soln);
392 
393  //Output x,y,...
394  for(unsigned i=0;i<DIM;i++)
395  {
396  outfile << x[i] << " ";
397  }
398 
399  //Output "exact solution"
400  for(unsigned i=0;i<exact_soln.size();i++)
401  {
402  outfile << exact_soln[i] << " ";
403  }
404 
405  outfile << std::endl;
406 
407  }
408 
409  // Write tecplot footer (e.g. FE connectivity lists)
410  write_tecplot_zone_footer(outfile,nplot);
411 
412 }
413 
414 //======================================================================
415 /// Output "exact" solution at a given time
416 /// Solution is provided via function pointer.
417 /// Plot at a given number of plot points.
418 /// Function prints as many components as are returned in solution Vector.
419 //=======================================================================
420 template<unsigned DIM>
421 void NavierStokesEquations<DIM>::output_fct(std::ostream &outfile,
422  const unsigned &nplot,
423  const double& time,
425 {
426 
427  //Vector of local coordinates
428  Vector<double> s(DIM);
429 
430  // Vector for coordintes
431  Vector<double> x(DIM);
432 
433  // Tecplot header info
434  outfile << tecplot_zone_string(nplot);
435 
436  // Exact solution Vector
437  Vector<double> exact_soln;
438 
439  // Loop over plot points
440  unsigned num_plot_points=nplot_points(nplot);
441  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
442  {
443 
444  // Get local coordinates of plot point
445  get_s_plot(iplot,nplot,s);
446 
447  // Get x position as Vector
448  interpolated_x(s,x);
449 
450  // Get exact solution at this point
451  (*exact_soln_pt)(time,x,exact_soln);
452 
453  //Output x,y,...
454  for(unsigned i=0;i<DIM;i++)
455  {
456  outfile << x[i] << " ";
457  }
458 
459  //Output "exact solution"
460  for(unsigned i=0;i<exact_soln.size();i++)
461  {
462  outfile << exact_soln[i] << " ";
463  }
464 
465  outfile << std::endl;
466  }
467 
468  // Write tecplot footer (e.g. FE connectivity lists)
469  write_tecplot_zone_footer(outfile,nplot);
470 
471 }
472 
473 //==============================================================
474 /// Output function: Velocities only
475 /// x,y,[z],u,v,[w]
476 /// in tecplot format at specified previous timestep (t=0: present;
477 /// t>0: previous timestep). Specified number of plot points in each
478 /// coordinate direction.
479 //==============================================================
480 template<unsigned DIM>
481 void NavierStokesEquations<DIM>::output_veloc(std::ostream& outfile,
482  const unsigned& nplot,
483  const unsigned& t)
484 {
485 
486  //Find number of nodes
487  unsigned n_node = nnode();
488 
489  //Local shape function
490  Shape psi(n_node);
491 
492  //Vectors of local coordinates and coords and velocities
493  Vector<double> s(DIM);
494  Vector<double> interpolated_x(DIM);
495  Vector<double> interpolated_u(DIM);
496 
497  // Tecplot header info
498  outfile << tecplot_zone_string(nplot);
499 
500  // Loop over plot points
501  unsigned num_plot_points=nplot_points(nplot);
502  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
503  {
504 
505  // Get local coordinates of plot point
506  get_s_plot(iplot,nplot,s);
507 
508  // Get shape functions
509  shape(s,psi);
510 
511  // Loop over directions
512  for(unsigned i=0;i<DIM;i++)
513  {
514  interpolated_x[i]=0.0;
515  interpolated_u[i]=0.0;
516  //Get the index at which velocity i is stored
517  unsigned u_nodal_index = u_index_nst(i);
518  //Loop over the local nodes and sum
519  for(unsigned l=0;l<n_node;l++)
520  {
521  interpolated_u[i] += nodal_value(t,l,u_nodal_index)*psi[l];
522  interpolated_x[i] += nodal_position(t,l,i)*psi[l];
523  }
524  }
525 
526  // Coordinates
527  for(unsigned i=0;i<DIM;i++)
528  {
529  outfile << interpolated_x[i] << " ";
530  }
531 
532  // Velocities
533  for(unsigned i=0;i<DIM;i++)
534  {
535  outfile << interpolated_u[i] << " ";
536  }
537 
538  outfile << std::endl;
539  }
540 
541  // Write tecplot footer (e.g. FE connectivity lists)
542  write_tecplot_zone_footer(outfile,nplot);
543 
544 }
545 
546 //==============================================================
547 /// Output function:
548 /// x,y,[z],u,v,[w],p
549 /// in tecplot format. Specified number of plot points in each
550 /// coordinate direction.
551 //==============================================================
552 template<unsigned DIM>
553 void NavierStokesEquations<DIM>::output(std::ostream &outfile,
554  const unsigned &nplot)
555 {
556 
557  //Vector of local coordinates
558  Vector<double> s(DIM);
559 
560  // Tecplot header info
561  outfile << tecplot_zone_string(nplot);
562 
563  // Loop over plot points
564  unsigned num_plot_points=nplot_points(nplot);
565  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
566  {
567 
568  // Get local coordinates of plot point
569  get_s_plot(iplot,nplot,s);
570 
571  // Coordinates
572  for(unsigned i=0;i<DIM;i++)
573  {
574  outfile << interpolated_x(s,i) << " ";
575  }
576 
577  // Velocities
578  for(unsigned i=0;i<DIM;i++)
579  {
580  outfile << interpolated_u_nst(s,i) << " ";
581  }
582 
583  // Pressure
584  outfile << interpolated_p_nst(s) << " ";
585 
586  outfile << std::endl;
587  }
588  outfile << std::endl;
589 
590  // Write tecplot footer (e.g. FE connectivity lists)
591  write_tecplot_zone_footer(outfile,nplot);
592 
593 }
594 
595 
596 //==============================================================
597 /// C-style output function:
598 /// x,y,[z],u,v,[w],p
599 /// in tecplot format. Specified number of plot points in each
600 /// coordinate direction.
601 //==============================================================
602 template<unsigned DIM>
604  const unsigned &nplot)
605 {
606 
607  //Vector of local coordinates
608  Vector<double> s(DIM);
609 
610  // Tecplot header info
611  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
612 
613  // Loop over plot points
614  unsigned num_plot_points=nplot_points(nplot);
615  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
616  {
617 
618  // Get local coordinates of plot point
619  get_s_plot(iplot,nplot,s);
620 
621  // Coordinates
622  for(unsigned i=0;i<DIM;i++)
623  {
624  fprintf(file_pt,"%g ",interpolated_x(s,i));
625  }
626 
627  // Velocities
628  for(unsigned i=0;i<DIM;i++)
629  {
630  fprintf(file_pt,"%g ",interpolated_u_nst(s,i));
631  }
632 
633  // Pressure
634  fprintf(file_pt,"%g \n",interpolated_p_nst(s));
635  }
636  fprintf(file_pt,"\n");
637 
638  // Write tecplot footer (e.g. FE connectivity lists)
639  write_tecplot_zone_footer(file_pt,nplot);
640 
641 }
642 
643 
644 //==============================================================
645 /// Full output function:
646 /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
647 /// in tecplot format. Specified number of plot points in each
648 /// coordinate direction
649 //==============================================================
650 template<unsigned DIM>
651 void NavierStokesEquations<DIM>::full_output(std::ostream &outfile,
652  const unsigned &nplot)
653 {
654 
655  //Vector of local coordinates
656  Vector<double> s(DIM);
657 
658  // Acceleration
659  Vector<double> dudt(DIM);
660 
661  // Mesh elocity
662  Vector<double> mesh_veloc(DIM,0.0);
663 
664  // Tecplot header info
665  outfile << tecplot_zone_string(nplot);
666 
667  //Find out how many nodes there are
668  unsigned n_node = nnode();
669 
670  //Set up memory for the shape functions
671  Shape psif(n_node);
672  DShape dpsifdx(n_node,DIM);
673 
674  // Loop over plot points
675  unsigned num_plot_points=nplot_points(nplot);
676  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
677  {
678 
679  // Get local coordinates of plot point
680  get_s_plot(iplot,nplot,s);
681 
682  //Call the derivatives of the shape and test functions
683  dshape_eulerian(s,psif,dpsifdx);
684 
685  //Allocate storage
686  Vector<double> mesh_veloc(DIM);
687  Vector<double> dudt(DIM);
688  Vector<double> dudt_ALE(DIM);
689  DenseMatrix<double> interpolated_dudx(DIM,DIM);
690 
691  //Initialise everything to zero
692  for(unsigned i=0;i<DIM;i++)
693  {
694  mesh_veloc[i]=0.0;
695  dudt[i]=0.0;
696  dudt_ALE[i]=0.0;
697  for(unsigned j=0;j<DIM;j++)
698  {
699  interpolated_dudx(i,j) = 0.0;
700  }
701  }
702 
703  //Calculate velocities and derivatives
704 
705 
706  //Loop over directions
707  for(unsigned i=0;i<DIM;i++)
708  {
709  //Get the index at which velocity i is stored
710  unsigned u_nodal_index = u_index_nst(i);
711  // Loop over nodes
712  for(unsigned l=0;l<n_node;l++)
713  {
714  dudt[i]+=du_dt_nst(l,u_nodal_index)*psif[l];
715  mesh_veloc[i]+=dnodal_position_dt(l,i)*psif[l];
716 
717  //Loop over derivative directions for velocity gradients
718  for(unsigned j=0;j<DIM;j++)
719  {
720  interpolated_dudx(i,j) += nodal_value(l,u_nodal_index)*dpsifdx(l,j);
721  }
722  }
723  }
724 
725 
726  // Get dudt in ALE form (incl mesh veloc)
727  for(unsigned i=0;i<DIM;i++)
728  {
729  dudt_ALE[i]=dudt[i];
730  for (unsigned k=0;k<DIM;k++)
731  {
732  dudt_ALE[i]-=mesh_veloc[k]*interpolated_dudx(i,k);
733  }
734  }
735 
736 
737  // Coordinates
738  for(unsigned i=0;i<DIM;i++)
739  {
740  outfile << interpolated_x(s,i) << " ";
741  }
742 
743  // Velocities
744  for(unsigned i=0;i<DIM;i++)
745  {
746  outfile << interpolated_u_nst(s,i) << " ";
747  }
748 
749  // Pressure
750  outfile << interpolated_p_nst(s) << " ";
751 
752  // Accelerations
753  for(unsigned i=0;i<DIM;i++)
754  {
755  outfile << dudt_ALE[i] << " ";
756  }
757 
758 // // Mesh velocity
759 // for(unsigned i=0;i<DIM;i++)
760 // {
761 // outfile << mesh_veloc[i] << " ";
762 // }
763 
764  // Dissipation
765  outfile << dissipation(s) << " ";
766 
767 
768  outfile << std::endl;
769 
770  }
771 
772  // Write tecplot footer (e.g. FE connectivity lists)
773  write_tecplot_zone_footer(outfile,nplot);
774 }
775 
776 
777 //==============================================================
778 /// Output function for vorticity.
779 /// x,y,[z],[omega_x,omega_y,[and/or omega_z]]
780 /// in tecplot format. Specified number of plot points in each
781 /// coordinate direction.
782 //==============================================================
783 template<unsigned DIM>
785  const unsigned &nplot)
786 {
787 
788  //Vector of local coordinates
789  Vector<double> s(DIM);
790 
791  // Create vorticity vector of the required size
792  Vector<double> vorticity;
793  if (DIM==2)
794  {
795  vorticity.resize(1);
796  }
797  else if (DIM==3)
798  {
799  vorticity.resize(3);
800  }
801  else
802  {
803  std::string error_message =
804  "Can't output vorticity in 1D - in fact, what do you mean\n";
805  error_message += "by the 1D Navier-Stokes equations?\n";
806 
807  throw OomphLibError(error_message,
808  OOMPH_CURRENT_FUNCTION,
809  OOMPH_EXCEPTION_LOCATION);
810  }
811 
812  // Tecplot header info
813  outfile << tecplot_zone_string(nplot);
814 
815  // Loop over plot points
816  unsigned num_plot_points=nplot_points(nplot);
817  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
818  {
819 
820  // Get local coordinates of plot point
821  get_s_plot(iplot,nplot,s);
822 
823  // Coordinates
824  for(unsigned i=0;i<DIM;i++)
825  {
826  outfile << interpolated_x(s,i) << " ";
827  }
828 
829  // Get vorticity
830  get_vorticity(s,vorticity);
831 
832  if (DIM==2)
833  {
834  outfile << vorticity[0];
835  }
836  else
837  {
838  outfile << vorticity[0] << " "
839  << vorticity[1] << " "
840  << vorticity[2] << " ";
841  }
842 
843  outfile << std::endl;
844  }
845  outfile << std::endl;
846 
847  // Write tecplot footer (e.g. FE connectivity lists)
848  write_tecplot_zone_footer(outfile,nplot);
849 
850 }
851 
852 
853 
854 //==============================================================
855 /// Return integral of dissipation over element
856 //==============================================================
857 template<unsigned DIM>
859 {
860 
861  // Initialise
862  double diss=0.0;
863 
864  //Set the value of n_intpt
865  unsigned n_intpt = integral_pt()->nweight();
866 
867  //Set the Vector to hold local coordinates
868  Vector<double> s(DIM);
869 
870  //Loop over the integration points
871  for(unsigned ipt=0;ipt<n_intpt;ipt++)
872  {
873 
874  //Assign values of s
875  for(unsigned i=0;i<DIM;i++)
876  {
877  s[i] = integral_pt()->knot(ipt,i);
878  }
879 
880  //Get the integral weight
881  double w = integral_pt()->weight(ipt);
882 
883  // Get Jacobian of mapping
884  double J = J_eulerian(s);
885 
886  // Get strain rate matrix
887  DenseMatrix<double> strainrate(DIM,DIM);
888  strain_rate(s,strainrate);
889 
890  // Initialise
891  double local_diss=0.0;
892  for(unsigned i=0;i<DIM;i++)
893  {
894  for(unsigned j=0;j<DIM;j++)
895  {
896  local_diss+=2.0*strainrate(i,j)*strainrate(i,j);
897  }
898  }
899 
900  diss+=local_diss*w*J;
901  }
902 
903  return diss;
904 
905 }
906 
907 //==============================================================
908 /// Compute traction (on the viscous scale) exerted onto
909 /// the fluid at local coordinate s. N has to be outer unit normal
910 /// to the fluid.
911 //==============================================================
912 template<unsigned DIM>
914  const Vector<double>& N,
915  Vector<double>& traction)
916 {
917 
918  // Get velocity gradients
919  DenseMatrix<double> strainrate(DIM,DIM);
920  strain_rate(s,strainrate);
921 
922  // Get pressure
923  double press=interpolated_p_nst(s);
924 
925  // Loop over traction components
926  for (unsigned i=0;i<DIM;i++)
927  {
928  traction[i]=-press*N[i];
929  for (unsigned j=0;j<DIM;j++)
930  {
931  traction[i]+=2.0*strainrate(i,j)*N[j];
932  }
933  }
934 
935 }
936 
937 //==============================================================
938 /// Compute traction (on the viscous scale) exerted onto
939 /// the fluid at local coordinate s, decomposed into pressure and
940 /// normal and tangential viscous components.
941 /// N has to be outer unit normal to the fluid.
942 //==============================================================
943 template<unsigned DIM>
945  const Vector<double>& N,
946  Vector<double>& traction_p,
947  Vector<double>& traction_visc_n,
948  Vector<double>& traction_visc_t)
949 {
950  Vector<double> traction_visc(DIM);
951 
952  // Get velocity gradients
953  DenseMatrix<double> strainrate(DIM,DIM);
954  strain_rate(s,strainrate);
955 
956  // Get pressure
957  double press=interpolated_p_nst(s);
958 
959  // Loop over traction components
960  for (unsigned i=0;i<DIM;i++)
961  {
962  // pressure
963  traction_p[i]=-press*N[i];
964  for (unsigned j=0;j<DIM;j++)
965  {
966  // viscous stress
967  traction_visc[i]+=2.0*strainrate(i,j)*N[j];
968  }
969  // Decompose viscous stress into normal and tangential components
970  traction_visc_n[i]=traction_visc[i]*N[i];
971  traction_visc_t[i]=traction_visc[i]-traction_visc_n[i];
972  }
973 
974 }
975 
976 //==============================================================
977 /// Return dissipation at local coordinate s
978 //==============================================================
979 template<unsigned DIM>
981 {
982  // Get strain rate matrix
983  DenseMatrix<double> strainrate(DIM,DIM);
984  strain_rate(s,strainrate);
985 
986  // Initialise
987  double local_diss=0.0;
988  for(unsigned i=0;i<DIM;i++)
989  {
990  for(unsigned j=0;j<DIM;j++)
991  {
992  local_diss+=2.0*strainrate(i,j)*strainrate(i,j);
993  }
994  }
995 
996  return local_diss;
997 }
998 
999 //==============================================================
1000 /// Get strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)
1001 //==============================================================
1002 template<unsigned DIM>
1004  DenseMatrix<double>& strainrate) const
1005  {
1006 
1007 #ifdef PARANOID
1008  if ((strainrate.ncol()!=DIM)||(strainrate.nrow()!=DIM))
1009  {
1010  std::ostringstream error_message;
1011  error_message << "The strain rate has incorrect dimensions "
1012  << strainrate.ncol() << " x "
1013  << strainrate.nrow() << " Not " << DIM << std::endl;
1014 
1015  throw OomphLibError(error_message.str(),
1016  OOMPH_CURRENT_FUNCTION,
1017  OOMPH_EXCEPTION_LOCATION);
1018  }
1019 #endif
1020 
1021  // Velocity gradient matrix
1022  DenseMatrix<double> dudx(DIM);
1023 
1024  //Find out how many nodes there are in the element
1025  unsigned n_node = nnode();
1026 
1027  //Set up memory for the shape and test functions
1028  Shape psi(n_node);
1029  DShape dpsidx(n_node,DIM);
1030 
1031  //Call the derivatives of the shape functions
1032  dshape_eulerian(s,psi,dpsidx);
1033 
1034  //Initialise to zero
1035  for(unsigned i=0;i<DIM;i++)
1036  {
1037  for(unsigned j=0;j<DIM;j++)
1038  {
1039  dudx(i,j) = 0.0;
1040  }
1041  }
1042 
1043  //Loop over veclocity components
1044  for(unsigned i=0;i<DIM;i++)
1045  {
1046  //Get the index at which the i-th velocity is stored
1047  unsigned u_nodal_index = u_index_nst(i);
1048  //Loop over derivative directions
1049  for(unsigned j=0;j<DIM;j++)
1050  {
1051  // Loop over nodes
1052  for(unsigned l=0;l<n_node;l++)
1053  {
1054  dudx(i,j) += nodal_value(l,u_nodal_index)*dpsidx(l,j);
1055  }
1056  }
1057  }
1058 
1059  //Loop over veclocity components
1060  for(unsigned i=0;i<DIM;i++)
1061  {
1062  //Loop over derivative directions
1063  for(unsigned j=0;j<DIM;j++)
1064  {
1065  strainrate(i,j)=0.5*(dudx(i,j)+dudx(j,i));
1066  }
1067  }
1068 
1069  }
1070 
1071 
1072 
1073 //==============================================================
1074 /// Compute 2D vorticity vector at local coordinate s (return in
1075 /// one and only component of vorticity vector
1076 //==============================================================
1077 template<>
1079  Vector<double>& vorticity) const
1080 {
1081 
1082 #ifdef PARANOID
1083  if (vorticity.size()!=1)
1084  {
1085  std::ostringstream error_message;
1086  error_message
1087  << "Computation of vorticity in 2D requires a 1D vector\n"
1088  << "which contains the only non-zero component of the\n"
1089  << "vorticity vector. You've passed a vector of size "
1090  << vorticity.size() << std::endl;
1091 
1092  throw OomphLibError(error_message.str(),
1093  OOMPH_CURRENT_FUNCTION,
1094  OOMPH_EXCEPTION_LOCATION);
1095  }
1096 #endif
1097 
1098  // Specify spatial dimension
1099  unsigned DIM=2;
1100 
1101  // Velocity gradient matrix
1102  DenseMatrix<double> dudx(DIM);
1103 
1104  //Find out how many nodes there are in the element
1105  unsigned n_node = nnode();
1106 
1107  //Set up memory for the shape and test functions
1108  Shape psi(n_node);
1109  DShape dpsidx(n_node,DIM);
1110 
1111  //Call the derivatives of the shape functions
1112  dshape_eulerian(s,psi,dpsidx);
1113 
1114  //Initialise to zero
1115  for(unsigned i=0;i<DIM;i++)
1116  {
1117  for(unsigned j=0;j<DIM;j++)
1118  {
1119  dudx(i,j) = 0.0;
1120  }
1121  }
1122 
1123  //Loop over veclocity components
1124  for(unsigned i=0;i<DIM;i++)
1125  {
1126  //Get the index at which the i-th velocity is stored
1127  unsigned u_nodal_index = u_index_nst(i);
1128  //Loop over derivative directions
1129  for(unsigned j=0;j<DIM;j++)
1130  {
1131  // Loop over nodes
1132  for(unsigned l=0;l<n_node;l++)
1133  {
1134  dudx(i,j) += nodal_value(l,u_nodal_index)*dpsidx(l,j);
1135  }
1136  }
1137  }
1138 
1139  // Z-component of vorticity
1140  vorticity[0]=dudx(1,0)-dudx(0,1);
1141 
1142 }
1143 
1144 
1145 
1146 
1147 //==============================================================
1148 /// Compute 3D vorticity vector at local coordinate s
1149 //==============================================================
1150 template<>
1152  Vector<double>& vorticity) const
1153 {
1154 
1155 #ifdef PARANOID
1156  if (vorticity.size()!=3)
1157  {
1158  std::ostringstream error_message;
1159  error_message
1160  << "Computation of vorticity in 3D requires a 3D vector\n"
1161  << "which contains the only non-zero component of the\n"
1162  << "vorticity vector. You've passed a vector of size "
1163  << vorticity.size() << std::endl;
1164 
1165  throw OomphLibError(error_message.str(),
1166  OOMPH_CURRENT_FUNCTION,
1167  OOMPH_EXCEPTION_LOCATION);
1168  }
1169 #endif
1170 
1171  // Specify spatial dimension
1172  unsigned DIM=3;
1173 
1174  // Velocity gradient matrix
1175  DenseMatrix<double> dudx(DIM);
1176 
1177  //Find out how many nodes there are in the element
1178  unsigned n_node = nnode();
1179 
1180  //Set up memory for the shape and test functions
1181  Shape psi(n_node);
1182  DShape dpsidx(n_node,DIM);
1183 
1184  //Call the derivatives of the shape functions
1185  dshape_eulerian(s,psi,dpsidx);
1186 
1187  //Initialise to zero
1188  for(unsigned i=0;i<DIM;i++)
1189  {
1190  for(unsigned j=0;j<DIM;j++)
1191  {
1192  dudx(i,j) = 0.0;
1193  }
1194  }
1195 
1196  //Loop over veclocity components
1197  for(unsigned i=0;i<DIM;i++)
1198  {
1199  //Get the index at which the i-th velocity component is stored
1200  unsigned u_nodal_index = u_index_nst(i);
1201  //Loop over derivative directions
1202  for(unsigned j=0;j<DIM;j++)
1203  {
1204  // Loop over nodes
1205  for(unsigned l=0;l<n_node;l++)
1206  {
1207  dudx(i,j) += nodal_value(l,u_nodal_index)*dpsidx(l,j);
1208  }
1209  }
1210  }
1211 
1212  vorticity[0]=dudx(2,1)-dudx(1,2);
1213  vorticity[1]=dudx(0,2)-dudx(2,0);
1214  vorticity[2]=dudx(1,0)-dudx(0,1);
1215 
1216 }
1217 
1218 
1219 //==============================================================
1220 /// \short Get integral of kinetic energy over element:
1221 /// Note that this is the "raw" kinetic energy in the sense
1222 /// that the density ratio has not been included. In problems
1223 /// with two or more fluids the user will have to remember
1224 /// to premultiply certain elements by the appropriate density
1225 /// ratio.
1226 //==============================================================
1227 template<unsigned DIM>
1229 {
1230 
1231  // Initialise
1232  double kin_en=0.0;
1233 
1234  //Set the value of n_intpt
1235  unsigned n_intpt = integral_pt()->nweight();
1236 
1237  //Set the Vector to hold local coordinates
1238  Vector<double> s(DIM);
1239 
1240  //Loop over the integration points
1241  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1242  {
1243  //Assign values of s
1244  for(unsigned i=0;i<DIM;i++)
1245  {
1246  s[i] = integral_pt()->knot(ipt,i);
1247  }
1248 
1249  //Get the integral weight
1250  double w = integral_pt()->weight(ipt);
1251 
1252  //Get Jacobian of mapping
1253  double J = J_eulerian(s);
1254 
1255  // Loop over directions
1256  double veloc_squared=0.0;
1257  for(unsigned i=0;i<DIM;i++)
1258  {
1259  veloc_squared+=interpolated_u_nst(s,i)*interpolated_u_nst(s,i);
1260  }
1261 
1262  kin_en+=0.5*veloc_squared*w*J;
1263  }
1264 
1265  return kin_en;
1266 
1267 }
1268 
1269 
1270 //==========================================================================
1271 /// \short Get integral of time derivative of kinetic energy over element:
1272 //==========================================================================
1273 template<unsigned DIM>
1275 {
1276  // Initialise
1277  double d_kin_en_dt=0.0;
1278 
1279  //Set the value of n_intpt
1280  const unsigned n_intpt = integral_pt()->nweight();
1281 
1282  //Set the Vector to hold local coordinates
1283  Vector<double> s(DIM);
1284 
1285  //Get the number of nodes
1286  const unsigned n_node = this->nnode();
1287 
1288  //Storage for the shape function
1289  Shape psi(n_node);
1290  DShape dpsidx(n_node,DIM);
1291 
1292  //Get the value at which the velocities are stored
1293  unsigned u_index[DIM];
1294  for(unsigned i=0;i<DIM;i++) {u_index[i] = this->u_index_nst(i);}
1295 
1296  //Loop over the integration points
1297  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1298  {
1299  //Get the jacobian of the mapping
1300  double J = dshape_eulerian_at_knot(ipt,psi,dpsidx);
1301 
1302  //Get the integral weight
1303  double w = integral_pt()->weight(ipt);
1304 
1305  //Now work out the velocity and the time derivative
1306  Vector<double> interpolated_u(DIM,0.0), interpolated_dudt(DIM,0.0);
1307 
1308  //Loop over the shape functions
1309  for(unsigned l=0;l<n_node;l++)
1310  {
1311  //Loop over the dimensions
1312  for(unsigned i=0;i<DIM;i++)
1313  {
1314  interpolated_u[i] += nodal_value(l,u_index[i])*psi(l);
1315  interpolated_dudt[i] += du_dt_nst(l,u_index[i])*psi(l);
1316  }
1317  }
1318 
1319  //Get mesh velocity bit
1320  if (!ALE_is_disabled)
1321  {
1322  Vector<double> mesh_velocity(DIM,0.0);
1323  DenseMatrix<double> interpolated_dudx(DIM,DIM,0.0);
1324 
1325  // Loop over nodes again
1326  for(unsigned l=0;l<n_node;l++)
1327  {
1328  for(unsigned i=0;i<DIM;i++)
1329  {
1330  mesh_velocity[i] += this->dnodal_position_dt(l,i)*psi(l);
1331 
1332  for(unsigned j=0;j<DIM;j++)
1333  {
1334  interpolated_dudx(i,j) +=
1335  this->nodal_value(l,u_index[i])*dpsidx(l,j);
1336  }
1337  }
1338  }
1339 
1340  //Subtract mesh velocity from du_dt
1341  for(unsigned i=0;i<DIM;i++)
1342  {
1343  for (unsigned k=0;k<DIM;k++)
1344  {
1345  interpolated_dudt[i] -= mesh_velocity[k]*interpolated_dudx(i,k);
1346  }
1347  }
1348  }
1349 
1350 
1351  // Loop over directions and add up u du/dt terms
1352  double sum=0.0;
1353  for(unsigned i=0;i<DIM;i++)
1354  { sum += interpolated_u[i]*interpolated_dudt[i];}
1355 
1356  d_kin_en_dt += sum*w*J;
1357  }
1358 
1359  return d_kin_en_dt;
1360 
1361 }
1362 
1363 
1364 //==============================================================
1365 /// Return pressure integrated over the element
1366 //==============================================================
1367 template<unsigned DIM>
1369 {
1370 
1371  // Initialise
1372  double press_int=0;
1373 
1374  //Set the value of n_intpt
1375  unsigned n_intpt = integral_pt()->nweight();
1376 
1377  //Set the Vector to hold local coordinates
1378  Vector<double> s(DIM);
1379 
1380  //Loop over the integration points
1381  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1382  {
1383 
1384  //Assign values of s
1385  for(unsigned i=0;i<DIM;i++)
1386  {
1387  s[i] = integral_pt()->knot(ipt,i);
1388  }
1389 
1390  //Get the integral weight
1391  double w = integral_pt()->weight(ipt);
1392 
1393  //Get Jacobian of mapping
1394  double J = J_eulerian(s);
1395 
1396  //Premultiply the weights and the Jacobian
1397  double W = w*J;
1398 
1399  // Get pressure
1400  double press=interpolated_p_nst(s);
1401 
1402  // Add
1403  press_int+=press*W;
1404 
1405  }
1406 
1407  return press_int;
1408 
1409 }
1410 
1411 
1412 //==============================================================
1413 /// Compute the residuals for the associated pressure advection
1414 /// diffusion problem. Used by the Fp preconditioner.
1415 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
1416 //==============================================================
1417 template<unsigned DIM>
1420  Vector<double> &residuals,
1421  DenseMatrix<double> &jacobian,
1422  unsigned flag)
1423 {
1424  // Return immediately if there are no dofs
1425  if (ndof()==0) return;
1426 
1427  //Find out how many nodes there are
1428  unsigned n_node = nnode();
1429 
1430  //Find out how many pressure dofs there are
1431  unsigned n_pres = npres_nst();
1432 
1433  //Find the indices at which the local velocities are stored
1434  unsigned u_nodal_index[DIM];
1435  for(unsigned i=0;i<DIM;i++) {u_nodal_index[i] = u_index_nst(i);}
1436 
1437  //Set up memory for the velocity shape fcts
1438  Shape psif(n_node);
1439  DShape dpsidx(n_node,DIM);
1440 
1441  //Set up memory for pressure shape and test functions
1442  Shape psip(n_pres), testp(n_pres);
1443  DShape dpsip(n_pres,DIM);
1444  DShape dtestp(n_pres,DIM);
1445 
1446  //Number of integration points
1447  unsigned n_intpt = integral_pt()->nweight();
1448 
1449  //Set the Vector to hold local coordinates
1450  Vector<double> s(DIM);
1451 
1452  //Get Physical Variables from Element
1453  //Reynolds number must be multiplied by the density ratio
1454  double scaled_re = re()*density_ratio();
1455 
1456  //Integers to store the local equations and unknowns
1457  int local_eqn=0, local_unknown=0;
1458 
1459  //Loop over the integration points
1460  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1461  {
1462  //Assign values of s
1463  for(unsigned i=0;i<DIM;i++) s[i] = integral_pt()->knot(ipt,i);
1464 
1465  //Get the integral weight
1466  double w = integral_pt()->weight(ipt);
1467 
1468  // Call the derivatives of the veloc shape functions
1469  // (Derivs not needed but they are free)
1470  double J = this->dshape_eulerian_at_knot(ipt,psif,dpsidx);
1471 
1472  //Call the pressure shape and test functions
1473  this->dpshape_and_dptest_eulerian_nst(s,psip,dpsip,testp,dtestp);
1474 
1475  //Premultiply the weights and the Jacobian
1476  double W = w*J;
1477 
1478  //Calculate local values of the pressure and velocity components
1479  //Allocate
1480  Vector<double> interpolated_u(DIM,0.0);
1481  Vector<double> interpolated_x(DIM,0.0);
1482  Vector<double> interpolated_dpdx(DIM,0.0);
1483 
1484  //Calculate pressure gradient
1485  for(unsigned l=0;l<n_pres;l++)
1486  {
1487  for (unsigned i=0;i<DIM;i++)
1488  {
1489  interpolated_dpdx[i] += p_nst(l)*dpsip(l,i);
1490  }
1491  }
1492 
1493  //Calculate velocities
1494 
1495  // Loop over nodes
1496  for(unsigned l=0;l<n_node;l++)
1497  {
1498  //Loop over directions
1499  for(unsigned i=0;i<DIM;i++)
1500  {
1501  //Get the nodal value
1502  double u_value = raw_nodal_value(l,u_nodal_index[i]);
1503  interpolated_u[i] += u_value*psif[l];
1504  interpolated_x[i] += raw_nodal_position(l,i)*psif[l];
1505  }
1506  }
1507 
1508  // Source function (for validaton only)
1509  double source=0.0;
1510  if (Press_adv_diff_source_fct_pt!=0)
1511  {
1512  source=Press_adv_diff_source_fct_pt(interpolated_x);
1513  }
1514 
1515  //Loop over the shape functions
1516  for(unsigned l=0;l<n_pres;l++)
1517  {
1518  local_eqn = p_local_eqn(l);
1519 
1520  //If not a boundary conditions
1521  if(local_eqn >= 0)
1522  {
1523  residuals[local_eqn] -= source*testp[l]*W;
1524  for(unsigned k=0;k<DIM;k++)
1525  {
1526  residuals[local_eqn] += interpolated_dpdx[k]*
1527  (scaled_re*interpolated_u[k]*testp[l]+dtestp(l,k))*W;
1528  }
1529 
1530  // Jacobian too?
1531  if(flag)
1532  {
1533  //Loop over the shape functions
1534  for(unsigned l2=0;l2<n_pres;l2++)
1535  {
1536  local_unknown = p_local_eqn(l2);
1537 
1538  //If not a boundary conditions
1539  if(local_unknown >= 0)
1540  {
1541  if ((int(eqn_number(local_eqn))!=Pinned_fp_pressure_eqn)&&
1542  (int(eqn_number(local_unknown))!=Pinned_fp_pressure_eqn))
1543  {
1544  for(unsigned k=0;k<DIM;k++)
1545  {
1546  jacobian(local_eqn,local_unknown)+=dtestp(l2,k)*
1547  (scaled_re*interpolated_u[k]*testp[l]+dtestp(l,k))*W;
1548  }
1549  }
1550  else
1551  {
1552  if ((int(eqn_number(local_eqn))==Pinned_fp_pressure_eqn)&&
1553  (int(eqn_number(local_unknown))==Pinned_fp_pressure_eqn))
1554  {
1555  jacobian(local_eqn,local_unknown)=1.0;
1556  }
1557  }
1558  }
1559  }
1560  } /*End of Jacobian calculation*/
1561  } //End of if not boundary condition
1562  }//End of loop over l
1563  }
1564 
1565  // Now add boundary contributions from Robin BCs
1566  unsigned nrobin=Pressure_advection_diffusion_robin_element_pt.size();
1567  for (unsigned e=0;e<nrobin;e++)
1568  {
1569  Pressure_advection_diffusion_robin_element_pt[e]->
1570  fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(
1571  residuals,jacobian,flag);
1572  }
1573 }
1574 
1575 
1576 //==============================================================
1577 /// Compute the residuals for the Navier--Stokes
1578 /// equations; flag=1(or 0): do (or don't) compute the
1579 /// Jacobian as well.
1580 //==============================================================
1581 template<unsigned DIM>
1584  DenseMatrix<double> &jacobian,
1585  DenseMatrix<double> &mass_matrix,
1586  unsigned flag)
1587 {
1588  // Return immediately if there are no dofs
1589  if (ndof()==0) return;
1590 
1591  //Find out how many nodes there are
1592  unsigned n_node = nnode();
1593 
1594  // Get continuous time from timestepper of first node
1595  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
1596 
1597  //Find out how many pressure dofs there are
1598  unsigned n_pres = npres_nst();
1599 
1600  //Find the indices at which the local velocities are stored
1601  unsigned u_nodal_index[DIM];
1602  for(unsigned i=0;i<DIM;i++) {u_nodal_index[i] = u_index_nst(i);}
1603 
1604  //Set up memory for the shape and test functions
1605  Shape psif(n_node), testf(n_node);
1606  DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
1607 
1608  //Set up memory for pressure shape and test functions
1609  Shape psip(n_pres), testp(n_pres);
1610 
1611  //Number of integration points
1612  unsigned n_intpt = integral_pt()->nweight();
1613 
1614  //Set the Vector to hold local coordinates
1615  Vector<double> s(DIM);
1616 
1617  //Get Physical Variables from Element
1618  //Reynolds number must be multiplied by the density ratio
1619  double scaled_re = re()*density_ratio();
1620  double scaled_re_st = re_st()*density_ratio();
1621  double scaled_re_inv_fr = re_invfr()*density_ratio();
1622  double visc_ratio = viscosity_ratio();
1623  Vector<double> G = g();
1624 
1625  //Integers to store the local equations and unknowns
1626  int local_eqn=0, local_unknown=0;
1627 
1628  //Loop over the integration points
1629  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1630  {
1631  //Assign values of s
1632  for(unsigned i=0;i<DIM;i++) s[i] = integral_pt()->knot(ipt,i);
1633  //Get the integral weight
1634  double w = integral_pt()->weight(ipt);
1635 
1636  //Call the derivatives of the shape and test functions
1637  double J =
1638  dshape_and_dtest_eulerian_at_knot_nst(ipt,psif,dpsifdx,testf,dtestfdx);
1639 
1640  //Call the pressure shape and test functions
1641  pshape_nst(s,psip,testp);
1642 
1643  //Premultiply the weights and the Jacobian
1644  double W = w*J;
1645 
1646  //Calculate local values of the pressure and velocity components
1647  //Allocate
1648  double interpolated_p=0.0;
1649  Vector<double> interpolated_u(DIM,0.0);
1650  Vector<double> interpolated_x(DIM,0.0);
1651  Vector<double> mesh_velocity(DIM,0.0);
1652  Vector<double> dudt(DIM,0.0);
1653  DenseMatrix<double> interpolated_dudx(DIM,DIM,0.0);
1654 
1655  //Calculate pressure
1656  for(unsigned l=0;l<n_pres;l++) interpolated_p += p_nst(l)*psip[l];
1657 
1658  //Calculate velocities and derivatives:
1659 
1660  // Loop over nodes
1661  for(unsigned l=0;l<n_node;l++)
1662  {
1663  //Loop over directions
1664  for(unsigned i=0;i<DIM;i++)
1665  {
1666  //Get the nodal value
1667  double u_value = raw_nodal_value(l,u_nodal_index[i]);
1668  interpolated_u[i] += u_value*psif[l];
1669  interpolated_x[i] += raw_nodal_position(l,i)*psif[l];
1670  dudt[i] += du_dt_nst(l,i)*psif[l];
1671 
1672  //Loop over derivative directions
1673  for(unsigned j=0;j<DIM;j++)
1674  {
1675  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
1676  }
1677  }
1678  }
1679 
1680  if (!ALE_is_disabled)
1681  {
1682  // Loop over nodes
1683  for(unsigned l=0;l<n_node;l++)
1684  {
1685  //Loop over directions
1686  for(unsigned i=0;i<DIM;i++)
1687  {
1688  mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif[l];
1689  }
1690  }
1691  }
1692 
1693  //Get the user-defined body force terms
1694  Vector<double> body_force(DIM);
1695  get_body_force_nst(time,ipt,s,interpolated_x,body_force);
1696 
1697  //Get the user-defined source function
1698  double source = get_source_nst(time,ipt,interpolated_x);
1699 
1700 
1701  //MOMENTUM EQUATIONS
1702  //------------------
1703 
1704  //Loop over the test functions
1705  for(unsigned l=0;l<n_node;l++)
1706  {
1707  //Loop over the velocity components
1708  for(unsigned i=0;i<DIM;i++)
1709  {
1710  /*IF it's not a boundary condition*/
1711  local_eqn = nodal_local_eqn(l,u_nodal_index[i]);
1712  if(local_eqn >= 0)
1713  {
1714  //Add the user-defined body force terms
1715  residuals[local_eqn] +=
1716  body_force[i]*testf[l]*W;
1717 
1718  //Add the gravitational body force term
1719  residuals[local_eqn] += scaled_re_inv_fr*testf[l]*G[i]*W;
1720 
1721  //Add the pressure gradient term
1722  residuals[local_eqn] += interpolated_p*dtestfdx(l,i)*W;
1723 
1724  //Add in the stress tensor terms
1725  //The viscosity ratio needs to go in here to ensure
1726  //continuity of normal stress is satisfied even in flows
1727  //with zero pressure gradient!
1728  for(unsigned k=0;k<DIM;k++)
1729  {
1730  residuals[local_eqn] -= visc_ratio*
1731  (interpolated_dudx(i,k) + Gamma[i]*interpolated_dudx(k,i))
1732  *dtestfdx(l,k)*W;
1733  }
1734 
1735  //Add in the inertial terms
1736  //du/dt term
1737  residuals[local_eqn] -= scaled_re_st*dudt[i]*testf[l]*W;
1738 
1739 
1740  //Convective terms, including mesh velocity
1741  for(unsigned k=0;k<DIM;k++)
1742  {
1743  double tmp=scaled_re*interpolated_u[k];
1744  if (!ALE_is_disabled) tmp-=scaled_re_st*mesh_velocity[k];
1745  residuals[local_eqn] -= tmp*interpolated_dudx(i,k)*testf[l]*W;
1746 
1747  }
1748 
1749  //CALCULATE THE JACOBIAN
1750  if(flag)
1751  {
1752  //Loop over the velocity shape functions again
1753  for(unsigned l2=0;l2<n_node;l2++)
1754  {
1755  //Loop over the velocity components again
1756  for(unsigned i2=0;i2<DIM;i2++)
1757  {
1758  //If at a non-zero degree of freedom add in the entry
1759  local_unknown = nodal_local_eqn(l2,u_nodal_index[i2]);
1760  if(local_unknown >= 0)
1761  {
1762  //Add contribution to Elemental Matrix
1763  jacobian(local_eqn,local_unknown)
1764  -= visc_ratio*Gamma[i]*dpsifdx(l2,i)*dtestfdx(l,i2)*W;
1765 
1766  //Extra component if i2 = i
1767  if(i2 == i)
1768  {
1769  /*Loop over velocity components*/
1770  for(unsigned k=0;k<DIM;k++)
1771  {
1772  jacobian(local_eqn,local_unknown)
1773  -= visc_ratio*dpsifdx(l2,k)*dtestfdx(l,k)*W;
1774  }
1775  }
1776 
1777  //Now add in the inertial terms
1778  jacobian(local_eqn,local_unknown)
1779  -= scaled_re*psif[l2]*interpolated_dudx(i,i2)*testf[l]*W;
1780 
1781  //Extra component if i2=i
1782  if(i2 == i)
1783  {
1784  //Add the mass matrix term (only diagonal entries)
1785  //Note that this is positive because the mass matrix
1786  //is taken to the other side of the equation when
1787  //formulating the generalised eigenproblem.
1788  if(flag==2)
1789  {
1790  mass_matrix(local_eqn,local_unknown) +=
1791  scaled_re_st*psif[l2]*testf[l]*W;
1792  }
1793 
1794  //du/dt term
1795  jacobian(local_eqn,local_unknown)
1796  -= scaled_re_st*
1797  node_pt(l2)->time_stepper_pt()->weight(1,0)*
1798  psif[l2]*testf[l]*W;
1799 
1800  //Loop over the velocity components
1801  for(unsigned k=0;k<DIM;k++)
1802  {
1803  double tmp=scaled_re*interpolated_u[k];
1804  if (!ALE_is_disabled) tmp-=scaled_re_st*mesh_velocity[k];
1805  jacobian(local_eqn,local_unknown) -=
1806  tmp*dpsifdx(l2,k)*testf[l]*W;
1807  }
1808  }
1809 
1810  }
1811  }
1812  }
1813 
1814  /*Now loop over pressure shape functions*/
1815  /*This is the contribution from pressure gradient*/
1816  for(unsigned l2=0;l2<n_pres;l2++)
1817  {
1818  /*If we are at a non-zero degree of freedom in the entry*/
1819  local_unknown = p_local_eqn(l2);
1820  if(local_unknown >= 0)
1821  {
1822  jacobian(local_eqn,local_unknown)
1823  += psip[l2]*dtestfdx(l,i)*W;
1824  }
1825  }
1826  } /*End of Jacobian calculation*/
1827 
1828  } //End of if not boundary condition statement
1829 
1830  } //End of loop over dimension
1831  } //End of loop over shape functions
1832 
1833 
1834 
1835  //CONTINUITY EQUATION
1836  //-------------------
1837 
1838  //Loop over the shape functions
1839  for(unsigned l=0;l<n_pres;l++)
1840  {
1841  local_eqn = p_local_eqn(l);
1842  //If not a boundary conditions
1843  if(local_eqn >= 0)
1844  {
1845 
1846  // Source term
1847  //residuals[local_eqn] -=source*testp[l]*W;
1848  double aux=-source;
1849 
1850  //Loop over velocity components
1851  for(unsigned k=0;k<DIM;k++)
1852  {
1853  //residuals[local_eqn] += interpolated_dudx(k,k)*testp[l]*W;
1854  aux += interpolated_dudx(k,k);
1855  }
1856 
1857  residuals[local_eqn]+=aux*testp[l]*W;
1858 
1859  /*CALCULATE THE JACOBIAN*/
1860  if(flag)
1861  {
1862  /*Loop over the velocity shape functions*/
1863  for(unsigned l2=0;l2<n_node;l2++)
1864  {
1865  /*Loop over velocity components*/
1866  for(unsigned i2=0;i2<DIM;i2++)
1867  {
1868  /*If we're at a non-zero degree of freedom add it in*/
1869  local_unknown = nodal_local_eqn(l2,u_nodal_index[i2]);
1870  if(local_unknown >= 0)
1871  {
1872  jacobian(local_eqn,local_unknown)
1873  += dpsifdx(l2,i2)*testp[l]*W;
1874  }
1875  } /*End of loop over i2*/
1876  } /*End of loop over l2*/
1877  } /*End of Jacobian calculation*/
1878 
1879  } //End of if not boundary condition
1880  } //End of loop over l
1881  }
1882 
1883 }
1884 
1885 //==============================================================
1886 /// Compute the derivatives of the residuals for the Navier--Stokes
1887 /// equations with respect to a parameter;
1888 /// flag=2 or 1 or 0: do (or don't) compute the
1889 /// Jacobian and mass matrix as well
1890 //==============================================================
1891 template<unsigned DIM>
1894  double* const &parameter_pt,
1895  Vector<double> &dres_dparam,
1896  DenseMatrix<double> &djac_dparam,
1897  DenseMatrix<double> &dmass_matrix_dparam,
1898  unsigned flag)
1899 {
1900  throw OomphLibError("Not yet implemented\n",
1901  OOMPH_CURRENT_FUNCTION,
1902  OOMPH_EXCEPTION_LOCATION);
1903 }
1904 
1905 //==================================================================
1906 /// \short Compute the hessian tensor vector products required to
1907 /// perform continuation of bifurcations analytically
1908 //================================================================
1909 template<unsigned DIM>
1912  Vector<double> const &Y,
1913  DenseMatrix<double> const &C,
1914  DenseMatrix<double> &product)
1915 {
1916  throw OomphLibError(
1917  "Not yet implemented\n",
1918  OOMPH_CURRENT_FUNCTION,
1919  OOMPH_EXCEPTION_LOCATION);
1920 }
1921 
1922 
1923 //======================================================================
1924 /// Compute derivatives of elemental residual vector with respect
1925 /// to nodal coordinates.
1926 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
1927 /// Overloads the FD-based version in the FE base class.
1928 //======================================================================
1929 template <unsigned DIM>
1932  dresidual_dnodal_coordinates)
1933 {
1934  // Return immediately if there are no dofs
1935  if(ndof()==0) { return; }
1936 
1937  // Determine number of nodes in element
1938  const unsigned n_node = nnode();
1939 
1940  // Get continuous time from timestepper of first node
1941  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
1942 
1943  // Determine number of pressure dofs in element
1944  const unsigned n_pres = npres_nst();
1945 
1946  // Find the indices at which the local velocities are stored
1947  unsigned u_nodal_index[DIM];
1948  for(unsigned i=0;i<DIM;i++) { u_nodal_index[i] = u_index_nst(i); }
1949 
1950  // Set up memory for the shape and test functions
1951  Shape psif(n_node), testf(n_node);
1952  DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
1953 
1954  // Set up memory for pressure shape and test functions
1955  Shape psip(n_pres), testp(n_pres);
1956 
1957  // Deriatives of shape fct derivatives w.r.t. nodal coords
1958  RankFourTensor<double> d_dpsifdx_dX(DIM,n_node,n_node,DIM);
1959  RankFourTensor<double> d_dtestfdx_dX(DIM,n_node,n_node,DIM);
1960 
1961  // Derivative of Jacobian of mapping w.r.t. to nodal coords
1962  DenseMatrix<double> dJ_dX(DIM,n_node);
1963 
1964  // Derivatives of derivative of u w.r.t. nodal coords
1965  RankFourTensor<double> d_dudx_dX(DIM,n_node,DIM,DIM);
1966 
1967  // Derivatives of nodal velocities w.r.t. nodal coords:
1968  // Assumption: Interaction only local via no-slip so
1969  // X_ij only affects U_ij.
1970  DenseMatrix<double> d_U_dX(DIM,n_node,0.0);
1971 
1972  // Determine the number of integration points
1973  const unsigned n_intpt = integral_pt()->nweight();
1974 
1975  // Vector to hold local coordinates
1976  Vector<double> s(DIM);
1977 
1978  // Get physical variables from element
1979  // (Reynolds number must be multiplied by the density ratio)
1980  double scaled_re = re()*density_ratio();
1981  double scaled_re_st = re_st()*density_ratio();
1982  double scaled_re_inv_fr = re_invfr()*density_ratio();
1983  double visc_ratio = viscosity_ratio();
1984  Vector<double> G = g();
1985 
1986  // FD step
1988 
1989  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1990  // Assumption: Interaction only local via no-slip so
1991  // X_ij only affects U_ij.
1992  bool element_has_node_with_aux_node_update_fct=false;
1993  for (unsigned q=0;q<n_node;q++)
1994  {
1995  Node* nod_pt=node_pt(q);
1996 
1997  // Only compute if there's a node-update fct involved
1998  if (nod_pt->has_auxiliary_node_update_fct_pt())
1999  {
2000  element_has_node_with_aux_node_update_fct=true;
2001 
2002  // Current nodal velocity
2003  Vector<double> u_ref(DIM);
2004  for (unsigned i=0;i<DIM;i++)
2005  {
2006  u_ref[i]=*(nod_pt->value_pt(u_nodal_index[i]));
2007  }
2008 
2009  // FD
2010  for (unsigned p=0;p<DIM;p++)
2011  {
2012  // Make backup
2013  double backup=nod_pt->x(p);
2014 
2015  // Do FD step. No node update required as we're
2016  // attacking the coordinate directly...
2017  nod_pt->x(p)+=eps_fd;
2018 
2019  // Do auxiliary node update (to apply no slip)
2021 
2022  // Evaluate
2023  d_U_dX(p,q)=(*(nod_pt->value_pt(u_nodal_index[p]))-u_ref[p])/eps_fd;
2024 
2025  // Reset
2026  nod_pt->x(p)=backup;
2027 
2028  // Do auxiliary node update (to apply no slip)
2030  }
2031  }
2032  }
2033 
2034  // Integer to store the local equation number
2035  int local_eqn=0;
2036 
2037  // Loop over the integration points
2038  for(unsigned ipt=0;ipt<n_intpt;ipt++)
2039  {
2040  // Assign values of s
2041  for(unsigned i=0;i<DIM;i++) { s[i] = integral_pt()->knot(ipt,i); }
2042 
2043  // Get the integral weight
2044  const double w = integral_pt()->weight(ipt);
2045 
2046  // Call the derivatives of the shape and test functions
2047  const double J = dshape_and_dtest_eulerian_at_knot_nst(ipt,psif,dpsifdx,
2048  d_dpsifdx_dX,
2049  testf,dtestfdx,
2050  d_dtestfdx_dX,dJ_dX);
2051 
2052  // Call the pressure shape and test functions
2053  pshape_nst(s,psip,testp);
2054 
2055  // Calculate local values of the pressure and velocity components
2056  // Allocate
2057  double interpolated_p=0.0;
2058  Vector<double> interpolated_u(DIM,0.0);
2059  Vector<double> interpolated_x(DIM,0.0);
2060  Vector<double> mesh_velocity(DIM,0.0);
2061  Vector<double> dudt(DIM,0.0);
2062  DenseMatrix<double> interpolated_dudx(DIM,DIM,0.0);
2063 
2064  // Calculate pressure
2065  for(unsigned l=0;l<n_pres;l++) { interpolated_p += p_nst(l)*psip[l]; }
2066 
2067  // Calculate velocities and derivatives:
2068 
2069  // Loop over nodes
2070  for(unsigned l=0;l<n_node;l++)
2071  {
2072  // Loop over directions
2073  for(unsigned i=0;i<DIM;i++)
2074  {
2075  // Get the nodal value
2076  double u_value = raw_nodal_value(l,u_nodal_index[i]);
2077  interpolated_u[i] += u_value*psif[l];
2078  interpolated_x[i] += raw_nodal_position(l,i)*psif[l];
2079  dudt[i] += du_dt_nst(l,i)*psif[l];
2080 
2081  // Loop over derivative directions
2082  for(unsigned j=0;j<DIM;j++)
2083  {
2084  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
2085  }
2086  }
2087  }
2088 
2089  if(!ALE_is_disabled)
2090  {
2091  // Loop over nodes
2092  for(unsigned l=0;l<n_node;l++)
2093  {
2094  // Loop over directions
2095  for(unsigned i=0;i<DIM;i++)
2096  {
2097  mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif[l];
2098  }
2099  }
2100  }
2101 
2102  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
2103  for(unsigned q=0;q<n_node;q++)
2104  {
2105  // Loop over coordinate directions
2106  for(unsigned p=0;p<DIM;p++)
2107  {
2108  for(unsigned i=0;i<DIM;i++)
2109  {
2110  for(unsigned k=0;k<DIM;k++)
2111  {
2112  double aux=0.0;
2113  for(unsigned j=0;j<n_node;j++)
2114  {
2115  aux +=
2116  raw_nodal_value(j,u_nodal_index[i])*d_dpsifdx_dX(p,q,j,k);
2117  }
2118  d_dudx_dX(p,q,i,k) = aux;
2119  }
2120  }
2121  }
2122  }
2123 
2124  // Get weight of actual nodal position/value in computation of mesh
2125  // velocity from positional/value time stepper
2126  const double pos_time_weight
2127  = node_pt(0)->position_time_stepper_pt()->weight(1,0);
2128  const double val_time_weight = node_pt(0)->time_stepper_pt()->weight(1,0);
2129 
2130  // Get the user-defined body force terms
2131  Vector<double> body_force(DIM);
2132  get_body_force_nst(time,ipt,s,interpolated_x,body_force);
2133 
2134  // Get the user-defined source function
2135  const double source = get_source_nst(time,ipt,interpolated_x);
2136 
2137  // Note: Can use raw values and avoid bypassing hanging information
2138  // because this is the non-refineable version!
2139 
2140  // Get gradient of body force function
2141  DenseMatrix<double> d_body_force_dx(DIM,DIM,0.0);
2142  get_body_force_gradient_nst(
2143  time,ipt,s,interpolated_x,d_body_force_dx);
2144 
2145  // Get gradient of source function
2146  Vector<double> source_gradient(DIM,0.0);
2147  get_source_gradient_nst(time,ipt,interpolated_x,source_gradient);
2148 
2149 
2150  // Assemble shape derivatives
2151  // --------------------------
2152 
2153  // MOMENTUM EQUATIONS
2154  // ------------------
2155 
2156  // Loop over the test functions
2157  for(unsigned l=0;l<n_node;l++)
2158  {
2159 
2160  // Loop over coordinate directions
2161  for(unsigned i=0;i<DIM;i++)
2162  {
2163  //Get the local equation
2164  local_eqn = nodal_local_eqn(l,u_nodal_index[i]);;
2165 
2166  // IF it's not a boundary condition
2167  if(local_eqn >= 0)
2168  {
2169  // Loop over coordinate directions
2170  for (unsigned p=0;p<DIM;p++)
2171  {
2172  // Loop over nodes
2173  for (unsigned q=0;q<n_node;q++)
2174  {
2175 
2176  // Residual x deriv of Jacobian
2177  //-----------------------------
2178 
2179  //Add the user-defined body force terms
2180  double sum = body_force[i]*testf[l];
2181 
2182  //Add the gravitational body force term
2183  sum += scaled_re_inv_fr*testf[l]*G[i];
2184 
2185  //Add the pressure gradient term
2186  sum += interpolated_p*dtestfdx(l,i);
2187 
2188  //Add in the stress tensor terms
2189  //The viscosity ratio needs to go in here to ensure
2190  //continuity of normal stress is satisfied even in flows
2191  //with zero pressure gradient!
2192  for(unsigned k=0;k<DIM;k++)
2193  {
2194  sum -= visc_ratio*
2195  (interpolated_dudx(i,k) + Gamma[i]*interpolated_dudx(k,i))
2196  *dtestfdx(l,k);
2197  }
2198 
2199  // Add in the inertial terms
2200 
2201  // du/dt term
2202  sum -= scaled_re_st*dudt[i]*testf[l];
2203 
2204  // Convective terms, including mesh velocity
2205  for(unsigned k=0;k<DIM;k++)
2206  {
2207  double tmp=scaled_re*interpolated_u[k];
2208  if (!ALE_is_disabled) { tmp-=scaled_re_st*mesh_velocity[k]; }
2209  sum -= tmp*interpolated_dudx(i,k)*testf[l];
2210  }
2211 
2212  // Multiply through by deriv of Jacobian and integration weight
2213  dresidual_dnodal_coordinates(local_eqn,p,q)+=sum*dJ_dX(p,q)*w;
2214 
2215  // Derivative of residual x Jacobian
2216  //----------------------------------
2217 
2218  // Body force
2219  sum=d_body_force_dx(i,p)*psif(q)*testf(l);
2220 
2221  // Pressure gradient term
2222  sum += interpolated_p*d_dtestfdx_dX(p,q,l,i);
2223 
2224  // Viscous term
2225  for (unsigned k=0;k<DIM;k++)
2226  {
2227  sum -= visc_ratio*(
2228  (interpolated_dudx(i,k) + Gamma[i]*interpolated_dudx(k,i))
2229  *d_dtestfdx_dX(p,q,l,k)+
2230  (d_dudx_dX(p,q,i,k) + Gamma[i]*d_dudx_dX(p,q,k,i))
2231  *dtestfdx(l,k));
2232  }
2233 
2234  // Convective terms, including mesh velocity
2235  for(unsigned k=0;k<DIM;k++)
2236  {
2237  double tmp=scaled_re*interpolated_u[k];
2238  if(!ALE_is_disabled) { tmp-=scaled_re_st*mesh_velocity[k]; }
2239  sum -= tmp*d_dudx_dX(p,q,i,k)*testf(l);
2240  }
2241  if(!ALE_is_disabled)
2242  {
2243  sum+=scaled_re_st*pos_time_weight*
2244  psif(q)*interpolated_dudx(i,p)*testf(l);
2245  }
2246 
2247  // Multiply through by Jacobian and integration weight
2248  dresidual_dnodal_coordinates(local_eqn,p,q)+=sum*J*w;
2249 
2250  // Derivs w.r.t. to nodal velocities
2251  // ---------------------------------
2252  if(element_has_node_with_aux_node_update_fct)
2253  {
2254  sum=-visc_ratio*Gamma[i]*dpsifdx(q,i)*dtestfdx(l,p)
2255  -scaled_re*psif(q)*interpolated_dudx(i,p)*testf(l);
2256  if (i==p)
2257  {
2258  sum-=scaled_re_st*val_time_weight*psif(q)*testf(l);
2259  for (unsigned k=0;k<DIM;k++)
2260  {
2261  sum-=visc_ratio*dpsifdx(q,k)*dtestfdx(l,k);
2262  double tmp=scaled_re*interpolated_u[k];
2263  if (!ALE_is_disabled) tmp-=scaled_re_st*mesh_velocity[k];
2264  sum-=tmp*dpsifdx(q,k)*testf(l);
2265  }
2266  }
2267  dresidual_dnodal_coordinates(local_eqn,p,q)+=
2268  sum*d_U_dX(p,q)*J*w;
2269  }
2270  }
2271  }
2272  }
2273  }
2274  } // End of loop over test functions
2275 
2276 
2277  // CONTINUITY EQUATION
2278  // -------------------
2279 
2280  //Loop over the shape functions
2281  for(unsigned l=0;l<n_pres;l++)
2282  {
2283  local_eqn = p_local_eqn(l);
2284 
2285  //If not a boundary conditions
2286  if(local_eqn >= 0)
2287  {
2288 
2289  // Loop over coordinate directions
2290  for (unsigned p=0;p<DIM;p++)
2291  {
2292  // Loop over nodes
2293  for (unsigned q=0;q<n_node;q++)
2294  {
2295 
2296  // Residual x deriv of Jacobian
2297  //-----------------------------
2298 
2299  // Source term
2300  double aux=-source;
2301 
2302  // Loop over velocity components
2303  for(unsigned k=0;k<DIM;k++)
2304  {
2305  aux += interpolated_dudx(k,k);
2306  }
2307 
2308  // Multiply through by deriv of Jacobian and integration weight
2309  dresidual_dnodal_coordinates(local_eqn,p,q)+=
2310  aux*dJ_dX(p,q)*testp[l]*w;
2311 
2312  // Derivative of residual x Jacobian
2313  //----------------------------------
2314 
2315  // Loop over velocity components
2316  aux=-source_gradient[p]*psif(q);
2317  for(unsigned k=0;k<DIM;k++)
2318  {
2319  aux += d_dudx_dX(p,q,k,k);
2320  }
2321  // Multiply through by Jacobian and integration weight
2322  dresidual_dnodal_coordinates(local_eqn,p,q)+=
2323  aux*testp[l]*J*w;
2324 
2325  // Derivs w.r.t. to nodal velocities
2326  //---------------------------------
2327  if(element_has_node_with_aux_node_update_fct)
2328  {
2329  aux=dpsifdx(q,p)*testp(l);
2330  dresidual_dnodal_coordinates(local_eqn,p,q)+=
2331  aux*d_U_dX(p,q)*J*w;
2332  }
2333  }
2334  }
2335  }
2336  }
2337  } // End of loop over integration points
2338 
2339 }
2340 
2341 
2342 
2343 
2344 
2345 
2346 
2347 
2348 //////////////////////////////////////////////////////////////////////
2349 //////////////////////////////////////////////////////////////////////
2350 //////////////////////////////////////////////////////////////////////
2351 
2352 ///2D Crouzeix-Raviart elements
2353 //Set the data for the number of Variables at each node
2354 template<>
2356 ={2,2,2,2,2,2,2,2,2};
2357 
2358 ///3D Crouzeix-Raviart elements
2359 //Set the data for the number of Variables at each node
2360 template<>
2362 ={3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3};
2363 
2364 
2365 //========================================================================
2366 /// Number of values (pinned or dofs) required at node n.
2367 //========================================================================
2368 template<unsigned DIM>
2369 unsigned QCrouzeixRaviartElement<DIM>::required_nvalue(const unsigned& n) const
2370  {return Initial_Nvalue[n];}
2371 
2372 
2373 //=========================================================================
2374 /// Add to the set \c paired_load_data pairs containing
2375 /// - the pointer to a Data object
2376 /// and
2377 /// - the index of the value in that Data object
2378 /// .
2379 /// for all values (pressures, velocities) that affect the
2380 /// load computed in the \c get_load(...) function.
2381 //=========================================================================
2382 template<unsigned DIM>
2384 identify_load_data(std::set<std::pair<Data*,unsigned> > &paired_load_data)
2385 {
2386  //Find the index at which the velocity is stored
2387  unsigned u_index[DIM];
2388  for(unsigned i=0;i<DIM;i++) {u_index[i] = this->u_index_nst(i);}
2389 
2390  //Loop over the nodes
2391  unsigned n_node = this->nnode();
2392  for(unsigned n=0;n<n_node;n++)
2393  {
2394  //Loop over the velocity components and add pointer to their data
2395  //and indices to the vectors
2396  for(unsigned i=0;i<DIM;i++)
2397  {
2398  paired_load_data.insert(std::make_pair(this->node_pt(n),u_index[i]));
2399  }
2400  }
2401 
2402  // Identify the pressure data
2403  identify_pressure_data(paired_load_data);
2404 
2405 }
2406 
2407 
2408 //=========================================================================
2409 /// Add to the set \c paired_pressue_data pairs containing
2410 /// - the pointer to a Data object
2411 /// and
2412 /// - the index of the value in that Data object
2413 /// .
2414 /// for all pressures values that affect the
2415 /// load computed in the \c get_load(...) function.
2416 //=========================================================================
2417 template<unsigned DIM>
2419 identify_pressure_data(std::set<std::pair<Data*,unsigned> > &paired_pressure_data)
2420 {
2421  //Loop over the internal data
2422  unsigned n_internal = this->ninternal_data();
2423  for(unsigned l=0;l<n_internal;l++)
2424  {
2425  unsigned nval=this->internal_data_pt(l)->nvalue();
2426  //Add internal data
2427  for (unsigned j=0;j<nval;j++)
2428  {
2429  paired_pressure_data.insert(std::make_pair(this->internal_data_pt(l),j));
2430  }
2431  }
2432 }
2433 
2434 
2435 //=============================================================================
2436 /// Create a list of pairs for all unknowns in this element,
2437 /// so that the first entry in each pair contains the global equation
2438 /// number of the unknown, while the second one contains the number
2439 /// of the DOF that this unknown is associated with.
2440 /// (Function can obviously only be called if the equation numbering
2441 /// scheme has been set up.)
2442 //=============================================================================
2443 template<unsigned DIM>
2445  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
2446 {
2447  // number of nodes
2448  unsigned n_node = this->nnode();
2449 
2450  // number of pressure values
2451  unsigned n_press = this->npres_nst();
2452 
2453  // temporary pair (used to store dof lookup prior to being added to list)
2454  std::pair<unsigned,unsigned> dof_lookup;
2455 
2456  // pressure dof number
2457  unsigned pressure_dof_number = DIM;
2458 
2459  // loop over the pressure values
2460  for (unsigned n = 0; n < n_press; n++)
2461  {
2462  // determine local eqn number
2463  int local_eqn_number = this->p_local_eqn(n);
2464 
2465  // ignore pinned values - far away degrees of freedom resulting
2466  // from hanging nodes can be ignored since these are be dealt
2467  // with by the element containing their master nodes
2468  if (local_eqn_number >= 0)
2469  {
2470  // store dof lookup in temporary pair: First entry in pair
2471  // is global equation number; second entry is dof type
2472  dof_lookup.first = this->eqn_number(local_eqn_number);
2473  dof_lookup.second = pressure_dof_number;
2474 
2475  // add to list
2476  dof_lookup_list.push_front(dof_lookup);
2477  }
2478  }
2479 
2480  // loop over the nodes
2481  for (unsigned n = 0; n < n_node; n++)
2482  {
2483  // find the number of values at this node
2484  unsigned nv = this->node_pt(n)->nvalue();
2485 
2486  //loop over these values
2487  for (unsigned v = 0; v < nv; v++)
2488  {
2489  // determine local eqn number
2490  int local_eqn_number = this->nodal_local_eqn(n, v);
2491 
2492  // ignore pinned values
2493  if (local_eqn_number >= 0)
2494  {
2495  // store dof lookup in temporary pair: First entry in pair
2496  // is global equation number; second entry is dof type
2497  dof_lookup.first = this->eqn_number(local_eqn_number);
2498  dof_lookup.second = v;
2499 
2500  // add to list
2501  dof_lookup_list.push_front(dof_lookup);
2502 
2503  }
2504  }
2505  }
2506 }
2507 
2508 
2509 
2510 
2511 ///////////////////////////////////////////////////////////////////////////
2512 ///////////////////////////////////////////////////////////////////////////
2513 ///////////////////////////////////////////////////////////////////////////
2514 
2515 
2516 
2517 //2D Taylor--Hood
2518 //Set the data for the number of Variables at each node
2519 template<>
2520 const unsigned QTaylorHoodElement<2>::Initial_Nvalue[9]={3,2,3,2,2,2,3,2,3};
2521 
2522 //Set the data for the pressure conversion array
2523 template<>
2524 const unsigned QTaylorHoodElement<2>::Pconv[4]={0,2,6,8};
2525 
2526 //3D Taylor--Hood
2527 //Set the data for the number of Variables at each node
2528 template<>
2529 const unsigned QTaylorHoodElement<3>::Initial_Nvalue[27]={4,3,4,3,3,3,4,3,4,3,3,3,3,3,3,3,3,3,4,3,4,3,3,3,4,3,4};
2530 
2531 //Set the data for the pressure conversion array
2532 template<>
2533 const unsigned QTaylorHoodElement<3>::Pconv[8]={0,2,6,8,18,20,24,26};
2534 
2535 //=========================================================================
2536 /// Add to the set \c paired_load_data pairs containing
2537 /// - the pointer to a Data object
2538 /// and
2539 /// - the index of the value in that Data object
2540 /// .
2541 /// for all values (pressures, velocities) that affect the
2542 /// load computed in the \c get_load(...) function.
2543 //=========================================================================
2544 template<unsigned DIM>
2546 identify_load_data(std::set<std::pair<Data*,unsigned> > &paired_load_data)
2547 {
2548  //Find the index at which the velocity is stored
2549  unsigned u_index[DIM];
2550  for(unsigned i=0;i<DIM;i++) {u_index[i] = this->u_index_nst(i);}
2551 
2552  //Loop over the nodes
2553  unsigned n_node = this->nnode();
2554  for(unsigned n=0;n<n_node;n++)
2555  {
2556  //Loop over the velocity components and add pointer to their data
2557  //and indices to the vectors
2558  for(unsigned i=0;i<DIM;i++)
2559  {
2560  paired_load_data.insert(std::make_pair(this->node_pt(n),u_index[i]));
2561  }
2562  }
2563 
2564  //Identify the pressure data
2565  this->identify_pressure_data(paired_load_data);
2566 
2567 }
2568 
2569 //=========================================================================
2570 /// Add to the set \c paired_pressure_data pairs containing
2571 /// - the pointer to a Data object
2572 /// and
2573 /// - the index of the value in that Data object
2574 /// .
2575 /// for pressure values that affect the
2576 /// load computed in the \c get_load(...) function.,
2577 //=========================================================================
2578 template<unsigned DIM>
2580 identify_pressure_data(std::set<std::pair<Data*,unsigned> > &paired_pressure_data)
2581 {
2582  //Find the index at which the pressure is stored
2583  unsigned p_index = static_cast<unsigned>(this->p_nodal_index_nst());
2584 
2585  //Loop over the pressure data
2586  unsigned n_pres= npres_nst();
2587  for(unsigned l=0;l<n_pres;l++)
2588  {
2589  //The DIMth entry in each nodal data is the pressure, which
2590  //affects the traction
2591  paired_pressure_data.insert(std::make_pair(this->node_pt(Pconv[l]),p_index));
2592  }
2593 }
2594 
2595 
2596 //============================================================================
2597 /// Create a list of pairs for all unknowns in this element,
2598 /// so the first entry in each pair contains the global equation
2599 /// number of the unknown, while the second one contains the number
2600 /// of the "DOF type" that this unknown is associated with.
2601 /// (Function can obviously only be called if the equation numbering
2602 /// scheme has been set up.)
2603 //============================================================================
2604 template<unsigned DIM>
2606  std::list<std::pair<unsigned long,
2607  unsigned> >& dof_lookup_list) const
2608 {
2609  // number of nodes
2610  unsigned n_node = this->nnode();
2611 
2612  // local eqn no for pressure unknown
2613 // unsigned p_index = this->p_nodal_index_nst();
2614 
2615  // temporary pair (used to store dof lookup prior to being added to list)
2616  std::pair<unsigned,unsigned> dof_lookup;
2617 
2618  // loop over the nodes
2619  for (unsigned n = 0; n < n_node; n++)
2620  {
2621  // find the number of values at this node
2622  unsigned nv = this->required_nvalue(n);
2623 
2624  //loop over these values
2625  for (unsigned v = 0; v < nv; v++)
2626  {
2627  // determine local eqn number
2628  int local_eqn_number = this->nodal_local_eqn(n, v);
2629 
2630  // ignore pinned values - far away degrees of freedom resulting
2631  // from hanging nodes can be ignored since these are be dealt
2632  // with by the element containing their master nodes
2633  if (local_eqn_number >= 0)
2634  {
2635  // store dof lookup in temporary pair: Global equation number
2636  // is the first entry in pair
2637  dof_lookup.first = this->eqn_number(local_eqn_number);
2638 
2639  // set dof numbers: Dof number is the second entry in pair
2640  dof_lookup.second = v;
2641 
2642  // add to list
2643  dof_lookup_list.push_front(dof_lookup);
2644  }
2645  }
2646  }
2647  }
2648 
2649 
2650 //====================================================================
2651 //// Force build of templates
2652 //====================================================================
2653 template class NavierStokesEquations<2>;
2654 template class QCrouzeixRaviartElement<2>;
2655 template class QTaylorHoodElement<2>;
2656 
2657 template class NavierStokesEquations<3>;
2658 template class QCrouzeixRaviartElement<3>;
2659 template class QTaylorHoodElement<3>;
2660 
2661 }
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
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_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. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element.
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...
void output_vorticity(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z], [omega_x,omega_y,[and/or omega_z]] in tecplot format. nplot points in each ...
cstr elem_len * i
Definition: cfortran.h:607
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
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:1509
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
char t
Definition: cfortran.h:572
double d_kin_energy_dt() const
Get integral of time derivative of kinetic energy over element.
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:322
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
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. Function prints as many components as are returned in solution Vector.
e
Definition: cfortran.h:575
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...
A Rank 4 Tensor class.
Definition: matrices.h:1625
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:504
bool has_auxiliary_node_update_fct_pt()
Boolean to indicate if node has a pointer to and auxiliary update function.
Definition: nodes.h:1500
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, unsigned flag)
Compute the derivatives of the residuals for the Navier–Stokes equations with respect to a parameter ...
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
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.
void output_veloc(std::ostream &outfile, const unsigned &nplot, const unsigned &t)
Output function: x,y,[z],u,v,[w] in tecplot format. nplot points in each coordinate direction at time...
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 get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. Overwrites default implementation in FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}.
A Rank 3 Tensor class.
Definition: matrices.h:1337
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void get_vorticity(const Vector< double > &s, Vector< double > &vorticity) const
Compute the vorticity vector at local coordinate s.
virtual void fill_in_generic_pressure_advection_diffusion_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
static char t char * s
Definition: cfortran.h:572
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:507
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
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 shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
Definition: shape.h:549
double pressure_integral() const
Integral of pressure over element.
double Default_Physical_Constant_Value
Default value for physical constants.
double kin_energy() const
Get integral of kinetic energy over element.
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
double dissipation() const
Return integral of dissipation over element.
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
virtual void fill_in_generic_residual_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations. Flag=1 (or 0): do (or don't) compute the Jacob...
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164