generalised_newtonian_axisym_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: 1194 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-27 08:44:43 +0100 (Fri, 27 May 2016) $
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
32 
33 
34 namespace oomph
35 {
36 
37 
38 /// Navier--Stokes equations static data
40 Gamma(2,1.0);
41 
42 //=================================================================
43 /// "Magic" negative number that indicates that the pressure is
44 /// not stored at a node. This cannot be -1 because that represents
45 /// the positional hanging scheme in the hanging_pt object of nodes
46 //=================================================================
49 
50 /// Navier--Stokes equations static data
53 
54  // Navier--Stokes equations static data
57 
58  /// Navier-Stokes equations default gravity vector
61 
62  /// By default disable the pre-multiplication of the pressure and continuity
63  /// equation by the viscosity ratio
66 
67 
68 //================================================================
69 /// Compute the diagonal of the velocity/pressure mass matrices.
70 /// If which one=0, both are computed, otherwise only the pressure
71 /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
72 /// LSC version of the preconditioner only needs that one)
73 /// NOTE: pressure versions isn't implemented yet because this
74 /// element has never been tried with Fp preconditoner.
75 //================================================================
78  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
79  const unsigned& which_one)
80  {
81 
82 #ifdef PARANOID
83  if ((which_one==0)||(which_one==1))
84  {
85  throw OomphLibError(
86  "Computation of diagonal of pressure mass matrix is not impmented yet.\n",
87  OOMPH_CURRENT_FUNCTION,
88  OOMPH_EXCEPTION_LOCATION);
89  }
90 #endif
91 
92  // Resize and initialise
93  veloc_mass_diag.assign(ndof(), 0.0);
94 
95  // find out how many nodes there are
96  const unsigned n_node = nnode();
97 
98  // find number of coordinates
99  const unsigned n_value = 3;
100 
101  // find the indices at which the local velocities are stored
102  Vector<unsigned> u_nodal_index(n_value);
103  for(unsigned i=0; i<n_value; i++)
104  {
105  u_nodal_index[i] = this->u_index_axi_nst(i);
106  }
107 
108  //Set up memory for test functions
109  Shape test(n_node);
110 
111  //Number of integration points
112  unsigned n_intpt = integral_pt()->nweight();
113 
114  //Integer to store the local equations no
115  int local_eqn=0;
116 
117  //Loop over the integration points
118  for(unsigned ipt=0; ipt<n_intpt; ipt++)
119  {
120 
121  //Get the integral weight
122  double w = integral_pt()->weight(ipt);
123 
124  //Get determinant of Jacobian of the mapping
125  double J = J_eulerian_at_knot(ipt);
126 
127  //Premultiply weights and Jacobian
128  double W = w*J;
129 
130  //Get the velocity test functions - there is no explicit
131  // function to give the test function so use shape
132  shape_at_knot(ipt,test);
133 
134  //Need to get the position to sort out the jacobian properly
135  double r = 0.0;
136  for(unsigned l=0;l<n_node;l++)
137  {
138  r += this->nodal_position(l,0)*test(l);
139  }
140 
141  //Multiply by the geometric factor
142  W *= r;
143 
144  //Loop over the veclocity test functions
145  for(unsigned l=0; l<n_node; l++)
146  {
147  //Loop over the velocity components
148  for(unsigned i=0; i<n_value; i++)
149  {
150  local_eqn = nodal_local_eqn(l,u_nodal_index[i]);
151 
152  //If not a boundary condition
153  if(local_eqn >= 0)
154  {
155  //Add the contribution
156  veloc_mass_diag[local_eqn] += test[l]*test[l] * W;
157  } //End of if not boundary condition statement
158  } //End of loop over dimension
159  } //End of loop over test functions
160 
161  }
162  }
163 
164 
165 //======================================================================
166 /// Validate against exact velocity solution at given time.
167 /// Solution is provided via function pointer.
168 /// Plot at a given number of plot points and compute L2 error
169 /// and L2 norm of velocity solution over element.
170 //=======================================================================
172 compute_error(std::ostream &outfile,
174  const double& time,
175  double& error, double& norm)
176 {
177  error=0.0;
178  norm=0.0;
179 
180  //Vector of local coordinates
181  Vector<double> s(2);
182 
183  // Vector for coordintes
184  Vector<double> x(2);
185 
186  //Set the value of Nintpt
187  unsigned Nintpt = integral_pt()->nweight();
188 
189  outfile << "ZONE" << std::endl;
190 
191  // Exact solution Vector (u,v,w,p)
192  Vector<double> exact_soln(4);
193 
194  //Loop over the integration points
195  for(unsigned ipt=0;ipt<Nintpt;ipt++)
196  {
197  //Assign values of s
198  for(unsigned i=0;i<2;i++) {s[i] = integral_pt()->knot(ipt,i);}
199 
200  //Get the integral weight
201  double w = integral_pt()->weight(ipt);
202 
203  // Get jacobian of mapping
204  double J = J_eulerian(s);
205 
206  // Get x position as Vector
207  interpolated_x(s,x);
208 
209  //Premultiply the weights and the Jacobian and r
210  double W = w*J*x[0];
211 
212  // Get exact solution at this point
213  (*exact_soln_pt)(time,x,exact_soln);
214 
215  // Velocity error
216  for(unsigned i=0;i<3;i++)
217  {
218  norm+=exact_soln[i]*exact_soln[i]*W;
219  error+=(exact_soln[i]-interpolated_u_axi_nst(s,i))*
220  (exact_soln[i]-interpolated_u_axi_nst(s,i))*W;
221  }
222 
223  //Output x,y,...,u_exact
224  for(unsigned i=0;i<2;i++) {outfile << x[i] << " ";}
225 
226  //Output x,y,[z],u_error,v_error,[w_error]
227  for(unsigned i=0;i<3;i++)
228  {outfile << exact_soln[i]-interpolated_u_axi_nst(s,i) << " ";}
229 
230  outfile << std::endl;
231  }
232 }
233 
234 //======================================================================
235 /// Validate against exact velocity solution
236 /// Solution is provided via function pointer.
237 /// Plot at a given number of plot points and compute L2 error
238 /// and L2 norm of velocity solution over element.
239 //=======================================================================
241 compute_error(std::ostream &outfile,
243  double& error, double& norm)
244 {
245  error=0.0;
246  norm=0.0;
247 
248  //Vector of local coordinates
249  Vector<double> s(2);
250 
251  // Vector for coordintes
252  Vector<double> x(2);
253 
254  //Set the value of Nintpt
255  unsigned Nintpt = integral_pt()->nweight();
256 
257  outfile << "ZONE" << std::endl;
258 
259  // Exact solution Vector (u,v,w,p)
260  Vector<double> exact_soln(4);
261 
262  //Loop over the integration points
263  for(unsigned ipt=0;ipt<Nintpt;ipt++)
264  {
265 
266  //Assign values of s
267  for(unsigned i=0;i<2;i++) {s[i] = integral_pt()->knot(ipt,i);}
268 
269  //Get the integral weight
270  double w = integral_pt()->weight(ipt);
271 
272  // Get jacobian of mapping
273  double J=J_eulerian(s);
274 
275  // Get x position as Vector
276  interpolated_x(s,x);
277 
278  //Premultiply the weights and the Jacobian and r
279  double W = w*J*x[0];
280 
281  // Get exact solution at this point
282  (*exact_soln_pt)(x,exact_soln);
283 
284  // Velocity error
285  for(unsigned i=0;i<3;i++)
286  {
287  norm+=exact_soln[i]*exact_soln[i]*W;
288  error+=(exact_soln[i]-interpolated_u_axi_nst(s,i))*
289  (exact_soln[i]-interpolated_u_axi_nst(s,i))*W;
290  }
291 
292  //Output x,y,...,u_exact
293  for(unsigned i=0;i<2;i++) {outfile << x[i] << " ";}
294 
295  //Output x,y,u_error,v_error,w_error
296  for(unsigned i=0;i<3;i++)
297  {outfile << exact_soln[i]-interpolated_u_axi_nst(s,i) << " ";}
298 
299  outfile << std::endl;
300  }
301 }
302 
303 //======================================================================
304 /// Output "exact" solution
305 /// Solution is provided via function pointer.
306 /// Plot at a given number of plot points.
307 /// Function prints as many components as are returned in solution Vector.
308 //=======================================================================
310 output_fct(std::ostream &outfile,
311  const unsigned &nplot,
313 {
314 
315  //Vector of local coordinates
316  Vector<double> s(2);
317 
318  // Vector for coordintes
319  Vector<double> x(2);
320 
321  // Tecplot header info
322  outfile << tecplot_zone_string(nplot);
323 
324  // Exact solution Vector
325  Vector<double> exact_soln;
326 
327  // Loop over plot points
328  unsigned num_plot_points=nplot_points(nplot);
329  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
330  {
331 
332  // Get local coordinates of plot point
333  get_s_plot(iplot,nplot,s);
334 
335  // Get x position as Vector
336  interpolated_x(s,x);
337 
338  // Get exact solution at this point
339  (*exact_soln_pt)(x,exact_soln);
340 
341  //Output x,y,...
342  for(unsigned i=0;i<2;i++)
343  {
344  outfile << x[i] << " ";
345  }
346 
347  //Output "exact solution"
348  for(unsigned i=0;i<exact_soln.size();i++)
349  {
350  outfile << exact_soln[i] << " ";
351  }
352 
353  outfile << std::endl;
354 
355  }
356 
357  // Write tecplot footer (e.g. FE connectivity lists)
358  write_tecplot_zone_footer(outfile,nplot);
359 
360 }
361 
362 //======================================================================
363 /// Output "exact" solution at a given time
364 /// Solution is provided via function pointer.
365 /// Plot at a given number of plot points.
366 /// Function prints as many components as are returned in solution Vector.
367 //=======================================================================
369 output_fct(std::ostream &outfile,
370  const unsigned &nplot,
371  const double& time,
373 {
374 
375  //Vector of local coordinates
376  Vector<double> s(2);
377 
378  // Vector for coordintes
379  Vector<double> x(2);
380 
381  // Tecplot header info
382  outfile << tecplot_zone_string(nplot);
383 
384  // Exact solution Vector
385  Vector<double> exact_soln;
386 
387  // Loop over plot points
388  unsigned num_plot_points=nplot_points(nplot);
389  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
390  {
391 
392  // Get local coordinates of plot point
393  get_s_plot(iplot,nplot,s);
394 
395  // Get x position as Vector
396  interpolated_x(s,x);
397 
398  // Get exact solution at this point
399  (*exact_soln_pt)(time,x,exact_soln);
400 
401  //Output x,y,...
402  for(unsigned i=0;i<2;i++)
403  {
404  outfile << x[i] << " ";
405  }
406 
407  //Output "exact solution"
408  for(unsigned i=0;i<exact_soln.size();i++)
409  {
410  outfile << exact_soln[i] << " ";
411  }
412 
413  outfile << std::endl;
414 
415  }
416 
417  // Write tecplot footer (e.g. FE connectivity lists)
418  write_tecplot_zone_footer(outfile,nplot);
419 
420 }
421 
422 //==============================================================
423 /// Output function: Velocities only
424 /// x,y,[z],u,v,[w]
425 /// in tecplot format at specified previous timestep (t=0: present;
426 /// t>0: previous timestep). Specified number of plot points in each
427 /// coordinate direction.
428 //==============================================================
430 output_veloc(std::ostream &outfile,
431  const unsigned &nplot,
432  const unsigned &t)
433 {
434  //Find number of nodes
435  unsigned n_node = nnode();
436 
437  //Local shape function
438  Shape psi(n_node);
439 
440  //Vectors of local coordinates and coords and velocities
441  Vector<double> s(2);
443  Vector<double> interpolated_u(3);
444 
445 
446  // Tecplot header info
447  outfile << tecplot_zone_string(nplot);
448 
449  // Loop over plot points
450  unsigned num_plot_points=nplot_points(nplot);
451  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
452  {
453 
454  // Get local coordinates of plot point
455  get_s_plot(iplot,nplot,s);
456 
457  // Get shape functions
458  shape(s,psi);
459 
460  // Loop over coordinate directions
461  for(unsigned i=0;i<2;i++)
462  {
463  interpolated_x[i]=0.0;
464  //Loop over the local nodes and sum
465  for(unsigned l=0;l<n_node;l++)
466  {interpolated_x[i] += nodal_position(t,l,i)*psi[l];}
467  }
468 
469  //Loop over the velocity components
470  for(unsigned i=0;i<3;i++)
471  {
472  //Get the index at which the velocity is stored
473  unsigned u_nodal_index = u_index_axi_nst(i);
474  interpolated_u[i]=0.0;
475  //Loop over the local nodes and sum
476  for(unsigned l=0;l<n_node;l++)
477  {interpolated_u[i] += nodal_value(t,l,u_nodal_index)*psi[l];}
478  }
479 
480  // Coordinates
481  for(unsigned i=0;i<2;i++) {outfile << interpolated_x[i] << " ";}
482 
483  // Velocities
484  for(unsigned i=0;i<3;i++) {outfile << interpolated_u[i] << " ";}
485 
486  outfile << std::endl;
487  }
488 
489  // Write tecplot footer (e.g. FE connectivity lists)
490  write_tecplot_zone_footer(outfile,nplot);
491 
492 }
493 
494 //==============================================================
495 /// Output function:
496 /// r,z,u,v,w,p
497 /// in tecplot format. Specified number of plot points in each
498 /// coordinate direction.
499 //==============================================================
501 output(std::ostream &outfile, const unsigned &nplot)
502 {
503  //Vector of local coordinates
504  Vector<double> s(2);
505 
506  // Tecplot header info
507  outfile << tecplot_zone_string(nplot);
508 
509  // Loop over plot points
510  unsigned num_plot_points=nplot_points(nplot);
511  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
512  {
513 
514  // Get local coordinates of plot point
515  get_s_plot(iplot,nplot,s);
516 
517  // Coordinates
518  for(unsigned i=0;i<2;i++) {outfile << interpolated_x(s,i) << " ";}
519 
520  // Velocities
521  for(unsigned i=0;i<3;i++) {outfile << interpolated_u_axi_nst(s,i) << " ";}
522 
523  // Pressure
524  outfile << interpolated_p_axi_nst(s) << " ";
525 
526  outfile << std::endl;
527  }
528  outfile << std::endl;
529 
530  // Write tecplot footer (e.g. FE connectivity lists)
531  write_tecplot_zone_footer(outfile,nplot);
532 
533 }
534 
535 
536 //==============================================================
537 /// Output function:
538 /// r,z,u,v,w,p
539 /// in tecplot format. Specified number of plot points in each
540 /// coordinate direction.
541 //==============================================================
543 output(FILE* file_pt, const unsigned &nplot)
544 {
545  //Vector of local coordinates
546  Vector<double> s(2);
547 
548  // Tecplot header info
549  fprintf(file_pt,"%s ",tecplot_zone_string(nplot).c_str());
550 
551  // Loop over plot points
552  unsigned num_plot_points=nplot_points(nplot);
553  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
554  {
555 
556  // Get local coordinates of plot point
557  get_s_plot(iplot,nplot,s);
558 
559  // Coordinates
560  for(unsigned i=0;i<2;i++)
561  {
562  //outfile << interpolated_x(s,i) << " ";
563  fprintf(file_pt,"%g ",interpolated_x(s,i));
564  }
565 
566  // Velocities
567  for(unsigned i=0;i<3;i++)
568  {
569  //outfile << interpolated_u(s,i) << " ";
570  fprintf(file_pt,"%g ",interpolated_u_axi_nst(s,i));
571  }
572 
573  // Pressure
574  //outfile << interpolated_p(s) << " ";
575  fprintf(file_pt,"%g ",interpolated_p_axi_nst(s));
576 
577  //outfile << std::endl;
578  fprintf(file_pt,"\n");
579  }
580  //outfile << std::endl;
581  fprintf(file_pt,"\n");
582 
583  // Write tecplot footer (e.g. FE connectivity lists)
584  write_tecplot_zone_footer(file_pt,nplot);
585 
586 }
587 
588 
589 
590 
591 //==============================================================
592 /// Return integral of dissipation over element
593 //==============================================================
595 dissipation() const
596 {
597  throw OomphLibError(
598  "Check the dissipation calculation for axisymmetric NSt",
599  OOMPH_CURRENT_FUNCTION,
600  OOMPH_EXCEPTION_LOCATION);
601 
602  // Initialise
603  double diss=0.0;
604 
605  //Set the value of Nintpt
606  unsigned Nintpt = integral_pt()->nweight();
607 
608  //Set the Vector to hold local coordinates
609  Vector<double> s(2);
610 
611  //Loop over the integration points
612  for(unsigned ipt=0;ipt<Nintpt;ipt++)
613  {
614 
615  //Assign values of s
616  for(unsigned i=0;i<2;i++)
617  {
618  s[i] = integral_pt()->knot(ipt,i);
619  }
620 
621  //Get the integral weight
622  double w = integral_pt()->weight(ipt);
623 
624  // Get Jacobian of mapping
625  double J = J_eulerian(s);
626 
627  // Get strain rate matrix
628  DenseMatrix<double> strainrate(3,3);
629  strain_rate(s,strainrate);
630 
631  // Initialise
632  double local_diss=0.0;
633  for(unsigned i=0;i<3;i++)
634  {
635  for(unsigned j=0;j<3;j++)
636  {
637  local_diss+=2.0*strainrate(i,j)*strainrate(i,j);
638  }
639  }
640 
641  diss+=local_diss*w*J;
642  }
643 
644  return diss;
645 
646 }
647 
648 //==============================================================
649 /// \short Compute traction (on the viscous scale) at local
650 /// coordinate s for outer unit normal N
651 //==============================================================
654  Vector<double>& traction)
655  const
656 {
657  //throw OomphLibError(
658  // "Check the traction calculation for axisymmetric NSt",
659  // OOMPH_CURRENT_FUNCTION,
660  // OOMPH_EXCEPTION_LOCATION);
661 
662  // Pad out normal vector if required
663  Vector<double> n_local(3,0.0);
664  n_local[0]=N[0];
665  n_local[1]=N[1];
666 
667 #ifdef PARANOID
668  if ((N.size()==3)&&(N[2]!=0.0))
669  {
670  throw OomphLibError(
671  "Unit normal passed into this fct should either be 2D (r,z) or have a zero component in the theta-direction",
672  OOMPH_CURRENT_FUNCTION,
673  OOMPH_EXCEPTION_LOCATION);
674  }
675 #endif
676 
677  // Get velocity gradients
678  DenseMatrix<double> strainrate(3,3,0.0);
679 
680  // Do we use the current or extrapolated strainrate to compute
681  // the second invariant?
683  {
684  strain_rate(s,strainrate);
685  }
686  else
687  {
688  extrapolated_strain_rate(s,strainrate);
689  }
690 
691  // Get the second invariant of the rate of strain tensor
693 
694  double visc=Constitutive_eqn_pt->viscosity(second_invariant);
695 
696  // Get pressure
697  double press=interpolated_p_axi_nst(s);
698 
699  // Loop over traction components
700  for (unsigned i=0;i<3;i++)
701  {
702  traction[i]=-press*n_local[i];
703  for (unsigned j=0;j<3;j++)
704  {
705  traction[i]+=visc*2.0*strainrate(i,j)*n_local[j];
706  }
707  }
708 }
709 
710 //==============================================================
711 /// Return dissipation at local coordinate s
712 //==============================================================
715 {
716  throw OomphLibError(
717  "Check the dissipation calculation for axisymmetric NSt",
718 OOMPH_CURRENT_FUNCTION,
719  OOMPH_EXCEPTION_LOCATION);
720 
721  // Get strain rate matrix
722  DenseMatrix<double> strainrate(3,3);
723  strain_rate(s,strainrate);
724 
725  // Initialise
726  double local_diss=0.0;
727  for(unsigned i=0;i<3;i++)
728  {
729  for(unsigned j=0;j<3;j++)
730  {
731  local_diss+=2.0*strainrate(i,j)*strainrate(i,j);
732  }
733  }
734 
735  return local_diss;
736 }
737 
738 //==============================================================
739 /// Get strain-rate tensor: \f$ e_{ij} \f$ where
740 /// \f$ i,j = r,z,\theta \f$ (in that order)
741 //==============================================================
744 {
745 
746 #ifdef PARANOID
747  if ((strainrate.ncol()!=3)||(strainrate.nrow()!=3))
748  {
749  std::ostringstream error_message;
750  error_message << "The strain rate has incorrect dimensions "
751  << strainrate.ncol() << " x "
752  << strainrate.nrow() << " Not 3" << std::endl;
753 
754  throw OomphLibError(error_message.str(),
755  OOMPH_CURRENT_FUNCTION,
756  OOMPH_EXCEPTION_LOCATION);
757  }
758 #endif
759 
760  //Find out how many nodes there are in the element
761  unsigned n_node = nnode();
762 
763  //Set up memory for the shape and test functions
764  Shape psi(n_node);
765  DShape dpsidx(n_node,2);
766 
767  //Call the derivatives of the shape functions
768  dshape_eulerian(s,psi,dpsidx);
769 
770  // Radius
771  double interpolated_r = 0.0;
772 
773  // Velocity components and their derivatives
774  double ur=0.0;
775  double durdr=0.0;
776  double durdz=0.0;
777  double uz=0.0;
778  double duzdr=0.0;
779  double duzdz=0.0;
780  double uphi=0.0;
781  double duphidr=0.0;
782  double duphidz=0.0;
783 
784  //Get the local storage for the indices
785  unsigned u_nodal_index[3];
786  for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}
787 
788  // Loop over nodes to assemble velocities and their derivatives
789  // w.r.t. to r and z (x_0 and x_1)
790  for(unsigned l=0;l<n_node;l++)
791  {
792  interpolated_r += nodal_position(l,0)*psi[l];
793 
794  ur += nodal_value(l,u_nodal_index[0])*psi[l];
795  uz += nodal_value(l,u_nodal_index[1])*psi[l];
796  uphi += nodal_value(l,u_nodal_index[2])*psi[l];
797 
798  durdr += nodal_value(l,u_nodal_index[0])*dpsidx(l,0);
799  durdz += nodal_value(l,u_nodal_index[0])*dpsidx(l,1);
800 
801  duzdr += nodal_value(l,u_nodal_index[1])*dpsidx(l,0);
802  duzdz += nodal_value(l,u_nodal_index[1])*dpsidx(l,1);
803 
804  duphidr += nodal_value(l,u_nodal_index[2])*dpsidx(l,0);
805  duphidz += nodal_value(l,u_nodal_index[2])*dpsidx(l,1);
806  }
807 
808 
809  // Assign strain rates without negative powers of the radius
810  // and zero those with:
811  strainrate(0,0)=durdr;
812  strainrate(0,1)=0.5*(durdz+duzdr);
813  strainrate(1,0)=strainrate(0,1);
814  strainrate(0,2)=0.0;
815  strainrate(2,0)=strainrate(0,2);
816  strainrate(1,1)=duzdz;
817  strainrate(1,2)=0.5*duphidz;
818  strainrate(2,1)=strainrate(1,2);
819  strainrate(2,2)=0.0;
820 
821 
822  // Overwrite the strain rates with negative powers of the radius
823  // unless we're at the origin
824  if (std::fabs(interpolated_r)>1.0e-16)
825  {
826  double inverse_radius=1.0/interpolated_r;
827  strainrate(0,2)=0.5*(duphidr-inverse_radius*uphi);
828  strainrate(2,0)=strainrate(0,2);
829  strainrate(2,2)=inverse_radius*ur;
830  }
831  else
832  {
833  // hierher L'Hospital? -- also for (2,0)
834  strainrate(2,2)=durdr;
835  }
836 
837 }
838 
839 //==============================================================
840 /// Get strain-rate tensor: \f$ e_{ij} \f$ where
841 /// \f$ i,j = r,z,\theta \f$ (in that order) at a specific time
842 //==============================================================
844 strain_rate(const unsigned& t, const Vector<double>& s,
845  DenseMatrix<double>& strainrate) const
846 {
847 
848 #ifdef PARANOID
849  if ((strainrate.ncol()!=3)||(strainrate.nrow()!=3))
850  {
851  std::ostringstream error_message;
852  error_message << "The strain rate has incorrect dimensions "
853  << strainrate.ncol() << " x "
854  << strainrate.nrow() << " Not 3" << std::endl;
855 
856  throw OomphLibError(error_message.str(),
857  "GeneralisedNewtonianAxisymmetricNavierStokeEquations::strain_rate()",
858  OOMPH_EXCEPTION_LOCATION);
859  }
860 #endif
861 
862  //Find out how many nodes there are in the element
863  unsigned n_node = nnode();
864 
865  //Set up memory for the shape and test functions
866  Shape psi(n_node);
867  DShape dpsidx(n_node,2);
868 
869  // Loop over all nodes to back up current positions and over-write them
870  // with the appropriate history values
871  DenseMatrix<double> backed_up_nodal_position(n_node,2);
872  for (unsigned j=0;j<n_node;j++)
873  {
874  backed_up_nodal_position(j,0)=node_pt(j)->x(0);
875  node_pt(j)->x(0)=node_pt(j)->x(t,0);
876  backed_up_nodal_position(j,1)=node_pt(j)->x(1);
877  node_pt(j)->x(1)=node_pt(j)->x(t,1);
878  }
879 
880  //Call the derivatives of the shape functions
881  dshape_eulerian(s,psi,dpsidx);
882 
883  // Radius
884  double interpolated_r = 0.0;
885  double interpolated_z = 0.0;
886 
887  // Velocity components and their derivatives
888  double ur=0.0;
889  double durdr=0.0;
890  double durdz=0.0;
891  double uz=0.0;
892  double duzdr=0.0;
893  double duzdz=0.0;
894  double uphi=0.0;
895  double duphidr=0.0;
896  double duphidz=0.0;
897 
898  //Get the local storage for the indices
899  unsigned u_nodal_index[3];
900  for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}
901 
902  // Loop over nodes to assemble velocities and their derivatives
903  // w.r.t. to r and z (x_0 and x_1)
904  for(unsigned l=0;l<n_node;l++)
905  {
906  interpolated_r += nodal_position(l,0)*psi[l];
907  interpolated_z += nodal_position(l,1)*psi[l];
908 
909  ur += nodal_value(t,l,u_nodal_index[0])*psi[l];
910  uz += nodal_value(t,l,u_nodal_index[1])*psi[l];
911  uphi += nodal_value(t,l,u_nodal_index[2])*psi[l];
912 
913  durdr += nodal_value(t,l,u_nodal_index[0])*dpsidx(l,0);
914  durdz += nodal_value(t,l,u_nodal_index[0])*dpsidx(l,1);
915 
916  duzdr += nodal_value(t,l,u_nodal_index[1])*dpsidx(l,0);
917  duzdz += nodal_value(t,l,u_nodal_index[1])*dpsidx(l,1);
918 
919  duphidr += nodal_value(t,l,u_nodal_index[2])*dpsidx(l,0);
920  duphidz += nodal_value(t,l,u_nodal_index[2])*dpsidx(l,1);
921  }
922 
923 
924  // Assign strain rates without negative powers of the radius
925  // and zero those with:
926  strainrate(0,0)=durdr;
927  strainrate(0,1)=0.5*(durdz+duzdr);
928  strainrate(1,0)=strainrate(0,1);
929  strainrate(0,2)=0.0;
930  strainrate(2,0)=strainrate(0,2);
931  strainrate(1,1)=duzdz;
932  strainrate(1,2)=0.5*duphidz;
933  strainrate(2,1)=strainrate(1,2);
934  strainrate(2,2)=0.0;
935 
936 
937  // Overwrite the strain rates with negative powers of the radius
938  // unless we're at the origin
939  if (std::fabs(interpolated_r)>1.0e-16)
940  {
941  double inverse_radius=1.0/interpolated_r;
942  strainrate(0,2)=0.5*(duphidr-inverse_radius*uphi);
943  strainrate(2,0)=strainrate(0,2);
944  strainrate(2,2)=inverse_radius*ur;
945  }
946  else
947  {
948  // hierher L'Hospital? --- also for strain(0,2)
949  strainrate(2,2)=durdr;
950  }
951 
952  // Reset current nodal positions
953  for (unsigned j=0;j<n_node;j++)
954  {
955  node_pt(j)->x(0)=backed_up_nodal_position(j,0);
956  node_pt(j)->x(1)=backed_up_nodal_position(j,1);
957  }
958 }
959 
960 //==============================================================
961 /// Get strain-rate tensor: \f$ e_{ij} \f$ where
962 /// \f$ i,j = r,z,\theta \f$ (in that order). Extrapolated
963 /// from history values.
964 //==============================================================
967  DenseMatrix<double>& strainrate) const
968 {
969 
970 #ifdef PARANOID
971  if ((strainrate.ncol()!=3)||(strainrate.nrow()!=3))
972  {
973  std::ostringstream error_message;
974  error_message << "The strain rate has incorrect dimensions "
975  << strainrate.ncol() << " x "
976  << strainrate.nrow() << " Not 3" << std::endl;
977 
978  throw OomphLibError(error_message.str(),
979  OOMPH_CURRENT_FUNCTION,
980  OOMPH_EXCEPTION_LOCATION);
981  }
982 #endif
983 
984  // Get strain rates at two previous times
985  DenseMatrix<double> strain_rate_minus_two(3,3);
986  strain_rate(2,s,strain_rate_minus_two);
987 
988  DenseMatrix<double> strain_rate_minus_one(3,3);
989  strain_rate(1,s,strain_rate_minus_one);
990 
991  // Get timestepper from first node
993 
994  // Current and previous timesteps
995  double dt_current=time_stepper_pt->time_pt()->dt(0);
996  double dt_prev=time_stepper_pt->time_pt()->dt(1);
997 
998  // Extrapolate
999  for (unsigned i=0;i<3;i++)
1000  {
1001  for (unsigned j=0;j<3;j++)
1002  {
1003  // Rate of changed based on previous two solutions
1004  double slope=(strain_rate_minus_one(i,j)-strain_rate_minus_two(i,j))/
1005  dt_prev;
1006 
1007  // Extrapolated value from previous computed one to current one
1008  strainrate(i,j)=strain_rate_minus_one(i,j)+slope*dt_current;
1009  }
1010  }
1011 }
1012 
1013 //==============================================================
1014 /// \short Get integral of kinetic energy over element:
1015 //==============================================================
1017 {
1018 
1019  throw OomphLibError(
1020  "Check the kinetic energy calculation for axisymmetric NSt",
1021  OOMPH_CURRENT_FUNCTION,
1022  OOMPH_EXCEPTION_LOCATION);
1023 
1024  // Initialise
1025  double kin_en=0.0;
1026 
1027  //Set the value of Nintpt
1028  unsigned Nintpt = integral_pt()->nweight();
1029 
1030  //Set the Vector to hold local coordinates
1031  Vector<double> s(2);
1032 
1033  //Loop over the integration points
1034  for(unsigned ipt=0;ipt<Nintpt;ipt++)
1035  {
1036  //Assign values of s
1037  for(unsigned i=0;i<2;i++) {s[i] = integral_pt()->knot(ipt,i);}
1038 
1039  //Get the integral weight
1040  double w = integral_pt()->weight(ipt);
1041 
1042  //Get Jacobian of mapping
1043  double J = J_eulerian(s);
1044 
1045  // Loop over directions
1046  double veloc_squared=0.0;
1047  for(unsigned i=0;i<3;i++)
1048  {
1049  veloc_squared+=interpolated_u_axi_nst(s,i)*interpolated_u_axi_nst(s,i);
1050  }
1051 
1052  kin_en+=0.5*veloc_squared*w*J*interpolated_x(s,0);
1053  }
1054 
1055  return kin_en;
1056 
1057 }
1058 
1059 //==============================================================
1060 /// Return pressure integrated over the element
1061 //==============================================================
1064 {
1065 
1066  // Initialise
1067  double press_int=0;
1068 
1069  //Set the value of Nintpt
1070  unsigned Nintpt = integral_pt()->nweight();
1071 
1072  //Set the Vector to hold local coordinates
1073  Vector<double> s(2);
1074 
1075  //Loop over the integration points
1076  for(unsigned ipt=0;ipt<Nintpt;ipt++)
1077  {
1078 
1079  //Assign values of s
1080  for(unsigned i=0;i<2;i++)
1081  {
1082  s[i] = integral_pt()->knot(ipt,i);
1083  }
1084 
1085  //Get the integral weight
1086  double w = integral_pt()->weight(ipt);
1087 
1088  //Get Jacobian of mapping
1089  double J = J_eulerian(s);
1090 
1091  //Premultiply the weights and the Jacobian
1092  double W = w*J*interpolated_x(s,0);
1093 
1094  // Get pressure
1095  double press=interpolated_p_axi_nst(s);
1096 
1097  // Add
1098  press_int+=press*W;
1099 
1100  }
1101 
1102  return press_int;
1103 
1104 }
1105 
1106 //==============================================================
1107 /// Get max. and min. strain rate invariant and visocosity
1108 /// over all integration points in element
1109 //==============================================================
1112  double& max_invariant,
1113  double& min_viscosity,
1114  double& max_viscosity) const
1115  {
1116  // Initialise
1117  min_invariant=DBL_MAX;
1118  max_invariant=-DBL_MAX;
1119  min_viscosity=DBL_MAX;
1120  max_viscosity=-DBL_MAX;
1121 
1122  //Number of integration points
1123  unsigned Nintpt = integral_pt()->nweight();
1124 
1125  //Set the Vector to hold local coordinates (two dimensions)
1126  Vector<double> s(2);
1127 
1128  //Loop over the integration points
1129  for(unsigned ipt=0;ipt<Nintpt;ipt++)
1130  {
1131  //Assign values of s
1132  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
1133 
1134  // The strainrate
1135  DenseMatrix<double> strainrate(3,3,0.0);
1136  strain_rate(s,strainrate);
1137 
1138  // Calculate the second invariant
1140  strainrate);
1141 
1142  // Get the viscosity according to the constitutive equation
1143  double viscosity=Constitutive_eqn_pt->viscosity(second_invariant);
1144 
1145  min_invariant=std::min(min_invariant,second_invariant);
1146  max_invariant=std::max(max_invariant,second_invariant);
1147  min_viscosity=std::min(min_viscosity,viscosity);
1148  max_viscosity=std::max(max_viscosity,viscosity);
1149  }
1150  }
1151 
1152 
1153 
1154 
1155 //==============================================================
1156 /// Compute the residuals for the Navier--Stokes
1157 /// equations; flag=1(or 0): do (or don't) compute the
1158 /// Jacobian as well.
1159 //==============================================================
1162  DenseMatrix<double> &jacobian,
1163  DenseMatrix<double> &mass_matrix,
1164  unsigned flag)
1165 {
1166  // Return immediately if there are no dofs
1167  if (ndof()==0) return;
1168 
1169  //Find out how many nodes there are
1170  unsigned n_node = nnode();
1171 
1172  // Get continuous time from timestepper of first node
1173  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
1174 
1175  //Find out how many pressure dofs there are
1176  unsigned n_pres = npres_axi_nst();
1177 
1178  //Get the nodal indices at which the velocity is stored
1179  unsigned u_nodal_index[3];
1180  for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}
1181 
1182  //Set up memory for the shape and test functions
1183  //Note that there are only two dimensions, r and z in this problem
1184  Shape psif(n_node), testf(n_node);
1185  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
1186 
1187  //Set up memory for pressure shape and test functions
1188  Shape psip(n_pres), testp(n_pres);
1189 
1190  //Number of integration points
1191  unsigned Nintpt = integral_pt()->nweight();
1192 
1193  //Set the Vector to hold local coordinates (two dimensions)
1194  Vector<double> s(2);
1195 
1196  //Get Physical Variables from Element
1197  //Reynolds number must be multiplied by the density ratio
1198  double scaled_re = re()*density_ratio();
1199  double scaled_re_st = re_st()*density_ratio();
1200  double scaled_re_inv_fr = re_invfr()*density_ratio();
1201  double scaled_re_inv_ro = re_invro()*density_ratio();
1202  // double visc_ratio = viscosity_ratio();
1203  Vector<double> G = g();
1204 
1205  //Integers used to store the local equation and unknown numbers
1206  int local_eqn=0, local_unknown=0;
1207 
1208  //Loop over the integration points
1209  for(unsigned ipt=0;ipt<Nintpt;ipt++)
1210  {
1211  //Assign values of s
1212  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
1213  //Get the integral weight
1214  double w = integral_pt()->weight(ipt);
1215 
1216  //Call the derivatives of the shape and test functions
1217  double J =
1218  dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,psif,dpsifdx,testf,dtestfdx);
1219 
1220  //Call the pressure shape and test functions
1221  pshape_axi_nst(s,psip,testp);
1222 
1223  //Premultiply the weights and the Jacobian
1224  double W = w*J;
1225 
1226  //Allocate storage for the position and the derivative of the
1227  //mesh positions wrt time
1229  Vector<double> mesh_velocity(2,0.0);
1230  //Allocate storage for the pressure, velocity components and their
1231  //time and spatial derivatives
1232  double interpolated_p=0.0;
1233  Vector<double> interpolated_u(3,0.0);
1234  Vector<double> dudt(3,0.0);
1235  DenseMatrix<double> interpolated_dudx(3,2,0.0);
1236 
1237  //Calculate pressure at integration point
1238  for(unsigned l=0;l<n_pres;l++) {interpolated_p += p_axi_nst(l)*psip[l];}
1239 
1240  //Calculate velocities and derivatives at integration point
1241 
1242  // Loop over nodes
1243  for(unsigned l=0;l<n_node;l++)
1244  {
1245  //Cache the shape function
1246  const double psif_ = psif(l);
1247  //Loop over the two coordinate directions
1248  for(unsigned i=0;i<2;i++)
1249  {
1250  interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
1251  }
1252 
1253  //Loop over the three velocity directions
1254  for(unsigned i=0;i<3;i++)
1255  {
1256  //Get the u_value
1257  const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
1258  interpolated_u[i] += u_value*psif_;
1259  dudt[i]+= du_dt_axi_nst(l,i)*psif_;
1260  //Loop over derivative directions
1261  for(unsigned j=0;j<2;j++)
1262  {interpolated_dudx(i,j) += u_value*dpsifdx(l,j);}
1263  }
1264  }
1265 
1266  //Get the mesh velocity if ALE is enabled
1267  if(!ALE_is_disabled)
1268  {
1269  // Loop over nodes
1270  for(unsigned l=0;l<n_node;l++)
1271  {
1272  //Loop over the two coordinate directions
1273  for(unsigned i=0;i<2;i++)
1274  {
1275  mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif(l);
1276  }
1277  }
1278  }
1279 
1280  // The strainrate used to compute the second invariant
1281  DenseMatrix<double> strainrate_to_compute_second_invariant(3,3,0.0);
1282 
1283  // the strainrate used to calculate the second invariant
1284  // can be either the current one or the one extrapolated from
1285  // previous velocity values
1287  {
1288  strain_rate(s,strainrate_to_compute_second_invariant);
1289  }
1290  else
1291  {
1292  extrapolated_strain_rate(ipt,strainrate_to_compute_second_invariant);
1293  }
1294 
1295  // Calculate the second invariant
1297  strainrate_to_compute_second_invariant);
1298 
1299  // Get the viscosity according to the constitutive equation
1300  double viscosity=Constitutive_eqn_pt->viscosity(second_invariant);
1301 
1302  //Get the user-defined body force terms
1303  Vector<double> body_force(3);
1304  get_body_force_axi_nst(time,ipt,s,interpolated_x,body_force);
1305 
1306  //Get the user-defined source function
1307  double source = get_source_fct(time,ipt,interpolated_x);
1308 
1309  //Get the user-defined viscosity function
1310  double visc_ratio;
1312  s,
1313  interpolated_x,
1314  visc_ratio);
1315 
1316  //r is the first position component
1317  double r = interpolated_x[0];
1318 
1319  // obacht set up vectors of the viscosity differentiated w.r.t.
1320  // the velocity components (radial, axial, azimuthal)
1321  Vector<double> dviscosity_dUr(n_node,0.0);
1322  Vector<double> dviscosity_dUz(n_node,0.0);
1323  Vector<double> dviscosity_dUphi(n_node,0.0);
1324 
1326  {
1327  // Calculate the derivate of the viscosity w.r.t. the second invariant
1328  double dviscosity_dsecond_invariant=
1329  Constitutive_eqn_pt->dviscosity_dinvariant(second_invariant);
1330 
1331  // FD step
1332  //double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
1333 
1334  // calculate a reference strainrate
1335  DenseMatrix<double> strainrate_ref(3,3,0.0);
1336  strain_rate(s,strainrate_ref);
1337 
1338  // pre-compute the derivative of the second invariant w.r.t. the
1339  // entries in the rate of strain tensor
1340  DenseMatrix<double> dinvariant_dstrainrate(3,3,0.0);
1341 
1342  // d I_2 / d epsilon_{r,r}
1343  dinvariant_dstrainrate(0,0)=strainrate_ref(1,1)+strainrate_ref(2,2);
1344  // d I_2 / d epsilon_{z,z}
1345  dinvariant_dstrainrate(1,1)=strainrate_ref(0,0)+strainrate_ref(2,2);
1346  // d I_2 / d epsilon_{phi,phi}
1347  dinvariant_dstrainrate(2,2)=strainrate_ref(0,0)+strainrate_ref(1,1);
1348  // d I_2 / d epsilon_{r,z}
1349  dinvariant_dstrainrate(0,1)=-strainrate_ref(1,0);
1350  // d I_2 / d epsilon_{z,r}
1351  dinvariant_dstrainrate(1,0)=-strainrate_ref(0,1);
1352  // d I_2 / d epsilon_{r,phi}
1353  dinvariant_dstrainrate(0,2)=-strainrate_ref(2,0);
1354  // d I_2 / d epsilon_{phi,r}
1355  dinvariant_dstrainrate(2,0)=-strainrate_ref(0,2);
1356  // d I_2 / d epsilon_{phi,z}
1357  dinvariant_dstrainrate(2,1)=-strainrate_ref(1,2);
1358  // d I_2 / d epsilon_{z,phi}
1359  dinvariant_dstrainrate(1,2)=-strainrate_ref(2,1);
1360 
1361  // loop over the nodes
1362  for(unsigned l=0;l<n_node;l++)
1363  {
1364  // Get pointer to l-th local node
1365  //Node* nod_pt = node_pt(l);
1366 
1367  // loop over the three velocity components
1368  for(unsigned i=0;i<3;i++)
1369  {
1370  // back up
1371  //double backup = nod_pt->value(u_nodal_index[i]);
1372 
1373  // do the FD step
1374  //nod_pt->set_value(u_nodal_index[i],
1375  // nod_pt->value(u_nodal_index[i])+eps_fd);
1376 
1377  // calculate updated strainrate
1378  //DenseMatrix<double> strainrate_fd(3,3,0.0);
1379  //strain_rate(s,strainrate_fd);
1380 
1381  // initialise the derivative of the second invariant w.r.t. the
1382  // unknown velocity U_{i,l}
1383  double dinvariant_dunknown=0.0;
1384 
1385  // loop over the entries of the rate of strain tensor
1386  for(unsigned m=0;m<3;m++)
1387  {
1388  for(unsigned n=0;n<3;n++)
1389  {
1390 
1391  // initialise the derivative of the strainrate w.r.t. the
1392  // unknown velocity U_{i,l}
1393  double dstrainrate_dunknown=0.0;
1394 
1395  // switch based on first index
1396  switch(m)
1397  {
1398  // epsilon_{r ...}
1399  case 0:
1400 
1401  // switch for second index
1402  switch(n)
1403  {
1404  // epsilon_{r r}
1405  case 0:
1406  if(i==0)
1407  {
1408  dstrainrate_dunknown=dpsifdx(l,0);
1409  }
1410  break;
1411 
1412  // epsilon_{r z}
1413  case 1:
1414  if(i==0)
1415  {
1416  dstrainrate_dunknown=0.5*dpsifdx(l,1);
1417  }
1418  else if(i==1)
1419  {
1420  dstrainrate_dunknown=0.5*dpsifdx(l,0);
1421  }
1422  break;
1423 
1424  // epsilon_{r phi}
1425  case 2:
1426  if(i==2)
1427  {
1428  dstrainrate_dunknown=0.5*dpsifdx(l,0)-0.5/r*psif[l];
1429  }
1430  break;
1431 
1432  default:
1433  std::ostringstream error_stream;
1434  error_stream << "Should never get here...";
1435  throw OomphLibError(
1436  error_stream.str(),
1437  OOMPH_CURRENT_FUNCTION,
1438  OOMPH_EXCEPTION_LOCATION);
1439 
1440  }
1441 
1442  break;
1443 
1444  // epsilon_{z ...}
1445  case 1:
1446 
1447  // switch for second index
1448  switch(n)
1449  {
1450  // epsilon_{z r}
1451  case 0:
1452  if(i==0)
1453  {
1454  dstrainrate_dunknown=0.5*dpsifdx(l,1);
1455  }
1456  else if(i==1)
1457  {
1458  dstrainrate_dunknown=0.5*dpsifdx(l,0);
1459  }
1460  break;
1461 
1462  // epsilon_{z z}
1463  case 1:
1464  if(i==1)
1465  {
1466  dstrainrate_dunknown=dpsifdx(l,1);
1467  }
1468  else
1469  {
1470  //dstrainrate_dunknown=0.0;
1471  }
1472  break;
1473 
1474  // epsilon_{z phi}
1475  case 2:
1476  if(i==2)
1477  {
1478  dstrainrate_dunknown=0.5*dpsifdx(l,1);
1479  }
1480  break;
1481 
1482  default:
1483  std::ostringstream error_stream;
1484  error_stream << "Should never get here...";
1485  throw OomphLibError(
1486  error_stream.str(),
1487  OOMPH_CURRENT_FUNCTION,
1488  OOMPH_EXCEPTION_LOCATION);
1489 
1490  }
1491 
1492  break;
1493 
1494  // epsilon_{phi ...}
1495  case 2:
1496 
1497  // switch for second index
1498  switch(n)
1499  {
1500  // epsilon_{phi r}
1501  case 0:
1502  if(i==2)
1503  {
1504  dstrainrate_dunknown=0.5*dpsifdx(l,0)-0.5/r*psif[l];
1505  }
1506  break;
1507 
1508  // epsilon_{phi z}
1509  case 1:
1510  if(i==2)
1511  {
1512  dstrainrate_dunknown=0.5*dpsifdx(l,1);
1513  }
1514  break;
1515 
1516  // epsilon_{phi phi}
1517  case 2:
1518  if(i==0)
1519  {
1520  dstrainrate_dunknown=1.0/r*psif[l];
1521  }
1522  break;
1523 
1524  default:
1525  std::ostringstream error_stream;
1526  error_stream << "Should never get here...";
1527  throw OomphLibError(
1528  error_stream.str(),
1529  OOMPH_CURRENT_FUNCTION,
1530  OOMPH_EXCEPTION_LOCATION);
1531 
1532  }
1533 
1534  break;
1535 
1536  default:
1537  std::ostringstream error_stream;
1538  error_stream << "Should never get here...";
1539  throw OomphLibError(
1540  error_stream.str(),
1541  OOMPH_CURRENT_FUNCTION,
1542  OOMPH_EXCEPTION_LOCATION);
1543 
1544  }
1545  // calculate the difference
1546  //double dstrainrate_dunknown =
1547  // (strainrate_fd(m,n)-strainrate_ref(m,n))/eps_fd;
1548 
1549  dinvariant_dunknown += dinvariant_dstrainrate(m,n)*
1550  dstrainrate_dunknown;
1551  }
1552  }
1553 
1554  // // get the invariant of the updated strainrate
1555  // double second_invariant_fd=
1556  // SecondInvariantHelper::second_invariant(strainrate_fd);
1557 
1558  // // calculate the difference
1559  // double dinvariant_dunknown =
1560  // (second_invariant_fd - second_invariant)/eps_fd;
1561 
1562  switch(i)
1563  {
1564  case 0:
1565  dviscosity_dUr[l] = dviscosity_dsecond_invariant*dinvariant_dunknown;
1566  break;
1567 
1568  case 1:
1569  dviscosity_dUz[l] = dviscosity_dsecond_invariant*dinvariant_dunknown;
1570  break;
1571 
1572  case 2:
1573  dviscosity_dUphi[l] = dviscosity_dsecond_invariant*
1574  dinvariant_dunknown;
1575  break;
1576 
1577  default:
1578  std::ostringstream error_stream;
1579  error_stream << "Should never get here...";
1580  throw OomphLibError(
1581  error_stream.str(),
1582  OOMPH_CURRENT_FUNCTION,
1583  OOMPH_EXCEPTION_LOCATION);
1584  }
1585 
1586  // Reset
1587  //nod_pt->set_value(u_nodal_index[i],backup);
1588  }
1589  }
1590  }
1591 
1592  //MOMENTUM EQUATIONS
1593  //------------------
1594 
1595  //Loop over the test functions
1596  for(unsigned l=0;l<n_node;l++)
1597  {
1598  //FIRST (RADIAL) MOMENTUM EQUATION
1599  local_eqn = nodal_local_eqn(l,u_nodal_index[0]);
1600  //If it's not a boundary condition
1601  if(local_eqn >= 0)
1602  {
1603  //Add the user-defined body force terms
1604  residuals[local_eqn] +=
1605  r*body_force[0]*testf[l]*W;
1606 
1607  //Add the gravitational body force term
1608  residuals[local_eqn] += r*scaled_re_inv_fr*testf[l]*G[0]*W;
1609 
1610  //Add the pressure gradient term
1611  // Potentially pre-multiply by viscosity ratio
1613  {
1614  residuals[local_eqn] +=
1615  visc_ratio*viscosity*interpolated_p*(testf[l] + r*dtestfdx(l,0))*W;
1616  }
1617  else
1618  {
1619  residuals[local_eqn] +=
1620  interpolated_p*(testf[l] + r*dtestfdx(l,0))*W;
1621  }
1622 
1623  //Add in the stress tensor terms
1624  //The viscosity ratio needs to go in here to ensure
1625  //continuity of normal stress is satisfied even in flows
1626  //with zero pressure gradient!
1627  // stress term 1
1628  residuals[local_eqn] -= visc_ratio*viscosity*
1629  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
1630 
1631  // stress term 2
1632  residuals[local_eqn] -= visc_ratio*viscosity*r*
1633  (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
1634  dtestfdx(l,1)*W;
1635 
1636  // stress term 3
1637  residuals[local_eqn] -=
1638  visc_ratio*viscosity*(1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
1639 
1640  //Add in the inertial terms
1641  //du/dt term
1642  residuals[local_eqn] -= scaled_re_st*r*dudt[0]*testf[l]*W;
1643 
1644  //Convective terms
1645  residuals[local_eqn] -=
1646  scaled_re*(r*interpolated_u[0]*interpolated_dudx(0,0)
1647  - interpolated_u[2]*interpolated_u[2]
1648  + r*interpolated_u[1]*interpolated_dudx(0,1))*testf[l]*W;
1649 
1650  //Mesh velocity terms
1651  if(!ALE_is_disabled)
1652  {
1653  for(unsigned k=0;k<2;k++)
1654  {
1655  residuals[local_eqn] +=
1656  scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(0,k)*testf[l]*W;
1657  }
1658  }
1659 
1660  //Add the Coriolis term
1661  residuals[local_eqn] +=
1662  2.0*r*scaled_re_inv_ro*interpolated_u[2]*testf[l]*W;
1663 
1664  //CALCULATE THE JACOBIAN
1665  if(flag)
1666  {
1667  //Loop over the velocity shape functions again
1668  for(unsigned l2=0;l2<n_node;l2++)
1669  {
1670  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
1671  //Radial velocity component
1672  if(local_unknown >= 0)
1673  {
1674  if(flag==2)
1675  {
1676  //Add the mass matrix
1677  mass_matrix(local_eqn,local_unknown) +=
1678  scaled_re_st*r*psif[l2]*testf[l]*W;
1679  }
1680 
1681  // stress term 1
1682  jacobian(local_eqn,local_unknown)
1683  -= visc_ratio*viscosity*r*(1.0+Gamma[0])
1684  *dpsifdx(l2,0)*dtestfdx(l,0)*W;
1685 
1687  {
1688  // stress term 1 contribution from generalised Newtonian model
1689  jacobian(local_eqn,local_unknown)
1690  -= visc_ratio*dviscosity_dUr[l2]*
1691  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
1692  }
1693 
1694  // stress term 2
1695  jacobian(local_eqn,local_unknown)
1696  -= visc_ratio*viscosity*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
1697 
1699  {
1700  // stress term 2 contribution from generalised Newtonian model
1701  jacobian(local_eqn,local_unknown)
1702  -= visc_ratio*dviscosity_dUr[l2]*
1703  r*(interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
1704  dtestfdx(l,1)*W;
1705  }
1706 
1707  // stress term 3
1708  jacobian(local_eqn,local_unknown)
1709  -= visc_ratio*viscosity*(1.0 + Gamma[0])*psif[l2]*testf[l]*W/r;
1710 
1712  {
1713  // stress term 3 contribution from generalised Newtonian model
1714  jacobian(local_eqn,local_unknown)
1715  -= visc_ratio*dviscosity_dUr[l2]*
1716  (1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
1717  }
1718 
1719  //Add in the inertial terms
1720  //du/dt term
1721  jacobian(local_eqn,local_unknown)
1722  -= scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1723  psif[l2]*testf[l]*W;
1724 
1725  //Convective terms
1726  jacobian(local_eqn,local_unknown) -=
1727  scaled_re*(r*psif[l2]*interpolated_dudx(0,0)
1728  + r*interpolated_u[0]*dpsifdx(l2,0)
1729  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
1730 
1731  //Mesh velocity terms
1732  if(!ALE_is_disabled)
1733  {
1734  for(unsigned k=0;k<2;k++)
1735  {
1736  jacobian(local_eqn,local_unknown) +=
1737  scaled_re_st*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
1738  }
1739  }
1740  }
1741 
1742 
1743  //Axial velocity component
1744  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
1745  if(local_unknown >= 0)
1746  {
1747 
1748  // no stress term 1
1749 
1751  {
1752  // stress term 1 contribution from generalised Newtonian model
1753  jacobian(local_eqn,local_unknown)
1754  -= visc_ratio*dviscosity_dUz[l2]*
1755  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
1756  }
1757 
1758  // stress term 2
1759  jacobian(local_eqn,local_unknown) -=
1760  visc_ratio*viscosity*r*Gamma[0]*dpsifdx(l2,0)*dtestfdx(l,1)*W;
1761 
1763  {
1764  // stress term 2 contribution from generalised Newtonian model
1765  jacobian(local_eqn,local_unknown) -=
1766  visc_ratio*dviscosity_dUz[l2]*
1767  r*(interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
1768  dtestfdx(l,1)*W;
1769  }
1770 
1771  // no stress term 3
1772 
1774  {
1775  // stress term 3 contribution from generalised Newtonian model
1776  jacobian(local_eqn,local_unknown)
1777  -= visc_ratio*dviscosity_dUz[l2]*
1778  (1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
1779  }
1780 
1781  //Convective terms
1782  jacobian(local_eqn,local_unknown) -=
1783  scaled_re*r*psif[l2]*interpolated_dudx(0,1)*testf[l]*W;
1784  }
1785 
1786  //Azimuthal velocity component
1787  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
1788  if(local_unknown >= 0)
1789  {
1790 
1791  // no stress term 1
1792 
1794  {
1795  // stress term 1 contribution from generalised Newtonian model
1796  jacobian(local_eqn,local_unknown)
1797  -= visc_ratio*dviscosity_dUphi[l2]*
1798  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
1799  }
1800 
1801  // no stress term 2
1802 
1804  {
1805  // stress term 2 contribution from generalised Newtonian model
1806  jacobian(local_eqn,local_unknown) -=
1807  visc_ratio*dviscosity_dUphi[l2]*
1808  r*(interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
1809  dtestfdx(l,1)*W;
1810  }
1811 
1812  // no stress term 3
1813 
1815  {
1816  // stress term 3 contribution from generalised Newtonian model
1817  jacobian(local_eqn,local_unknown)
1818  -= visc_ratio*dviscosity_dUphi[l2]*
1819  (1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
1820  }
1821 
1822  //Convective terms
1823  jacobian(local_eqn,local_unknown) -=
1824  - scaled_re*2.0*interpolated_u[2]*psif[l2]*testf[l]*W;
1825 
1826  //Coriolis terms
1827  jacobian(local_eqn,local_unknown) +=
1828  2.0*r*scaled_re_inv_ro*psif[l2]*testf[l]*W;
1829  }
1830  }
1831 
1832  /*Now loop over pressure shape functions*/
1833  /*This is the contribution from pressure gradient*/
1834  // Potentially pre-multiply by viscosity ratio
1835  for(unsigned l2=0;l2<n_pres;l2++)
1836  {
1837  local_unknown = p_local_eqn(l2);
1838  /*If we are at a non-zero degree of freedom in the entry*/
1839  if(local_unknown >= 0)
1840  {
1842  {
1843  jacobian(local_eqn,local_unknown)
1844  += visc_ratio*viscosity*psip[l2]*(testf[l] + r*dtestfdx(l,0))*W;
1845  }
1846  else
1847  {
1848  jacobian(local_eqn,local_unknown)
1849  += psip[l2]*(testf[l] + r*dtestfdx(l,0))*W;
1850  }
1851  }
1852  }
1853  } /*End of Jacobian calculation*/
1854 
1855  } //End of if not boundary condition statement
1856 
1857  //SECOND (AXIAL) MOMENTUM EQUATION
1858  local_eqn = nodal_local_eqn(l,u_nodal_index[1]);
1859  //If it's not a boundary condition
1860  if(local_eqn >= 0)
1861  {
1862  //Add the user-defined body force terms
1863  //Remember to multiply by the density ratio!
1864  residuals[local_eqn] +=
1865  r*body_force[1]*testf[l]*W;
1866 
1867  //Add the gravitational body force term
1868  residuals[local_eqn] += r*scaled_re_inv_fr*testf[l]*G[1]*W;
1869 
1870  //Add the pressure gradient term
1871  // Potentially pre-multiply by viscosity ratio
1873  {
1874  residuals[local_eqn] += r*visc_ratio*viscosity*interpolated_p*
1875  dtestfdx(l,1)*W;
1876  }
1877  else
1878  {
1879  residuals[local_eqn] += r*interpolated_p*dtestfdx(l,1)*W;
1880  }
1881 
1882  //Add in the stress tensor terms
1883  //The viscosity ratio needs to go in here to ensure
1884  //continuity of normal stress is satisfied even in flows
1885  //with zero pressure gradient!
1886  // stress term 1
1887  residuals[local_eqn] -= visc_ratio*viscosity*
1888  r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
1889  *dtestfdx(l,0)*W;
1890 
1891  // stress term 2
1892  residuals[local_eqn] -= visc_ratio*viscosity*r*
1893  (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
1894 
1895  //Add in the inertial terms
1896  //du/dt term
1897  residuals[local_eqn] -= scaled_re_st*r*dudt[1]*testf[l]*W;
1898 
1899  //Convective terms
1900  residuals[local_eqn] -=
1901  scaled_re*(r*interpolated_u[0]*interpolated_dudx(1,0)
1902  + r*interpolated_u[1]*interpolated_dudx(1,1))*testf[l]*W;
1903 
1904  //Mesh velocity terms
1905  if(!ALE_is_disabled)
1906  {
1907  for(unsigned k=0;k<2;k++)
1908  {
1909  residuals[local_eqn] +=
1910  scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(1,k)*testf[l]*W;
1911  }
1912  }
1913 
1914  //CALCULATE THE JACOBIAN
1915  if(flag)
1916  {
1917  //Loop over the velocity shape functions again
1918  for(unsigned l2=0;l2<n_node;l2++)
1919  {
1920  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
1921  //Radial velocity component
1922  if(local_unknown >= 0)
1923  {
1924  //Add in the stress tensor terms
1925  //The viscosity ratio needs to go in here to ensure
1926  //continuity of normal stress is satisfied even in flows
1927  //with zero pressure gradient!
1928  // stress term 1
1929  jacobian(local_eqn,local_unknown) -=
1930  visc_ratio*viscosity*r*Gamma[1]*dpsifdx(l2,1)*dtestfdx(l,0)*W;
1931 
1933  {
1934  // stress term 1 contribution from generalised Newtonian model
1935  jacobian(local_eqn,local_unknown) -=
1936  visc_ratio*dviscosity_dUr[l2]*
1937  r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
1938  *dtestfdx(l,0)*W;
1939  }
1940 
1941  // no stress term 2
1942 
1944  {
1945  // stress term 2 contribution from generalised Newtonian model
1946  jacobian(local_eqn,local_unknown) -=
1947  visc_ratio*dviscosity_dUr[l2]*
1948  r*(1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
1949  }
1950 
1951  //Convective terms
1952  jacobian(local_eqn,local_unknown) -=
1953  scaled_re*r*psif[l2]*interpolated_dudx(1,0)*testf[l]*W;
1954  }
1955 
1956  //Axial velocity component
1957  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
1958  if(local_unknown >= 0)
1959  {
1960  if(flag==2)
1961  {
1962  //Add the mass matrix
1963  mass_matrix(local_eqn,local_unknown) +=
1964  scaled_re_st*r*psif[l2]*testf[l]*W;
1965  }
1966 
1967  //Add in the stress tensor terms
1968  //The viscosity ratio needs to go in here to ensure
1969  //continuity of normal stress is satisfied even in flows
1970  //with zero pressure gradient!
1971  // stress term 1
1972  jacobian(local_eqn,local_unknown) -=
1973  visc_ratio*viscosity*r*dpsifdx(l2,0)*dtestfdx(l,0)*W;
1974 
1976  {
1977  // stress term 1 contribution from generalised Newtonian model
1978  jacobian(local_eqn,local_unknown) -=
1979  visc_ratio*dviscosity_dUz[l2]*
1980  r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
1981  *dtestfdx(l,0)*W;
1982  }
1983 
1984  // stress term 2
1985  jacobian(local_eqn,local_unknown) -=
1986  visc_ratio*viscosity*r*(1.0 + Gamma[1])*dpsifdx(l2,1)*
1987  dtestfdx(l,1)*W;
1988 
1990  {
1991  // stress term 2 contribution from generalised Newtonian model
1992  jacobian(local_eqn,local_unknown) -=
1993  visc_ratio*dviscosity_dUz[l2]*
1994  r*(1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
1995  }
1996 
1997  //Add in the inertial terms
1998  //du/dt term
1999  jacobian(local_eqn,local_unknown) -=
2000  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
2001  psif[l2]*testf[l]*W;
2002 
2003  //Convective terms
2004  jacobian(local_eqn,local_unknown) -=
2005  scaled_re*(r*interpolated_u[0]*dpsifdx(l2,0)
2006  + r*psif[l2]*interpolated_dudx(1,1)
2007  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
2008 
2009 
2010  //Mesh velocity terms
2011  if(!ALE_is_disabled)
2012  {
2013  for(unsigned k=0;k<2;k++)
2014  {
2015  jacobian(local_eqn,local_unknown) +=
2016  scaled_re_st*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
2017  }
2018  }
2019  }
2020 
2021  //There are no azimithal terms in the axial momentum equation
2022  // edit: for the generalised Newtonian elements there are...
2023  //Azimuthal velocity component
2024  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
2025  if(local_unknown >= 0)
2026  {
2028  {
2029  // stress term 1 contribution from generalised Newtonian model
2030  jacobian(local_eqn,local_unknown) -=
2031  visc_ratio*dviscosity_dUphi[l2]*
2032  r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
2033  *dtestfdx(l,0)*W;
2034 
2035  // stress term 2 contribution from generalised Newtonian model
2036  jacobian(local_eqn,local_unknown) -=
2037  visc_ratio*dviscosity_dUphi[l2]*
2038  r*(1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
2039  }
2040  }
2041  } //End of loop over velocity shape functions
2042 
2043  /*Now loop over pressure shape functions*/
2044  /*This is the contribution from pressure gradient*/
2045  // Potentially pre-multiply by viscosity ratio
2046  for(unsigned l2=0;l2<n_pres;l2++)
2047  {
2048  local_unknown = p_local_eqn(l2);
2049  /*If we are at a non-zero degree of freedom in the entry*/
2050  if(local_unknown >= 0)
2051  {
2053  {
2054  jacobian(local_eqn,local_unknown)
2055  += r*visc_ratio*viscosity*psip[l2]*dtestfdx(l,1)*W;
2056  }
2057  else
2058  {
2059  jacobian(local_eqn,local_unknown)
2060  += r*psip[l2]*dtestfdx(l,1)*W;
2061  }
2062  }
2063  }
2064  } /*End of Jacobian calculation*/
2065 
2066  } //End of AXIAL MOMENTUM EQUATION
2067 
2068  //THIRD (AZIMUTHAL) MOMENTUM EQUATION
2069  local_eqn = nodal_local_eqn(l,u_nodal_index[2]);
2070  if(local_eqn >= 0)
2071  {
2072  //Add the user-defined body force terms
2073  //Remember to multiply by the density ratio!
2074  residuals[local_eqn] +=
2075  r*body_force[2]*testf[l]*W;
2076 
2077  //Add the gravitational body force term
2078  residuals[local_eqn] += r*scaled_re_inv_fr*testf[l]*G[2]*W;
2079 
2080  //There is NO pressure gradient term
2081 
2082  //Add in the stress tensor terms
2083  //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  // stress term 1
2087  residuals[local_eqn] -= visc_ratio*viscosity*
2088  (r*interpolated_dudx(2,0) -
2089  Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
2090 
2091  // stress term 2
2092  residuals[local_eqn] -= visc_ratio*viscosity*r*
2093  interpolated_dudx(2,1)*dtestfdx(l,1)*W;
2094 
2095  // stress term 3
2096  residuals[local_eqn] -= visc_ratio*viscosity*
2097  ((interpolated_u[2]/r) -Gamma[0]*interpolated_dudx(2,0))*testf[l]*W;
2098 
2099 
2100  //Add in the inertial terms
2101  //du/dt term
2102  residuals[local_eqn] -= scaled_re_st*r*dudt[2]*testf[l]*W;
2103 
2104  //Convective terms
2105  residuals[local_eqn] -=
2106  scaled_re*(r*interpolated_u[0]*interpolated_dudx(2,0)
2107  + interpolated_u[0]*interpolated_u[2]
2108  + r*interpolated_u[1]*interpolated_dudx(2,1))*testf[l]*W;
2109 
2110  //Mesh velocity terms
2111  if(!ALE_is_disabled)
2112  {
2113  for(unsigned k=0;k<2;k++)
2114  {
2115  residuals[local_eqn] +=
2116  scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(2,k)*testf[l]*W;
2117  }
2118  }
2119 
2120  //Add the Coriolis term
2121  residuals[local_eqn] -=
2122  2.0*r*scaled_re_inv_ro*interpolated_u[0]*testf[l]*W;
2123 
2124  //CALCULATE THE JACOBIAN
2125  if(flag)
2126  {
2127  //Loop over the velocity shape functions again
2128  for(unsigned l2=0;l2<n_node;l2++)
2129  {
2130  //Radial velocity component
2131  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
2132  if(local_unknown >= 0)
2133  {
2135  {
2136  // stress term 1 contribution from generalised Newtonian model
2137  jacobian(local_eqn,local_unknown) -=
2138  visc_ratio*dviscosity_dUr[l2]*
2139  (r*interpolated_dudx(2,0) -
2140  Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
2141 
2142  // stress term 2 contribution from generalised Newtonian model
2143  jacobian(local_eqn,local_unknown) -=
2144  visc_ratio*dviscosity_dUr[l2]*
2145  r*interpolated_dudx(2,1)*dtestfdx(l,1)*W;
2146 
2147  // stress term 3 contribution from generalised Newtonian model
2148  jacobian(local_eqn,local_unknown) -=
2149  visc_ratio*dviscosity_dUr[l2]*
2150  ((interpolated_u[2]/r) -Gamma[0]*interpolated_dudx(2,0))*
2151  testf[l]*W;
2152  }
2153 
2154  //Convective terms
2155  jacobian(local_eqn,local_unknown) -=
2156  scaled_re*(r*psif[l2]*interpolated_dudx(2,0)
2157  + psif[l2]*interpolated_u[2])*testf[l]*W;
2158 
2159  //Coriolis term
2160  jacobian(local_eqn,local_unknown) -=
2161  2.0*r*scaled_re_inv_ro*psif[l2]*testf[l]*W;
2162  }
2163 
2164  //Axial velocity component
2165  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
2166  if(local_unknown >= 0)
2167  {
2169  {
2170  // stress term 1 contribution from generalised Newtonian model
2171  jacobian(local_eqn,local_unknown) -=
2172  visc_ratio*dviscosity_dUz[l2]*
2173  (r*interpolated_dudx(2,0) -
2174  Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
2175 
2176  // stress term 2 contribution from generalised Newtonian model
2177  jacobian(local_eqn,local_unknown) -=
2178  visc_ratio*dviscosity_dUz[l2]*
2179  r*interpolated_dudx(2,1)*dtestfdx(l,1)*W;
2180 
2181  // stress term 3 contribution from generalised Newtonian model
2182  jacobian(local_eqn,local_unknown) -=
2183  visc_ratio*dviscosity_dUz[l2]*
2184  ((interpolated_u[2]/r) -Gamma[0]*interpolated_dudx(2,0))*
2185  testf[l]*W;
2186  }
2187 
2188  //Convective terms
2189  jacobian(local_eqn,local_unknown) -=
2190  scaled_re*r*psif[l2]*interpolated_dudx(2,1)*testf[l]*W;
2191  }
2192 
2193  //Azimuthal velocity component
2194  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
2195  if(local_unknown >= 0)
2196  {
2197  if(flag==2)
2198  {
2199  //Add the mass matrix
2200  mass_matrix(local_eqn,local_unknown) +=
2201  scaled_re_st*r*psif[l2]*testf[l]*W;
2202  }
2203 
2204  //Add in the stress tensor terms
2205  //The viscosity ratio needs to go in here to ensure
2206  //continuity of normal stress is satisfied even in flows
2207  //with zero pressure gradient!
2208  // stress term 1
2209  jacobian(local_eqn,local_unknown) -=
2210  visc_ratio*viscosity*(r*dpsifdx(l2,0) -
2211  Gamma[0]*psif[l2])*dtestfdx(l,0)*W;
2212 
2214  {
2215  // stress term 1 contribution from generalised Newtonian model
2216  jacobian(local_eqn,local_unknown) -=
2217  visc_ratio*dviscosity_dUphi[l2]*
2218  (r*interpolated_dudx(2,0) -
2219  Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
2220  }
2221 
2222  // stress term 2
2223  jacobian(local_eqn,local_unknown) -=
2224  visc_ratio*viscosity*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
2225 
2227  {
2228  // stress term 2 contribution from generalised Newtonian model
2229  jacobian(local_eqn,local_unknown) -=
2230  visc_ratio*dviscosity_dUphi[l2]*
2231  r*interpolated_dudx(2,1)*dtestfdx(l,1)*W;
2232  }
2233 
2234  // stress term 3
2235  jacobian(local_eqn,local_unknown) -=
2236  visc_ratio*viscosity*((psif[l2]/r) - Gamma[0]*dpsifdx(l2,0))
2237  *testf[l]*W;
2238 
2240  {
2241  // stress term 3 contribution from generalised Newtonian model
2242  jacobian(local_eqn,local_unknown) -=
2243  visc_ratio*dviscosity_dUphi[l2]*
2244  ((interpolated_u[2]/r) -Gamma[0]*interpolated_dudx(2,0))*
2245  testf[l]*W;
2246  }
2247 
2248  //Add in the inertial terms
2249  //du/dt term
2250  jacobian(local_eqn,local_unknown) -=
2251  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
2252  psif[l2]*testf[l]*W;
2253 
2254  //Convective terms
2255  jacobian(local_eqn,local_unknown) -=
2256  scaled_re*(r*interpolated_u[0]*dpsifdx(l2,0)
2257  + interpolated_u[0]*psif[l2]
2258  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
2259 
2260  //Mesh velocity terms
2261  if(!ALE_is_disabled)
2262  {
2263  for(unsigned k=0;k<2;k++)
2264  {
2265  jacobian(local_eqn,local_unknown) +=
2266  scaled_re_st*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
2267  }
2268  }
2269  }
2270  }
2271 
2272  //There are no pressure terms
2273  } //End of Jacobian
2274 
2275  } //End of AZIMUTHAL EQUATION
2276 
2277  } //End of loop over shape functions
2278 
2279 
2280  //CONTINUITY EQUATION
2281  //-------------------
2282 
2283  //Loop over the shape functions
2284  for(unsigned l=0;l<n_pres;l++)
2285  {
2286  local_eqn = p_local_eqn(l);
2287  //If not a boundary conditions
2288  if(local_eqn >= 0)
2289  {
2290  // Source term
2291  residuals[local_eqn] -= r*source*testp[l]*W;
2292 
2293  //Gradient terms
2294  // Potentially pre-multiply by viscosity ratio
2296  {
2297  residuals[local_eqn] +=
2298  visc_ratio*viscosity*(interpolated_u[0] + r*interpolated_dudx(0,0)
2299  + r*interpolated_dudx(1,1))*testp[l]*W;
2300  }
2301  else
2302  {
2303  residuals[local_eqn] +=
2304  (interpolated_u[0] + r*interpolated_dudx(0,0)
2305  + r*interpolated_dudx(1,1))*testp[l]*W;
2306  }
2307 
2308  /*CALCULATE THE JACOBIAN*/
2309  if(flag)
2310  {
2311  /*Loop over the velocity shape functions*/
2312  for(unsigned l2=0;l2<n_node;l2++)
2313  {
2314  //Radial velocity component
2315  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
2316  if(local_unknown >= 0)
2317  {
2319  {
2320  jacobian(local_eqn,local_unknown) +=
2321  visc_ratio*viscosity*(psif[l2] + r*dpsifdx(l2,0))*testp[l]*W;
2322  }
2323  else
2324  {
2325  jacobian(local_eqn,local_unknown) +=
2326  (psif[l2] + r*dpsifdx(l2,0))*testp[l]*W;
2327  }
2328  }
2329 
2330  //Axial velocity component
2331  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
2332  if(local_unknown >= 0)
2333  {
2335  {
2336  jacobian(local_eqn,local_unknown) +=
2337  r*visc_ratio*viscosity*dpsifdx(l2,1)*testp[l]*W;
2338  }
2339  else
2340  {
2341  jacobian(local_eqn,local_unknown) +=
2342  r*dpsifdx(l2,1)*testp[l]*W;
2343  }
2344  }
2345  } /*End of loop over l2*/
2346  } /*End of Jacobian calculation*/
2347 
2348  } //End of if not boundary condition
2349 
2350  } //End of loop over l
2351  }
2352 
2353 }
2354 
2355 
2356 
2357 //======================================================================
2358 /// Compute derivatives of elemental residual vector with respect
2359 /// to nodal coordinates.
2360 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
2361 /// Overloads the FD-based version in the FiniteElement base class.
2362 //======================================================================
2365  RankThreeTensor<double>& dresidual_dnodal_coordinates)
2366 {
2367 
2368  throw OomphLibError(
2369  "This has not been checked for generalised Newtonian elements!",
2370  OOMPH_CURRENT_FUNCTION,
2371  OOMPH_EXCEPTION_LOCATION);
2372 
2373  // Return immediately if there are no dofs
2374  if(ndof()==0) { return; }
2375 
2376  // Determine number of nodes in element
2377  const unsigned n_node = nnode();
2378 
2379  // Get continuous time from timestepper of first node
2380  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
2381 
2382  // Determine number of pressure dofs in element
2383  const unsigned n_pres = npres_axi_nst();
2384 
2385  // Find the indices at which the local velocities are stored
2386  unsigned u_nodal_index[3];
2387  for(unsigned i=0;i<3;i++) { u_nodal_index[i] = u_index_axi_nst(i); }
2388 
2389  // Set up memory for the shape and test functions
2390  // Note that there are only two dimensions, r and z, in this problem
2391  Shape psif(n_node), testf(n_node);
2392  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
2393 
2394  // Set up memory for pressure shape and test functions
2395  Shape psip(n_pres), testp(n_pres);
2396 
2397  // Deriatives of shape fct derivatives w.r.t. nodal coords
2398  RankFourTensor<double> d_dpsifdx_dX(2,n_node,n_node,2);
2399  RankFourTensor<double> d_dtestfdx_dX(2,n_node,n_node,2);
2400 
2401  // Derivative of Jacobian of mapping w.r.t. to nodal coords
2402  DenseMatrix<double> dJ_dX(2,n_node);
2403 
2404  // Derivatives of derivative of u w.r.t. nodal coords
2405  // Note that the entry d_dudx_dX(p,q,i,k) contains the derivative w.r.t.
2406  // nodal coordinate X_pq of du_i/dx_k. Since there are three velocity
2407  // components, i can take on the values 0, 1 and 2, while k and p only
2408  // take on the values 0 and 1 since there are only two spatial dimensions.
2409  RankFourTensor<double> d_dudx_dX(2,n_node,3,2);
2410 
2411  // Derivatives of nodal velocities w.r.t. nodal coords:
2412  // Assumption: Interaction only local via no-slip so
2413  // X_pq only affects U_iq.
2414  // Note that the entry d_U_dX(p,q,i) contains the derivative w.r.t. nodal
2415  // coordinate X_pq of nodal value U_iq. In other words this entry is the
2416  // derivative of the i-th velocity component w.r.t. the p-th spatial
2417  // coordinate, taken at the q-th local node.
2418  RankThreeTensor<double> d_U_dX(2,n_node,3,0.0);
2419 
2420  // Determine the number of integration points
2421  const unsigned n_intpt = integral_pt()->nweight();
2422 
2423  // Vector to hold local coordinates (two dimensions)
2424  Vector<double> s(2);
2425 
2426  // Get physical variables from element
2427  // (Reynolds number must be multiplied by the density ratio)
2428  const double scaled_re = re()*density_ratio();
2429  const double scaled_re_st = re_st()*density_ratio();
2430  const double scaled_re_inv_fr = re_invfr()*density_ratio();
2431  const double scaled_re_inv_ro = re_invro()*density_ratio();
2432  const double visc_ratio = viscosity_ratio();
2433  const Vector<double> G = g();
2434 
2435  // FD step
2437 
2438  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
2439  // Assumption: Interaction only local via no-slip so
2440  // X_ij only affects U_ij.
2441  bool element_has_node_with_aux_node_update_fct=false;
2442  for(unsigned q=0;q<n_node;q++)
2443  {
2444  // Get pointer to q-th local node
2445  Node* nod_pt = node_pt(q);
2446 
2447  // Only compute if there's a node-update fct involved
2448  if(nod_pt->has_auxiliary_node_update_fct_pt())
2449  {
2450  element_has_node_with_aux_node_update_fct=true;
2451 
2452  // This functionality has not been tested so issue a warning
2453  // to this effect
2454  std::ostringstream warning_stream;
2455  warning_stream << "\nThe functionality to evaluate the additional"
2456  << "\ncontribution to the deriv of the residual eqn"
2457  << "\nw.r.t. the nodal coordinates which comes about"
2458  << "\nif a node's values are updated using an auxiliary"
2459  << "\nnode update function has NOT been tested for"
2460  << "\naxisymmetric Navier-Stokes elements. Use at your"
2461  << "\nown risk"
2462  << std::endl;
2464  warning_stream.str(),
2465  "GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates",
2466  OOMPH_EXCEPTION_LOCATION);
2467 
2468  // Current nodal velocity
2469  Vector<double> u_ref(3);
2470  for(unsigned i=0;i<3;i++)
2471  {
2472  u_ref[i]=*(nod_pt->value_pt(u_nodal_index[i]));
2473  }
2474 
2475  // FD
2476  for(unsigned p=0;p<2;p++)
2477  {
2478  // Make backup
2479  double backup=nod_pt->x(p);
2480 
2481  // Do FD step. No node update required as we're
2482  // attacking the coordinate directly...
2483  nod_pt->x(p) += eps_fd;
2484 
2485  // Do auxiliary node update (to apply no slip)
2487 
2488  // Loop over velocity components
2489  for(unsigned i=0;i<3;i++)
2490  {
2491  // Evaluate
2492  d_U_dX(p,q,i)=(*(nod_pt->value_pt(u_nodal_index[i]))-u_ref[i])/eps_fd;
2493  }
2494 
2495  // Reset
2496  nod_pt->x(p)=backup;
2497 
2498  // Do auxiliary node update (to apply no slip)
2500  }
2501  }
2502  }
2503 
2504  // Integer to store the local equation number
2505  int local_eqn=0;
2506 
2507  // Loop over the integration points
2508  for(unsigned ipt=0;ipt<n_intpt;ipt++)
2509  {
2510  // Assign values of s
2511  for(unsigned i=0;i<2;i++) { s[i] = integral_pt()->knot(ipt,i); }
2512 
2513  // Get the integral weight
2514  const double w = integral_pt()->weight(ipt);
2515 
2516  // Call the derivatives of the shape and test functions
2518  ipt,psif,dpsifdx,d_dpsifdx_dX,testf,dtestfdx,d_dtestfdx_dX,dJ_dX);
2519 
2520  // Call the pressure shape and test functions
2521  pshape_axi_nst(s,psip,testp);
2522 
2523  // Allocate storage for the position and the derivative of the
2524  // mesh positions w.r.t. time
2526  Vector<double> mesh_velocity(2,0.0);
2527 
2528  // Allocate storage for the pressure, velocity components and their
2529  // time and spatial derivatives
2530  double interpolated_p=0.0;
2531  Vector<double> interpolated_u(3,0.0);
2532  Vector<double> dudt(3,0.0);
2533  DenseMatrix<double> interpolated_dudx(3,2,0.0);
2534 
2535  // Calculate pressure at integration point
2536  for(unsigned l=0;l<n_pres;l++) { interpolated_p += p_axi_nst(l)*psip[l]; }
2537 
2538  // Calculate velocities and derivatives at integration point
2539  // ---------------------------------------------------------
2540 
2541  // Loop over nodes
2542  for(unsigned l=0;l<n_node;l++)
2543  {
2544  // Cache the shape function
2545  const double psif_ = psif(l);
2546 
2547  // Loop over the two coordinate directions
2548  for(unsigned i=0;i<2;i++)
2549  {
2550  interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
2551  }
2552 
2553  // Loop over the three velocity directions
2554  for(unsigned i=0;i<3;i++)
2555  {
2556  // Get the nodal value
2557  const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
2558  interpolated_u[i] += u_value*psif_;
2559  dudt[i] += du_dt_axi_nst(l,i)*psif_;
2560 
2561  // Loop over derivative directions
2562  for(unsigned j=0;j<2;j++)
2563  {
2564  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
2565  }
2566  }
2567  }
2568 
2569  // Get the mesh velocity if ALE is enabled
2570  if(!ALE_is_disabled)
2571  {
2572  // Loop over nodes
2573  for(unsigned l=0;l<n_node;l++)
2574  {
2575  // Loop over the two coordinate directions
2576  for(unsigned i=0;i<2;i++)
2577  {
2578  mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif[l];
2579  }
2580  }
2581  }
2582 
2583  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
2584  for(unsigned q=0;q<n_node;q++)
2585  {
2586  // Loop over the two coordinate directions
2587  for(unsigned p=0;p<2;p++)
2588  {
2589  // Loop over the three velocity components
2590  for(unsigned i=0;i<3;i++)
2591  {
2592  // Loop over the two coordinate directions
2593  for(unsigned k=0;k<2;k++)
2594  {
2595  double aux = 0.0;
2596 
2597  // Loop over nodes and add contribution
2598  for(unsigned j=0;j<n_node;j++)
2599  {
2600  aux +=
2601  this->raw_nodal_value(j,u_nodal_index[i])*d_dpsifdx_dX(p,q,j,k);
2602  }
2603  d_dudx_dX(p,q,i,k) = aux;
2604  }
2605  }
2606  }
2607  }
2608 
2609  // Get weight of actual nodal position/value in computation of mesh
2610  // velocity from positional/value time stepper
2611  const double pos_time_weight
2612  = node_pt(0)->position_time_stepper_pt()->weight(1,0);
2613  const double val_time_weight = node_pt(0)->time_stepper_pt()->weight(1,0);
2614 
2615  // Get the user-defined body force terms
2616  Vector<double> body_force(3);
2617  get_body_force_axi_nst(time,ipt,s,interpolated_x,body_force);
2618 
2619  // Get the user-defined source function
2620  const double source = get_source_fct(time,ipt,interpolated_x);
2621 
2622  // Note: Can use raw values and avoid bypassing hanging information
2623  // because this is the non-refineable version!
2624 
2625  // Get gradient of body force function
2626  DenseMatrix<double> d_body_force_dx(3,2,0.0);
2627  get_body_force_gradient_axi_nst(time,ipt,s,interpolated_x,
2628  d_body_force_dx);
2629 
2630  // Get gradient of source function
2631  Vector<double> source_gradient(2,0.0);
2632  get_source_fct_gradient(time,ipt,interpolated_x,source_gradient);
2633 
2634  // Cache r (the first position component)
2635  const double r = interpolated_x[0];
2636 
2637  // Assemble shape derivatives
2638  // --------------------------
2639 
2640  // ==================
2641  // MOMENTUM EQUATIONS
2642  // ==================
2643 
2644  // Loop over the test functions
2645  for(unsigned l=0;l<n_node;l++)
2646  {
2647  // Cache the test function
2648  const double testf_ = testf[l];
2649 
2650  // --------------------------------
2651  // FIRST (RADIAL) MOMENTUM EQUATION
2652  // --------------------------------
2653  local_eqn = nodal_local_eqn(l,u_nodal_index[0]);;
2654 
2655  // If it's not a boundary condition
2656  if(local_eqn >= 0)
2657  {
2658  // Loop over the two coordinate directions
2659  for(unsigned p=0;p<2;p++)
2660  {
2661  // Loop over nodes
2662  for(unsigned q=0;q<n_node;q++)
2663  {
2664 
2665  // Residual x deriv of Jacobian
2666  // ----------------------------
2667 
2668  // Add the user-defined body force terms
2669  double sum = r*body_force[0]*testf_;
2670 
2671  // Add the gravitational body force term
2672  sum += r*scaled_re_inv_fr*testf_*G[0];
2673 
2674  // Add the pressure gradient term
2675  sum += interpolated_p*(testf_ + r*dtestfdx(l,0));
2676 
2677  // Add the stress tensor terms
2678  // The viscosity ratio needs to go in here to ensure
2679  // continuity of normal stress is satisfied even in flows
2680  // with zero pressure gradient!
2681  sum -= visc_ratio*
2682  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0);
2683 
2684  sum -= visc_ratio*r*
2685  (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
2686  dtestfdx(l,1);
2687 
2688  sum -= visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]*testf_/r;
2689 
2690  // Add the du/dt term
2691  sum -= scaled_re_st*r*dudt[0]*testf_;
2692 
2693  // Add the convective terms
2694  sum -=
2695  scaled_re*(r*interpolated_u[0]*interpolated_dudx(0,0)
2696  - interpolated_u[2]*interpolated_u[2]
2697  + r*interpolated_u[1]*interpolated_dudx(0,1))*testf_;
2698 
2699  // Add the mesh velocity terms
2700  if(!ALE_is_disabled)
2701  {
2702  for(unsigned k=0;k<2;k++)
2703  {
2704  sum += scaled_re_st*r*mesh_velocity[k]
2705  *interpolated_dudx(0,k)*testf_;
2706  }
2707  }
2708 
2709  // Add the Coriolis term
2710  sum += 2.0*r*scaled_re_inv_ro*interpolated_u[2]*testf_;
2711 
2712  // Multiply through by deriv of Jacobian and integration weight
2713  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*dJ_dX(p,q)*w;
2714 
2715  // Derivative of residual x Jacobian
2716  // ---------------------------------
2717 
2718  // Body force
2719  sum = r*d_body_force_dx(0,p)*psif[q]*testf_;
2720  if(p==0) { sum += body_force[0]*psif[q]*testf_; }
2721 
2722  // Gravitational body force
2723  if(p==0) { sum += scaled_re_inv_fr*G[0]*psif[q]*testf_; }
2724 
2725  // Pressure gradient term
2726  sum += r*interpolated_p*d_dtestfdx_dX(p,q,l,0);
2727  if(p==0) { sum += interpolated_p*psif[q]*dtestfdx(l,0); }
2728 
2729  // Viscous terms
2730  sum -= r*visc_ratio*(
2731  (1.0+Gamma[0])*(
2732  d_dudx_dX(p,q,0,0)*dtestfdx(l,0)
2733  + interpolated_dudx(0,0)*d_dtestfdx_dX(p,q,l,0))
2734  + (d_dudx_dX(p,q,0,1)
2735  + Gamma[0]*d_dudx_dX(p,q,1,0))*dtestfdx(l,1)
2736  + (interpolated_dudx(0,1)
2737  + Gamma[0]*interpolated_dudx(1,0))
2738  *d_dtestfdx_dX(p,q,l,1));
2739  if(p==0)
2740  {
2741  sum -= visc_ratio*(
2742  (1.0+Gamma[0])*(
2743  interpolated_dudx(0,0)*psif[q]*dtestfdx(l,0)
2744  - interpolated_u[0]*psif[q]*testf_/(r*r))
2745  + (interpolated_dudx(0,1)
2746  + Gamma[0]*interpolated_dudx(1,0))*psif[q]*dtestfdx(l,1));
2747  }
2748 
2749  // Convective terms, including mesh velocity
2750  for(unsigned k=0;k<2;k++)
2751  {
2752  double tmp = scaled_re*interpolated_u[k];
2753  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
2754  sum -= r*tmp*d_dudx_dX(p,q,0,k)*testf_;
2755  if(p==0) { sum -= tmp*interpolated_dudx(0,k)*psif[q]*testf_; }
2756  }
2757  if(!ALE_is_disabled)
2758  {
2759  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(0,p)
2760  *psif[q]*testf_;
2761  }
2762 
2763  // du/dt term
2764  if(p==0) { sum -= scaled_re_st*dudt[0]*psif[q]*testf_; }
2765 
2766  // Coriolis term
2767  if(p==0)
2768  {
2769  sum += 2.0*scaled_re_inv_ro*interpolated_u[2]*psif[q]*testf_;
2770  }
2771 
2772  // Multiply through by Jacobian and integration weight
2773  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
2774 
2775  // Derivs w.r.t. nodal velocities
2776  // ------------------------------
2777  if(element_has_node_with_aux_node_update_fct)
2778  {
2779  // Contribution from deriv of first velocity component
2780  double tmp = scaled_re_st*r*val_time_weight*psif[q]*testf_;
2781  tmp += r*scaled_re*interpolated_dudx(0,0)*psif[q]*testf_;
2782  for(unsigned k=0;k<2;k++)
2783  {
2784  double tmp2 = scaled_re*interpolated_u[k];
2785  if(!ALE_is_disabled) { tmp2 -= scaled_re_st*mesh_velocity[k]; }
2786  tmp += r*tmp2*dpsifdx(q,k)*testf_;
2787  }
2788 
2789  tmp += r*visc_ratio*(1.0+Gamma[0])*dpsifdx(q,0)*dtestfdx(l,0);
2790  tmp += r*visc_ratio*dpsifdx(q,1)*dtestfdx(l,1);
2791  tmp += visc_ratio*(1.0+Gamma[0])*psif[q]*testf_/r;
2792 
2793  // Multiply through by dU_0q/dX_pq
2794  sum = -d_U_dX(p,q,0)*tmp;
2795 
2796  // Contribution from deriv of second velocity component
2797  tmp = scaled_re*r*interpolated_dudx(0,1)*psif[q]*testf_;
2798  tmp += r*visc_ratio*Gamma[0]*dpsifdx(q,0)*dtestfdx(l,1);
2799 
2800  // Multiply through by dU_1q/dX_pq
2801  sum -= d_U_dX(p,q,1)*tmp;
2802 
2803  // Contribution from deriv of third velocity component
2804  tmp = 2.0*(r*scaled_re_inv_ro
2805  + scaled_re*interpolated_u[2])*psif[q]*testf_;
2806 
2807  // Multiply through by dU_2q/dX_pq
2808  sum += d_U_dX(p,q,2)*tmp;
2809 
2810  // Multiply through by Jacobian and integration weight
2811  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
2812  }
2813  } // End of loop over q
2814  } // End of loop over p
2815  } // End of if it's not a boundary condition
2816 
2817  // --------------------------------
2818  // SECOND (AXIAL) MOMENTUM EQUATION
2819  // --------------------------------
2820  local_eqn = nodal_local_eqn(l,u_nodal_index[1]);;
2821 
2822  // If it's not a boundary condition
2823  if(local_eqn >= 0)
2824  {
2825  // Loop over the two coordinate directions
2826  for(unsigned p=0;p<2;p++)
2827  {
2828  // Loop over nodes
2829  for(unsigned q=0;q<n_node;q++)
2830  {
2831 
2832  // Residual x deriv of Jacobian
2833  // ----------------------------
2834 
2835  // Add the user-defined body force terms
2836  double sum = r*body_force[1]*testf_;
2837 
2838  // Add the gravitational body force term
2839  sum += r*scaled_re_inv_fr*testf_*G[1];
2840 
2841  // Add the pressure gradient term
2842  sum += r*interpolated_p*dtestfdx(l,1);
2843 
2844  // Add the stress tensor terms
2845  // The viscosity ratio needs to go in here to ensure
2846  // continuity of normal stress is satisfied even in flows
2847  // with zero pressure gradient!
2848  sum -= visc_ratio*
2849  r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
2850  *dtestfdx(l,0);
2851 
2852  sum -= visc_ratio*r*
2853  (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1);
2854 
2855  // Add the du/dt term
2856  sum -= scaled_re_st*r*dudt[1]*testf_;
2857 
2858  // Add the convective terms
2859  sum -=
2860  scaled_re*(r*interpolated_u[0]*interpolated_dudx(1,0)
2861  + r*interpolated_u[1]*interpolated_dudx(1,1))*testf_;
2862 
2863  // Add the mesh velocity terms
2864  if(!ALE_is_disabled)
2865  {
2866  for(unsigned k=0;k<2;k++)
2867  {
2868  sum +=
2869  scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(1,k)*testf_;
2870  }
2871  }
2872 
2873  // Multiply through by deriv of Jacobian and integration weight
2874  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*dJ_dX(p,q)*w;
2875 
2876  // Derivative of residual x Jacobian
2877  // ---------------------------------
2878 
2879  // Body force
2880  sum = r*d_body_force_dx(1,p)*psif[q]*testf_;
2881  if(p==0) { sum += body_force[1]*psif[q]*testf_; }
2882 
2883  // Gravitational body force
2884  if(p==0) { sum += scaled_re_inv_fr*G[1]*psif[q]*testf_; }
2885 
2886  // Pressure gradient term
2887  sum += r*interpolated_p*d_dtestfdx_dX(p,q,l,1);
2888  if(p==0) { sum += interpolated_p*psif[q]*dtestfdx(l,1); }
2889 
2890  // Viscous terms
2891  sum -= r*visc_ratio*(
2892  (d_dudx_dX(p,q,1,0) + Gamma[1]*d_dudx_dX(p,q,0,1))*dtestfdx(l,0)
2893  + (interpolated_dudx(1,0)
2894  + Gamma[1]*interpolated_dudx(0,1))*d_dtestfdx_dX(p,q,l,0)
2895  + (1.0+Gamma[1])*d_dudx_dX(p,q,1,1)*dtestfdx(l,1)
2896  + (1.0+Gamma[1])*interpolated_dudx(1,1)*d_dtestfdx_dX(p,q,l,1));
2897  if(p==0)
2898  {
2899  sum -= visc_ratio*(
2900  (interpolated_dudx(1,0)
2901  + Gamma[1]*interpolated_dudx(0,1))*psif[q]*dtestfdx(l,0)
2902  + (1.0+Gamma[1])*interpolated_dudx(1,1)*psif[q]*dtestfdx(l,1));
2903  }
2904 
2905  // Convective terms, including mesh velocity
2906  for(unsigned k=0;k<2;k++)
2907  {
2908  double tmp = scaled_re*interpolated_u[k];
2909  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
2910  sum -= r*tmp*d_dudx_dX(p,q,1,k)*testf_;
2911  if(p==0) { sum -= tmp*interpolated_dudx(1,k)*psif[q]*testf_; }
2912  }
2913  if(!ALE_is_disabled)
2914  {
2915  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(1,p)
2916  *psif[q]*testf_;
2917  }
2918 
2919  // du/dt term
2920  if(p==0) { sum -= scaled_re_st*dudt[1]*psif[q]*testf_; }
2921 
2922  // Multiply through by Jacobian and integration weight
2923  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
2924 
2925  // Derivs w.r.t. to nodal velocities
2926  // ---------------------------------
2927  if(element_has_node_with_aux_node_update_fct)
2928  {
2929  // Contribution from deriv of first velocity component
2930  double tmp = scaled_re*r*interpolated_dudx(1,0)*psif[q]*testf_;
2931  tmp += r*visc_ratio*Gamma[1]*dpsifdx(q,1)*dtestfdx(l,0);
2932 
2933  // Multiply through by dU_0q/dX_pq
2934  sum = -d_U_dX(p,q,0)*tmp;
2935 
2936  // Contribution from deriv of second velocity component
2937  tmp = scaled_re_st*r*val_time_weight*psif[q]*testf_;
2938  tmp += r*scaled_re*interpolated_dudx(1,1)*psif[q]*testf_;
2939  for(unsigned k=0;k<2;k++)
2940  {
2941  double tmp2 = scaled_re*interpolated_u[k];
2942  if(!ALE_is_disabled) { tmp2 -= scaled_re_st*mesh_velocity[k]; }
2943  tmp += r*tmp2*dpsifdx(q,k)*testf_;
2944  }
2945  tmp += r*visc_ratio*(dpsifdx(q,0)*dtestfdx(l,0)
2946  + (1.0+Gamma[1])*dpsifdx(q,1)*dtestfdx(l,1));
2947 
2948  // Multiply through by dU_1q/dX_pq
2949  sum -= d_U_dX(p,q,1)*tmp;
2950 
2951  // Multiply through by Jacobian and integration weight
2952  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
2953  }
2954  } // End of loop over q
2955  } // End of loop over p
2956  } // End of if it's not a boundary condition
2957 
2958  // -----------------------------------
2959  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
2960  // -----------------------------------
2961  local_eqn = nodal_local_eqn(l,u_nodal_index[2]);;
2962 
2963  // If it's not a boundary condition
2964  if(local_eqn >= 0)
2965  {
2966  // Loop over the two coordinate directions
2967  for(unsigned p=0;p<2;p++)
2968  {
2969  // Loop over nodes
2970  for(unsigned q=0;q<n_node;q++)
2971  {
2972 
2973  // Residual x deriv of Jacobian
2974  // ----------------------------
2975 
2976  // Add the user-defined body force terms
2977  double sum = r*body_force[2]*testf_;
2978 
2979  // Add the gravitational body force term
2980  sum += r*scaled_re_inv_fr*testf_*G[2];
2981 
2982  // There is NO pressure gradient term
2983 
2984  // Add in the stress tensor terms
2985  // The viscosity ratio needs to go in here to ensure
2986  // continuity of normal stress is satisfied even in flows
2987  // with zero pressure gradient!
2988  sum -= visc_ratio*(r*interpolated_dudx(2,0) -
2989  Gamma[0]*interpolated_u[2])*dtestfdx(l,0);
2990 
2991  sum -= visc_ratio*r*interpolated_dudx(2,1)*dtestfdx(l,1);
2992 
2993  sum -= visc_ratio*((interpolated_u[2]/r)
2994  - Gamma[0]*interpolated_dudx(2,0))*testf_;
2995 
2996  // Add the du/dt term
2997  sum -= scaled_re_st*r*dudt[2]*testf_;
2998 
2999  // Add the convective terms
3000  sum -=
3001  scaled_re*(r*interpolated_u[0]*interpolated_dudx(2,0)
3002  + interpolated_u[0]*interpolated_u[2]
3003  + r*interpolated_u[1]*interpolated_dudx(2,1))*testf_;
3004 
3005  // Add the mesh velocity terms
3006  if(!ALE_is_disabled)
3007  {
3008  for(unsigned k=0;k<2;k++)
3009  {
3010  sum +=
3011  scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(2,k)*testf_;
3012  }
3013  }
3014 
3015  // Add the Coriolis term
3016  sum -= 2.0*r*scaled_re_inv_ro*interpolated_u[0]*testf_;
3017 
3018  // Multiply through by deriv of Jacobian and integration weight
3019  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*dJ_dX(p,q)*w;
3020 
3021  // Derivative of residual x Jacobian
3022  // ---------------------------------
3023 
3024  // Body force
3025  sum = r*d_body_force_dx(2,p)*psif[q]*testf_;
3026  if(p==0) { sum += body_force[2]*psif[q]*testf_; }
3027 
3028  // Gravitational body force
3029  if(p==0) { sum += scaled_re_inv_fr*G[2]*psif[q]*testf_; }
3030 
3031  // There is NO pressure gradient term
3032 
3033  // Viscous terms
3034  sum -= visc_ratio*r*(
3035  d_dudx_dX(p,q,2,0)*dtestfdx(l,0)
3036  + interpolated_dudx(2,0)*d_dtestfdx_dX(p,q,l,0)
3037  + d_dudx_dX(p,q,2,1)*dtestfdx(l,1)
3038  + interpolated_dudx(2,1)*d_dtestfdx_dX(p,q,l,1));
3039 
3040  sum += visc_ratio*Gamma[0]*d_dudx_dX(p,q,2,0)*testf_;
3041 
3042  if(p==0)
3043  {
3044  sum -= visc_ratio*(interpolated_dudx(2,0)*psif[q]*dtestfdx(l,0)
3045  + interpolated_dudx(2,1)*psif[q]*dtestfdx(l,1)
3046  + interpolated_u[2]*psif[q]*testf_/(r*r));
3047  }
3048 
3049  // Convective terms, including mesh velocity
3050  for(unsigned k=0;k<2;k++)
3051  {
3052  double tmp = scaled_re*interpolated_u[k];
3053  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
3054  sum -= r*tmp*d_dudx_dX(p,q,2,k)*testf_;
3055  if(p==0) { sum -= tmp*interpolated_dudx(2,k)*psif[q]*testf_; }
3056  }
3057  if(!ALE_is_disabled)
3058  {
3059  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(2,p)
3060  *psif[q]*testf_;
3061  }
3062 
3063  // du/dt term
3064  if(p==0) { sum -= scaled_re_st*dudt[2]*psif[q]*testf_; }
3065 
3066  // Coriolis term
3067  if(p==0)
3068  {
3069  sum -= 2.0*scaled_re_inv_ro*interpolated_u[0]*psif[q]*testf_;
3070  }
3071 
3072  // Multiply through by Jacobian and integration weight
3073  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
3074 
3075  // Derivs w.r.t. to nodal velocities
3076  // ---------------------------------
3077  if(element_has_node_with_aux_node_update_fct)
3078  {
3079  // Contribution from deriv of first velocity component
3080  double tmp = (2.0*r*scaled_re_inv_ro
3081  + r*scaled_re*interpolated_dudx(2,0)
3082  - scaled_re*interpolated_u[2])*psif[q]*testf_;
3083 
3084  // Multiply through by dU_0q/dX_pq
3085  sum = -d_U_dX(p,q,0)*tmp;
3086 
3087  // Contribution from deriv of second velocity component
3088  // Multiply through by dU_1q/dX_pq
3089  sum -= d_U_dX(p,q,1)*scaled_re*r*interpolated_dudx(2,1)
3090  *psif[q]*testf_;
3091 
3092  // Contribution from deriv of third velocity component
3093  tmp = scaled_re_st*r*val_time_weight*psif[q]*testf_;
3094  tmp -= scaled_re*interpolated_u[0]*psif[q]*testf_;
3095  for(unsigned k=0;k<2;k++)
3096  {
3097  double tmp2 = scaled_re*interpolated_u[k];
3098  if(!ALE_is_disabled) { tmp2 -= scaled_re_st*mesh_velocity[k]; }
3099  tmp += r*tmp2*dpsifdx(q,k)*testf_;
3100  }
3101  tmp += visc_ratio*(r*dpsifdx(q,0)
3102  - Gamma[0]*psif[q])*dtestfdx(l,0);
3103  tmp += r*visc_ratio*dpsifdx(q,1)*dtestfdx(l,1);
3104  tmp += visc_ratio*((psif[q]/r) - Gamma[0]*dpsifdx(q,0))*testf_;
3105 
3106  // Multiply through by dU_2q/dX_pq
3107  sum -= d_U_dX(p,q,2)*tmp;
3108 
3109  // Multiply through by Jacobian and integration weight
3110  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
3111  }
3112  } // End of loop over q
3113  } // End of loop over p
3114  } // End of if it's not a boundary condition
3115 
3116  } // End of loop over test functions
3117 
3118 
3119  // ===================
3120  // CONTINUITY EQUATION
3121  // ===================
3122 
3123  // Loop over the shape functions
3124  for(unsigned l=0;l<n_pres;l++)
3125  {
3126  local_eqn = p_local_eqn(l);
3127 
3128  // Cache the test function
3129  const double testp_ = testp[l];
3130 
3131  // If not a boundary conditions
3132  if(local_eqn >= 0)
3133  {
3134 
3135  // Loop over the two coordinate directions
3136  for(unsigned p=0;p<2;p++)
3137  {
3138  // Loop over nodes
3139  for(unsigned q=0;q<n_node;q++)
3140  {
3141 
3142  // Residual x deriv of Jacobian
3143  //-----------------------------
3144 
3145  // Source term
3146  double aux = -r*source;
3147 
3148  // Gradient terms
3149  aux += (interpolated_u[0]
3150  + r*interpolated_dudx(0,0)
3151  + r*interpolated_dudx(1,1));
3152 
3153  // Multiply through by deriv of Jacobian and integration weight
3154  dresidual_dnodal_coordinates(local_eqn,p,q) +=
3155  aux*dJ_dX(p,q)*testp_*w;
3156 
3157  // Derivative of residual x Jacobian
3158  // ---------------------------------
3159 
3160  // Gradient of source function
3161  aux = -r*source_gradient[p]*psif[q];
3162  if(p==0) { aux -= source*psif[q]; }
3163 
3164  // Gradient terms
3165  aux += r*(d_dudx_dX(p,q,0,0) + d_dudx_dX(p,q,1,1));
3166  if(p==0)
3167  {
3168  aux += (interpolated_dudx(0,0) + interpolated_dudx(1,1))*psif[q];
3169  }
3170 
3171  // Derivs w.r.t. nodal velocities
3172  // ------------------------------
3173  if(element_has_node_with_aux_node_update_fct)
3174  {
3175  aux += d_U_dX(p,q,0)*(psif[q] + r*dpsifdx(q,0));
3176  aux += d_U_dX(p,q,1)*r*dpsifdx(q,1);
3177  }
3178 
3179  // Multiply through by Jacobian, test fct and integration weight
3180  dresidual_dnodal_coordinates(local_eqn,p,q) += aux*testp_*J*w;
3181 
3182  }
3183  }
3184  }
3185  } // End of loop over shape functions for continuity eqn
3186 
3187  } // End of loop over integration points
3188 }
3189 
3190 
3191 
3192 
3193 
3194 //==============================================================
3195 /// Compute the residuals for the Navier--Stokes
3196 /// equations; flag=1(or 0): do (or don't) compute the
3197 /// Jacobian as well.
3198 //==============================================================
3201  double* const &parameter_pt,
3202  Vector<double> &dres_dparam,
3203  DenseMatrix<double> &djac_dparam,
3204  DenseMatrix<double> &dmass_matrix_dparam,
3205  unsigned flag)
3206 {
3207  //Die if the parameter is not the Reynolds number
3208  if(parameter_pt!=this->re_pt())
3209  {
3210  std::ostringstream error_stream;
3211  error_stream <<
3212  "Cannot compute analytic jacobian for parameter addressed by "
3213  << parameter_pt << "\n";
3214  error_stream << "Can only compute derivatives wrt Re ("
3215  << Re_pt << ")\n";
3216  throw OomphLibError(
3217  error_stream.str(),
3218  OOMPH_CURRENT_FUNCTION,
3219  OOMPH_EXCEPTION_LOCATION);
3220  }
3221 
3222  //Which parameters are we differentiating with respect to
3223  bool diff_re = false;
3224  bool diff_re_st = false;
3225  bool diff_re_inv_fr = false;
3226  bool diff_re_inv_ro = false;
3227 
3228  //Set the boolean flags according to the parameter pointer
3229  if(parameter_pt==this->re_pt()) {diff_re = true;}
3230  if(parameter_pt==this->re_st_pt()) {diff_re_st = true;}
3231  if(parameter_pt==this->re_invfr_pt()) {diff_re_inv_fr = true;}
3232  if(parameter_pt==this->re_invro_pt()) {diff_re_inv_ro = true;}
3233 
3234 
3235  //Find out how many nodes there are
3236  unsigned n_node = nnode();
3237 
3238  //Find out how many pressure dofs there are
3239  unsigned n_pres = npres_axi_nst();
3240 
3241  //Get the nodal indices at which the velocity is stored
3242  unsigned u_nodal_index[3];
3243  for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}
3244 
3245  //Set up memory for the shape and test functions
3246  //Note that there are only two dimensions, r and z in this problem
3247  Shape psif(n_node), testf(n_node);
3248  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
3249 
3250  //Set up memory for pressure shape and test functions
3251  Shape psip(n_pres), testp(n_pres);
3252 
3253  //Number of integration points
3254  unsigned Nintpt = integral_pt()->nweight();
3255 
3256  //Set the Vector to hold local coordinates (two dimensions)
3257  Vector<double> s(2);
3258 
3259  //Get Physical Variables from Element
3260  //Reynolds number must be multiplied by the density ratio
3261  //double scaled_re = re()*density_ratio();
3262  //double scaled_re_st = re_st()*density_ratio();
3263  //double scaled_re_inv_fr = re_invfr()*density_ratio();
3264  //double scaled_re_inv_ro = re_invro()*density_ratio();
3265  double dens_ratio = this->density_ratio();
3266  // double visc_ratio = viscosity_ratio();
3267  Vector<double> G = g();
3268 
3269  //Integers used to store the local equation and unknown numbers
3270  int local_eqn=0, local_unknown=0;
3271 
3272  //Loop over the integration points
3273  for(unsigned ipt=0;ipt<Nintpt;ipt++)
3274  {
3275  //Assign values of s
3276  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
3277  //Get the integral weight
3278  double w = integral_pt()->weight(ipt);
3279 
3280  //Call the derivatives of the shape and test functions
3281  double J =
3282  dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,psif,dpsifdx,testf,dtestfdx);
3283 
3284  //Call the pressure shape and test functions
3285  pshape_axi_nst(s,psip,testp);
3286 
3287  //Premultiply the weights and the Jacobian
3288  double W = w*J;
3289 
3290  //Allocate storage for the position and the derivative of the
3291  //mesh positions wrt time
3293  Vector<double> mesh_velocity(2,0.0);
3294  //Allocate storage for the pressure, velocity components and their
3295  //time and spatial derivatives
3296  double interpolated_p=0.0;
3297  Vector<double> interpolated_u(3,0.0);
3298  Vector<double> dudt(3,0.0);
3299  DenseMatrix<double> interpolated_dudx(3,2,0.0);
3300 
3301  //Calculate pressure at integration point
3302  for(unsigned l=0;l<n_pres;l++) {interpolated_p += p_axi_nst(l)*psip[l];}
3303 
3304  //Calculate velocities and derivatives at integration point
3305 
3306  // Loop over nodes
3307  for(unsigned l=0;l<n_node;l++)
3308  {
3309  //Cache the shape function
3310  const double psif_ = psif(l);
3311  //Loop over the two coordinate directions
3312  for(unsigned i=0;i<2;i++)
3313  {
3314  interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
3315  }
3316  //mesh_velocity[i] += dnodal_position_dt(l,i)*psif[l];
3317 
3318  //Loop over the three velocity directions
3319  for(unsigned i=0;i<3;i++)
3320  {
3321  //Get the u_value
3322  const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
3323  interpolated_u[i] += u_value*psif_;
3324  dudt[i]+= du_dt_axi_nst(l,i)*psif_;
3325  //Loop over derivative directions
3326  for(unsigned j=0;j<2;j++)
3327  {interpolated_dudx(i,j) += u_value*dpsifdx(l,j);}
3328  }
3329  }
3330 
3331  //Get the mesh velocity if ALE is enabled
3332  if(!ALE_is_disabled)
3333  {
3334  // Loop over nodes
3335  for(unsigned l=0;l<n_node;l++)
3336  {
3337  //Loop over the two coordinate directions
3338  for(unsigned i=0;i<2;i++)
3339  {
3340  mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif(l);
3341  }
3342  }
3343  }
3344 
3345 
3346  //Get the user-defined body force terms
3347  //Vector<double> body_force(3);
3348  //get_body_force(time(),ipt,interpolated_x,body_force);
3349 
3350  //Get the user-defined source function
3351  //double source = get_source_fct(time(),ipt,interpolated_x);
3352 
3353  //Get the user-defined viscosity function
3354  double visc_ratio;
3356  s,
3357  interpolated_x,
3358  visc_ratio);
3359 
3360  //r is the first position component
3361  double r = interpolated_x[0];
3362 
3363 
3364  //MOMENTUM EQUATIONS
3365  //------------------
3366 
3367  //Loop over the test functions
3368  for(unsigned l=0;l<n_node;l++)
3369  {
3370 
3371  //FIRST (RADIAL) MOMENTUM EQUATION
3372  local_eqn = nodal_local_eqn(l,u_nodal_index[0]);
3373  //If it's not a boundary condition
3374  if(local_eqn >= 0)
3375  {
3376  //No user-defined body force terms
3377  //dres_dparam[local_eqn] +=
3378  // r*body_force[0]*testf[l]*W;
3379 
3380  //Add the gravitational body force term if the reynolds number
3381  //is equal to re_inv_fr
3382  if(diff_re_inv_fr)
3383  {
3384  dres_dparam[local_eqn] += r*dens_ratio*testf[l]*G[0]*W;
3385  }
3386 
3387  //No pressure gradient term
3388  //residuals[local_eqn] +=
3389  // interpolated_p*(testf[l] + r*dtestfdx(l,0))*W;
3390 
3391  //No in the stress tensor terms
3392  //The viscosity ratio needs to go in here to ensure
3393  //continuity of normal stress is satisfied even in flows
3394  //with zero pressure gradient!
3395  //residuals[local_eqn] -= visc_ratio*
3396  // r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
3397 
3398  //residuals[local_eqn] -= visc_ratio*r*
3399  // (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
3400  // dtestfdx(l,1)*W;
3401 
3402  //residuals[local_eqn] -=
3403  // visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
3404 
3405  //Add in the inertial terms
3406  //du/dt term
3407  if(diff_re_st)
3408  {
3409  dres_dparam[local_eqn] -= dens_ratio*r*dudt[0]*testf[l]*W;
3410  }
3411 
3412  //Convective terms
3413  if(diff_re)
3414  {
3415  dres_dparam[local_eqn] -=
3416  dens_ratio*(r*interpolated_u[0]*interpolated_dudx(0,0)
3417  - interpolated_u[2]*interpolated_u[2]
3418  + r*interpolated_u[1]*interpolated_dudx(0,1))*testf[l]*W;
3419  }
3420 
3421  //Mesh velocity terms
3422  if(!ALE_is_disabled)
3423  {
3424  if(diff_re_st)
3425  {
3426  for(unsigned k=0;k<2;k++)
3427  {
3428  dres_dparam[local_eqn] +=
3429  dens_ratio*r*mesh_velocity[k]*interpolated_dudx(0,k)*testf[l]*W;
3430  }
3431  }
3432  }
3433 
3434  //Add the Coriolis term
3435  if(diff_re_inv_ro)
3436  {
3437  dres_dparam[local_eqn] +=
3438  2.0*r*dens_ratio*interpolated_u[2]*testf[l]*W;
3439  }
3440 
3441  //CALCULATE THE JACOBIAN
3442  if(flag)
3443  {
3444  //Loop over the velocity shape functions again
3445  for(unsigned l2=0;l2<n_node;l2++)
3446  {
3447  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
3448  //Radial velocity component
3449  if(local_unknown >= 0)
3450  {
3451  if(flag==2)
3452  {
3453  if(diff_re_st)
3454  {
3455  //Add the mass matrix
3456  dmass_matrix_dparam(local_eqn,local_unknown) +=
3457  dens_ratio*r*psif[l2]*testf[l]*W;
3458  }
3459  }
3460 
3461  //Add contribution to the Jacobian matrix
3462  //jacobian(local_eqn,local_unknown)
3463  // -= visc_ratio*r*(1.0+Gamma[0])
3464  //*dpsifdx(l2,0)*dtestfdx(l,0)*W;
3465 
3466  //jacobian(local_eqn,local_unknown)
3467  // -= visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
3468 
3469  //jacobian(local_eqn,local_unknown)
3470  // -= visc_ratio*(1.0 + Gamma[0])*psif[l2]*testf[l]*W/r;
3471 
3472  //Add in the inertial terms
3473  //du/dt term
3474  if(diff_re_st)
3475  {
3476  djac_dparam(local_eqn,local_unknown)
3477  -= dens_ratio*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
3478  psif[l2]*testf[l]*W;
3479  }
3480 
3481  //Convective terms
3482  if(diff_re)
3483  {
3484  djac_dparam(local_eqn,local_unknown) -=
3485  dens_ratio*(r*psif[l2]*interpolated_dudx(0,0)
3486  + r*interpolated_u[0]*dpsifdx(l2,0)
3487  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
3488  }
3489 
3490  //Mesh velocity terms
3491  if(!ALE_is_disabled)
3492  {
3493  for(unsigned k=0;k<2;k++)
3494  {
3495  if(diff_re_st)
3496  {
3497  djac_dparam(local_eqn,local_unknown) +=
3498  dens_ratio*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
3499  }
3500  }
3501  }
3502  }
3503 
3504  //Axial velocity component
3505  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
3506  if(local_unknown >= 0)
3507  {
3508  //jacobian(local_eqn,local_unknown) -=
3509  // visc_ratio*r*Gamma[0]*dpsifdx(l2,0)*dtestfdx(l,1)*W;
3510 
3511  //Convective terms
3512  if(diff_re)
3513  {
3514  djac_dparam(local_eqn,local_unknown) -=
3515  dens_ratio*r*psif[l2]*interpolated_dudx(0,1)*testf[l]*W;
3516  }
3517  }
3518 
3519  //Azimuthal velocity component
3520  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
3521  if(local_unknown >= 0)
3522  {
3523  //Convective terms
3524  if(diff_re)
3525  {
3526  djac_dparam(local_eqn,local_unknown) -=
3527  - dens_ratio*2.0*interpolated_u[2]*psif[l2]*testf[l]*W;
3528  }
3529  //Coriolis terms
3530  if(diff_re_inv_ro)
3531  {
3532  djac_dparam(local_eqn,local_unknown) +=
3533  2.0*r*dens_ratio*psif[l2]*testf[l]*W;
3534  }
3535  }
3536  }
3537 
3538  /*Now loop over pressure shape functions*/
3539  /*This is the contribution from pressure gradient*/
3540  //for(unsigned l2=0;l2<n_pres;l2++)
3541  // {
3542  // local_unknown = p_local_eqn(l2);
3543  // /*If we are at a non-zero degree of freedom in the entry*/
3544  // if(local_unknown >= 0)
3545  // {
3546  // jacobian(local_eqn,local_unknown)
3547  // += psip[l2]*(testf[l] + r*dtestfdx(l,0))*W;
3548  // }
3549  // }
3550  } /*End of Jacobian calculation*/
3551 
3552  } //End of if not boundary condition statement
3553 
3554  //SECOND (AXIAL) MOMENTUM EQUATION
3555  local_eqn = nodal_local_eqn(l,u_nodal_index[1]);
3556  //If it's not a boundary condition
3557  if(local_eqn >= 0)
3558  {
3559  //Add the user-defined body force terms
3560  //Remember to multiply by the density ratio!
3561  //residuals[local_eqn] +=
3562  // r*body_force[1]*testf[l]*W;
3563 
3564  //Add the gravitational body force term
3565  if(diff_re_inv_fr)
3566  {
3567  dres_dparam[local_eqn] += r*dens_ratio*testf[l]*G[1]*W;
3568  }
3569 
3570  //Add the pressure gradient term
3571  //residuals[local_eqn] += r*interpolated_p*dtestfdx(l,1)*W;
3572 
3573  //Add in the stress tensor terms
3574  //The viscosity ratio needs to go in here to ensure
3575  //continuity of normal stress is satisfied even in flows
3576  //with zero pressure gradient!
3577  //residuals[local_eqn] -= visc_ratio*
3578  // r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
3579  // *dtestfdx(l,0)*W;
3580 
3581  //residuals[local_eqn] -= visc_ratio*r*
3582  // (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
3583 
3584  //Add in the inertial terms
3585  //du/dt term
3586  if(diff_re_st)
3587  {
3588  dres_dparam[local_eqn] -= dens_ratio*r*dudt[1]*testf[l]*W;
3589  }
3590 
3591  //Convective terms
3592  if(diff_re)
3593  {
3594  dres_dparam[local_eqn] -=
3595  dens_ratio*(r*interpolated_u[0]*interpolated_dudx(1,0)
3596  + r*interpolated_u[1]*interpolated_dudx(1,1))*testf[l]*W;
3597  }
3598 
3599  //Mesh velocity terms
3600  if(!ALE_is_disabled)
3601  {
3602  if(diff_re_st)
3603  {
3604  for(unsigned k=0;k<2;k++)
3605  {
3606  dres_dparam[local_eqn] +=
3607  dens_ratio*r*mesh_velocity[k]*interpolated_dudx(1,k)*testf[l]*W;
3608  }
3609  }
3610  }
3611 
3612  //CALCULATE THE JACOBIAN
3613  if(flag)
3614  {
3615  //Loop over the velocity shape functions again
3616  for(unsigned l2=0;l2<n_node;l2++)
3617  {
3618  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
3619  //Radial velocity component
3620  if(local_unknown >= 0)
3621  {
3622  //Add in the stress tensor terms
3623  //The viscosity ratio needs to go in here to ensure
3624  //continuity of normal stress is satisfied even in flows
3625  //with zero pressure gradient!
3626  //jacobian(local_eqn,local_unknown) -=
3627  //visc_ratio*r*Gamma[1]*dpsifdx(l2,1)*dtestfdx(l,0)*W;
3628 
3629  //Convective terms
3630  if(diff_re)
3631  {
3632  djac_dparam(local_eqn,local_unknown) -=
3633  dens_ratio*r*psif[l2]*interpolated_dudx(1,0)*testf[l]*W;
3634  }
3635  }
3636 
3637  //Axial velocity component
3638  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
3639  if(local_unknown >= 0)
3640  {
3641  if(flag==2)
3642  {
3643  if(diff_re_st)
3644  {
3645  //Add the mass matrix
3646  dmass_matrix_dparam(local_eqn,local_unknown) +=
3647  dens_ratio*r*psif[l2]*testf[l]*W;
3648  }
3649  }
3650 
3651 
3652  //Add in the stress tensor terms
3653  //The viscosity ratio needs to go in here to ensure
3654  //continuity of normal stress is satisfied even in flows
3655  //with zero pressure gradient!
3656  //jacobian(local_eqn,local_unknown) -=
3657  //visc_ratio*r*dpsifdx(l2,0)*dtestfdx(l,0)*W;
3658 
3659  //jacobian(local_eqn,local_unknown) -=
3660  // visc_ratio*r*(1.0 + Gamma[1])*dpsifdx(l2,1)*
3661  // dtestfdx(l,1)*W;
3662 
3663  //Add in the inertial terms
3664  //du/dt term
3665  if(diff_re_st)
3666  {
3667  djac_dparam(local_eqn,local_unknown) -=
3668  dens_ratio*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
3669  psif[l2]*testf[l]*W;
3670  }
3671 
3672  //Convective terms
3673  if(diff_re)
3674  {
3675  djac_dparam(local_eqn,local_unknown) -=
3676  dens_ratio*(r*interpolated_u[0]*dpsifdx(l2,0)
3677  + r*psif[l2]*interpolated_dudx(1,1)
3678  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
3679  }
3680 
3681  //Mesh velocity terms
3682  if(!ALE_is_disabled)
3683  {
3684  for(unsigned k=0;k<2;k++)
3685  {
3686  if(diff_re_st)
3687  {
3688  djac_dparam(local_eqn,local_unknown) +=
3689  dens_ratio*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
3690  }
3691  }
3692  }
3693  }
3694 
3695  //There are no azimithal terms in the axial momentum equation
3696  } //End of loop over velocity shape functions
3697 
3698  } /*End of Jacobian calculation*/
3699 
3700  } //End of AXIAL MOMENTUM EQUATION
3701 
3702  //THIRD (AZIMUTHAL) MOMENTUM EQUATION
3703  local_eqn = nodal_local_eqn(l,u_nodal_index[2]);
3704  if(local_eqn >= 0)
3705  {
3706  //Add the user-defined body force terms
3707  //Remember to multiply by the density ratio!
3708  //residuals[local_eqn] +=
3709  // r*body_force[2]*testf[l]*W;
3710 
3711  //Add the gravitational body force term
3712  if(diff_re_inv_fr)
3713  {
3714  dres_dparam[local_eqn] += r*dens_ratio*testf[l]*G[2]*W;
3715  }
3716 
3717  //There is NO pressure gradient term
3718 
3719  //Add in the stress tensor terms
3720  //The viscosity ratio needs to go in here to ensure
3721  //continuity of normal stress is satisfied even in flows
3722  //with zero pressure gradient!
3723  //residuals[local_eqn] -= visc_ratio*
3724  // (r*interpolated_dudx(2,0) -
3725  // Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
3726 
3727  //residuals[local_eqn] -= visc_ratio*r*
3728  // interpolated_dudx(2,1)*dtestfdx(l,1)*W;
3729 
3730  //residuals[local_eqn] -= visc_ratio*
3731  // ((interpolated_u[2]/r) - Gamma[0]*interpolated_dudx(2,0))*testf[l]*W;
3732 
3733 
3734  //Add in the inertial terms
3735  //du/dt term
3736  if(diff_re_st)
3737  {
3738  dres_dparam[local_eqn] -= dens_ratio*r*dudt[2]*testf[l]*W;
3739  }
3740 
3741  //Convective terms
3742  if(diff_re)
3743  {
3744  dres_dparam[local_eqn] -=
3745  dens_ratio*(r*interpolated_u[0]*interpolated_dudx(2,0)
3746  + interpolated_u[0]*interpolated_u[2]
3747  + r*interpolated_u[1]*interpolated_dudx(2,1))*testf[l]*W;
3748  }
3749 
3750  //Mesh velocity terms
3751  if(!ALE_is_disabled)
3752  {
3753  if(diff_re_st)
3754  {
3755  for(unsigned k=0;k<2;k++)
3756  {
3757  dres_dparam[local_eqn] +=
3758  dens_ratio*r*mesh_velocity[k]*interpolated_dudx(2,k)*testf[l]*W;
3759  }
3760  }
3761  }
3762 
3763  //Add the Coriolis term
3764  if(diff_re_inv_ro)
3765  {
3766  dres_dparam[local_eqn] -=
3767  2.0*r*dens_ratio*interpolated_u[0]*testf[l]*W;
3768  }
3769 
3770  //CALCULATE THE JACOBIAN
3771  if(flag)
3772  {
3773  //Loop over the velocity shape functions again
3774  for(unsigned l2=0;l2<n_node;l2++)
3775  {
3776  //Radial velocity component
3777  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
3778  if(local_unknown >= 0)
3779  {
3780  //Convective terms
3781  if(diff_re)
3782  {
3783  djac_dparam(local_eqn,local_unknown) -=
3784  dens_ratio*(r*psif[l2]*interpolated_dudx(2,0)
3785  + psif[l2]*interpolated_u[2])*testf[l]*W;
3786  }
3787 
3788  //Coriolis term
3789  if(diff_re_inv_ro)
3790  {
3791  djac_dparam(local_eqn,local_unknown) -=
3792  2.0*r*dens_ratio*psif[l2]*testf[l]*W;
3793  }
3794  }
3795 
3796  //Axial velocity component
3797  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
3798  if(local_unknown >= 0)
3799  {
3800  //Convective terms
3801  if(diff_re)
3802  {
3803  djac_dparam(local_eqn,local_unknown) -=
3804  dens_ratio*r*psif[l2]*interpolated_dudx(2,1)*testf[l]*W;
3805  }
3806  }
3807 
3808  //Azimuthal velocity component
3809  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
3810  if(local_unknown >= 0)
3811  {
3812  if(flag==2)
3813  {
3814  //Add the mass matrix
3815  if(diff_re_st)
3816  {
3817  dmass_matrix_dparam(local_eqn,local_unknown) +=
3818  dens_ratio*r*psif[l2]*testf[l]*W;
3819  }
3820  }
3821 
3822  //Add in the stress tensor terms
3823  //The viscosity ratio needs to go in here to ensure
3824  //continuity of normal stress is satisfied even in flows
3825  //with zero pressure gradient!
3826  //jacobian(local_eqn,local_unknown) -=
3827  // visc_ratio*(r*dpsifdx(l2,0) -
3828  // Gamma[0]*psif[l2])*dtestfdx(l,0)*W;
3829 
3830  //jacobian(local_eqn,local_unknown) -=
3831  // visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
3832 
3833  //jacobian(local_eqn,local_unknown) -=
3834  // visc_ratio*((psif[l2]/r) - Gamma[0]*dpsifdx(l2,0))
3835  // *testf[l]*W;
3836 
3837  //Add in the inertial terms
3838  //du/dt term
3839  if(diff_re_st)
3840  {
3841  djac_dparam(local_eqn,local_unknown) -=
3842  dens_ratio*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
3843  psif[l2]*testf[l]*W;
3844  }
3845 
3846  //Convective terms
3847  if(diff_re)
3848  {
3849  djac_dparam(local_eqn,local_unknown) -=
3850  dens_ratio*(r*interpolated_u[0]*dpsifdx(l2,0)
3851  + interpolated_u[0]*psif[l2]
3852  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
3853  }
3854 
3855  //Mesh velocity terms
3856  if(!ALE_is_disabled)
3857  {
3858  if(diff_re_st)
3859  {
3860  for(unsigned k=0;k<2;k++)
3861  {
3862  djac_dparam(local_eqn,local_unknown) +=
3863  dens_ratio*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
3864  }
3865  }
3866  }
3867  }
3868  }
3869 
3870  //There are no pressure terms
3871  } //End of Jacobian
3872 
3873  } //End of AZIMUTHAL EQUATION
3874 
3875  } //End of loop over shape functions
3876 
3877 
3878  //CONTINUITY EQUATION NO PARAMETERS
3879  //-------------------
3880  }
3881 
3882 }
3883 
3884 //=========================================================================
3885 /// \short Compute the hessian tensor vector products required to
3886 /// perform continuation of bifurcations analytically
3887 //=========================================================================
3890  Vector<double> const &Y,
3891  DenseMatrix<double> const &C,
3892  DenseMatrix<double> &product)
3893 {
3894  //Find out how many nodes there are
3895  unsigned n_node = nnode();
3896 
3897  //Get the nodal indices at which the velocity is stored
3898  unsigned u_nodal_index[3];
3899  for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}
3900 
3901  //Set up memory for the shape and test functions
3902  //Note that there are only two dimensions, r and z in this problem
3903  Shape psif(n_node), testf(n_node);
3904  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
3905 
3906  //Number of integration points
3907  unsigned Nintpt = integral_pt()->nweight();
3908 
3909  //Set the Vector to hold local coordinates (two dimensions)
3910  Vector<double> s(2);
3911 
3912  //Get Physical Variables from Element
3913  //Reynolds number must be multiplied by the density ratio
3914  double scaled_re = re()*density_ratio();
3915  // double visc_ratio = viscosity_ratio();
3916  Vector<double> G = g();
3917 
3918  //Integers used to store the local equation and unknown numbers
3919  int local_eqn=0, local_unknown=0, local_freedom=0;
3920 
3921  //How may dofs are there
3922  const unsigned n_dof = this->ndof();
3923 
3924  //Create a local matrix eigenvector product contribution
3925  DenseMatrix<double> jac_y(n_dof,n_dof,0.0);
3926 
3927  //Loop over the integration points
3928  for(unsigned ipt=0;ipt<Nintpt;ipt++)
3929  {
3930  //Assign values of s
3931  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
3932  //Get the integral weight
3933  double w = integral_pt()->weight(ipt);
3934 
3935  //Call the derivatives of the shape and test functions
3936  double J =
3937  dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,psif,dpsifdx,testf,dtestfdx);
3938 
3939  //Premultiply the weights and the Jacobian
3940  double W = w*J;
3941 
3942  //Allocate storage for the position and the derivative of the
3943  //mesh positions wrt time
3945  //Vector<double> mesh_velocity(2,0.0);
3946  //Allocate storage for the pressure, velocity components and their
3947  //time and spatial derivatives
3948  Vector<double> interpolated_u(3,0.0);
3949  //Vector<double> dudt(3,0.0);
3950  DenseMatrix<double> interpolated_dudx(3,2,0.0);
3951 
3952  //Calculate velocities and derivatives at integration point
3953 
3954  // Loop over nodes
3955  for(unsigned l=0;l<n_node;l++)
3956  {
3957  //Cache the shape function
3958  const double psif_ = psif(l);
3959  //Loop over the two coordinate directions
3960  for(unsigned i=0;i<2;i++)
3961  {
3962  interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
3963  }
3964 
3965  //Loop over the three velocity directions
3966  for(unsigned i=0;i<3;i++)
3967  {
3968  //Get the u_value
3969  const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
3970  interpolated_u[i] += u_value*psif_;
3971  //dudt[i]+= du_dt_axi_nst(l,i)*psif_;
3972  //Loop over derivative directions
3973  for(unsigned j=0;j<2;j++)
3974  {interpolated_dudx(i,j) += u_value*dpsifdx(l,j);}
3975  }
3976  }
3977 
3978  //Get the mesh velocity if ALE is enabled
3979  if(!ALE_is_disabled)
3980  {
3981  throw OomphLibError("Moving nodes not implemented\n",
3982  OOMPH_CURRENT_FUNCTION,
3983  OOMPH_EXCEPTION_LOCATION);
3984  }
3985 
3986  //r is the first position component
3987  double r = interpolated_x[0];
3988 
3989 
3990  //MOMENTUM EQUATIONS
3991  //------------------
3992 
3993  //Loop over the test functions
3994  for(unsigned l=0;l<n_node;l++)
3995  {
3996 
3997  //FIRST (RADIAL) MOMENTUM EQUATION
3998  local_eqn = nodal_local_eqn(l,u_nodal_index[0]);
3999  //If it's not a boundary condition
4000  if(local_eqn >= 0)
4001  {
4002  //Loop over the velocity shape functions yet again
4003  for(unsigned l3=0;l3<n_node;l3++)
4004  {
4005  //Derivative of jacobian terms with respect to radial velocity
4006  local_freedom = nodal_local_eqn(l3,u_nodal_index[0]);
4007  if(local_freedom >= 0)
4008  {
4009  //Storage for the sums
4010  double temp = 0.0;
4011 
4012  //Loop over the velocity shape functions again
4013  for(unsigned l2=0;l2<n_node;l2++)
4014  {
4015  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
4016  //Radial velocity component
4017  if(local_unknown >= 0)
4018  {
4019  //Remains of convective terms
4020  temp -=scaled_re*(r*psif[l2]*dpsifdx(l3,0)
4021  + r*psif[l3]*dpsifdx(l2,0))*
4022  Y[local_unknown]*testf[l]*W;
4023  }
4024 
4025  //Axial velocity component
4026  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
4027  if(local_unknown >= 0)
4028  {
4029  //Convective terms
4030  temp -=
4031  scaled_re*r*psif[l2]*dpsifdx(l3,1)*Y[local_unknown]*testf[l]*W;
4032  }
4033  }
4034  //Add the summed contribution to the product matrix
4035  jac_y(local_eqn,local_freedom) += temp;
4036  } //End of derivative wrt radial coordinate
4037 
4038 
4039  //Derivative of jacobian terms with respect to axial velocity
4040  local_freedom = nodal_local_eqn(l3,u_nodal_index[1]);
4041  if(local_freedom >= 0)
4042  {
4043  double temp = 0.0;
4044 
4045  //Loop over the velocity shape functions again
4046  for(unsigned l2=0;l2<n_node;l2++)
4047  {
4048  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
4049  //Radial velocity component
4050  if(local_unknown >= 0)
4051  {
4052  //Remains of convective terms
4053  temp -= scaled_re*(r*psif[l3]*dpsifdx(l2,1))*
4054  Y[local_unknown]*testf[l]*W;
4055  }
4056  }
4057  //Add the summed contribution to the product matrix
4058  jac_y(local_eqn,local_freedom) += temp;
4059  } //End of derivative wrt axial coordinate
4060 
4061  //Derivative of jacobian terms with respect to azimuthal velocity
4062  local_freedom = nodal_local_eqn(l3,u_nodal_index[2]);
4063  if(local_freedom >= 0)
4064  {
4065  double temp = 0.0;
4066 
4067  //Loop over the velocity shape functions again
4068  for(unsigned l2=0;l2<n_node;l2++)
4069  {
4070  //Azimuthal velocity component
4071  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
4072  if(local_unknown >= 0)
4073  {
4074  //Convective terms
4075  temp -=
4076  - scaled_re*2.0*psif[l3]*psif[l2]*Y[local_unknown]*testf[l]*W;
4077  }
4078  }
4079  //Add the summed contibution
4080  jac_y(local_eqn,local_freedom) += temp;
4081 
4082  } //End of if not boundary condition statement
4083  } //End of loop over freedoms
4084  } //End of RADIAL MOMENTUM EQUATION
4085 
4086 
4087  //SECOND (AXIAL) MOMENTUM EQUATION
4088  local_eqn = nodal_local_eqn(l,u_nodal_index[1]);
4089  //If it's not a boundary condition
4090  if(local_eqn >= 0)
4091  {
4092  //Loop over the velocity shape functions yet again
4093  for(unsigned l3=0;l3<n_node;l3++)
4094  {
4095  //Derivative of jacobian terms with respect to radial velocity
4096  local_freedom = nodal_local_eqn(l3,u_nodal_index[0]);
4097  if(local_freedom >= 0)
4098  {
4099  double temp = 0.0;
4100 
4101  //Loop over the velocity shape functions again
4102  for(unsigned l2=0;l2<n_node;l2++)
4103  {
4104  //Axial velocity component
4105  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
4106  if(local_unknown >= 0)
4107  {
4108  //Convective terms
4109  temp -=
4110  scaled_re*(r*psif[l3]*dpsifdx(l2,0))*
4111  Y[local_unknown]*testf[l]*W;
4112  }
4113  }
4114  jac_y(local_eqn,local_freedom) += temp;
4115 
4116  //There are no azimithal terms in the axial momentum equation
4117  } //End of loop over velocity shape functions
4118 
4119 
4120  //Derivative of jacobian terms with respect to axial velocity
4121  local_freedom = nodal_local_eqn(l3,u_nodal_index[1]);
4122  if(local_freedom >= 0)
4123  {
4124  double temp = 0.0;
4125 
4126  //Loop over the velocity shape functions again
4127  for(unsigned l2=0;l2<n_node;l2++)
4128  {
4129  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
4130  //Radial velocity component
4131  if(local_unknown >= 0)
4132  {
4133  //Convective terms
4134  temp -=
4135  scaled_re*r*psif[l2]*dpsifdx(l3,0)*Y[local_unknown]*testf[l]*W;
4136  }
4137 
4138  //Axial velocity component
4139  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
4140  if(local_unknown >= 0)
4141  {
4142  //Convective terms
4143  temp -=
4144  scaled_re*(
4145  r*psif[l2]*dpsifdx(l3,1)
4146  + r*psif[l3]*dpsifdx(l2,1))*Y[local_unknown]*testf[l]*W;
4147  }
4148 
4149  //There are no azimithal terms in the axial momentum equation
4150  } //End of loop over velocity shape functions
4151 
4152  //Add summed contributiont to jacobian product matrix
4153  jac_y(local_eqn,local_freedom) += temp;
4154  }
4155  } //End of loop over local freedoms
4156 
4157  } //End of AXIAL MOMENTUM EQUATION
4158 
4159  //THIRD (AZIMUTHAL) MOMENTUM EQUATION
4160  local_eqn = nodal_local_eqn(l,u_nodal_index[2]);
4161  if(local_eqn >= 0)
4162  {
4163  //Loop over the velocity shape functions yet again
4164  for(unsigned l3=0;l3<n_node;l3++)
4165  {
4166  //Derivative of jacobian terms with respect to radial velocity
4167  local_freedom = nodal_local_eqn(l3,u_nodal_index[0]);
4168  if(local_freedom >= 0)
4169  {
4170  double temp = 0.0;
4171 
4172  //Loop over the velocity shape functions again
4173  for(unsigned l2=0;l2<n_node;l2++)
4174  {
4175  //Azimuthal velocity component
4176  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
4177  if(local_unknown >= 0)
4178  {
4179  //Convective terms
4180  temp -=
4181  scaled_re*(r*psif[l3]*dpsifdx(l2,0)
4182  + psif[l3]*psif[l2])*Y[local_unknown]*testf[l]*W;
4183  }
4184  }
4185  //Add the summed contribution to the jacobian eigenvector sum
4186  jac_y(local_eqn,local_freedom) += temp;
4187  }
4188 
4189  //Derivative of jacobian terms with respect to axial velocity
4190  local_freedom = nodal_local_eqn(l3,u_nodal_index[1]);
4191  if(local_freedom >= 0)
4192  {
4193  double temp = 0.0;
4194 
4195  //Loop over the velocity shape functions again
4196  for(unsigned l2=0;l2<n_node;l2++)
4197  {
4198  //Azimuthal velocity component
4199  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
4200  if(local_unknown >= 0)
4201  {
4202  //Convective terms
4203  temp -=
4204  scaled_re*(r*psif[l3]*dpsifdx(l2,1))*
4205  Y[local_unknown]*testf[l]*W;
4206  }
4207  }
4208  //Add the summed contribution to the jacobian eigenvector sum
4209  jac_y(local_eqn,local_freedom) += temp;
4210  }
4211 
4212 
4213  //Derivative of jacobian terms with respect to azimuthal velocity
4214  local_freedom = nodal_local_eqn(l3,u_nodal_index[2]);
4215  if(local_freedom >= 0)
4216  {
4217  double temp = 0.0;
4218 
4219 
4220  //Loop over the velocity shape functions again
4221  for(unsigned l2=0;l2<n_node;l2++)
4222  {
4223  //Radial velocity component
4224  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
4225  if(local_unknown >= 0)
4226  {
4227  //Convective terms
4228  temp -=
4229  scaled_re*(r*psif[l2]*dpsifdx(l3,0)
4230  + psif[l2]*psif[l3])*Y[local_unknown]*testf[l]*W;
4231  }
4232 
4233  //Axial velocity component
4234  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
4235  if(local_unknown >= 0)
4236  {
4237  //Convective terms
4238  temp -=
4239  scaled_re*r*psif[l2]*dpsifdx(l3,1)*Y[local_unknown]*testf[l]*W;
4240  }
4241  }
4242  //Add the fully summed contribution
4243  jac_y(local_eqn,local_freedom) += temp;
4244  }
4245  } //End of loop over freedoms
4246 
4247  //There are no pressure terms
4248  } //End of AZIMUTHAL EQUATION
4249 
4250  } //End of loop over shape functions
4251  }
4252 
4253  //We have now assembled the matrix (d J_{ij} Y_j)/d u_{k}
4254  //and simply need to sum over the vectors
4255  const unsigned n_vec = C.nrow();
4256  for(unsigned i=0;i<n_dof;i++)
4257  {
4258  for(unsigned k=0;k<n_dof;k++)
4259  {
4260  //Cache the value of the hessian y product
4261  const double j_y = jac_y(i,k);
4262  //Loop over the possible vectors
4263  for(unsigned v=0;v<n_vec;v++)
4264  {
4265  product(v,i) += j_y*C(v,k);
4266  }
4267  }
4268  }
4269 }
4270 
4271 
4272 //////////////////////////////////////////////////////////////////////
4273 //////////////////////////////////////////////////////////////////////
4274 //////////////////////////////////////////////////////////////////////
4275 
4276 //=============================================================================
4277 /// Create a list of pairs for all unknowns in this element,
4278 /// so that the first entry in each pair contains the global equation
4279 /// number of the unknown, while the second one contains the number
4280 /// of the "DOF type" that this unknown is associated with.
4281 /// (Function can obviously only be called if the equation numbering
4282 /// scheme has been set up.)
4283 //=============================================================================
4286  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
4287 {
4288  // number of nodes
4289  unsigned n_node = this->nnode();
4290 
4291  // number of pressure values
4292  unsigned n_press = this->npres_axi_nst();
4293 
4294  // temporary pair (used to store dof lookup prior to being added to list)
4295  std::pair<unsigned,unsigned> dof_lookup;
4296 
4297  // pressure dof number (is this really OK?)
4298  unsigned pressure_dof_number = 3;
4299 
4300  // loop over the pressure values
4301  for (unsigned n = 0; n < n_press; n++)
4302  {
4303  // determine local eqn number
4304  int local_eqn_number = this->p_local_eqn(n);
4305 
4306  // ignore pinned values - far away degrees of freedom resulting
4307  // from hanging nodes can be ignored since these are be dealt
4308  // with by the element containing their master nodes
4309  if (local_eqn_number >= 0)
4310  {
4311  // store dof lookup in temporary pair: First entry in pair
4312  // is global equation number; second entry is dof type
4313  dof_lookup.first = this->eqn_number(local_eqn_number);
4314  dof_lookup.second = pressure_dof_number;
4315 
4316  // add to list
4317  dof_lookup_list.push_front(dof_lookup);
4318  }
4319  }
4320 
4321  // loop over the nodes
4322  for (unsigned n = 0; n < n_node; n++)
4323  {
4324  // find the number of values at this node
4325  unsigned nv = this->node_pt(n)->nvalue();
4326 
4327  //loop over these values
4328  for (unsigned v = 0; v < nv; v++)
4329  {
4330  // determine local eqn number
4331  int local_eqn_number = this->nodal_local_eqn(n, v);
4332 
4333  // ignore pinned values
4334  if (local_eqn_number >= 0)
4335  {
4336  // store dof lookup in temporary pair: First entry in pair
4337  // is global equation number; second entry is dof type
4338  dof_lookup.first = this->eqn_number(local_eqn_number);
4339  dof_lookup.second = v;
4340 
4341  // add to list
4342  dof_lookup_list.push_front(dof_lookup);
4343 
4344  }
4345  }
4346  }
4347 }
4348 
4349 
4350 
4351 ///GeneralisedNewtonianAxisymmetric Crouzeix-Raviart elements
4352 //Set the data for the number of Variables at each node
4354 Initial_Nvalue[9]={3,3,3,3,3,3,3,3,3};
4355 
4356 //========================================================================
4357 /// Number of values (pinned or dofs) required at node n.
4358 //========================================================================
4360 required_nvalue(const unsigned &n) const {return Initial_Nvalue[n];}
4361 
4362 //========================================================================
4363 /// Compute traction at local coordinate s for outer unit normal N
4364 //========================================================================
4367  Vector<double>& traction)
4368  const
4369 {
4371 }
4372 
4373 ///////////////////////////////////////////////////////////////////////////
4374 ///////////////////////////////////////////////////////////////////////////
4375 ///////////////////////////////////////////////////////////////////////////
4376 
4377 
4378 
4379 //============================================================================
4380 /// Create a list of pairs for all unknowns in this element,
4381 /// so the first entry in each pair contains the global equation
4382 /// number of the unknown, while the second one contains the number
4383 /// of the "DOF type" that this unknown is associated with.
4384 /// (Function can obviously only be called if the equation numbering
4385 /// scheme has been set up.)
4386 //============================================================================
4389  std::list<std::pair<unsigned long,
4390  unsigned> >& dof_lookup_list) const
4391 {
4392  // number of nodes
4393  unsigned n_node = this->nnode();
4394 
4395  // temporary pair (used to store dof lookup prior to being added to list)
4396  std::pair<unsigned,unsigned> dof_lookup;
4397 
4398  // loop over the nodes
4399  for (unsigned n = 0; n < n_node; n++)
4400  {
4401  // find the number of values at this node
4402  unsigned nv = this->node_pt(n)->nvalue();
4403 
4404  //loop over these values
4405  for (unsigned v = 0; v < nv; v++)
4406  {
4407  // determine local eqn number
4408  int local_eqn_number = this->nodal_local_eqn(n, v);
4409 
4410  // ignore pinned values - far away degrees of freedom resulting
4411  // from hanging nodes can be ignored since these are be dealt
4412  // with by the element containing their master nodes
4413  if (local_eqn_number >= 0)
4414  {
4415  // store dof lookup in temporary pair: Global equation number
4416  // is the first entry in pair
4417  dof_lookup.first = this->eqn_number(local_eqn_number);
4418 
4419  // set dof numbers: Dof number is the second entry in pair
4420  dof_lookup.second = v;
4421 
4422  // add to list
4423  dof_lookup_list.push_front(dof_lookup);
4424  }
4425  }
4426  }
4427 }
4428 
4429 
4430 
4431 //GeneralisedNewtonianAxisymmetric Taylor--Hood
4432 //Set the data for the number of Variables at each node
4434 Initial_Nvalue[9]={4,3,4,3,3,3,4,3,4};
4435 
4436 //Set the data for the pressure conversion array
4438 Pconv[4]={0,2,6,8};
4439 
4440 
4441 //========================================================================
4442 /// Compute traction at local coordinate s for outer unit normal N
4443 //========================================================================
4446  Vector<double>& traction)
4447  const
4448 {
4450 }
4451 
4452 }
int p_local_eqn(const unsigned &n) const
Overload the access function for the pressure's local equation numbers.
virtual void extrapolated_strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Extrapolated strain-rate tensor: where (in that order) based on the previous time steps evaluated a...
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
virtual void fill_in_generic_dresidual_contribution_axi_nst(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, unsigned flag)
Compute the derivative of residuals for the Navier–Stokes equations; with respect to a parameeter fla...
GeneralisedNewtonianConstitutiveEquation< 3 > * Constitutive_eqn_pt
Pointer to the generalised Newtonian constitutive equation.
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
virtual double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)=0
cstr elem_len * i
Definition: cfortran.h:607
virtual double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point...
Definition: elements.cc:4049
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
virtual unsigned u_index_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.
int local_eqn_number(const unsigned long &ieqn_global) const
Return the local equation number corresponding to the ieqn_global-th global equation number...
Definition: elements.h:731
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
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. This function computes these terms analytically and overwrites the default implementation in the FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}.
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
virtual void fill_in_generic_residual_contribution_axi_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=2 or 1 or 0: compute the Jacobian and/or ...
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:3981
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:557
char t
Definition: cfortran.h:572
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
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
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
Compute traction at local coordinate s for outer unit normal N.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation. NOTE: Moved to cc file because of a possible compiler bug in gcc (yes, really!). The move to the cc file avoids inlining which appears to cause problems (only) when compiled with gcc and -O3; offensive "illegal read" is in optimised-out section of code and data that is allegedly illegal is readily readable (by other means) just before this function is called so I can't really see how we could possibly be responsible for this...
Definition: elements.cc:1657
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...
e
Definition: cfortran.h:575
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
virtual double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and derivatives w.r.t. global coords at ipt-th integration point Return J...
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
A Rank 4 Tensor class.
Definition: matrices.h:1625
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
Compute traction at local coordinate s for outer unit normal N.
const Vector< double > & g() const
Vector of gravitational components.
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:504
virtual void get_body_force_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Calculate the body force fct at a given time and Eulerian position.
bool has_auxiliary_node_update_fct_pt()
Boolean to indicate if node has a pointer to and auxiliary update function.
Definition: nodes.h:1500
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
void traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
Compute traction (on the viscous scale) at local coordinate s for outer unit normal N...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
A Rank 3 Tensor class.
Definition: matrices.h:1337
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
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.
static char t char * s
Definition: cfortran.h:572
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
double & dt(const unsigned &t=0)
Definition: timesteppers.h:137
double get_source_fct(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
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...
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
virtual void get_source_fct_gradient(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
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 shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3158
void interpolated_u_axi_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
virtual double p_axi_nst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes...
double interpolated_p_axi_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
virtual int p_local_eqn(const unsigned &n) const =0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
virtual void pshape_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
double raw_dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n. Do not use the hanging node repre...
Definition: elements.h:2170
virtual void get_viscosity_ratio_axisym_nst(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, double &visc_ratio)
Calculate the viscosity ratio relative to the viscosity used in the definition of the Reynolds number...
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2215
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2446
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...
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed...
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
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...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
virtual void get_body_force_gradient_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
Get gradient of body force term at (Eulerian) position x. Computed via function pointer (if set) or b...
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 std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
virtual unsigned npres_axi_nst() const =0
Function to return number of pressure degrees of freedom.
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.
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
virtual double viscosity(const double &second_invariant_of_rate_of_strain_tensor)=0
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
double du_dt_axi_nst(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
double second_invariant(const DenseMatrix< double > &tensor)
Compute second invariant of tensor.
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164
void max_and_min_invariant_and_viscosity(double &min_invariant, double &max_invariant, double &min_viscosity, double &max_viscosity) const
Get max. and min. strain rate invariant and visocosity over all integration points in element...
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.