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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 // Non-inline functions for NS elements
32 
33 
34 namespace oomph
35 {
36 
37 /// Navier--Stokes equations static data
38 Vector<double> AxisymmetricNavierStokesEquations::Gamma(2,1.0);
39 
40 //=================================================================
41 /// "Magic" negative number that indicates that the pressure is
42 /// not stored at a node. This cannot be -1 because that represents
43 /// the positional hanging scheme in the hanging_pt object of nodes
44 //=================================================================
46 
47 /// Navier--Stokes equations static data
50 
51 // Navier--Stokes equations static data
54 
55 /// Navier-Stokes equations default gravity vector
58 
59 
60 //================================================================
61 /// Compute the diagonal of the velocity/pressure mass matrices.
62 /// If which one=0, both are computed, otherwise only the pressure
63 /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
64 /// LSC version of the preconditioner only needs that one)
65 /// NOTE: pressure versions isn't implemented yet because this
66 /// element has never been tried with Fp preconditoner.
67 //================================================================
70  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
71  const unsigned& which_one)
72  {
73 
74 #ifdef PARANOID
75  if ((which_one==0)||(which_one==1))
76  {
77  throw OomphLibError(
78  "Computation of diagonal of pressure mass matrix is not impmented yet.\n",
79  OOMPH_CURRENT_FUNCTION,
80  OOMPH_EXCEPTION_LOCATION);
81  }
82 #endif
83 
84  // Resize and initialise
85  veloc_mass_diag.assign(ndof(), 0.0);
86 
87  // find out how many nodes there are
88  const unsigned n_node = nnode();
89 
90  // find number of coordinates
91  const unsigned n_value = 3;
92 
93  // find the indices at which the local velocities are stored
94  Vector<unsigned> u_nodal_index(n_value);
95  for(unsigned i=0; i<n_value; i++)
96  {
97  u_nodal_index[i] = this->u_index_axi_nst(i);
98  }
99 
100  //Set up memory for test functions
101  Shape test(n_node);
102 
103  //Number of integration points
104  unsigned n_intpt = integral_pt()->nweight();
105 
106  //Integer to store the local equations no
107  int local_eqn=0;
108 
109  //Loop over the integration points
110  for(unsigned ipt=0; ipt<n_intpt; ipt++)
111  {
112 
113  //Get the integral weight
114  double w = integral_pt()->weight(ipt);
115 
116  //Get determinant of Jacobian of the mapping
117  double J = J_eulerian_at_knot(ipt);
118 
119  //Premultiply weights and Jacobian
120  double W = w*J;
121 
122  //Get the velocity test functions - there is no explicit
123  // function to give the test function so use shape
124  shape_at_knot(ipt,test);
125 
126  //Need to get the position to sort out the jacobian properly
127  double r = 0.0;
128  for(unsigned l=0;l<n_node;l++)
129  {
130  r += this->nodal_position(l,0)*test(l);
131  }
132 
133  //Multiply by the geometric factor
134  W *= r;
135 
136  //Loop over the veclocity test functions
137  for(unsigned l=0; l<n_node; l++)
138  {
139  //Loop over the velocity components
140  for(unsigned i=0; i<n_value; i++)
141  {
142  local_eqn = nodal_local_eqn(l,u_nodal_index[i]);
143 
144  //If not a boundary condition
145  if(local_eqn >= 0)
146  {
147  //Add the contribution
148  veloc_mass_diag[local_eqn] += test[l]*test[l] * W;
149  } //End of if not boundary condition statement
150  } //End of loop over dimension
151  } //End of loop over test functions
152 
153  }
154  }
155 
156 
157 //======================================================================
158 /// Validate against exact velocity solution at given time.
159 /// Solution is provided via function pointer.
160 /// Plot at a given number of plot points and compute L2 error
161 /// and L2 norm of velocity solution over element.
162 //=======================================================================
164 compute_error(std::ostream &outfile,
166  const double& time,
167  double& error, double& norm)
168 {
169  error=0.0;
170  norm=0.0;
171 
172  //Vector of local coordinates
173  Vector<double> s(2);
174 
175  // Vector for coordintes
176  Vector<double> x(2);
177 
178  //Set the value of Nintpt
179  unsigned Nintpt = integral_pt()->nweight();
180 
181  outfile << "ZONE" << std::endl;
182 
183  // Exact solution Vector (u,v,w,p)
184  Vector<double> exact_soln(4);
185 
186  //Loop over the integration points
187  for(unsigned ipt=0;ipt<Nintpt;ipt++)
188  {
189  //Assign values of s
190  for(unsigned i=0;i<2;i++) {s[i] = integral_pt()->knot(ipt,i);}
191 
192  //Get the integral weight
193  double w = integral_pt()->weight(ipt);
194 
195  // Get jacobian of mapping
196  double J = J_eulerian(s);
197 
198  // Get x position as Vector
199  interpolated_x(s,x);
200 
201  //Premultiply the weights and the Jacobian and r
202  double W = w*J*x[0];
203 
204  // Get exact solution at this point
205  (*exact_soln_pt)(time,x,exact_soln);
206 
207  // Velocity error
208  for(unsigned i=0;i<3;i++)
209  {
210  norm+=exact_soln[i]*exact_soln[i]*W;
211  error+=(exact_soln[i]-interpolated_u_axi_nst(s,i))*
212  (exact_soln[i]-interpolated_u_axi_nst(s,i))*W;
213  }
214 
215  //Output x,y,...,u_exact
216  for(unsigned i=0;i<2;i++) {outfile << x[i] << " ";}
217 
218  //Output x,y,[z],u_error,v_error,[w_error]
219  for(unsigned i=0;i<3;i++)
220  {outfile << exact_soln[i]-interpolated_u_axi_nst(s,i) << " ";}
221 
222  outfile << std::endl;
223  }
224 }
225 
226 //======================================================================
227 /// Validate against exact velocity solution
228 /// Solution is provided via function pointer.
229 /// Plot at a given number of plot points and compute L2 error
230 /// and L2 norm of velocity solution over element.
231 //=======================================================================
233 compute_error(std::ostream &outfile,
235  double& error, double& norm)
236 {
237  error=0.0;
238  norm=0.0;
239 
240  //Vector of local coordinates
241  Vector<double> s(2);
242 
243  // Vector for coordintes
244  Vector<double> x(2);
245 
246  //Set the value of Nintpt
247  unsigned Nintpt = integral_pt()->nweight();
248 
249  outfile << "ZONE" << std::endl;
250 
251  // Exact solution Vector (u,v,w,p)
252  Vector<double> exact_soln(4);
253 
254  //Loop over the integration points
255  for(unsigned ipt=0;ipt<Nintpt;ipt++)
256  {
257 
258  //Assign values of s
259  for(unsigned i=0;i<2;i++) {s[i] = integral_pt()->knot(ipt,i);}
260 
261  //Get the integral weight
262  double w = integral_pt()->weight(ipt);
263 
264  // Get jacobian of mapping
265  double J=J_eulerian(s);
266 
267  // Get x position as Vector
268  interpolated_x(s,x);
269 
270  //Premultiply the weights and the Jacobian and r
271  double W = w*J*x[0];
272 
273  // Get exact solution at this point
274  (*exact_soln_pt)(x,exact_soln);
275 
276  // Velocity error
277  for(unsigned i=0;i<3;i++)
278  {
279  norm+=exact_soln[i]*exact_soln[i]*W;
280  error+=(exact_soln[i]-interpolated_u_axi_nst(s,i))*
281  (exact_soln[i]-interpolated_u_axi_nst(s,i))*W;
282  }
283 
284  //Output x,y,...,u_exact
285  for(unsigned i=0;i<2;i++) {outfile << x[i] << " ";}
286 
287  //Output x,y,u_error,v_error,w_error
288  for(unsigned i=0;i<3;i++)
289  {outfile << exact_soln[i]-interpolated_u_axi_nst(s,i) << " ";}
290 
291  outfile << std::endl;
292  }
293 }
294 
295 //======================================================================
296 /// Output "exact" solution
297 /// Solution is provided via function pointer.
298 /// Plot at a given number of plot points.
299 /// Function prints as many components as are returned in solution Vector.
300 //=======================================================================
302 output_fct(std::ostream &outfile,
303  const unsigned &nplot,
305 {
306 
307  //Vector of local coordinates
308  Vector<double> s(2);
309 
310  // Vector for coordintes
311  Vector<double> x(2);
312 
313  // Tecplot header info
314  outfile << tecplot_zone_string(nplot);
315 
316  // Exact solution Vector
317  Vector<double> exact_soln;
318 
319  // Loop over plot points
320  unsigned num_plot_points=nplot_points(nplot);
321  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
322  {
323 
324  // Get local coordinates of plot point
325  get_s_plot(iplot,nplot,s);
326 
327  // Get x position as Vector
328  interpolated_x(s,x);
329 
330  // Get exact solution at this point
331  (*exact_soln_pt)(x,exact_soln);
332 
333  //Output x,y,...
334  for(unsigned i=0;i<2;i++)
335  {
336  outfile << x[i] << " ";
337  }
338 
339  //Output "exact solution"
340  for(unsigned i=0;i<exact_soln.size();i++)
341  {
342  outfile << exact_soln[i] << " ";
343  }
344 
345  outfile << std::endl;
346 
347  }
348 
349  // Write tecplot footer (e.g. FE connectivity lists)
350  write_tecplot_zone_footer(outfile,nplot);
351 
352 }
353 
354 //======================================================================
355 /// Output "exact" solution at a given time
356 /// Solution is provided via function pointer.
357 /// Plot at a given number of plot points.
358 /// Function prints as many components as are returned in solution Vector.
359 //=======================================================================
361 output_fct(std::ostream &outfile,
362  const unsigned &nplot,
363  const double& time,
365 {
366 
367  //Vector of local coordinates
368  Vector<double> s(2);
369 
370  // Vector for coordintes
371  Vector<double> x(2);
372 
373  // Tecplot header info
374  outfile << tecplot_zone_string(nplot);
375 
376  // Exact solution Vector
377  Vector<double> exact_soln;
378 
379  // Loop over plot points
380  unsigned num_plot_points=nplot_points(nplot);
381  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
382  {
383 
384  // Get local coordinates of plot point
385  get_s_plot(iplot,nplot,s);
386 
387  // Get x position as Vector
388  interpolated_x(s,x);
389 
390  // Get exact solution at this point
391  (*exact_soln_pt)(time,x,exact_soln);
392 
393  //Output x,y,...
394  for(unsigned i=0;i<2;i++)
395  {
396  outfile << x[i] << " ";
397  }
398 
399  //Output "exact solution"
400  for(unsigned i=0;i<exact_soln.size();i++)
401  {
402  outfile << exact_soln[i] << " ";
403  }
404 
405  outfile << std::endl;
406 
407  }
408 
409  // Write tecplot footer (e.g. FE connectivity lists)
410  write_tecplot_zone_footer(outfile,nplot);
411 
412 }
413 
414 //==============================================================
415 /// Output function: Velocities only
416 /// x,y,[z],u,v,[w]
417 /// in tecplot format at specified previous timestep (t=0: present;
418 /// t>0: previous timestep). Specified number of plot points in each
419 /// coordinate direction.
420 //==============================================================
422 output_veloc(std::ostream &outfile,
423  const unsigned &nplot,
424  const unsigned &t)
425 {
426  //Find number of nodes
427  unsigned n_node = nnode();
428 
429  //Local shape function
430  Shape psi(n_node);
431 
432  //Vectors of local coordinates and coords and velocities
433  Vector<double> s(2);
435  Vector<double> interpolated_u(3);
436 
437 
438  // Tecplot header info
439  outfile << tecplot_zone_string(nplot);
440 
441  // Loop over plot points
442  unsigned num_plot_points=nplot_points(nplot);
443  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
444  {
445 
446  // Get local coordinates of plot point
447  get_s_plot(iplot,nplot,s);
448 
449  // Get shape functions
450  shape(s,psi);
451 
452  // Loop over coordinate directions
453  for(unsigned i=0;i<2;i++)
454  {
455  interpolated_x[i]=0.0;
456  //Loop over the local nodes and sum
457  for(unsigned l=0;l<n_node;l++)
458  {interpolated_x[i] += nodal_position(t,l,i)*psi[l];}
459  }
460 
461  //Loop over the velocity components
462  for(unsigned i=0;i<3;i++)
463  {
464  //Get the index at which the velocity is stored
465  unsigned u_nodal_index = u_index_axi_nst(i);
466  interpolated_u[i]=0.0;
467  //Loop over the local nodes and sum
468  for(unsigned l=0;l<n_node;l++)
469  {interpolated_u[i] += nodal_value(t,l,u_nodal_index)*psi[l];}
470  }
471 
472  // Coordinates
473  for(unsigned i=0;i<2;i++) {outfile << interpolated_x[i] << " ";}
474 
475  // Velocities
476  for(unsigned i=0;i<3;i++) {outfile << interpolated_u[i] << " ";}
477 
478  outfile << std::endl;
479  }
480 
481  // Write tecplot footer (e.g. FE connectivity lists)
482  write_tecplot_zone_footer(outfile,nplot);
483 
484 }
485 
486 //==============================================================
487 /// Output function:
488 /// r,z,u,v,w,p
489 /// in tecplot format. Specified number of plot points in each
490 /// coordinate direction.
491 //==============================================================
492 void AxisymmetricNavierStokesEquations::output(std::ostream &outfile,
493  const unsigned &nplot)
494 {
495  //Vector of local coordinates
496  Vector<double> s(2);
497 
498  // Tecplot header info
499  outfile << tecplot_zone_string(nplot);
500 
501  // Loop over plot points
502  unsigned num_plot_points=nplot_points(nplot);
503  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
504  {
505 
506  // Get local coordinates of plot point
507  get_s_plot(iplot,nplot,s);
508 
509  // Coordinates
510  for(unsigned i=0;i<2;i++) {outfile << interpolated_x(s,i) << " ";}
511 
512  // Velocities
513  for(unsigned i=0;i<3;i++) {outfile << interpolated_u_axi_nst(s,i) << " ";}
514 
515  // Pressure
516  outfile << interpolated_p_axi_nst(s) << " ";
517 
518  outfile << std::endl;
519  }
520  outfile << std::endl;
521 
522  // Write tecplot footer (e.g. FE connectivity lists)
523  write_tecplot_zone_footer(outfile,nplot);
524 
525 }
526 
527 
528 //==============================================================
529 /// Output function:
530 /// r,z,u,v,w,p
531 /// in tecplot format. Specified number of plot points in each
532 /// coordinate direction.
533 //==============================================================
535  const unsigned &nplot)
536 {
537  //Vector of local coordinates
538  Vector<double> s(2);
539 
540  // Tecplot header info
541  fprintf(file_pt,"%s ",tecplot_zone_string(nplot).c_str());
542 
543  // Loop over plot points
544  unsigned num_plot_points=nplot_points(nplot);
545  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
546  {
547 
548  // Get local coordinates of plot point
549  get_s_plot(iplot,nplot,s);
550 
551  // Coordinates
552  for(unsigned i=0;i<2;i++)
553  {
554  //outfile << interpolated_x(s,i) << " ";
555  fprintf(file_pt,"%g ",interpolated_x(s,i));
556  }
557 
558  // Velocities
559  for(unsigned i=0;i<3;i++)
560  {
561  //outfile << interpolated_u(s,i) << " ";
562  fprintf(file_pt,"%g ",interpolated_u_axi_nst(s,i));
563  }
564 
565  // Pressure
566  //outfile << interpolated_p(s) << " ";
567  fprintf(file_pt,"%g ",interpolated_p_axi_nst(s));
568 
569  //outfile << std::endl;
570  fprintf(file_pt,"\n");
571  }
572  //outfile << std::endl;
573  fprintf(file_pt,"\n");
574 
575  // Write tecplot footer (e.g. FE connectivity lists)
576  write_tecplot_zone_footer(file_pt,nplot);
577 
578 }
579 
580 
581 
582 
583 //==============================================================
584 /// Return integral of dissipation over element
585 //==============================================================
587 {
588  throw OomphLibError(
589  "Check the dissipation calculation for axisymmetric NSt",
590  OOMPH_CURRENT_FUNCTION,
591  OOMPH_EXCEPTION_LOCATION);
592 
593  // Initialise
594  double diss=0.0;
595 
596  //Set the value of Nintpt
597  unsigned Nintpt = integral_pt()->nweight();
598 
599  //Set the Vector to hold local coordinates
600  Vector<double> s(2);
601 
602  //Loop over the integration points
603  for(unsigned ipt=0;ipt<Nintpt;ipt++)
604  {
605 
606  //Assign values of s
607  for(unsigned i=0;i<2;i++)
608  {
609  s[i] = integral_pt()->knot(ipt,i);
610  }
611 
612  //Get the integral weight
613  double w = integral_pt()->weight(ipt);
614 
615  // Get Jacobian of mapping
616  double J = J_eulerian(s);
617 
618  // Get strain rate matrix
619  DenseMatrix<double> strainrate(3,3);
620  strain_rate(s,strainrate);
621 
622  // Initialise
623  double local_diss=0.0;
624  for(unsigned i=0;i<3;i++)
625  {
626  for(unsigned j=0;j<3;j++)
627  {
628  local_diss+=2.0*strainrate(i,j)*strainrate(i,j);
629  }
630  }
631 
632  diss+=local_diss*w*J;
633  }
634 
635  return diss;
636 
637 }
638 
639 //==============================================================
640 /// \short Compute traction (on the viscous scale) at local
641 /// coordinate s for outer unit normal N
642 //==============================================================
644  const Vector<double>& N,
645  Vector<double>& traction)
646  const
647 {
648  //throw OomphLibError(
649  // "Check the traction calculation for axisymmetric NSt",
650  // OOMPH_CURRENT_FUNCTION,
651  // OOMPH_EXCEPTION_LOCATION);
652 
653  // Pad out normal vector if required
654  Vector<double> n_local(3,0.0);
655  n_local[0]=N[0];
656  n_local[1]=N[1];
657 
658 #ifdef PARANOID
659  if ((N.size()==3)&&(N[2]!=0.0))
660  {
661  throw OomphLibError(
662  "Unit normal passed into this fct should either be 2D (r,z) or have a zero component in the theta-direction",
663  OOMPH_CURRENT_FUNCTION,
664  OOMPH_EXCEPTION_LOCATION);
665  }
666 #endif
667 
668  // Get velocity gradients
669  DenseMatrix<double> strainrate(3,3);
670  strain_rate(s,strainrate);
671 
672  // Get pressure
673  double press=interpolated_p_axi_nst(s);
674 
675  // Loop over traction components
676  for (unsigned i=0;i<3;i++)
677  {
678  traction[i]=-press*n_local[i];
679  for (unsigned j=0;j<3;j++)
680  {
681  traction[i]+=2.0*strainrate(i,j)*n_local[j];
682  }
683  }
684 }
685 
686 //==============================================================
687 /// Return dissipation at local coordinate s
688 //==============================================================
691 {
692  throw OomphLibError(
693  "Check the dissipation calculation for axisymmetric NSt",
694 OOMPH_CURRENT_FUNCTION,
695  OOMPH_EXCEPTION_LOCATION);
696 
697  // Get strain rate matrix
698  DenseMatrix<double> strainrate(3,3);
699  strain_rate(s,strainrate);
700 
701  // Initialise
702  double local_diss=0.0;
703  for(unsigned i=0;i<3;i++)
704  {
705  for(unsigned j=0;j<3;j++)
706  {
707  local_diss+=2.0*strainrate(i,j)*strainrate(i,j);
708  }
709  }
710 
711  return local_diss;
712 }
713 
714 //==============================================================
715 /// Get strain-rate tensor: \f$ e_{ij} \f$ where
716 /// \f$ i,j = r,z,\theta \f$ (in that order)
717 //==============================================================
720 {
721 
722 #ifdef PARANOID
723  if ((strainrate.ncol()!=3)||(strainrate.nrow()!=3))
724  {
725  std::ostringstream error_message;
726  error_message << "The strain rate has incorrect dimensions "
727  << strainrate.ncol() << " x "
728  << strainrate.nrow() << " Not 3" << std::endl;
729 
730  throw OomphLibError(error_message.str(),
731 OOMPH_CURRENT_FUNCTION,
732  OOMPH_EXCEPTION_LOCATION);
733  }
734 #endif
735 
736  //Find out how many nodes there are in the element
737  unsigned n_node = nnode();
738 
739  //Set up memory for the shape and test functions
740  Shape psi(n_node);
741  DShape dpsidx(n_node,2);
742 
743  //Call the derivatives of the shape functions
744  dshape_eulerian(s,psi,dpsidx);
745 
746  // Radius
747  double interpolated_r = 0.0;
748 
749  // Velocity components and their derivatives
750  double ur=0.0;
751  double durdr=0.0;
752  double durdz=0.0;
753  double uz=0.0;
754  double duzdr=0.0;
755  double duzdz=0.0;
756  double uphi=0.0;
757  double duphidr=0.0;
758  double duphidz=0.0;
759 
760  //Get the local storage for the indices
761  unsigned u_nodal_index[3];
762  for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}
763 
764  // Loop over nodes to assemble velocities and their derivatives
765  // w.r.t. to r and z (x_0 and x_1)
766  for(unsigned l=0;l<n_node;l++)
767  {
768  interpolated_r += nodal_position(l,0)*psi[l];
769 
770  ur += nodal_value(l,u_nodal_index[0])*psi[l];
771  uz += nodal_value(l,u_nodal_index[1])*psi[l];
772  uphi += nodal_value(l,u_nodal_index[2])*psi[l];
773 
774  durdr += nodal_value(l,u_nodal_index[0])*dpsidx(l,0);
775  durdz += nodal_value(l,u_nodal_index[0])*dpsidx(l,1);
776 
777  duzdr += nodal_value(l,u_nodal_index[1])*dpsidx(l,0);
778  duzdz += nodal_value(l,u_nodal_index[1])*dpsidx(l,1);
779 
780  duphidr += nodal_value(l,u_nodal_index[2])*dpsidx(l,0);
781  duphidz += nodal_value(l,u_nodal_index[2])*dpsidx(l,1);
782  }
783 
784 
785  // Assign strain rates without negative powers of the radius
786  // and zero those with:
787  strainrate(0,0)=durdr;
788  strainrate(0,1)=0.5*(durdz+duzdr);
789  strainrate(1,0)=strainrate(0,1);
790  strainrate(0,2)=0.0;
791  strainrate(2,0)=strainrate(0,2);
792  strainrate(1,1)=duzdz;
793  strainrate(1,2)=0.5*duphidz;
794  strainrate(2,1)=strainrate(1,2);
795  strainrate(2,2)=0.0;
796 
797 
798  // Overwrite the strain rates with negative powers of the radius
799  // unless we're at the origin
800  if (std::fabs(interpolated_r)>1.0e-16)
801  {
802  double inverse_radius=1.0/interpolated_r;
803  strainrate(0,2)=0.5*(duphidr-inverse_radius*uphi);
804  strainrate(2,0)=strainrate(0,2);
805  strainrate(2,2)=inverse_radius*ur;
806  }
807 
808 
809 }
810 
811 
812 //==============================================================
813 /// \short Get integral of kinetic energy over element:
814 //==============================================================
816 {
817 
818  throw OomphLibError(
819  "Check the kinetic energy calculation for axisymmetric NSt",
820 OOMPH_CURRENT_FUNCTION,
821  OOMPH_EXCEPTION_LOCATION);
822 
823  // Initialise
824  double kin_en=0.0;
825 
826  //Set the value of Nintpt
827  unsigned Nintpt = integral_pt()->nweight();
828 
829  //Set the Vector to hold local coordinates
830  Vector<double> s(2);
831 
832  //Loop over the integration points
833  for(unsigned ipt=0;ipt<Nintpt;ipt++)
834  {
835  //Assign values of s
836  for(unsigned i=0;i<2;i++) {s[i] = integral_pt()->knot(ipt,i);}
837 
838  //Get the integral weight
839  double w = integral_pt()->weight(ipt);
840 
841  //Get Jacobian of mapping
842  double J = J_eulerian(s);
843 
844  // Loop over directions
845  double veloc_squared=0.0;
846  for(unsigned i=0;i<3;i++)
847  {
848  veloc_squared+=interpolated_u_axi_nst(s,i)*interpolated_u_axi_nst(s,i);
849  }
850 
851  kin_en+=0.5*veloc_squared*w*J*interpolated_x(s,0);
852  }
853 
854  return kin_en;
855 
856 }
857 
858 //==============================================================
859 /// Return pressure integrated over the element
860 //==============================================================
862 {
863 
864  // Initialise
865  double press_int=0;
866 
867  //Set the value of Nintpt
868  unsigned Nintpt = integral_pt()->nweight();
869 
870  //Set the Vector to hold local coordinates
871  Vector<double> s(2);
872 
873  //Loop over the integration points
874  for(unsigned ipt=0;ipt<Nintpt;ipt++)
875  {
876 
877  //Assign values of s
878  for(unsigned i=0;i<2;i++)
879  {
880  s[i] = integral_pt()->knot(ipt,i);
881  }
882 
883  //Get the integral weight
884  double w = integral_pt()->weight(ipt);
885 
886  //Get Jacobian of mapping
887  double J = J_eulerian(s);
888 
889  //Premultiply the weights and the Jacobian
890  double W = w*J*interpolated_x(s,0);
891 
892  // Get pressure
893  double press=interpolated_p_axi_nst(s);
894 
895  // Add
896  press_int+=press*W;
897 
898  }
899 
900  return press_int;
901 
902 }
903 
904 //==============================================================
905 /// Compute the residuals for the Navier--Stokes
906 /// equations; flag=1(or 0): do (or don't) compute the
907 /// Jacobian as well.
908 //==============================================================
911  DenseMatrix<double> &jacobian,
912  DenseMatrix<double> &mass_matrix,
913  unsigned flag)
914 {
915  //Find out how many nodes there are
916  unsigned n_node = nnode();
917 
918  // Get continuous time from timestepper of first node
919  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
920 
921  //Find out how many pressure dofs there are
922  unsigned n_pres = npres_axi_nst();
923 
924  //Get the nodal indices at which the velocity is stored
925  unsigned u_nodal_index[3];
926  for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}
927 
928  //Set up memory for the shape and test functions
929  //Note that there are only two dimensions, r and z in this problem
930  Shape psif(n_node), testf(n_node);
931  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
932 
933  //Set up memory for pressure shape and test functions
934  Shape psip(n_pres), testp(n_pres);
935 
936  //Number of integration points
937  unsigned Nintpt = integral_pt()->nweight();
938 
939  //Set the Vector to hold local coordinates (two dimensions)
940  Vector<double> s(2);
941 
942  //Get Physical Variables from Element
943  //Reynolds number must be multiplied by the density ratio
944  double scaled_re = re()*density_ratio();
945  double scaled_re_st = re_st()*density_ratio();
946  double scaled_re_inv_fr = re_invfr()*density_ratio();
947  double scaled_re_inv_ro = re_invro()*density_ratio();
948  // double visc_ratio = viscosity_ratio();
949  Vector<double> G = g();
950 
951  //Integers used to store the local equation and unknown numbers
952  int local_eqn=0, local_unknown=0;
953 
954  //Loop over the integration points
955  for(unsigned ipt=0;ipt<Nintpt;ipt++)
956  {
957  //Assign values of s
958  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
959  //Get the integral weight
960  double w = integral_pt()->weight(ipt);
961 
962  //Call the derivatives of the shape and test functions
963  double J =
964  dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,psif,dpsifdx,testf,dtestfdx);
965 
966  //Call the pressure shape and test functions
967  pshape_axi_nst(s,psip,testp);
968 
969  //Premultiply the weights and the Jacobian
970  double W = w*J;
971 
972  //Allocate storage for the position and the derivative of the
973  //mesh positions wrt time
975  Vector<double> mesh_velocity(2,0.0);
976  //Allocate storage for the pressure, velocity components and their
977  //time and spatial derivatives
978  double interpolated_p=0.0;
979  Vector<double> interpolated_u(3,0.0);
980  Vector<double> dudt(3,0.0);
981  DenseMatrix<double> interpolated_dudx(3,2,0.0);
982 
983  //Calculate pressure at integration point
984  for(unsigned l=0;l<n_pres;l++) {interpolated_p += p_axi_nst(l)*psip[l];}
985 
986  //Calculate velocities and derivatives at integration point
987 
988  // Loop over nodes
989  for(unsigned l=0;l<n_node;l++)
990  {
991  //Cache the shape function
992  const double psif_ = psif(l);
993  //Loop over the two coordinate directions
994  for(unsigned i=0;i<2;i++)
995  {
996  interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
997  }
998  //mesh_velocity[i] += dnodal_position_dt(l,i)*psif[l];
999 
1000  //Loop over the three velocity directions
1001  for(unsigned i=0;i<3;i++)
1002  {
1003  //Get the u_value
1004  const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
1005  interpolated_u[i] += u_value*psif_;
1006  dudt[i]+= du_dt_axi_nst(l,i)*psif_;
1007  //Loop over derivative directions
1008  for(unsigned j=0;j<2;j++)
1009  {interpolated_dudx(i,j) += u_value*dpsifdx(l,j);}
1010  }
1011  }
1012 
1013  //Get the mesh velocity if ALE is enabled
1014  if(!ALE_is_disabled)
1015  {
1016  // Loop over nodes
1017  for(unsigned l=0;l<n_node;l++)
1018  {
1019  //Loop over the two coordinate directions
1020  for(unsigned i=0;i<2;i++)
1021  {
1022  mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif(l);
1023  }
1024  }
1025  }
1026 
1027 
1028  //Get the user-defined body force terms
1029  Vector<double> body_force(3);
1030  get_body_force_axi_nst(time,ipt,s,interpolated_x,body_force);
1031 
1032  //Get the user-defined source function
1033  double source = get_source_fct(time,ipt,interpolated_x);
1034 
1035  //Get the user-defined viscosity function
1036  double visc_ratio;
1038  s,
1039  interpolated_x,
1040  visc_ratio);
1041 
1042  //r is the first position component
1043  double r = interpolated_x[0];
1044 
1045  //MOMENTUM EQUATIONS
1046  //------------------
1047 
1048  //Loop over the test functions
1049  for(unsigned l=0;l<n_node;l++)
1050  {
1051  //FIRST (RADIAL) MOMENTUM EQUATION
1052  local_eqn = nodal_local_eqn(l,u_nodal_index[0]);
1053  //If it's not a boundary condition
1054  if(local_eqn >= 0)
1055  {
1056  //Add the user-defined body force terms
1057  residuals[local_eqn] +=
1058  r*body_force[0]*testf[l]*W;
1059 
1060  //Add the gravitational body force term
1061  residuals[local_eqn] += r*scaled_re_inv_fr*testf[l]*G[0]*W;
1062 
1063  //Add the pressure gradient term
1064  residuals[local_eqn] +=
1065  interpolated_p*(testf[l] + r*dtestfdx(l,0))*W;
1066 
1067  //Add in the stress tensor terms
1068  //The viscosity ratio needs to go in here to ensure
1069  //continuity of normal stress is satisfied even in flows
1070  //with zero pressure gradient!
1071  residuals[local_eqn] -= visc_ratio*
1072  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
1073 
1074  residuals[local_eqn] -= visc_ratio*r*
1075  (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
1076  dtestfdx(l,1)*W;
1077 
1078  residuals[local_eqn] -=
1079  visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
1080 
1081  //Add in the inertial terms
1082  //du/dt term
1083  residuals[local_eqn] -= scaled_re_st*r*dudt[0]*testf[l]*W;
1084 
1085  //Convective terms
1086  residuals[local_eqn] -=
1087  scaled_re*(r*interpolated_u[0]*interpolated_dudx(0,0)
1088  - interpolated_u[2]*interpolated_u[2]
1089  + r*interpolated_u[1]*interpolated_dudx(0,1))*testf[l]*W;
1090 
1091  //Mesh velocity terms
1092  if(!ALE_is_disabled)
1093  {
1094  for(unsigned k=0;k<2;k++)
1095  {
1096  residuals[local_eqn] +=
1097  scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(0,k)*testf[l]*W;
1098  }
1099  }
1100 
1101  //Add the Coriolis term
1102  residuals[local_eqn] +=
1103  2.0*r*scaled_re_inv_ro*interpolated_u[2]*testf[l]*W;
1104 
1105  //CALCULATE THE JACOBIAN
1106  if(flag)
1107  {
1108  //Loop over the velocity shape functions again
1109  for(unsigned l2=0;l2<n_node;l2++)
1110  {
1111  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
1112  //Radial velocity component
1113  if(local_unknown >= 0)
1114  {
1115  if(flag==2)
1116  {
1117  //Add the mass matrix
1118  mass_matrix(local_eqn,local_unknown) +=
1119  scaled_re_st*r*psif[l2]*testf[l]*W;
1120  }
1121 
1122  //Add contribution to the Jacobian matrix
1123  jacobian(local_eqn,local_unknown)
1124  -= visc_ratio*r*(1.0+Gamma[0])
1125  *dpsifdx(l2,0)*dtestfdx(l,0)*W;
1126 
1127  jacobian(local_eqn,local_unknown)
1128  -= visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
1129 
1130  jacobian(local_eqn,local_unknown)
1131  -= visc_ratio*(1.0 + Gamma[0])*psif[l2]*testf[l]*W/r;
1132 
1133  //Add in the inertial terms
1134  //du/dt term
1135  jacobian(local_eqn,local_unknown)
1136  -= scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1137  psif[l2]*testf[l]*W;
1138 
1139  //Convective terms
1140  jacobian(local_eqn,local_unknown) -=
1141  scaled_re*(r*psif[l2]*interpolated_dudx(0,0)
1142  + r*interpolated_u[0]*dpsifdx(l2,0)
1143  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
1144 
1145  //Mesh velocity terms
1146  if(!ALE_is_disabled)
1147  {
1148  for(unsigned k=0;k<2;k++)
1149  {
1150  jacobian(local_eqn,local_unknown) +=
1151  scaled_re_st*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
1152  }
1153  }
1154  }
1155 
1156 
1157  //Axial velocity component
1158  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
1159  if(local_unknown >= 0)
1160  {
1161  jacobian(local_eqn,local_unknown) -=
1162  visc_ratio*r*Gamma[0]*dpsifdx(l2,0)*dtestfdx(l,1)*W;
1163 
1164  //Convective terms
1165  jacobian(local_eqn,local_unknown) -=
1166  scaled_re*r*psif[l2]*interpolated_dudx(0,1)*testf[l]*W;
1167  }
1168 
1169  //Azimuthal velocity component
1170  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
1171  if(local_unknown >= 0)
1172  {
1173  //Convective terms
1174  jacobian(local_eqn,local_unknown) -=
1175  - scaled_re*2.0*interpolated_u[2]*psif[l2]*testf[l]*W;
1176 
1177  //Coriolis terms
1178  jacobian(local_eqn,local_unknown) +=
1179  2.0*r*scaled_re_inv_ro*psif[l2]*testf[l]*W;
1180  }
1181  }
1182 
1183  /*Now loop over pressure shape functions*/
1184  /*This is the contribution from pressure gradient*/
1185  for(unsigned l2=0;l2<n_pres;l2++)
1186  {
1187  local_unknown = p_local_eqn(l2);
1188  /*If we are at a non-zero degree of freedom in the entry*/
1189  if(local_unknown >= 0)
1190  {
1191  jacobian(local_eqn,local_unknown)
1192  += psip[l2]*(testf[l] + r*dtestfdx(l,0))*W;
1193  }
1194  }
1195  } /*End of Jacobian calculation*/
1196 
1197  } //End of if not boundary condition statement
1198 
1199  //SECOND (AXIAL) MOMENTUM EQUATION
1200  local_eqn = nodal_local_eqn(l,u_nodal_index[1]);
1201  //If it's not a boundary condition
1202  if(local_eqn >= 0)
1203  {
1204  //Add the user-defined body force terms
1205  //Remember to multiply by the density ratio!
1206  residuals[local_eqn] +=
1207  r*body_force[1]*testf[l]*W;
1208 
1209  //Add the gravitational body force term
1210  residuals[local_eqn] += r*scaled_re_inv_fr*testf[l]*G[1]*W;
1211 
1212  //Add the pressure gradient term
1213  residuals[local_eqn] += r*interpolated_p*dtestfdx(l,1)*W;
1214 
1215  //Add in the stress tensor terms
1216  //The viscosity ratio needs to go in here to ensure
1217  //continuity of normal stress is satisfied even in flows
1218  //with zero pressure gradient!
1219  residuals[local_eqn] -= visc_ratio*
1220  r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
1221  *dtestfdx(l,0)*W;
1222 
1223  residuals[local_eqn] -= visc_ratio*r*
1224  (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
1225 
1226  //Add in the inertial terms
1227  //du/dt term
1228  residuals[local_eqn] -= scaled_re_st*r*dudt[1]*testf[l]*W;
1229 
1230  //Convective terms
1231  residuals[local_eqn] -=
1232  scaled_re*(r*interpolated_u[0]*interpolated_dudx(1,0)
1233  + r*interpolated_u[1]*interpolated_dudx(1,1))*testf[l]*W;
1234 
1235  //Mesh velocity terms
1236  if(!ALE_is_disabled)
1237  {
1238  for(unsigned k=0;k<2;k++)
1239  {
1240  residuals[local_eqn] +=
1241  scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(1,k)*testf[l]*W;
1242  }
1243  }
1244 
1245  //CALCULATE THE JACOBIAN
1246  if(flag)
1247  {
1248  //Loop over the velocity shape functions again
1249  for(unsigned l2=0;l2<n_node;l2++)
1250  {
1251  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
1252  //Radial velocity component
1253  if(local_unknown >= 0)
1254  {
1255  //Add in the stress tensor terms
1256  //The viscosity ratio needs to go in here to ensure
1257  //continuity of normal stress is satisfied even in flows
1258  //with zero pressure gradient!
1259  jacobian(local_eqn,local_unknown) -=
1260  visc_ratio*r*Gamma[1]*dpsifdx(l2,1)*dtestfdx(l,0)*W;
1261 
1262  //Convective terms
1263  jacobian(local_eqn,local_unknown) -=
1264  scaled_re*r*psif[l2]*interpolated_dudx(1,0)*testf[l]*W;
1265  }
1266 
1267  //Axial velocity component
1268  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
1269  if(local_unknown >= 0)
1270  {
1271  if(flag==2)
1272  {
1273  //Add the mass matrix
1274  mass_matrix(local_eqn,local_unknown) +=
1275  scaled_re_st*r*psif[l2]*testf[l]*W;
1276  }
1277 
1278 
1279  //Add in the stress tensor terms
1280  //The viscosity ratio needs to go in here to ensure
1281  //continuity of normal stress is satisfied even in flows
1282  //with zero pressure gradient!
1283  jacobian(local_eqn,local_unknown) -=
1284  visc_ratio*r*dpsifdx(l2,0)*dtestfdx(l,0)*W;
1285 
1286  jacobian(local_eqn,local_unknown) -=
1287  visc_ratio*r*(1.0 + Gamma[1])*dpsifdx(l2,1)*
1288  dtestfdx(l,1)*W;
1289 
1290  //Add in the inertial terms
1291  //du/dt term
1292  jacobian(local_eqn,local_unknown) -=
1293  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1294  psif[l2]*testf[l]*W;
1295 
1296  //Convective terms
1297  jacobian(local_eqn,local_unknown) -=
1298  scaled_re*(r*interpolated_u[0]*dpsifdx(l2,0)
1299  + r*psif[l2]*interpolated_dudx(1,1)
1300  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
1301 
1302 
1303  //Mesh velocity terms
1304  if(!ALE_is_disabled)
1305  {
1306  for(unsigned k=0;k<2;k++)
1307  {
1308  jacobian(local_eqn,local_unknown) +=
1309  scaled_re_st*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
1310  }
1311  }
1312  }
1313 
1314  //There are no azimithal terms in the axial momentum equation
1315  } //End of loop over velocity shape functions
1316 
1317  /*Now loop over pressure shape functions*/
1318  /*This is the contribution from pressure gradient*/
1319  for(unsigned l2=0;l2<n_pres;l2++)
1320  {
1321  local_unknown = p_local_eqn(l2);
1322  /*If we are at a non-zero degree of freedom in the entry*/
1323  if(local_unknown >= 0)
1324  {
1325  jacobian(local_eqn,local_unknown)
1326  += r*psip[l2]*dtestfdx(l,1)*W;
1327  }
1328  }
1329  } /*End of Jacobian calculation*/
1330 
1331  } //End of AXIAL MOMENTUM EQUATION
1332 
1333  //THIRD (AZIMUTHAL) MOMENTUM EQUATION
1334  local_eqn = nodal_local_eqn(l,u_nodal_index[2]);
1335  if(local_eqn >= 0)
1336  {
1337  //Add the user-defined body force terms
1338  //Remember to multiply by the density ratio!
1339  residuals[local_eqn] +=
1340  r*body_force[2]*testf[l]*W;
1341 
1342  //Add the gravitational body force term
1343  residuals[local_eqn] += r*scaled_re_inv_fr*testf[l]*G[2]*W;
1344 
1345  //There is NO pressure gradient term
1346 
1347  //Add in the stress tensor terms
1348  //The viscosity ratio needs to go in here to ensure
1349  //continuity of normal stress is satisfied even in flows
1350  //with zero pressure gradient!
1351  residuals[local_eqn] -= visc_ratio*
1352  (r*interpolated_dudx(2,0) -
1353  Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
1354 
1355  residuals[local_eqn] -= visc_ratio*r*
1356  interpolated_dudx(2,1)*dtestfdx(l,1)*W;
1357 
1358  residuals[local_eqn] -= visc_ratio*
1359  ((interpolated_u[2]/r) - Gamma[0]*interpolated_dudx(2,0))*testf[l]*W;
1360 
1361 
1362  //Add in the inertial terms
1363  //du/dt term
1364  residuals[local_eqn] -= scaled_re_st*r*dudt[2]*testf[l]*W;
1365 
1366  //Convective terms
1367  residuals[local_eqn] -=
1368  scaled_re*(r*interpolated_u[0]*interpolated_dudx(2,0)
1369  + interpolated_u[0]*interpolated_u[2]
1370  + r*interpolated_u[1]*interpolated_dudx(2,1))*testf[l]*W;
1371 
1372  //Mesh velocity terms
1373  if(!ALE_is_disabled)
1374  {
1375  for(unsigned k=0;k<2;k++)
1376  {
1377  residuals[local_eqn] +=
1378  scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(2,k)*testf[l]*W;
1379  }
1380  }
1381 
1382  //Add the Coriolis term
1383  residuals[local_eqn] -=
1384  2.0*r*scaled_re_inv_ro*interpolated_u[0]*testf[l]*W;
1385 
1386  //CALCULATE THE JACOBIAN
1387  if(flag)
1388  {
1389  //Loop over the velocity shape functions again
1390  for(unsigned l2=0;l2<n_node;l2++)
1391  {
1392  //Radial velocity component
1393  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
1394  if(local_unknown >= 0)
1395  {
1396  //Convective terms
1397  jacobian(local_eqn,local_unknown) -=
1398  scaled_re*(r*psif[l2]*interpolated_dudx(2,0)
1399  + psif[l2]*interpolated_u[2])*testf[l]*W;
1400 
1401  //Coriolis term
1402  jacobian(local_eqn,local_unknown) -=
1403  2.0*r*scaled_re_inv_ro*psif[l2]*testf[l]*W;
1404  }
1405 
1406  //Axial velocity component
1407  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
1408  if(local_unknown >= 0)
1409  {
1410  //Convective terms
1411  jacobian(local_eqn,local_unknown) -=
1412  scaled_re*r*psif[l2]*interpolated_dudx(2,1)*testf[l]*W;
1413  }
1414 
1415  //Azimuthal velocity component
1416  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
1417  if(local_unknown >= 0)
1418  {
1419  if(flag==2)
1420  {
1421  //Add the mass matrix
1422  mass_matrix(local_eqn,local_unknown) +=
1423  scaled_re_st*r*psif[l2]*testf[l]*W;
1424  }
1425 
1426  //Add in the stress tensor terms
1427  //The viscosity ratio needs to go in here to ensure
1428  //continuity of normal stress is satisfied even in flows
1429  //with zero pressure gradient!
1430  jacobian(local_eqn,local_unknown) -=
1431  visc_ratio*(r*dpsifdx(l2,0) -
1432  Gamma[0]*psif[l2])*dtestfdx(l,0)*W;
1433 
1434  jacobian(local_eqn,local_unknown) -=
1435  visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
1436 
1437  jacobian(local_eqn,local_unknown) -=
1438  visc_ratio*((psif[l2]/r) - Gamma[0]*dpsifdx(l2,0))
1439  *testf[l]*W;
1440 
1441  //Add in the inertial terms
1442  //du/dt term
1443  jacobian(local_eqn,local_unknown) -=
1444  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1445  psif[l2]*testf[l]*W;
1446 
1447  //Convective terms
1448  jacobian(local_eqn,local_unknown) -=
1449  scaled_re*(r*interpolated_u[0]*dpsifdx(l2,0)
1450  + interpolated_u[0]*psif[l2]
1451  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
1452 
1453  //Mesh velocity terms
1454  if(!ALE_is_disabled)
1455  {
1456  for(unsigned k=0;k<2;k++)
1457  {
1458  jacobian(local_eqn,local_unknown) +=
1459  scaled_re_st*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
1460  }
1461  }
1462  }
1463  }
1464 
1465  //There are no pressure terms
1466  } //End of Jacobian
1467 
1468  } //End of AZIMUTHAL EQUATION
1469 
1470  } //End of loop over shape functions
1471 
1472 
1473  //CONTINUITY EQUATION
1474  //-------------------
1475 
1476  //Loop over the shape functions
1477  for(unsigned l=0;l<n_pres;l++)
1478  {
1479  local_eqn = p_local_eqn(l);
1480  //If not a boundary conditions
1481  if(local_eqn >= 0)
1482  {
1483  // Source term
1484  residuals[local_eqn] -= r*source*testp[l]*W;
1485 
1486  //Gradient terms
1487  residuals[local_eqn] +=
1488  (interpolated_u[0] + r*interpolated_dudx(0,0)
1489  + r*interpolated_dudx(1,1))*testp[l]*W;
1490 
1491  /*CALCULATE THE JACOBIAN*/
1492  if(flag)
1493  {
1494  /*Loop over the velocity shape functions*/
1495  for(unsigned l2=0;l2<n_node;l2++)
1496  {
1497  //Radial velocity component
1498  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
1499  if(local_unknown >= 0)
1500  {
1501  jacobian(local_eqn,local_unknown) +=
1502  (psif[l2] + r*dpsifdx(l2,0))*testp[l]*W;
1503 
1504  }
1505 
1506  //Axial velocity component
1507  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
1508  if(local_unknown >= 0)
1509  {
1510  jacobian(local_eqn,local_unknown) +=
1511  r*dpsifdx(l2,1)*testp[l]*W;
1512  }
1513  } /*End of loop over l2*/
1514  } /*End of Jacobian calculation*/
1515 
1516  } //End of if not boundary condition
1517 
1518  } //End of loop over l
1519  }
1520 
1521 }
1522 
1523 
1524 
1525 //======================================================================
1526 /// Compute derivatives of elemental residual vector with respect
1527 /// to nodal coordinates.
1528 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
1529 /// Overloads the FD-based version in the FiniteElement base class.
1530 //======================================================================
1533  RankThreeTensor<double>& dresidual_dnodal_coordinates)
1534 {
1535  // Return immediately if there are no dofs
1536  if(ndof()==0) { return; }
1537 
1538  // Get continuous time from timestepper of first node
1539  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
1540 
1541  // Determine number of nodes in element
1542  const unsigned n_node = nnode();
1543 
1544  // Determine number of pressure dofs in element
1545  const unsigned n_pres = npres_axi_nst();
1546 
1547  // Find the indices at which the local velocities are stored
1548  unsigned u_nodal_index[3];
1549  for(unsigned i=0;i<3;i++) { u_nodal_index[i] = u_index_axi_nst(i); }
1550 
1551  // Set up memory for the shape and test functions
1552  // Note that there are only two dimensions, r and z, in this problem
1553  Shape psif(n_node), testf(n_node);
1554  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
1555 
1556  // Set up memory for pressure shape and test functions
1557  Shape psip(n_pres), testp(n_pres);
1558 
1559  // Deriatives of shape fct derivatives w.r.t. nodal coords
1560  RankFourTensor<double> d_dpsifdx_dX(2,n_node,n_node,2);
1561  RankFourTensor<double> d_dtestfdx_dX(2,n_node,n_node,2);
1562 
1563  // Derivative of Jacobian of mapping w.r.t. to nodal coords
1564  DenseMatrix<double> dJ_dX(2,n_node);
1565 
1566  // Derivatives of derivative of u w.r.t. nodal coords
1567  // Note that the entry d_dudx_dX(p,q,i,k) contains the derivative w.r.t.
1568  // nodal coordinate X_pq of du_i/dx_k. Since there are three velocity
1569  // components, i can take on the values 0, 1 and 2, while k and p only
1570  // take on the values 0 and 1 since there are only two spatial dimensions.
1571  RankFourTensor<double> d_dudx_dX(2,n_node,3,2);
1572 
1573  // Derivatives of nodal velocities w.r.t. nodal coords:
1574  // Assumption: Interaction only local via no-slip so
1575  // X_pq only affects U_iq.
1576  // Note that the entry d_U_dX(p,q,i) contains the derivative w.r.t. nodal
1577  // coordinate X_pq of nodal value U_iq. In other words this entry is the
1578  // derivative of the i-th velocity component w.r.t. the p-th spatial
1579  // coordinate, taken at the q-th local node.
1580  RankThreeTensor<double> d_U_dX(2,n_node,3,0.0);
1581 
1582  // Determine the number of integration points
1583  const unsigned n_intpt = integral_pt()->nweight();
1584 
1585  // Vector to hold local coordinates (two dimensions)
1586  Vector<double> s(2);
1587 
1588  // Get physical variables from element
1589  // (Reynolds number must be multiplied by the density ratio)
1590  const double scaled_re = re()*density_ratio();
1591  const double scaled_re_st = re_st()*density_ratio();
1592  const double scaled_re_inv_fr = re_invfr()*density_ratio();
1593  const double scaled_re_inv_ro = re_invro()*density_ratio();
1594  const double visc_ratio = viscosity_ratio();
1595  const Vector<double> G = g();
1596 
1597  // FD step
1599 
1600  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1601  // Assumption: Interaction only local via no-slip so
1602  // X_ij only affects U_ij.
1603  bool element_has_node_with_aux_node_update_fct=false;
1604  for(unsigned q=0;q<n_node;q++)
1605  {
1606  // Get pointer to q-th local node
1607  Node* nod_pt = node_pt(q);
1608 
1609  // Only compute if there's a node-update fct involved
1610  if(nod_pt->has_auxiliary_node_update_fct_pt())
1611  {
1612  element_has_node_with_aux_node_update_fct=true;
1613 
1614  // This functionality has not been tested so issue a warning
1615  // to this effect
1616  std::ostringstream warning_stream;
1617  warning_stream << "\nThe functionality to evaluate the additional"
1618  << "\ncontribution to the deriv of the residual eqn"
1619  << "\nw.r.t. the nodal coordinates which comes about"
1620  << "\nif a node's values are updated using an auxiliary"
1621  << "\nnode update function has NOT been tested for"
1622  << "\naxisymmetric Navier-Stokes elements. Use at your"
1623  << "\nown risk"
1624  << std::endl;
1626  warning_stream.str(),
1627  "AxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates",
1628  OOMPH_EXCEPTION_LOCATION);
1629 
1630  // Current nodal velocity
1631  Vector<double> u_ref(3);
1632  for(unsigned i=0;i<3;i++)
1633  {
1634  u_ref[i]=*(nod_pt->value_pt(u_nodal_index[i]));
1635  }
1636 
1637  // FD
1638  for(unsigned p=0;p<2;p++)
1639  {
1640  // Make backup
1641  double backup=nod_pt->x(p);
1642 
1643  // Do FD step. No node update required as we're
1644  // attacking the coordinate directly...
1645  nod_pt->x(p) += eps_fd;
1646 
1647  // Do auxiliary node update (to apply no slip)
1649 
1650  // Loop over velocity components
1651  for(unsigned i=0;i<3;i++)
1652  {
1653  // Evaluate
1654  d_U_dX(p,q,i)=(*(nod_pt->value_pt(u_nodal_index[i]))-u_ref[i])/eps_fd;
1655  }
1656 
1657  // Reset
1658  nod_pt->x(p)=backup;
1659 
1660  // Do auxiliary node update (to apply no slip)
1662  }
1663  }
1664  }
1665 
1666  // Integer to store the local equation number
1667  int local_eqn=0;
1668 
1669  // Loop over the integration points
1670  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1671  {
1672  // Assign values of s
1673  for(unsigned i=0;i<2;i++) { s[i] = integral_pt()->knot(ipt,i); }
1674 
1675  // Get the integral weight
1676  const double w = integral_pt()->weight(ipt);
1677 
1678  // Call the derivatives of the shape and test functions
1680  ipt,psif,dpsifdx,d_dpsifdx_dX,testf,dtestfdx,d_dtestfdx_dX,dJ_dX);
1681 
1682  // Call the pressure shape and test functions
1683  pshape_axi_nst(s,psip,testp);
1684 
1685  // Allocate storage for the position and the derivative of the
1686  // mesh positions w.r.t. time
1688  Vector<double> mesh_velocity(2,0.0);
1689 
1690  // Allocate storage for the pressure, velocity components and their
1691  // time and spatial derivatives
1692  double interpolated_p=0.0;
1693  Vector<double> interpolated_u(3,0.0);
1694  Vector<double> dudt(3,0.0);
1695  DenseMatrix<double> interpolated_dudx(3,2,0.0);
1696 
1697  // Calculate pressure at integration point
1698  for(unsigned l=0;l<n_pres;l++) { interpolated_p += p_axi_nst(l)*psip[l]; }
1699 
1700  // Calculate velocities and derivatives at integration point
1701  // ---------------------------------------------------------
1702 
1703  // Loop over nodes
1704  for(unsigned l=0;l<n_node;l++)
1705  {
1706  // Cache the shape function
1707  const double psif_ = psif(l);
1708 
1709  // Loop over the two coordinate directions
1710  for(unsigned i=0;i<2;i++)
1711  {
1712  interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
1713  }
1714 
1715  // Loop over the three velocity directions
1716  for(unsigned i=0;i<3;i++)
1717  {
1718  // Get the nodal value
1719  const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
1720  interpolated_u[i] += u_value*psif_;
1721  dudt[i] += du_dt_axi_nst(l,i)*psif_;
1722 
1723  // Loop over derivative directions
1724  for(unsigned j=0;j<2;j++)
1725  {
1726  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
1727  }
1728  }
1729  }
1730 
1731  // Get the mesh velocity if ALE is enabled
1732  if(!ALE_is_disabled)
1733  {
1734  // Loop over nodes
1735  for(unsigned l=0;l<n_node;l++)
1736  {
1737  // Loop over the two coordinate directions
1738  for(unsigned i=0;i<2;i++)
1739  {
1740  mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif[l];
1741  }
1742  }
1743  }
1744 
1745  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1746  for(unsigned q=0;q<n_node;q++)
1747  {
1748  // Loop over the two coordinate directions
1749  for(unsigned p=0;p<2;p++)
1750  {
1751  // Loop over the three velocity components
1752  for(unsigned i=0;i<3;i++)
1753  {
1754  // Loop over the two coordinate directions
1755  for(unsigned k=0;k<2;k++)
1756  {
1757  double aux = 0.0;
1758 
1759  // Loop over nodes and add contribution
1760  for(unsigned j=0;j<n_node;j++)
1761  {
1762  aux +=
1763  this->raw_nodal_value(j,u_nodal_index[i])*d_dpsifdx_dX(p,q,j,k);
1764  }
1765  d_dudx_dX(p,q,i,k) = aux;
1766  }
1767  }
1768  }
1769  }
1770 
1771  // Get weight of actual nodal position/value in computation of mesh
1772  // velocity from positional/value time stepper
1773  const double pos_time_weight
1774  = node_pt(0)->position_time_stepper_pt()->weight(1,0);
1775  const double val_time_weight = node_pt(0)->time_stepper_pt()->weight(1,0);
1776 
1777  // Get the user-defined body force terms
1778  Vector<double> body_force(3);
1779  get_body_force_axi_nst(time,ipt,s,interpolated_x,body_force);
1780 
1781  // Get the user-defined source function
1782  const double source = get_source_fct(time,ipt,interpolated_x);
1783 
1784  // Note: Can use raw values and avoid bypassing hanging information
1785  // because this is the non-refineable version!
1786 
1787  // Get gradient of body force function
1788  DenseMatrix<double> d_body_force_dx(3,2,0.0);
1789  get_body_force_gradient_axi_nst(time,ipt,s,interpolated_x,
1790  d_body_force_dx);
1791 
1792  // Get gradient of source function
1793  Vector<double> source_gradient(2,0.0);
1794  get_source_fct_gradient(time,ipt,interpolated_x,source_gradient);
1795 
1796  // Cache r (the first position component)
1797  const double r = interpolated_x[0];
1798 
1799  // Assemble shape derivatives
1800  // --------------------------
1801 
1802  // ==================
1803  // MOMENTUM EQUATIONS
1804  // ==================
1805 
1806  // Loop over the test functions
1807  for(unsigned l=0;l<n_node;l++)
1808  {
1809  // Cache the test function
1810  const double testf_ = testf[l];
1811 
1812  // --------------------------------
1813  // FIRST (RADIAL) MOMENTUM EQUATION
1814  // --------------------------------
1815  local_eqn = nodal_local_eqn(l,u_nodal_index[0]);;
1816 
1817  // If it's not a boundary condition
1818  if(local_eqn >= 0)
1819  {
1820  // Loop over the two coordinate directions
1821  for(unsigned p=0;p<2;p++)
1822  {
1823  // Loop over nodes
1824  for(unsigned q=0;q<n_node;q++)
1825  {
1826 
1827  // Residual x deriv of Jacobian
1828  // ----------------------------
1829 
1830  // Add the user-defined body force terms
1831  double sum = r*body_force[0]*testf_;
1832 
1833  // Add the gravitational body force term
1834  sum += r*scaled_re_inv_fr*testf_*G[0];
1835 
1836  // Add the pressure gradient term
1837  sum += interpolated_p*(testf_ + r*dtestfdx(l,0));
1838 
1839  // Add the stress tensor terms
1840  // The viscosity ratio needs to go in here to ensure
1841  // continuity of normal stress is satisfied even in flows
1842  // with zero pressure gradient!
1843  sum -= visc_ratio*
1844  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0);
1845 
1846  sum -= visc_ratio*r*
1847  (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
1848  dtestfdx(l,1);
1849 
1850  sum -= visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]*testf_/r;
1851 
1852  // Add the du/dt term
1853  sum -= scaled_re_st*r*dudt[0]*testf_;
1854 
1855  // Add the convective terms
1856  sum -=
1857  scaled_re*(r*interpolated_u[0]*interpolated_dudx(0,0)
1858  - interpolated_u[2]*interpolated_u[2]
1859  + r*interpolated_u[1]*interpolated_dudx(0,1))*testf_;
1860 
1861  // Add the mesh velocity terms
1862  if(!ALE_is_disabled)
1863  {
1864  for(unsigned k=0;k<2;k++)
1865  {
1866  sum += scaled_re_st*r*mesh_velocity[k]
1867  *interpolated_dudx(0,k)*testf_;
1868  }
1869  }
1870 
1871  // Add the Coriolis term
1872  sum += 2.0*r*scaled_re_inv_ro*interpolated_u[2]*testf_;
1873 
1874  // Multiply through by deriv of Jacobian and integration weight
1875  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*dJ_dX(p,q)*w;
1876 
1877  // Derivative of residual x Jacobian
1878  // ---------------------------------
1879 
1880  // Body force
1881  sum = r*d_body_force_dx(0,p)*psif[q]*testf_;
1882  if(p==0) { sum += body_force[0]*psif[q]*testf_; }
1883 
1884  // Gravitational body force
1885  if(p==0) { sum += scaled_re_inv_fr*G[0]*psif[q]*testf_; }
1886 
1887  // Pressure gradient term
1888  sum += r*interpolated_p*d_dtestfdx_dX(p,q,l,0);
1889  if(p==0) { sum += interpolated_p*psif[q]*dtestfdx(l,0); }
1890 
1891  // Viscous terms
1892  sum -= r*visc_ratio*(
1893  (1.0+Gamma[0])*(
1894  d_dudx_dX(p,q,0,0)*dtestfdx(l,0)
1895  + interpolated_dudx(0,0)*d_dtestfdx_dX(p,q,l,0))
1896  + (d_dudx_dX(p,q,0,1)
1897  + Gamma[0]*d_dudx_dX(p,q,1,0))*dtestfdx(l,1)
1898  + (interpolated_dudx(0,1)
1899  + Gamma[0]*interpolated_dudx(1,0))
1900  *d_dtestfdx_dX(p,q,l,1));
1901  if(p==0)
1902  {
1903  sum -= visc_ratio*(
1904  (1.0+Gamma[0])*(
1905  interpolated_dudx(0,0)*psif[q]*dtestfdx(l,0)
1906  - interpolated_u[0]*psif[q]*testf_/(r*r))
1907  + (interpolated_dudx(0,1)
1908  + Gamma[0]*interpolated_dudx(1,0))*psif[q]*dtestfdx(l,1));
1909  }
1910 
1911  // Convective terms, including mesh velocity
1912  for(unsigned k=0;k<2;k++)
1913  {
1914  double tmp = scaled_re*interpolated_u[k];
1915  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
1916  sum -= r*tmp*d_dudx_dX(p,q,0,k)*testf_;
1917  if(p==0) { sum -= tmp*interpolated_dudx(0,k)*psif[q]*testf_; }
1918  }
1919  if(!ALE_is_disabled)
1920  {
1921  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(0,p)
1922  *psif[q]*testf_;
1923  }
1924 
1925  // du/dt term
1926  if(p==0) { sum -= scaled_re_st*dudt[0]*psif[q]*testf_; }
1927 
1928  // Coriolis term
1929  if(p==0)
1930  {
1931  sum += 2.0*scaled_re_inv_ro*interpolated_u[2]*psif[q]*testf_;
1932  }
1933 
1934  // Multiply through by Jacobian and integration weight
1935  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
1936 
1937  // Derivs w.r.t. nodal velocities
1938  // ------------------------------
1939  if(element_has_node_with_aux_node_update_fct)
1940  {
1941  // Contribution from deriv of first velocity component
1942  double tmp = scaled_re_st*r*val_time_weight*psif[q]*testf_;
1943  tmp += r*scaled_re*interpolated_dudx(0,0)*psif[q]*testf_;
1944  for(unsigned k=0;k<2;k++)
1945  {
1946  double tmp2 = scaled_re*interpolated_u[k];
1947  if(!ALE_is_disabled) { tmp2 -= scaled_re_st*mesh_velocity[k]; }
1948  tmp += r*tmp2*dpsifdx(q,k)*testf_;
1949  }
1950 
1951  tmp += r*visc_ratio*(1.0+Gamma[0])*dpsifdx(q,0)*dtestfdx(l,0);
1952  tmp += r*visc_ratio*dpsifdx(q,1)*dtestfdx(l,1);
1953  tmp += visc_ratio*(1.0+Gamma[0])*psif[q]*testf_/r;
1954 
1955  // Multiply through by dU_0q/dX_pq
1956  sum = -d_U_dX(p,q,0)*tmp;
1957 
1958  // Contribution from deriv of second velocity component
1959  tmp = scaled_re*r*interpolated_dudx(0,1)*psif[q]*testf_;
1960  tmp += r*visc_ratio*Gamma[0]*dpsifdx(q,0)*dtestfdx(l,1);
1961 
1962  // Multiply through by dU_1q/dX_pq
1963  sum -= d_U_dX(p,q,1)*tmp;
1964 
1965  // Contribution from deriv of third velocity component
1966  tmp = 2.0*(r*scaled_re_inv_ro
1967  + scaled_re*interpolated_u[2])*psif[q]*testf_;
1968 
1969  // Multiply through by dU_2q/dX_pq
1970  sum += d_U_dX(p,q,2)*tmp;
1971 
1972  // Multiply through by Jacobian and integration weight
1973  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
1974  }
1975  } // End of loop over q
1976  } // End of loop over p
1977  } // End of if it's not a boundary condition
1978 
1979  // --------------------------------
1980  // SECOND (AXIAL) MOMENTUM EQUATION
1981  // --------------------------------
1982  local_eqn = nodal_local_eqn(l,u_nodal_index[1]);;
1983 
1984  // If it's not a boundary condition
1985  if(local_eqn >= 0)
1986  {
1987  // Loop over the two coordinate directions
1988  for(unsigned p=0;p<2;p++)
1989  {
1990  // Loop over nodes
1991  for(unsigned q=0;q<n_node;q++)
1992  {
1993 
1994  // Residual x deriv of Jacobian
1995  // ----------------------------
1996 
1997  // Add the user-defined body force terms
1998  double sum = r*body_force[1]*testf_;
1999 
2000  // Add the gravitational body force term
2001  sum += r*scaled_re_inv_fr*testf_*G[1];
2002 
2003  // Add the pressure gradient term
2004  sum += r*interpolated_p*dtestfdx(l,1);
2005 
2006  // Add the stress tensor terms
2007  // The viscosity ratio needs to go in here to ensure
2008  // continuity of normal stress is satisfied even in flows
2009  // with zero pressure gradient!
2010  sum -= visc_ratio*
2011  r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
2012  *dtestfdx(l,0);
2013 
2014  sum -= visc_ratio*r*
2015  (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1);
2016 
2017  // Add the du/dt term
2018  sum -= scaled_re_st*r*dudt[1]*testf_;
2019 
2020  // Add the convective terms
2021  sum -=
2022  scaled_re*(r*interpolated_u[0]*interpolated_dudx(1,0)
2023  + r*interpolated_u[1]*interpolated_dudx(1,1))*testf_;
2024 
2025  // Add the mesh velocity terms
2026  if(!ALE_is_disabled)
2027  {
2028  for(unsigned k=0;k<2;k++)
2029  {
2030  sum +=
2031  scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(1,k)*testf_;
2032  }
2033  }
2034 
2035  // Multiply through by deriv of Jacobian and integration weight
2036  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*dJ_dX(p,q)*w;
2037 
2038  // Derivative of residual x Jacobian
2039  // ---------------------------------
2040 
2041  // Body force
2042  sum = r*d_body_force_dx(1,p)*psif[q]*testf_;
2043  if(p==0) { sum += body_force[1]*psif[q]*testf_; }
2044 
2045  // Gravitational body force
2046  if(p==0) { sum += scaled_re_inv_fr*G[1]*psif[q]*testf_; }
2047 
2048  // Pressure gradient term
2049  sum += r*interpolated_p*d_dtestfdx_dX(p,q,l,1);
2050  if(p==0) { sum += interpolated_p*psif[q]*dtestfdx(l,1); }
2051 
2052  // Viscous terms
2053  sum -= r*visc_ratio*(
2054  (d_dudx_dX(p,q,1,0) + Gamma[1]*d_dudx_dX(p,q,0,1))*dtestfdx(l,0)
2055  + (interpolated_dudx(1,0)
2056  + Gamma[1]*interpolated_dudx(0,1))*d_dtestfdx_dX(p,q,l,0)
2057  + (1.0+Gamma[1])*d_dudx_dX(p,q,1,1)*dtestfdx(l,1)
2058  + (1.0+Gamma[1])*interpolated_dudx(1,1)*d_dtestfdx_dX(p,q,l,1));
2059  if(p==0)
2060  {
2061  sum -= visc_ratio*(
2062  (interpolated_dudx(1,0)
2063  + Gamma[1]*interpolated_dudx(0,1))*psif[q]*dtestfdx(l,0)
2064  + (1.0+Gamma[1])*interpolated_dudx(1,1)*psif[q]*dtestfdx(l,1));
2065  }
2066 
2067  // Convective terms, including mesh velocity
2068  for(unsigned k=0;k<2;k++)
2069  {
2070  double tmp = scaled_re*interpolated_u[k];
2071  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
2072  sum -= r*tmp*d_dudx_dX(p,q,1,k)*testf_;
2073  if(p==0) { sum -= tmp*interpolated_dudx(1,k)*psif[q]*testf_; }
2074  }
2075  if(!ALE_is_disabled)
2076  {
2077  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(1,p)
2078  *psif[q]*testf_;
2079  }
2080 
2081  // du/dt term
2082  if(p==0) { sum -= scaled_re_st*dudt[1]*psif[q]*testf_; }
2083 
2084  // Multiply through by Jacobian and integration weight
2085  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
2086 
2087  // Derivs w.r.t. to nodal velocities
2088  // ---------------------------------
2089  if(element_has_node_with_aux_node_update_fct)
2090  {
2091  // Contribution from deriv of first velocity component
2092  double tmp = scaled_re*r*interpolated_dudx(1,0)*psif[q]*testf_;
2093  tmp += r*visc_ratio*Gamma[1]*dpsifdx(q,1)*dtestfdx(l,0);
2094 
2095  // Multiply through by dU_0q/dX_pq
2096  sum = -d_U_dX(p,q,0)*tmp;
2097 
2098  // Contribution from deriv of second velocity component
2099  tmp = scaled_re_st*r*val_time_weight*psif[q]*testf_;
2100  tmp += r*scaled_re*interpolated_dudx(1,1)*psif[q]*testf_;
2101  for(unsigned k=0;k<2;k++)
2102  {
2103  double tmp2 = scaled_re*interpolated_u[k];
2104  if(!ALE_is_disabled) { tmp2 -= scaled_re_st*mesh_velocity[k]; }
2105  tmp += r*tmp2*dpsifdx(q,k)*testf_;
2106  }
2107  tmp += r*visc_ratio*(dpsifdx(q,0)*dtestfdx(l,0)
2108  + (1.0+Gamma[1])*dpsifdx(q,1)*dtestfdx(l,1));
2109 
2110  // Multiply through by dU_1q/dX_pq
2111  sum -= d_U_dX(p,q,1)*tmp;
2112 
2113  // Multiply through by Jacobian and integration weight
2114  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
2115  }
2116  } // End of loop over q
2117  } // End of loop over p
2118  } // End of if it's not a boundary condition
2119 
2120  // -----------------------------------
2121  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
2122  // -----------------------------------
2123  local_eqn = nodal_local_eqn(l,u_nodal_index[2]);;
2124 
2125  // If it's not a boundary condition
2126  if(local_eqn >= 0)
2127  {
2128  // Loop over the two coordinate directions
2129  for(unsigned p=0;p<2;p++)
2130  {
2131  // Loop over nodes
2132  for(unsigned q=0;q<n_node;q++)
2133  {
2134 
2135  // Residual x deriv of Jacobian
2136  // ----------------------------
2137 
2138  // Add the user-defined body force terms
2139  double sum = r*body_force[2]*testf_;
2140 
2141  // Add the gravitational body force term
2142  sum += r*scaled_re_inv_fr*testf_*G[2];
2143 
2144  // There is NO pressure gradient term
2145 
2146  // Add in the stress tensor terms
2147  // The viscosity ratio needs to go in here to ensure
2148  // continuity of normal stress is satisfied even in flows
2149  // with zero pressure gradient!
2150  sum -= visc_ratio*(r*interpolated_dudx(2,0) -
2151  Gamma[0]*interpolated_u[2])*dtestfdx(l,0);
2152 
2153  sum -= visc_ratio*r*interpolated_dudx(2,1)*dtestfdx(l,1);
2154 
2155  sum -= visc_ratio*((interpolated_u[2]/r)
2156  - Gamma[0]*interpolated_dudx(2,0))*testf_;
2157 
2158  // Add the du/dt term
2159  sum -= scaled_re_st*r*dudt[2]*testf_;
2160 
2161  // Add the convective terms
2162  sum -=
2163  scaled_re*(r*interpolated_u[0]*interpolated_dudx(2,0)
2164  + interpolated_u[0]*interpolated_u[2]
2165  + r*interpolated_u[1]*interpolated_dudx(2,1))*testf_;
2166 
2167  // Add the mesh velocity terms
2168  if(!ALE_is_disabled)
2169  {
2170  for(unsigned k=0;k<2;k++)
2171  {
2172  sum +=
2173  scaled_re_st*r*mesh_velocity[k]*interpolated_dudx(2,k)*testf_;
2174  }
2175  }
2176 
2177  // Add the Coriolis term
2178  sum -= 2.0*r*scaled_re_inv_ro*interpolated_u[0]*testf_;
2179 
2180  // Multiply through by deriv of Jacobian and integration weight
2181  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*dJ_dX(p,q)*w;
2182 
2183  // Derivative of residual x Jacobian
2184  // ---------------------------------
2185 
2186  // Body force
2187  sum = r*d_body_force_dx(2,p)*psif[q]*testf_;
2188  if(p==0) { sum += body_force[2]*psif[q]*testf_; }
2189 
2190  // Gravitational body force
2191  if(p==0) { sum += scaled_re_inv_fr*G[2]*psif[q]*testf_; }
2192 
2193  // There is NO pressure gradient term
2194 
2195  // Viscous terms
2196  sum -= visc_ratio*r*(
2197  d_dudx_dX(p,q,2,0)*dtestfdx(l,0)
2198  + interpolated_dudx(2,0)*d_dtestfdx_dX(p,q,l,0)
2199  + d_dudx_dX(p,q,2,1)*dtestfdx(l,1)
2200  + interpolated_dudx(2,1)*d_dtestfdx_dX(p,q,l,1));
2201 
2202  sum += visc_ratio*Gamma[0]*d_dudx_dX(p,q,2,0)*testf_;
2203 
2204  if(p==0)
2205  {
2206  sum -= visc_ratio*(interpolated_dudx(2,0)*psif[q]*dtestfdx(l,0)
2207  + interpolated_dudx(2,1)*psif[q]*dtestfdx(l,1)
2208  + interpolated_u[2]*psif[q]*testf_/(r*r));
2209  }
2210 
2211  // Convective terms, including mesh velocity
2212  for(unsigned k=0;k<2;k++)
2213  {
2214  double tmp = scaled_re*interpolated_u[k];
2215  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
2216  sum -= r*tmp*d_dudx_dX(p,q,2,k)*testf_;
2217  if(p==0) { sum -= tmp*interpolated_dudx(2,k)*psif[q]*testf_; }
2218  }
2219  if(!ALE_is_disabled)
2220  {
2221  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(2,p)
2222  *psif[q]*testf_;
2223  }
2224 
2225  // du/dt term
2226  if(p==0) { sum -= scaled_re_st*dudt[2]*psif[q]*testf_; }
2227 
2228  // Coriolis term
2229  if(p==0)
2230  {
2231  sum -= 2.0*scaled_re_inv_ro*interpolated_u[0]*psif[q]*testf_;
2232  }
2233 
2234  // Multiply through by Jacobian and integration weight
2235  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
2236 
2237  // Derivs w.r.t. to nodal velocities
2238  // ---------------------------------
2239  if(element_has_node_with_aux_node_update_fct)
2240  {
2241  // Contribution from deriv of first velocity component
2242  double tmp = (2.0*r*scaled_re_inv_ro
2243  + r*scaled_re*interpolated_dudx(2,0)
2244  - scaled_re*interpolated_u[2])*psif[q]*testf_;
2245 
2246  // Multiply through by dU_0q/dX_pq
2247  sum = -d_U_dX(p,q,0)*tmp;
2248 
2249  // Contribution from deriv of second velocity component
2250  // Multiply through by dU_1q/dX_pq
2251  sum -= d_U_dX(p,q,1)*scaled_re*r*interpolated_dudx(2,1)
2252  *psif[q]*testf_;
2253 
2254  // Contribution from deriv of third velocity component
2255  tmp = scaled_re_st*r*val_time_weight*psif[q]*testf_;
2256  tmp -= scaled_re*interpolated_u[0]*psif[q]*testf_;
2257  for(unsigned k=0;k<2;k++)
2258  {
2259  double tmp2 = scaled_re*interpolated_u[k];
2260  if(!ALE_is_disabled) { tmp2 -= scaled_re_st*mesh_velocity[k]; }
2261  tmp += r*tmp2*dpsifdx(q,k)*testf_;
2262  }
2263  tmp += visc_ratio*(r*dpsifdx(q,0)
2264  - Gamma[0]*psif[q])*dtestfdx(l,0);
2265  tmp += r*visc_ratio*dpsifdx(q,1)*dtestfdx(l,1);
2266  tmp += visc_ratio*((psif[q]/r) - Gamma[0]*dpsifdx(q,0))*testf_;
2267 
2268  // Multiply through by dU_2q/dX_pq
2269  sum -= d_U_dX(p,q,2)*tmp;
2270 
2271  // Multiply through by Jacobian and integration weight
2272  dresidual_dnodal_coordinates(local_eqn,p,q) += sum*J*w;
2273  }
2274  } // End of loop over q
2275  } // End of loop over p
2276  } // End of if it's not a boundary condition
2277 
2278  } // End of loop over test functions
2279 
2280 
2281  // ===================
2282  // CONTINUITY EQUATION
2283  // ===================
2284 
2285  // Loop over the shape functions
2286  for(unsigned l=0;l<n_pres;l++)
2287  {
2288  local_eqn = p_local_eqn(l);
2289 
2290  // Cache the test function
2291  const double testp_ = testp[l];
2292 
2293  // If not a boundary conditions
2294  if(local_eqn >= 0)
2295  {
2296 
2297  // Loop over the two coordinate directions
2298  for(unsigned p=0;p<2;p++)
2299  {
2300  // Loop over nodes
2301  for(unsigned q=0;q<n_node;q++)
2302  {
2303 
2304  // Residual x deriv of Jacobian
2305  //-----------------------------
2306 
2307  // Source term
2308  double aux = -r*source;
2309 
2310  // Gradient terms
2311  aux += (interpolated_u[0]
2312  + r*interpolated_dudx(0,0)
2313  + r*interpolated_dudx(1,1));
2314 
2315  // Multiply through by deriv of Jacobian and integration weight
2316  dresidual_dnodal_coordinates(local_eqn,p,q) +=
2317  aux*dJ_dX(p,q)*testp_*w;
2318 
2319  // Derivative of residual x Jacobian
2320  // ---------------------------------
2321 
2322  // Gradient of source function
2323  aux = -r*source_gradient[p]*psif[q];
2324  if(p==0) { aux -= source*psif[q]; }
2325 
2326  // Gradient terms
2327  aux += r*(d_dudx_dX(p,q,0,0) + d_dudx_dX(p,q,1,1));
2328  if(p==0)
2329  {
2330  aux += (interpolated_dudx(0,0) + interpolated_dudx(1,1))*psif[q];
2331  }
2332 
2333  // Derivs w.r.t. nodal velocities
2334  // ------------------------------
2335  if(element_has_node_with_aux_node_update_fct)
2336  {
2337  aux += d_U_dX(p,q,0)*(psif[q] + r*dpsifdx(q,0));
2338  aux += d_U_dX(p,q,1)*r*dpsifdx(q,1);
2339  }
2340 
2341  // Multiply through by Jacobian, test fct and integration weight
2342  dresidual_dnodal_coordinates(local_eqn,p,q) += aux*testp_*J*w;
2343 
2344  }
2345  }
2346  }
2347  } // End of loop over shape functions for continuity eqn
2348 
2349  } // End of loop over integration points
2350 }
2351 
2352 
2353 
2354 
2355 
2356 //==============================================================
2357 /// Compute the residuals for the Navier--Stokes
2358 /// equations; flag=1(or 0): do (or don't) compute the
2359 /// Jacobian as well.
2360 //==============================================================
2363  double* const &parameter_pt,
2364  Vector<double> &dres_dparam,
2365  DenseMatrix<double> &djac_dparam,
2366  DenseMatrix<double> &dmass_matrix_dparam,
2367  unsigned flag)
2368 {
2369  //Die if the parameter is not the Reynolds number
2370  if(parameter_pt!=this->re_pt())
2371  {
2372  std::ostringstream error_stream;
2373  error_stream <<
2374  "Cannot compute analytic jacobian for parameter addressed by "
2375  << parameter_pt << "\n";
2376  error_stream << "Can only compute derivatives wrt Re ("
2377  << Re_pt << ")\n";
2378  throw OomphLibError(
2379  error_stream.str(),
2380  OOMPH_CURRENT_FUNCTION,
2381  OOMPH_EXCEPTION_LOCATION);
2382  }
2383 
2384  //Which parameters are we differentiating with respect to
2385  bool diff_re = false;
2386  bool diff_re_st = false;
2387  bool diff_re_inv_fr = false;
2388  bool diff_re_inv_ro = false;
2389 
2390  //Set the boolean flags according to the parameter pointer
2391  if(parameter_pt==this->re_pt()) {diff_re = true;}
2392  if(parameter_pt==this->re_st_pt()) {diff_re_st = true;}
2393  if(parameter_pt==this->re_invfr_pt()) {diff_re_inv_fr = true;}
2394  if(parameter_pt==this->re_invro_pt()) {diff_re_inv_ro = true;}
2395 
2396 
2397  //Find out how many nodes there are
2398  unsigned n_node = nnode();
2399 
2400  //Find out how many pressure dofs there are
2401  unsigned n_pres = npres_axi_nst();
2402 
2403  //Get the nodal indices at which the velocity is stored
2404  unsigned u_nodal_index[3];
2405  for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}
2406 
2407  //Set up memory for the shape and test functions
2408  //Note that there are only two dimensions, r and z in this problem
2409  Shape psif(n_node), testf(n_node);
2410  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
2411 
2412  //Set up memory for pressure shape and test functions
2413  Shape psip(n_pres), testp(n_pres);
2414 
2415  //Number of integration points
2416  unsigned Nintpt = integral_pt()->nweight();
2417 
2418  //Set the Vector to hold local coordinates (two dimensions)
2419  Vector<double> s(2);
2420 
2421  //Get Physical Variables from Element
2422  //Reynolds number must be multiplied by the density ratio
2423  //double scaled_re = re()*density_ratio();
2424  //double scaled_re_st = re_st()*density_ratio();
2425  //double scaled_re_inv_fr = re_invfr()*density_ratio();
2426  //double scaled_re_inv_ro = re_invro()*density_ratio();
2427  double dens_ratio = this->density_ratio();
2428  // double visc_ratio = viscosity_ratio();
2429  Vector<double> G = g();
2430 
2431  //Integers used to store the local equation and unknown numbers
2432  int local_eqn=0, local_unknown=0;
2433 
2434  //Loop over the integration points
2435  for(unsigned ipt=0;ipt<Nintpt;ipt++)
2436  {
2437  //Assign values of s
2438  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
2439  //Get the integral weight
2440  double w = integral_pt()->weight(ipt);
2441 
2442  //Call the derivatives of the shape and test functions
2443  double J =
2444  dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,psif,dpsifdx,testf,dtestfdx);
2445 
2446  //Call the pressure shape and test functions
2447  pshape_axi_nst(s,psip,testp);
2448 
2449  //Premultiply the weights and the Jacobian
2450  double W = w*J;
2451 
2452  //Allocate storage for the position and the derivative of the
2453  //mesh positions wrt time
2455  Vector<double> mesh_velocity(2,0.0);
2456  //Allocate storage for the pressure, velocity components and their
2457  //time and spatial derivatives
2458  double interpolated_p=0.0;
2459  Vector<double> interpolated_u(3,0.0);
2460  Vector<double> dudt(3,0.0);
2461  DenseMatrix<double> interpolated_dudx(3,2,0.0);
2462 
2463  //Calculate pressure at integration point
2464  for(unsigned l=0;l<n_pres;l++) {interpolated_p += p_axi_nst(l)*psip[l];}
2465 
2466  //Calculate velocities and derivatives at integration point
2467 
2468  // Loop over nodes
2469  for(unsigned l=0;l<n_node;l++)
2470  {
2471  //Cache the shape function
2472  const double psif_ = psif(l);
2473  //Loop over the two coordinate directions
2474  for(unsigned i=0;i<2;i++)
2475  {
2476  interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
2477  }
2478  //mesh_velocity[i] += dnodal_position_dt(l,i)*psif[l];
2479 
2480  //Loop over the three velocity directions
2481  for(unsigned i=0;i<3;i++)
2482  {
2483  //Get the u_value
2484  const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
2485  interpolated_u[i] += u_value*psif_;
2486  dudt[i]+= du_dt_axi_nst(l,i)*psif_;
2487  //Loop over derivative directions
2488  for(unsigned j=0;j<2;j++)
2489  {interpolated_dudx(i,j) += u_value*dpsifdx(l,j);}
2490  }
2491  }
2492 
2493  //Get the mesh velocity if ALE is enabled
2494  if(!ALE_is_disabled)
2495  {
2496  // Loop over nodes
2497  for(unsigned l=0;l<n_node;l++)
2498  {
2499  //Loop over the two coordinate directions
2500  for(unsigned i=0;i<2;i++)
2501  {
2502  mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif(l);
2503  }
2504  }
2505  }
2506 
2507 
2508  //Get the user-defined body force terms
2509  //Vector<double> body_force(3);
2510  //get_body_force(time(),ipt,interpolated_x,body_force);
2511 
2512  //Get the user-defined source function
2513  //double source = get_source_fct(time(),ipt,interpolated_x);
2514 
2515  //Get the user-defined viscosity function
2516  double visc_ratio;
2518  s,
2519  interpolated_x,
2520  visc_ratio);
2521 
2522  //r is the first position component
2523  double r = interpolated_x[0];
2524 
2525 
2526  //MOMENTUM EQUATIONS
2527  //------------------
2528 
2529  //Loop over the test functions
2530  for(unsigned l=0;l<n_node;l++)
2531  {
2532 
2533  //FIRST (RADIAL) MOMENTUM EQUATION
2534  local_eqn = nodal_local_eqn(l,u_nodal_index[0]);
2535  //If it's not a boundary condition
2536  if(local_eqn >= 0)
2537  {
2538  //No user-defined body force terms
2539  //dres_dparam[local_eqn] +=
2540  // r*body_force[0]*testf[l]*W;
2541 
2542  //Add the gravitational body force term if the reynolds number
2543  //is equal to re_inv_fr
2544  if(diff_re_inv_fr)
2545  {
2546  dres_dparam[local_eqn] += r*dens_ratio*testf[l]*G[0]*W;
2547  }
2548 
2549  //No pressure gradient term
2550  //residuals[local_eqn] +=
2551  // interpolated_p*(testf[l] + r*dtestfdx(l,0))*W;
2552 
2553  //No in the stress tensor terms
2554  //The viscosity ratio needs to go in here to ensure
2555  //continuity of normal stress is satisfied even in flows
2556  //with zero pressure gradient!
2557  //residuals[local_eqn] -= visc_ratio*
2558  // r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
2559 
2560  //residuals[local_eqn] -= visc_ratio*r*
2561  // (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
2562  // dtestfdx(l,1)*W;
2563 
2564  //residuals[local_eqn] -=
2565  // visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
2566 
2567  //Add in the inertial terms
2568  //du/dt term
2569  if(diff_re_st)
2570  {
2571  dres_dparam[local_eqn] -= dens_ratio*r*dudt[0]*testf[l]*W;
2572  }
2573 
2574  //Convective terms
2575  if(diff_re)
2576  {
2577  dres_dparam[local_eqn] -=
2578  dens_ratio*(r*interpolated_u[0]*interpolated_dudx(0,0)
2579  - interpolated_u[2]*interpolated_u[2]
2580  + r*interpolated_u[1]*interpolated_dudx(0,1))*testf[l]*W;
2581  }
2582 
2583  //Mesh velocity terms
2584  if(!ALE_is_disabled)
2585  {
2586  if(diff_re_st)
2587  {
2588  for(unsigned k=0;k<2;k++)
2589  {
2590  dres_dparam[local_eqn] +=
2591  dens_ratio*r*mesh_velocity[k]*interpolated_dudx(0,k)*testf[l]*W;
2592  }
2593  }
2594  }
2595 
2596  //Add the Coriolis term
2597  if(diff_re_inv_ro)
2598  {
2599  dres_dparam[local_eqn] +=
2600  2.0*r*dens_ratio*interpolated_u[2]*testf[l]*W;
2601  }
2602 
2603  //CALCULATE THE JACOBIAN
2604  if(flag)
2605  {
2606  //Loop over the velocity shape functions again
2607  for(unsigned l2=0;l2<n_node;l2++)
2608  {
2609  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
2610  //Radial velocity component
2611  if(local_unknown >= 0)
2612  {
2613  if(flag==2)
2614  {
2615  if(diff_re_st)
2616  {
2617  //Add the mass matrix
2618  dmass_matrix_dparam(local_eqn,local_unknown) +=
2619  dens_ratio*r*psif[l2]*testf[l]*W;
2620  }
2621  }
2622 
2623  //Add contribution to the Jacobian matrix
2624  //jacobian(local_eqn,local_unknown)
2625  // -= visc_ratio*r*(1.0+Gamma[0])
2626  //*dpsifdx(l2,0)*dtestfdx(l,0)*W;
2627 
2628  //jacobian(local_eqn,local_unknown)
2629  // -= visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
2630 
2631  //jacobian(local_eqn,local_unknown)
2632  // -= visc_ratio*(1.0 + Gamma[0])*psif[l2]*testf[l]*W/r;
2633 
2634  //Add in the inertial terms
2635  //du/dt term
2636  if(diff_re_st)
2637  {
2638  djac_dparam(local_eqn,local_unknown)
2639  -= dens_ratio*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
2640  psif[l2]*testf[l]*W;
2641  }
2642 
2643  //Convective terms
2644  if(diff_re)
2645  {
2646  djac_dparam(local_eqn,local_unknown) -=
2647  dens_ratio*(r*psif[l2]*interpolated_dudx(0,0)
2648  + r*interpolated_u[0]*dpsifdx(l2,0)
2649  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
2650  }
2651 
2652  //Mesh velocity terms
2653  if(!ALE_is_disabled)
2654  {
2655  for(unsigned k=0;k<2;k++)
2656  {
2657  if(diff_re_st)
2658  {
2659  djac_dparam(local_eqn,local_unknown) +=
2660  dens_ratio*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
2661  }
2662  }
2663  }
2664  }
2665 
2666  //Axial velocity component
2667  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
2668  if(local_unknown >= 0)
2669  {
2670  //jacobian(local_eqn,local_unknown) -=
2671  // visc_ratio*r*Gamma[0]*dpsifdx(l2,0)*dtestfdx(l,1)*W;
2672 
2673  //Convective terms
2674  if(diff_re)
2675  {
2676  djac_dparam(local_eqn,local_unknown) -=
2677  dens_ratio*r*psif[l2]*interpolated_dudx(0,1)*testf[l]*W;
2678  }
2679  }
2680 
2681  //Azimuthal velocity component
2682  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
2683  if(local_unknown >= 0)
2684  {
2685  //Convective terms
2686  if(diff_re)
2687  {
2688  djac_dparam(local_eqn,local_unknown) -=
2689  - dens_ratio*2.0*interpolated_u[2]*psif[l2]*testf[l]*W;
2690  }
2691  //Coriolis terms
2692  if(diff_re_inv_ro)
2693  {
2694  djac_dparam(local_eqn,local_unknown) +=
2695  2.0*r*dens_ratio*psif[l2]*testf[l]*W;
2696  }
2697  }
2698  }
2699 
2700  /*Now loop over pressure shape functions*/
2701  /*This is the contribution from pressure gradient*/
2702  //for(unsigned l2=0;l2<n_pres;l2++)
2703  // {
2704  // local_unknown = p_local_eqn(l2);
2705  // /*If we are at a non-zero degree of freedom in the entry*/
2706  // if(local_unknown >= 0)
2707  // {
2708  // jacobian(local_eqn,local_unknown)
2709  // += psip[l2]*(testf[l] + r*dtestfdx(l,0))*W;
2710  // }
2711  // }
2712  } /*End of Jacobian calculation*/
2713 
2714  } //End of if not boundary condition statement
2715 
2716  //SECOND (AXIAL) MOMENTUM EQUATION
2717  local_eqn = nodal_local_eqn(l,u_nodal_index[1]);
2718  //If it's not a boundary condition
2719  if(local_eqn >= 0)
2720  {
2721  //Add the user-defined body force terms
2722  //Remember to multiply by the density ratio!
2723  //residuals[local_eqn] +=
2724  // r*body_force[1]*testf[l]*W;
2725 
2726  //Add the gravitational body force term
2727  if(diff_re_inv_fr)
2728  {
2729  dres_dparam[local_eqn] += r*dens_ratio*testf[l]*G[1]*W;
2730  }
2731 
2732  //Add the pressure gradient term
2733  //residuals[local_eqn] += r*interpolated_p*dtestfdx(l,1)*W;
2734 
2735  //Add in the stress tensor terms
2736  //The viscosity ratio needs to go in here to ensure
2737  //continuity of normal stress is satisfied even in flows
2738  //with zero pressure gradient!
2739  //residuals[local_eqn] -= visc_ratio*
2740  // r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
2741  // *dtestfdx(l,0)*W;
2742 
2743  //residuals[local_eqn] -= visc_ratio*r*
2744  // (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
2745 
2746  //Add in the inertial terms
2747  //du/dt term
2748  if(diff_re_st)
2749  {
2750  dres_dparam[local_eqn] -= dens_ratio*r*dudt[1]*testf[l]*W;
2751  }
2752 
2753  //Convective terms
2754  if(diff_re)
2755  {
2756  dres_dparam[local_eqn] -=
2757  dens_ratio*(r*interpolated_u[0]*interpolated_dudx(1,0)
2758  + r*interpolated_u[1]*interpolated_dudx(1,1))*testf[l]*W;
2759  }
2760 
2761  //Mesh velocity terms
2762  if(!ALE_is_disabled)
2763  {
2764  if(diff_re_st)
2765  {
2766  for(unsigned k=0;k<2;k++)
2767  {
2768  dres_dparam[local_eqn] +=
2769  dens_ratio*r*mesh_velocity[k]*interpolated_dudx(1,k)*testf[l]*W;
2770  }
2771  }
2772  }
2773 
2774  //CALCULATE THE JACOBIAN
2775  if(flag)
2776  {
2777  //Loop over the velocity shape functions again
2778  for(unsigned l2=0;l2<n_node;l2++)
2779  {
2780  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
2781  //Radial velocity component
2782  if(local_unknown >= 0)
2783  {
2784  //Add in the stress tensor terms
2785  //The viscosity ratio needs to go in here to ensure
2786  //continuity of normal stress is satisfied even in flows
2787  //with zero pressure gradient!
2788  //jacobian(local_eqn,local_unknown) -=
2789  //visc_ratio*r*Gamma[1]*dpsifdx(l2,1)*dtestfdx(l,0)*W;
2790 
2791  //Convective terms
2792  if(diff_re)
2793  {
2794  djac_dparam(local_eqn,local_unknown) -=
2795  dens_ratio*r*psif[l2]*interpolated_dudx(1,0)*testf[l]*W;
2796  }
2797  }
2798 
2799  //Axial velocity component
2800  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
2801  if(local_unknown >= 0)
2802  {
2803  if(flag==2)
2804  {
2805  if(diff_re_st)
2806  {
2807  //Add the mass matrix
2808  dmass_matrix_dparam(local_eqn,local_unknown) +=
2809  dens_ratio*r*psif[l2]*testf[l]*W;
2810  }
2811  }
2812 
2813 
2814  //Add in the stress tensor terms
2815  //The viscosity ratio needs to go in here to ensure
2816  //continuity of normal stress is satisfied even in flows
2817  //with zero pressure gradient!
2818  //jacobian(local_eqn,local_unknown) -=
2819  //visc_ratio*r*dpsifdx(l2,0)*dtestfdx(l,0)*W;
2820 
2821  //jacobian(local_eqn,local_unknown) -=
2822  // visc_ratio*r*(1.0 + Gamma[1])*dpsifdx(l2,1)*
2823  // dtestfdx(l,1)*W;
2824 
2825  //Add in the inertial terms
2826  //du/dt term
2827  if(diff_re_st)
2828  {
2829  djac_dparam(local_eqn,local_unknown) -=
2830  dens_ratio*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
2831  psif[l2]*testf[l]*W;
2832  }
2833 
2834  //Convective terms
2835  if(diff_re)
2836  {
2837  djac_dparam(local_eqn,local_unknown) -=
2838  dens_ratio*(r*interpolated_u[0]*dpsifdx(l2,0)
2839  + r*psif[l2]*interpolated_dudx(1,1)
2840  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
2841  }
2842 
2843  //Mesh velocity terms
2844  if(!ALE_is_disabled)
2845  {
2846  for(unsigned k=0;k<2;k++)
2847  {
2848  if(diff_re_st)
2849  {
2850  djac_dparam(local_eqn,local_unknown) +=
2851  dens_ratio*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
2852  }
2853  }
2854  }
2855  }
2856 
2857  //There are no azimithal terms in the axial momentum equation
2858  } //End of loop over velocity shape functions
2859 
2860  } /*End of Jacobian calculation*/
2861 
2862  } //End of AXIAL MOMENTUM EQUATION
2863 
2864  //THIRD (AZIMUTHAL) MOMENTUM EQUATION
2865  local_eqn = nodal_local_eqn(l,u_nodal_index[2]);
2866  if(local_eqn >= 0)
2867  {
2868  //Add the user-defined body force terms
2869  //Remember to multiply by the density ratio!
2870  //residuals[local_eqn] +=
2871  // r*body_force[2]*testf[l]*W;
2872 
2873  //Add the gravitational body force term
2874  if(diff_re_inv_fr)
2875  {
2876  dres_dparam[local_eqn] += r*dens_ratio*testf[l]*G[2]*W;
2877  }
2878 
2879  //There is NO pressure gradient term
2880 
2881  //Add in the stress tensor terms
2882  //The viscosity ratio needs to go in here to ensure
2883  //continuity of normal stress is satisfied even in flows
2884  //with zero pressure gradient!
2885  //residuals[local_eqn] -= visc_ratio*
2886  // (r*interpolated_dudx(2,0) -
2887  // Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
2888 
2889  //residuals[local_eqn] -= visc_ratio*r*
2890  // interpolated_dudx(2,1)*dtestfdx(l,1)*W;
2891 
2892  //residuals[local_eqn] -= visc_ratio*
2893  // ((interpolated_u[2]/r) - Gamma[0]*interpolated_dudx(2,0))*testf[l]*W;
2894 
2895 
2896  //Add in the inertial terms
2897  //du/dt term
2898  if(diff_re_st)
2899  {
2900  dres_dparam[local_eqn] -= dens_ratio*r*dudt[2]*testf[l]*W;
2901  }
2902 
2903  //Convective terms
2904  if(diff_re)
2905  {
2906  dres_dparam[local_eqn] -=
2907  dens_ratio*(r*interpolated_u[0]*interpolated_dudx(2,0)
2908  + interpolated_u[0]*interpolated_u[2]
2909  + r*interpolated_u[1]*interpolated_dudx(2,1))*testf[l]*W;
2910  }
2911 
2912  //Mesh velocity terms
2913  if(!ALE_is_disabled)
2914  {
2915  if(diff_re_st)
2916  {
2917  for(unsigned k=0;k<2;k++)
2918  {
2919  dres_dparam[local_eqn] +=
2920  dens_ratio*r*mesh_velocity[k]*interpolated_dudx(2,k)*testf[l]*W;
2921  }
2922  }
2923  }
2924 
2925  //Add the Coriolis term
2926  if(diff_re_inv_ro)
2927  {
2928  dres_dparam[local_eqn] -=
2929  2.0*r*dens_ratio*interpolated_u[0]*testf[l]*W;
2930  }
2931 
2932  //CALCULATE THE JACOBIAN
2933  if(flag)
2934  {
2935  //Loop over the velocity shape functions again
2936  for(unsigned l2=0;l2<n_node;l2++)
2937  {
2938  //Radial velocity component
2939  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
2940  if(local_unknown >= 0)
2941  {
2942  //Convective terms
2943  if(diff_re)
2944  {
2945  djac_dparam(local_eqn,local_unknown) -=
2946  dens_ratio*(r*psif[l2]*interpolated_dudx(2,0)
2947  + psif[l2]*interpolated_u[2])*testf[l]*W;
2948  }
2949 
2950  //Coriolis term
2951  if(diff_re_inv_ro)
2952  {
2953  djac_dparam(local_eqn,local_unknown) -=
2954  2.0*r*dens_ratio*psif[l2]*testf[l]*W;
2955  }
2956  }
2957 
2958  //Axial velocity component
2959  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
2960  if(local_unknown >= 0)
2961  {
2962  //Convective terms
2963  if(diff_re)
2964  {
2965  djac_dparam(local_eqn,local_unknown) -=
2966  dens_ratio*r*psif[l2]*interpolated_dudx(2,1)*testf[l]*W;
2967  }
2968  }
2969 
2970  //Azimuthal velocity component
2971  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
2972  if(local_unknown >= 0)
2973  {
2974  if(flag==2)
2975  {
2976  //Add the mass matrix
2977  if(diff_re_st)
2978  {
2979  dmass_matrix_dparam(local_eqn,local_unknown) +=
2980  dens_ratio*r*psif[l2]*testf[l]*W;
2981  }
2982  }
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  //jacobian(local_eqn,local_unknown) -=
2989  // visc_ratio*(r*dpsifdx(l2,0) -
2990  // Gamma[0]*psif[l2])*dtestfdx(l,0)*W;
2991 
2992  //jacobian(local_eqn,local_unknown) -=
2993  // visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
2994 
2995  //jacobian(local_eqn,local_unknown) -=
2996  // visc_ratio*((psif[l2]/r) - Gamma[0]*dpsifdx(l2,0))
2997  // *testf[l]*W;
2998 
2999  //Add in the inertial terms
3000  //du/dt term
3001  if(diff_re_st)
3002  {
3003  djac_dparam(local_eqn,local_unknown) -=
3004  dens_ratio*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
3005  psif[l2]*testf[l]*W;
3006  }
3007 
3008  //Convective terms
3009  if(diff_re)
3010  {
3011  djac_dparam(local_eqn,local_unknown) -=
3012  dens_ratio*(r*interpolated_u[0]*dpsifdx(l2,0)
3013  + interpolated_u[0]*psif[l2]
3014  + r*interpolated_u[1]*dpsifdx(l2,1))*testf[l]*W;
3015  }
3016 
3017  //Mesh velocity terms
3018  if(!ALE_is_disabled)
3019  {
3020  if(diff_re_st)
3021  {
3022  for(unsigned k=0;k<2;k++)
3023  {
3024  djac_dparam(local_eqn,local_unknown) +=
3025  dens_ratio*r*mesh_velocity[k]*dpsifdx(l2,k)*testf[l]*W;
3026  }
3027  }
3028  }
3029  }
3030  }
3031 
3032  //There are no pressure terms
3033  } //End of Jacobian
3034 
3035  } //End of AZIMUTHAL EQUATION
3036 
3037  } //End of loop over shape functions
3038 
3039 
3040  //CONTINUITY EQUATION NO PARAMETERS
3041  //-------------------
3042  }
3043 
3044 }
3045 
3046 //=========================================================================
3047 /// \short Compute the hessian tensor vector products required to
3048 /// perform continuation of bifurcations analytically
3049 //=========================================================================
3052  Vector<double> const &Y,
3053  DenseMatrix<double> const &C,
3054  DenseMatrix<double> &product)
3055 {
3056  //Find out how many nodes there are
3057  unsigned n_node = nnode();
3058 
3059  //Get the nodal indices at which the velocity is stored
3060  unsigned u_nodal_index[3];
3061  for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}
3062 
3063  //Set up memory for the shape and test functions
3064  //Note that there are only two dimensions, r and z in this problem
3065  Shape psif(n_node), testf(n_node);
3066  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
3067 
3068  //Number of integration points
3069  unsigned Nintpt = integral_pt()->nweight();
3070 
3071  //Set the Vector to hold local coordinates (two dimensions)
3072  Vector<double> s(2);
3073 
3074  //Get Physical Variables from Element
3075  //Reynolds number must be multiplied by the density ratio
3076  double scaled_re = re()*density_ratio();
3077  // double visc_ratio = viscosity_ratio();
3078  Vector<double> G = g();
3079 
3080  //Integers used to store the local equation and unknown numbers
3081  int local_eqn=0, local_unknown=0, local_freedom=0;
3082 
3083  //How may dofs are there
3084  const unsigned n_dof = this->ndof();
3085 
3086  //Create a local matrix eigenvector product contribution
3087  DenseMatrix<double> jac_y(n_dof,n_dof,0.0);
3088 
3089  //Loop over the integration points
3090  for(unsigned ipt=0;ipt<Nintpt;ipt++)
3091  {
3092  //Assign values of s
3093  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
3094  //Get the integral weight
3095  double w = integral_pt()->weight(ipt);
3096 
3097  //Call the derivatives of the shape and test functions
3098  double J =
3099  dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,psif,dpsifdx,testf,dtestfdx);
3100 
3101  //Premultiply the weights and the Jacobian
3102  double W = w*J;
3103 
3104  //Allocate storage for the position and the derivative of the
3105  //mesh positions wrt time
3107  //Vector<double> mesh_velocity(2,0.0);
3108  //Allocate storage for the pressure, velocity components and their
3109  //time and spatial derivatives
3110  Vector<double> interpolated_u(3,0.0);
3111  //Vector<double> dudt(3,0.0);
3112  DenseMatrix<double> interpolated_dudx(3,2,0.0);
3113 
3114  //Calculate velocities and derivatives at integration point
3115 
3116  // Loop over nodes
3117  for(unsigned l=0;l<n_node;l++)
3118  {
3119  //Cache the shape function
3120  const double psif_ = psif(l);
3121  //Loop over the two coordinate directions
3122  for(unsigned i=0;i<2;i++)
3123  {
3124  interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
3125  }
3126 
3127  //Loop over the three velocity directions
3128  for(unsigned i=0;i<3;i++)
3129  {
3130  //Get the u_value
3131  const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
3132  interpolated_u[i] += u_value*psif_;
3133  //dudt[i]+= du_dt_axi_nst(l,i)*psif_;
3134  //Loop over derivative directions
3135  for(unsigned j=0;j<2;j++)
3136  {interpolated_dudx(i,j) += u_value*dpsifdx(l,j);}
3137  }
3138  }
3139 
3140  //Get the mesh velocity if ALE is enabled
3141  if(!ALE_is_disabled)
3142  {
3143  throw OomphLibError("Moving nodes not implemented\n",
3144  OOMPH_CURRENT_FUNCTION,
3145  OOMPH_EXCEPTION_LOCATION);
3146  }
3147 
3148  //r is the first position component
3149  double r = interpolated_x[0];
3150 
3151 
3152  //MOMENTUM EQUATIONS
3153  //------------------
3154 
3155  //Loop over the test functions
3156  for(unsigned l=0;l<n_node;l++)
3157  {
3158 
3159  //FIRST (RADIAL) MOMENTUM EQUATION
3160  local_eqn = nodal_local_eqn(l,u_nodal_index[0]);
3161  //If it's not a boundary condition
3162  if(local_eqn >= 0)
3163  {
3164  //Loop over the velocity shape functions yet again
3165  for(unsigned l3=0;l3<n_node;l3++)
3166  {
3167  //Derivative of jacobian terms with respect to radial velocity
3168  local_freedom = nodal_local_eqn(l3,u_nodal_index[0]);
3169  if(local_freedom >= 0)
3170  {
3171  //Storage for the sums
3172  double temp = 0.0;
3173 
3174  //Loop over the velocity shape functions again
3175  for(unsigned l2=0;l2<n_node;l2++)
3176  {
3177  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
3178  //Radial velocity component
3179  if(local_unknown >= 0)
3180  {
3181  //Remains of convective terms
3182  temp -=scaled_re*(r*psif[l2]*dpsifdx(l3,0)
3183  + r*psif[l3]*dpsifdx(l2,0))*
3184  Y[local_unknown]*testf[l]*W;
3185  }
3186 
3187  //Axial velocity component
3188  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
3189  if(local_unknown >= 0)
3190  {
3191  //Convective terms
3192  temp -=
3193  scaled_re*r*psif[l2]*dpsifdx(l3,1)*Y[local_unknown]*testf[l]*W;
3194  }
3195  }
3196  //Add the summed contribution to the product matrix
3197  jac_y(local_eqn,local_freedom) += temp;
3198  } //End of derivative wrt radial coordinate
3199 
3200 
3201  //Derivative of jacobian terms with respect to axial velocity
3202  local_freedom = nodal_local_eqn(l3,u_nodal_index[1]);
3203  if(local_freedom >= 0)
3204  {
3205  double temp = 0.0;
3206 
3207  //Loop over the velocity shape functions again
3208  for(unsigned l2=0;l2<n_node;l2++)
3209  {
3210  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
3211  //Radial velocity component
3212  if(local_unknown >= 0)
3213  {
3214  //Remains of convective terms
3215  temp -= scaled_re*(r*psif[l3]*dpsifdx(l2,1))*
3216  Y[local_unknown]*testf[l]*W;
3217  }
3218  }
3219  //Add the summed contribution to the product matrix
3220  jac_y(local_eqn,local_freedom) += temp;
3221  } //End of derivative wrt axial coordinate
3222 
3223  //Derivative of jacobian terms with respect to azimuthal velocity
3224  local_freedom = nodal_local_eqn(l3,u_nodal_index[2]);
3225  if(local_freedom >= 0)
3226  {
3227  double temp = 0.0;
3228 
3229  //Loop over the velocity shape functions again
3230  for(unsigned l2=0;l2<n_node;l2++)
3231  {
3232  //Azimuthal velocity component
3233  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
3234  if(local_unknown >= 0)
3235  {
3236  //Convective terms
3237  temp -=
3238  - scaled_re*2.0*psif[l3]*psif[l2]*Y[local_unknown]*testf[l]*W;
3239  }
3240  }
3241  //Add the summed contibution
3242  jac_y(local_eqn,local_freedom) += temp;
3243 
3244  } //End of if not boundary condition statement
3245  } //End of loop over freedoms
3246  } //End of RADIAL MOMENTUM EQUATION
3247 
3248 
3249  //SECOND (AXIAL) MOMENTUM EQUATION
3250  local_eqn = nodal_local_eqn(l,u_nodal_index[1]);
3251  //If it's not a boundary condition
3252  if(local_eqn >= 0)
3253  {
3254  //Loop over the velocity shape functions yet again
3255  for(unsigned l3=0;l3<n_node;l3++)
3256  {
3257  //Derivative of jacobian terms with respect to radial velocity
3258  local_freedom = nodal_local_eqn(l3,u_nodal_index[0]);
3259  if(local_freedom >= 0)
3260  {
3261  double temp = 0.0;
3262 
3263  //Loop over the velocity shape functions again
3264  for(unsigned l2=0;l2<n_node;l2++)
3265  {
3266  //Axial velocity component
3267  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
3268  if(local_unknown >= 0)
3269  {
3270  //Convective terms
3271  temp -=
3272  scaled_re*(r*psif[l3]*dpsifdx(l2,0))*
3273  Y[local_unknown]*testf[l]*W;
3274  }
3275  }
3276  jac_y(local_eqn,local_freedom) += temp;
3277 
3278  //There are no azimithal terms in the axial momentum equation
3279  } //End of loop over velocity shape functions
3280 
3281 
3282  //Derivative of jacobian terms with respect to axial velocity
3283  local_freedom = nodal_local_eqn(l3,u_nodal_index[1]);
3284  if(local_freedom >= 0)
3285  {
3286  double temp = 0.0;
3287 
3288  //Loop over the velocity shape functions again
3289  for(unsigned l2=0;l2<n_node;l2++)
3290  {
3291  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
3292  //Radial velocity component
3293  if(local_unknown >= 0)
3294  {
3295  //Convective terms
3296  temp -=
3297  scaled_re*r*psif[l2]*dpsifdx(l3,0)*Y[local_unknown]*testf[l]*W;
3298  }
3299 
3300  //Axial velocity component
3301  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
3302  if(local_unknown >= 0)
3303  {
3304  //Convective terms
3305  temp -=
3306  scaled_re*(
3307  r*psif[l2]*dpsifdx(l3,1)
3308  + r*psif[l3]*dpsifdx(l2,1))*Y[local_unknown]*testf[l]*W;
3309  }
3310 
3311  //There are no azimithal terms in the axial momentum equation
3312  } //End of loop over velocity shape functions
3313 
3314  //Add summed contributiont to jacobian product matrix
3315  jac_y(local_eqn,local_freedom) += temp;
3316  }
3317  } //End of loop over local freedoms
3318 
3319  } //End of AXIAL MOMENTUM EQUATION
3320 
3321  //THIRD (AZIMUTHAL) MOMENTUM EQUATION
3322  local_eqn = nodal_local_eqn(l,u_nodal_index[2]);
3323  if(local_eqn >= 0)
3324  {
3325  //Loop over the velocity shape functions yet again
3326  for(unsigned l3=0;l3<n_node;l3++)
3327  {
3328  //Derivative of jacobian terms with respect to radial velocity
3329  local_freedom = nodal_local_eqn(l3,u_nodal_index[0]);
3330  if(local_freedom >= 0)
3331  {
3332  double temp = 0.0;
3333 
3334  //Loop over the velocity shape functions again
3335  for(unsigned l2=0;l2<n_node;l2++)
3336  {
3337  //Azimuthal velocity component
3338  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
3339  if(local_unknown >= 0)
3340  {
3341  //Convective terms
3342  temp -=
3343  scaled_re*(r*psif[l3]*dpsifdx(l2,0)
3344  + psif[l3]*psif[l2])*Y[local_unknown]*testf[l]*W;
3345  }
3346  }
3347  //Add the summed contribution to the jacobian eigenvector sum
3348  jac_y(local_eqn,local_freedom) += temp;
3349  }
3350 
3351  //Derivative of jacobian terms with respect to axial velocity
3352  local_freedom = nodal_local_eqn(l3,u_nodal_index[1]);
3353  if(local_freedom >= 0)
3354  {
3355  double temp = 0.0;
3356 
3357  //Loop over the velocity shape functions again
3358  for(unsigned l2=0;l2<n_node;l2++)
3359  {
3360  //Azimuthal velocity component
3361  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
3362  if(local_unknown >= 0)
3363  {
3364  //Convective terms
3365  temp -=
3366  scaled_re*(r*psif[l3]*dpsifdx(l2,1))*
3367  Y[local_unknown]*testf[l]*W;
3368  }
3369  }
3370  //Add the summed contribution to the jacobian eigenvector sum
3371  jac_y(local_eqn,local_freedom) += temp;
3372  }
3373 
3374 
3375  //Derivative of jacobian terms with respect to azimuthal velocity
3376  local_freedom = nodal_local_eqn(l3,u_nodal_index[2]);
3377  if(local_freedom >= 0)
3378  {
3379  double temp = 0.0;
3380 
3381 
3382  //Loop over the velocity shape functions again
3383  for(unsigned l2=0;l2<n_node;l2++)
3384  {
3385  //Radial velocity component
3386  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
3387  if(local_unknown >= 0)
3388  {
3389  //Convective terms
3390  temp -=
3391  scaled_re*(r*psif[l2]*dpsifdx(l3,0)
3392  + psif[l2]*psif[l3])*Y[local_unknown]*testf[l]*W;
3393  }
3394 
3395  //Axial velocity component
3396  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
3397  if(local_unknown >= 0)
3398  {
3399  //Convective terms
3400  temp -=
3401  scaled_re*r*psif[l2]*dpsifdx(l3,1)*Y[local_unknown]*testf[l]*W;
3402  }
3403  }
3404  //Add the fully summed contribution
3405  jac_y(local_eqn,local_freedom) += temp;
3406  }
3407  } //End of loop over freedoms
3408 
3409  //There are no pressure terms
3410  } //End of AZIMUTHAL EQUATION
3411 
3412  } //End of loop over shape functions
3413  }
3414 
3415  //We have now assembled the matrix (d J_{ij} Y_j)/d u_{k}
3416  //and simply need to sum over the vectors
3417  const unsigned n_vec = C.nrow();
3418  for(unsigned i=0;i<n_dof;i++)
3419  {
3420  for(unsigned k=0;k<n_dof;k++)
3421  {
3422  //Cache the value of the hessian y product
3423  const double j_y = jac_y(i,k);
3424  //Loop over the possible vectors
3425  for(unsigned v=0;v<n_vec;v++)
3426  {
3427  product(v,i) += j_y*C(v,k);
3428  }
3429  }
3430  }
3431 }
3432 
3433 
3434 //////////////////////////////////////////////////////////////////////
3435 //////////////////////////////////////////////////////////////////////
3436 //////////////////////////////////////////////////////////////////////
3437 
3438 //=============================================================================
3439 /// Create a list of pairs for all unknowns in this element,
3440 /// so that the first entry in each pair contains the global equation
3441 /// number of the unknown, while the second one contains the number
3442 /// of the "DOF type" that this unknown is associated with.
3443 /// (Function can obviously only be called if the equation numbering
3444 /// scheme has been set up.)
3445 //=============================================================================
3447  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
3448 {
3449  // number of nodes
3450  unsigned n_node = this->nnode();
3451 
3452  // number of pressure values
3453  unsigned n_press = this->npres_axi_nst();
3454 
3455  // temporary pair (used to store dof lookup prior to being added to list)
3456  std::pair<unsigned,unsigned> dof_lookup;
3457 
3458  // pressure dof number (is this really OK?)
3459  unsigned pressure_dof_number = 3;
3460 
3461  // loop over the pressure values
3462  for (unsigned n = 0; n < n_press; n++)
3463  {
3464  // determine local eqn number
3465  int local_eqn_number = this->p_local_eqn(n);
3466 
3467  // ignore pinned values - far away degrees of freedom resulting
3468  // from hanging nodes can be ignored since these are be dealt
3469  // with by the element containing their master nodes
3470  if (local_eqn_number >= 0)
3471  {
3472  // store dof lookup in temporary pair: First entry in pair
3473  // is global equation number; second entry is dof type
3474  dof_lookup.first = this->eqn_number(local_eqn_number);
3475  dof_lookup.second = pressure_dof_number;
3476 
3477  // add to list
3478  dof_lookup_list.push_front(dof_lookup);
3479  }
3480  }
3481 
3482  // loop over the nodes
3483  for (unsigned n = 0; n < n_node; n++)
3484  {
3485  // find the number of values at this node
3486  unsigned nv = this->node_pt(n)->nvalue();
3487 
3488  //loop over these values
3489  for (unsigned v = 0; v < nv; v++)
3490  {
3491  // determine local eqn number
3492  int local_eqn_number = this->nodal_local_eqn(n, v);
3493 
3494  // ignore pinned values
3495  if (local_eqn_number >= 0)
3496  {
3497  // store dof lookup in temporary pair: First entry in pair
3498  // is global equation number; second entry is dof type
3499  dof_lookup.first = this->eqn_number(local_eqn_number);
3500  dof_lookup.second = v;
3501 
3502  // add to list
3503  dof_lookup_list.push_front(dof_lookup);
3504 
3505  }
3506  }
3507  }
3508 }
3509 
3510 
3511 
3512 ///Axisymmetric Crouzeix-Raviart elements
3513 //Set the data for the number of Variables at each node
3515 Initial_Nvalue[9]={3,3,3,3,3,3,3,3,3};
3516 
3517 //========================================================================
3518 /// Number of values (pinned or dofs) required at node n.
3519 //========================================================================
3521 required_nvalue(const unsigned &n) const {return Initial_Nvalue[n];}
3522 
3523 //========================================================================
3524 /// Compute traction at local coordinate s for outer unit normal N
3525 //========================================================================
3527  const Vector<double>& N,
3528  Vector<double>& traction)
3529  const
3530 {
3532 }
3533 
3534 ///////////////////////////////////////////////////////////////////////////
3535 ///////////////////////////////////////////////////////////////////////////
3536 ///////////////////////////////////////////////////////////////////////////
3537 
3538 
3539 
3540 //============================================================================
3541 /// Create a list of pairs for all unknowns in this element,
3542 /// so the first entry in each pair contains the global equation
3543 /// number of the unknown, while the second one contains the number
3544 /// of the "DOF type" that this unknown is associated with.
3545 /// (Function can obviously only be called if the equation numbering
3546 /// scheme has been set up.)
3547 //============================================================================
3549  std::list<std::pair<unsigned long,
3550  unsigned> >& dof_lookup_list) const
3551 {
3552  // number of nodes
3553  unsigned n_node = this->nnode();
3554 
3555  // temporary pair (used to store dof lookup prior to being added to list)
3556  std::pair<unsigned,unsigned> dof_lookup;
3557 
3558  // loop over the nodes
3559  for (unsigned n = 0; n < n_node; n++)
3560  {
3561  // find the number of values at this node
3562  unsigned nv = this->node_pt(n)->nvalue();
3563 
3564  //loop over these values
3565  for (unsigned v = 0; v < nv; v++)
3566  {
3567  // determine local eqn number
3568  int local_eqn_number = this->nodal_local_eqn(n, v);
3569 
3570  // ignore pinned values - far away degrees of freedom resulting
3571  // from hanging nodes can be ignored since these are be dealt
3572  // with by the element containing their master nodes
3573  if (local_eqn_number >= 0)
3574  {
3575  // store dof lookup in temporary pair: Global equation number
3576  // is the first entry in pair
3577  dof_lookup.first = this->eqn_number(local_eqn_number);
3578 
3579  // set dof numbers: Dof number is the second entry in pair
3580  dof_lookup.second = v;
3581 
3582  // add to list
3583  dof_lookup_list.push_front(dof_lookup);
3584  }
3585  }
3586  }
3587 }
3588 
3589 
3590 
3591 //Axisymmetric Taylor--Hood
3592 //Set the data for the number of Variables at each node
3593 const unsigned AxisymmetricQTaylorHoodElement::
3594 Initial_Nvalue[9]={4,3,4,3,3,3,4,3,4};
3595 
3596 //Set the data for the pressure conversion array
3597 const unsigned AxisymmetricQTaylorHoodElement::Pconv[4]={0,2,6,8};
3598 
3599 
3600 //========================================================================
3601 /// Compute traction at local coordinate s for outer unit normal N
3602 //========================================================================
3604  const Vector<double>& N,
3605  Vector<double>& traction)
3606  const
3607 {
3609 }
3610 
3611 }
virtual void pshape_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
double *& re_invro_pt()
Pointer to global inverse Froude number.
virtual unsigned u_index_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.
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...
double *& re_invfr_pt()
Pointer to global inverse Froude number.
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.
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...
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 int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
double get_source_fct(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
const double & re_invfr() const
Global inverse Froude number.
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 ...
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
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...
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 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
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
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
unsigned npres_axi_nst() const
Return number of pressure values.
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
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
double kin_energy() const
Get integral of kinetic energy over element.
void interpolated_u_axi_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
virtual void get_source_fct_gradient(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
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
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.
double interpolated_p_axi_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
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 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
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...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed...
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
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers...
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...
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
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 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...
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
e
Definition: cfortran.h:575
int p_local_eqn(const unsigned &n) const
Overload the access function for the pressure's local equation numbers.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
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
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:504
bool has_auxiliary_node_update_fct_pt()
Boolean to indicate if node has a pointer to and auxiliary update function.
Definition: nodes.h:1500
virtual void 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.
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
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...
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 & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
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
virtual unsigned npres_axi_nst() const =0
Function to return number of pressure degrees of freedom.
A Rank 3 Tensor class.
Definition: matrices.h:1337
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...
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.
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
double pressure_integral() const
Integral of pressure 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 * Re_pt
Pointer to global Reynolds number.
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}.
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(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
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
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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
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
const double & re() const
Reynolds number.
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
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
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.
double *& re_pt()
Pointer to Reynolds number.
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
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
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed...
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
const Vector< double > & g() const
Vector of gravitational components.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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
double dissipation() const
Return integral of dissipation over element.
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
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
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...
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
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...
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164
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 double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)