polar_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//====================================================================
31 
32 namespace oomph
33 {
34 
35 //////////////////////////////////////////////////////////////////////////
36 //======================================================================//
37 /// Start of what would've been navier_stokes_elements.cc //
38 //======================================================================//
39 //////////////////////////////////////////////////////////////////////////
40 
41 /// Navier--Stokes equations static data
42 Vector<double> PolarNavierStokesEquations::Gamma(2,0.0);
43 
44 //=================================================================
45 /// "Magic" negative number that indicates that the pressure is
46 /// not stored at a node. This cannot be -1 because that represents
47 /// the positional hanging scheme in the hanging_pt object of nodes
48 //=================================================================
50 
51 /// Navier--Stokes equations static data
53 
54 /// Navier--Stokes equations static data
56 
57 /// Navier-Stokes equations default gravity vector
59 
60 //======================================================================
61 /// Validate against exact velocity solution at given time.
62 /// Solution is provided via function pointer.
63 /// Plot at a given number of plot points and compute L2 error
64 /// and L2 norm of velocity solution over element.
65 //=======================================================================
67 compute_error(std::ostream &outfile,
69  const double& time,
70  double& error, double& norm)
71 {
72  error=0.0;
73  norm=0.0;
74 
75  //Vector of local coordinates
76  Vector<double> s(2);
77 
78  // Vector for coordintes
79  Vector<double> x(2);
80 
81  //Set the value of n_intpt
82  unsigned n_intpt = integral_pt()->nweight();
83 
84  outfile << "ZONE" << std::endl;
85 
86  // Exact solution Vector (u,v,[w],p)
87  Vector<double> exact_soln(2+1);
88 
89  //Loop over the integration points
90  for(unsigned ipt=0;ipt<n_intpt;ipt++)
91  {
92 
93  //Assign values of s
94  for(unsigned i=0;i<2;i++)
95  {
96  s[i] = integral_pt()->knot(ipt,i);
97  }
98 
99  //Get the integral weight
100  double w = integral_pt()->weight(ipt);
101 
102  // Get jacobian of mapping
103  double J=J_eulerian(s);
104 
105  //Premultiply the weights and the Jacobian
106  double W = w*J;
107 
108  // Get x position as Vector
109  interpolated_x(s,x);
110 
111  // Get exact solution at this point
112  (*exact_soln_pt)(time,x,exact_soln);
113 
114  // Velocity error
115  for(unsigned i=0;i<2;i++)
116  {
117  norm+=exact_soln[i]*exact_soln[i]*W;
118  error+=(exact_soln[i]-interpolated_u_pnst(s,i))*
119  (exact_soln[i]-interpolated_u_pnst(s,i))*W;
120  }
121 
122  //Output x,y,...,u_exact
123  for(unsigned i=0;i<2;i++)
124  {
125  outfile << x[i] << " ";
126  }
127 
128  //Output x,y,[z],u_error,v_error,[w_error]
129  for(unsigned i=0;i<2;i++)
130  {
131  outfile << exact_soln[i]-interpolated_u_pnst(s,i) << " ";
132  }
133 
134  outfile << std::endl;
135 
136  }
137 }
138 
139 //======================================================================
140 /// Validate against exact velocity solution
141 /// Solution is provided via function pointer.
142 /// Plot at a given number of plot points and compute L2 error
143 /// and L2 norm of velocity solution over element.
144 //=======================================================================
145 void PolarNavierStokesEquations::compute_error(std::ostream &outfile,
147  double& error, double& norm)
148 {
149 
150  error=0.0;
151  norm=0.0;
152 
153  //Vector of local coordinates
154  Vector<double> s(2);
155 
156  // Vector for coordintes
157  Vector<double> x(2);
158 
159  //Set the value of n_intpt
160  unsigned n_intpt = integral_pt()->nweight();
161 
162  outfile << "ZONE" << std::endl;
163 
164  // Exact solution Vector (u,v,[w],p)
165  Vector<double> exact_soln(2+1);
166 
167  //Loop over the integration points
168  for(unsigned ipt=0;ipt<n_intpt;ipt++)
169  {
170 
171  //Assign values of s
172  for(unsigned i=0;i<2;i++)
173  {
174  s[i] = integral_pt()->knot(ipt,i);
175  }
176 
177  //Get the integral weight
178  double w = integral_pt()->weight(ipt);
179 
180  // Get jacobian of mapping
181  double J=J_eulerian(s);
182 
183  //Premultiply the weights and the Jacobian
184  double W = w*J;
185 
186  // Get x position as Vector
187  interpolated_x(s,x);
188 
189  // Get exact solution at this point
190  (*exact_soln_pt)(x,exact_soln);
191 
192  // Velocity error
193  for(unsigned i=0;i<2;i++)
194  {
195  norm+=exact_soln[i]*exact_soln[i]*W;
196  error+=(exact_soln[i]-interpolated_u_pnst(s,i))*
197  (exact_soln[i]-interpolated_u_pnst(s,i))*W;
198  }
199 
200  //Output x,y,...,u_exact
201  for(unsigned i=0;i<2;i++)
202  {
203  outfile << x[i] << " ";
204  }
205 
206  //Output x,y,[z],u_error,v_error,[w_error]
207  for(unsigned i=0;i<2;i++)
208  {
209  outfile << exact_soln[i]-interpolated_u_pnst(s,i) << " ";
210  }
211 
212  outfile << std::endl;
213 
214  }
215 }
216 
217 //======================================================================
218 /// Output "exact" solution
219 /// Solution is provided via function pointer.
220 /// Plot at a given number of plot points.
221 /// Function prints as many components as are returned in solution Vector.
222 //=======================================================================
223 void PolarNavierStokesEquations::output_fct(std::ostream &outfile,
224  const unsigned &nplot,
226 {
227 
228  //Vector of local coordinates
229  Vector<double> s(2);
230 
231  // Vector for coordintes
232  Vector<double> x(2);
233 
234  // Tecplot header info
235  outfile << tecplot_zone_string(nplot);
236 
237  // Exact solution Vector
238  Vector<double> exact_soln;
239 
240  // Loop over plot points
241  unsigned num_plot_points=nplot_points(nplot);
242  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
243  {
244 
245  // Get local coordinates of plot point
246  get_s_plot(iplot,nplot,s);
247 
248  // Get x position as Vector
249  interpolated_x(s,x);
250 
251  // Get exact solution at this point
252  (*exact_soln_pt)(x,exact_soln);
253 
254  //Output x,y,...
255  for(unsigned i=0;i<2;i++)
256  {
257  outfile << x[i] << " ";
258  }
259 
260  //Output "exact solution"
261  for(unsigned i=0;i<exact_soln.size();i++)
262  {
263  outfile << exact_soln[i] << " ";
264  }
265 
266  outfile << std::endl;
267  }
268 
269 // Write tecplot footer (e.g. FE connectivity lists)
270  write_tecplot_zone_footer(outfile,nplot);
271 }
272 
273 //======================================================================
274 /// Output "exact" solution at a given time
275 /// Solution is provided via function pointer.
276 /// Plot at a given number of plot points.
277 /// Function prints as many components as are returned in solution Vector.
278 //=======================================================================
279 void PolarNavierStokesEquations::output_fct(std::ostream &outfile,
280  const unsigned &nplot,
281  const double& time,
283 {
284 
285  //Vector of local coordinates
286  Vector<double> s(2);
287 
288  // Vector for coordintes
289  Vector<double> x(2);
290 
291  // Tecplot header info
292  outfile << tecplot_zone_string(nplot);
293 
294  // Exact solution Vector
295  Vector<double> exact_soln;
296 
297  // Loop over plot points
298  unsigned num_plot_points=nplot_points(nplot);
299  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
300  {
301 
302  // Get local coordinates of plot point
303  get_s_plot(iplot,nplot,s);
304 
305  // Get x position as Vector
306  interpolated_x(s,x);
307 
308  // Get exact solution at this point
309  (*exact_soln_pt)(time,x,exact_soln);
310 
311  //Output x,y,...
312  for(unsigned i=0;i<2;i++)
313  {
314  outfile << x[i] << " ";
315  }
316 
317  //Output "exact solution"
318  for(unsigned i=0;i<exact_soln.size();i++)
319  {
320  outfile << exact_soln[i] << " ";
321  }
322 
323  outfile << std::endl;
324 
325  }
326 
327  // Write tecplot footer (e.g. FE connectivity lists)
328  write_tecplot_zone_footer(outfile,nplot);
329 
330 
331 }
332 
333 //==============================================================
334 /// Output function: Velocities only
335 /// x,y,[z],u,v,[w]
336 /// in tecplot format at specified previous timestep (t=0: present;
337 /// t>0: previous timestep). Specified number of plot points in each
338 /// coordinate direction.
339 //==============================================================
340 void PolarNavierStokesEquations::output_veloc(std::ostream& outfile,
341  const unsigned& nplot,
342  const unsigned& t)
343 {
344 
345  //Find number of nodes
346  unsigned n_node = nnode();
347 
348  //Local shape function
349  Shape psi(n_node);
350 
351  //Vectors of local coordinates and coords and velocities
352  Vector<double> s(2);
354  Vector<double> interpolated_u(2);
355 
356  // Tecplot header info
357  outfile << tecplot_zone_string(nplot);
358 
359  // Loop over plot points
360  unsigned num_plot_points=nplot_points(nplot);
361  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
362  {
363 
364  // Get local coordinates of plot point
365  get_s_plot(iplot,nplot,s);
366 
367  // Get shape functions
368  shape(s,psi);
369 
370  // Loop over directions
371  for(unsigned i=0;i<2;i++)
372  {
373  interpolated_x[i]=0.0;
374  interpolated_u[i]=0.0;
375  //Loop over the local nodes and sum
376  for(unsigned l=0;l<n_node;l++)
377  {
378  interpolated_u[i] += u_pnst(t,l,i)*psi[l];
379  interpolated_x[i] += nodal_position(t,l,i)*psi[l];
380  }
381  }
382 
383  // Coordinates
384  for(unsigned i=0;i<2;i++)
385  {
386  outfile << interpolated_x[i] << " ";
387  }
388 
389  // Velocities
390  for(unsigned i=0;i<2;i++)
391  {
392  outfile << interpolated_u[i] << " ";
393  }
394 
395  outfile << std::endl;
396  }
397 // Write tecplot footer (e.g. FE connectivity lists)
398  write_tecplot_zone_footer(outfile,nplot);
399 
400 
401 }
402 
403 // keyword: primary
404 //==============================================================
405 /// Output function:
406 /// x,y,[z],u,v,[w],p
407 /// in tecplot format. Specified number of plot points in each
408 /// coordinate direction.
409 //==============================================================
410 void PolarNavierStokesEquations::output(std::ostream &outfile,
411  const unsigned &nplot)
412 {
413 
414  //Vector of local coordinates
415  Vector<double> s(2);
416 
417  // Tecplot header info
418  outfile << tecplot_zone_string(nplot);
419 
420  // Loop over plot points
421  unsigned num_plot_points=nplot_points(nplot);
422  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
423  {
424 
425  // Get local coordinates of plot point
426  get_s_plot(iplot,nplot,s);
427 
428  //Work out global physical coordinate
429  const double Alpha = alpha();
430  double r = interpolated_x(s,0);
431  double phi = interpolated_x(s,1);
432  double theta = Alpha*phi;
433 
434  // Coordinates
435  outfile << r*cos(theta) << " " << r*sin(theta) << " ";
436 
437  // Velocities
438  outfile << interpolated_u_pnst(s,0)*cos(theta) - interpolated_u_pnst(s,1)*sin(theta)
439  << " ";
440  outfile << interpolated_u_pnst(s,0)*sin(theta) + interpolated_u_pnst(s,1)*cos(theta)
441  << " ";
442 
443  // Pressure
444  outfile << interpolated_p_pnst(s) << " ";
445 
446  // Radial and Azimuthal velocities
447  outfile << interpolated_u_pnst(s,0) << " " << interpolated_u_pnst(s,1) << " ";
448  //comment start here
449  /*
450  double similarity_solution,dsimilarity_solution;
451  get_similarity_solution(theta,Alpha,similarity_solution,dsimilarity_solution);
452  // Similarity solution:
453  outfile << similarity_solution/(Alpha*r) << " ";
454  // Error from similarity solution:
455  outfile << interpolated_u_pnst(s,0) - similarity_solution/(Alpha*r) << " ";
456  */
457  /*
458  //Work out the Stokes flow similarity solution
459  double mult = (Alpha/(sin(2.*Alpha)-2.*Alpha*cos(2.*Alpha)));
460  double similarity_solution = mult*(cos(2.*Alpha*phi)-cos(2.*Alpha));
461  // Similarity solution:
462  outfile << similarity_solution/(Alpha*r) << " ";
463  // Error from similarity solution:
464  outfile << interpolated_u_pnst(s,0) - similarity_solution/(Alpha*r) << " ";
465  //comment end here
466  */
467  outfile << 0 << " " << 1 << " ";
468 
469  // r and theta for better plotting
470  outfile << r << " " << phi << " ";
471 
472  outfile << std::endl;
473  }
474  outfile << std::endl;
475 
476 }
477 
478 /*
479 //Full_convergence_checks output:
480 //==============================================================
481 /// Output the exact similarity solution and error
482 /// Output function in tecplot format.
483 /// Specified number of plot points in each
484 /// coordinate direction.
485 //==============================================================
486 void PolarNavierStokesEquations::output(std::ostream &outfile,
487  const unsigned &nplot)
488 {
489 
490  //Vector of local coordinates
491  Vector<double> s(2);
492 
493  // Tecplot header info
494  outfile << tecplot_zone_string(nplot);
495 
496  // Loop over plot points
497  unsigned num_plot_points=nplot_points(nplot);
498  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
499  {
500 
501  // Get local coordinates of plot point
502  get_s_plot(iplot,nplot,s);
503 
504  //Work out global physical coordinate
505  const double Alpha = alpha();
506  double r = interpolated_x(s,0);
507  double phi = interpolated_x(s,1);
508  double theta = Alpha*phi;
509 
510  // Coordinates
511  outfile << r*cos(theta) << " " << r*sin(theta) << " ";
512 
513  // Velocities
514  outfile << interpolated_u_pnst(s,0)*cos(theta) - interpolated_u_pnst(s,1)*sin(theta)
515  << " ";
516  outfile << interpolated_u_pnst(s,0)*sin(theta) + interpolated_u_pnst(s,1)*cos(theta)
517  << " ";
518 
519  // Pressure
520  outfile << interpolated_p_pnst(s) << " ";
521 
522  // Radial and Azimuthal velocities
523  outfile << interpolated_u_pnst(s,0) << " " << interpolated_u_pnst(s,1) << " ";
524 
525  // I need to add exact pressure, radial velocity and radial velocity derivatives here
526  // Plus the error from these
527  double m=2.*Alpha;
528  double mu=(1./(sin(m)-m*cos(m)));
529  double exact_u=(mu/r)*(cos(m*phi)-cos(m));
530  double exact_p=(2.*mu)*((cos(m*phi)/(r*r))-cos(m));
531  double exact_dudr=-(exact_u/r);
532  double exact_dudphi=-mu*m*sin(m*phi)/r;
533 
534  // If we don't have Stokes flow then we need to overwrite the Stokes solution by
535  // reading in the correct similarity solution from file
536  const double Re = re();
537  if(Re>1.e-8)
538  {
539  double similarity_solution,dsimilarity_solution;
540  get_similarity_solution(theta,Alpha,similarity_solution,dsimilarity_solution);
541  double A = Global_Physical_Variables::P2/Alpha;
542 
543  exact_u = similarity_solution/(r*Alpha);
544  exact_p = 2.*(exact_u/r)-(A/(2.*Alpha*Alpha))*((1./(r*r))-1.);
545  exact_dudr = -(exact_u/r);
546  exact_dudphi = dsimilarity_solution/(r*Alpha);
547  }
548 
549  // exact pressure and error
550  outfile << exact_p << " " << (interpolated_p_pnst(s)-exact_p) << " ";
551 
552  // r and theta for better plotting
553  outfile << r << " " << phi << " ";
554 
555  // exact and oomph du/dr and du/dphi?
556  outfile << exact_u << " " << (interpolated_u_pnst(s,0)-exact_u) << " ";
557  outfile << exact_dudr << " " << (interpolated_dudx_pnst(s,0,0)-exact_dudr) << " ";
558  outfile << exact_dudphi << " " << (interpolated_dudx_pnst(s,0,1)-exact_dudphi) << " ";
559 
560  outfile << std::endl;
561  }
562  outfile << std::endl;
563 }
564 */
565 //==============================================================
566 /// C-style output function:
567 /// x,y,[z],u,v,[w],p
568 /// in tecplot format. Specified number of plot points in each
569 /// coordinate direction.
570 //==============================================================
572  const unsigned &nplot)
573 {
574 
575  //Vector of local coordinates
576  Vector<double> s(2);
577 
578  // Tecplot header info
579  //outfile << tecplot_zone_string(nplot);
580  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
581 
582  // Loop over plot points
583  unsigned num_plot_points=nplot_points(nplot);
584  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
585  {
586 
587  // Get local coordinates of plot point
588  get_s_plot(iplot,nplot,s);
589 
590  // Coordinates
591  for(unsigned i=0;i<2;i++)
592  {
593  //outfile << interpolated_x(s,i) << " ";
594  fprintf(file_pt,"%g ",interpolated_x(s,i));
595  }
596 
597  // Velocities
598  for(unsigned i=0;i<2;i++)
599  {
600  //outfile << interpolated_u_pnst(s,i) << " ";
601  fprintf(file_pt,"%g ",interpolated_u_pnst(s,i));
602  }
603 
604  // Pressure
605  //outfile << interpolated_p_pnst(s) << " ";
606  fprintf(file_pt,"%g \n",interpolated_p_pnst(s));
607  }
608  fprintf(file_pt,"\n");
609 
610 }
611 
612 
613 //==============================================================
614 /// Full output function:
615 /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
616 /// in tecplot format. Specified number of plot points in each
617 /// coordinate direction
618 //==============================================================
619 void PolarNavierStokesEquations::full_output(std::ostream &outfile,
620  const unsigned &nplot)
621 {
622 
623  //Vector of local coordinates
624  Vector<double> s(2);
625 
626  // Acceleration
627  Vector<double> dudt(2);
628 
629  // Mesh elocity
630  Vector<double> mesh_veloc(2);
631 
632  // Tecplot header info
633  outfile << tecplot_zone_string(nplot);
634 
635  //Find out how many nodes there are
636  unsigned n_node = nnode();
637 
638  //Set up memory for the shape functions
639  Shape psif(n_node);
640  DShape dpsifdx(n_node,2);
641 
642  // Loop over plot points
643  unsigned num_plot_points=nplot_points(nplot);
644  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
645  {
646 
647  // Get local coordinates of plot point
648  get_s_plot(iplot,nplot,s);
649 
650  //Call the derivatives of the shape and test functions
651  dshape_eulerian(s,psif,dpsifdx);
652 
653  //Allocate storage
654  Vector<double> mesh_veloc(2);
655  Vector<double> dudt(2);
656  Vector<double> dudt_ALE(2);
657  DenseMatrix<double> interpolated_dudx(2,2);
658  //DenseDoubleMatrix interpolated_dudx(2,2);
659 
660  //Initialise everything to zero
661  for(unsigned i=0;i<2;i++)
662  {
663  mesh_veloc[i]=0.0;
664  dudt[i]=0.0;
665  dudt_ALE[i]=0.0;
666  for(unsigned j=0;j<2;j++)
667  {
668  interpolated_dudx(i,j) = 0.0;
669  }
670  }
671 
672  //Calculate velocities and derivatives
673 
674  // Loop over nodes
675  for(unsigned l=0;l<n_node;l++)
676  {
677  //Loop over directions
678  for(unsigned i=0;i<2;i++)
679  {
680 
681  dudt[i]+=du_dt_pnst(l,i)*psif[l];
682  mesh_veloc[i]+=dnodal_position_dt(l,i)*psif[l];
683 
684  //Loop over derivative directions for velocity gradients
685  for(unsigned j=0;j<2;j++)
686  {
687  interpolated_dudx(i,j) += u_pnst(l,i)*dpsifdx(l,j);
688  }
689 
690  }
691  }
692 
693 
694  // Get dudt in ALE form (incl mesh veloc)
695  for(unsigned i=0;i<2;i++)
696  {
697  dudt_ALE[i]=dudt[i];
698  for (unsigned k=0;k<2;k++)
699  {
700  dudt_ALE[i]-=mesh_veloc[k]*interpolated_dudx(i,k);
701  }
702  }
703 
704 
705  // Coordinates
706  for(unsigned i=0;i<2;i++)
707  {
708  outfile << interpolated_x(s,i) << " ";
709  }
710 
711  // Velocities
712  for(unsigned i=0;i<2;i++)
713  {
714  outfile << interpolated_u_pnst(s,i) << " ";
715  }
716 
717  // Pressure
718  outfile << interpolated_p_pnst(s) << " ";
719 
720  // Accelerations
721  for(unsigned i=0;i<2;i++)
722  {
723  outfile << dudt_ALE[i] << " ";
724  }
725 
726  // Dissipation
727  outfile << dissipation(s) << " ";
728 
729 
730  outfile << std::endl;
731 
732  }
733 
734 }
735 
736 //==============================================================
737 /// Return integral of dissipation over element
738 //==============================================================
740 {
741 
742  // Initialise
743  double diss=0.0;
744 
745  //Set the value of n_intpt
746  unsigned n_intpt = integral_pt()->nweight();
747 
748  //Set the Vector to hold local coordinates
749  Vector<double> s(2);
750 
751  //Loop over the integration points
752  for(unsigned ipt=0;ipt<n_intpt;ipt++)
753  {
754 
755  //Assign values of s
756  for(unsigned i=0;i<2;i++)
757  {
758  s[i] = integral_pt()->knot(ipt,i);
759  }
760 
761  //Get the integral weight
762  double w = integral_pt()->weight(ipt);
763 
764  // Get Jacobian of mapping
765  double J = J_eulerian(s);
766 
767  // Get strain rate matrix
768  DenseMatrix<double> strainrate(2,2);
769  //DenseDoubleMatrix strainrate(2,2);
770  strain_rate(s,strainrate);
771 
772  // Initialise
773  double local_diss=0.0;
774  for(unsigned i=0;i<2;i++)
775  {
776  for(unsigned j=0;j<2;j++)
777  {
778  local_diss+=2.0*strainrate(i,j)*strainrate(i,j);
779  }
780  }
781 
782  diss+=local_diss*w*J;
783  }
784 
785  return diss;
786 
787 }
788 
789 //==============================================================
790 /// \short Compute traction (on the viscous scale) at local
791 /// coordinate s for outer unit normal N
792 //==============================================================
794  const Vector<double>& N,
795  Vector<double>& traction)
796 {
797 
798  // Get velocity gradients
799  DenseMatrix<double> strainrate(2,2);
800  //DenseDoubleMatrix strainrate(2,2);
801  strain_rate(s,strainrate);
802 
803  // Get pressure
804  double press=interpolated_p_pnst(s);
805 
806  // Loop over traction components
807  for (unsigned i=0;i<2;i++)
808  {
809  traction[i]=-press*N[i];
810  for (unsigned j=0;j<2;j++)
811  {
812  traction[i]+=2.0*strainrate(i,j)*N[j];
813  }
814  }
815 }
816 
817 //==============================================================
818 /// Return dissipation at local coordinate s
819 //==============================================================
821 {
822  // Get strain rate matrix
823  DenseMatrix<double> strainrate(2,2);
824  //DenseDoubleMatrix strainrate(2,2);
825  strain_rate(s,strainrate);
826 
827  // Initialise
828  double local_diss=0.0;
829  for(unsigned i=0;i<2;i++)
830  {
831  for(unsigned j=0;j<2;j++)
832  {
833  local_diss+=2.0*strainrate(i,j)*strainrate(i,j);
834  }
835  }
836 
837  return local_diss;
838 }
839 
840 //==============================================================
841 /// Get strain-rate tensor:
842 /// Slightly more complicated in polar coordinates (see eg. Aris)
843 //==============================================================
845  DenseMatrix<double>& strainrate) const
846  {
847 
848 #ifdef PARANOID
849  if ((strainrate.ncol()!=2)||(strainrate.nrow()!=2))
850  {
851  std::ostringstream error_stream;
852  error_stream << "Wrong size " << strainrate.ncol() << " "
853  << strainrate.nrow() << std::endl;
854  throw OomphLibError(error_stream.str(),
855  OOMPH_CURRENT_FUNCTION,
856  OOMPH_EXCEPTION_LOCATION);
857  }
858 #endif
859 
860  // Velocity gradient matrix
861  DenseMatrix<double> dudx(2);
862 
863  //Find out how many nodes there are in the element
864  unsigned n_node = nnode();
865 
866  //Set up memory for the shape and test functions
867  Shape psif(n_node);
868  DShape dpsifdx(n_node,2);
869 
870  //Call the derivatives of the shape functions
871  dshape_eulerian(s,psif,dpsifdx);
872 
873  //Find the indices at which the local velocities are stored
874  unsigned u_nodal_index[2];
875  for(unsigned i=0;i<2;i++) {u_nodal_index[i] = u_index_pnst(i);}
876 
877  //Calculate local values of the velocity components
878  //Allocate
879  Vector<double> interpolated_u(2);
881  DenseMatrix<double> interpolated_dudx(2,2);
882 
883  //Initialise to zero
884  for(unsigned i=0;i<2;i++)
885  {
886  interpolated_u[i] = 0.0;
887  interpolated_x[i] = 0.0;
888  for(unsigned j=0;j<2;j++)
889  {
890  interpolated_dudx(i,j) = 0.0;
891  }
892  }
893 
894  //Calculate velocities and derivatives:
895  // Loop over nodes
896  for(unsigned l=0;l<n_node;l++)
897  {
898  //Loop over directions
899  for(unsigned i=0;i<2;i++)
900  {
901  //Get the nodal value
902  double u_value = this->nodal_value(l,u_nodal_index[i]);
903  interpolated_u[i] += u_value*psif[l];
904  interpolated_x[i] += this->nodal_position(l,i)*psif[l];
905 
906  //Loop over derivative directions
907  for(unsigned j=0;j<2;j++)
908  {
909  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
910  }
911  }
912  }
913 
914  const double Alpha = alpha();
915  //Can no longer loop over strain components
916  strainrate(0,0)=interpolated_dudx(0,0);
917  strainrate(1,1)=(1./(Alpha*interpolated_x[0]))*interpolated_dudx(1,1)+(interpolated_u[0]/interpolated_x[0]);
918  strainrate(1,0)=0.5*(interpolated_dudx(1,0)-(interpolated_u[1]/interpolated_x[0])+(1./(Alpha*interpolated_x[0]))*interpolated_dudx(0,1));
919  strainrate(0,1)=strainrate(1,0);
920 
921  }
922 
923 //==============================================================
924 /// Return polar strain-rate tensor multiplied by r
925 /// Slightly more complicated in polar coordinates (see eg. Aris)
926 //==============================================================
928  DenseMatrix<double>& strainrate) const
929  {
930 
931 #ifdef PARANOID
932  if ((strainrate.ncol()!=2)||(strainrate.nrow()!=2))
933  {
934  std::ostringstream error_stream;
935  error_stream << "Wrong size " << strainrate.ncol() << " "
936  << strainrate.nrow() << std::endl;
937  throw OomphLibError(error_stream.str(),
938  OOMPH_CURRENT_FUNCTION,
939  OOMPH_EXCEPTION_LOCATION);
940  }
941 #endif
942 
943  // Velocity gradient matrix
944  DenseMatrix<double> dudx(2);
945 
946  //Find out how many nodes there are in the element
947  unsigned n_node = nnode();
948 
949  //Set up memory for the shape and test functions
950  Shape psif(n_node);
951  DShape dpsifdx(n_node,2);
952 
953  //Call the derivatives of the shape functions
954  dshape_eulerian(s,psif,dpsifdx);
955 
956  //Find the indices at which the local velocities are stored
957  unsigned u_nodal_index[2];
958  for(unsigned i=0;i<2;i++) {u_nodal_index[i] = u_index_pnst(i);}
959 
960  //Calculate local values of the velocity components
961  //Allocate
962  Vector<double> interpolated_u(2);
964  DenseMatrix<double> interpolated_dudx(2,2);
965 
966  //Initialise to zero
967  for(unsigned i=0;i<2;i++)
968  {
969  interpolated_u[i] = 0.0;
970  interpolated_x[i] = 0.0;
971  for(unsigned j=0;j<2;j++)
972  {
973  interpolated_dudx(i,j) = 0.0;
974  }
975  }
976 
977  //Calculate velocities and derivatives:
978  // Loop over nodes
979  for(unsigned l=0;l<n_node;l++)
980  {
981  //Loop over directions
982  for(unsigned i=0;i<2;i++)
983  {
984  //Get the nodal value
985  double u_value = this->nodal_value(l,u_nodal_index[i]);
986  interpolated_u[i] += u_value*psif[l];
987  interpolated_x[i] += this->nodal_position(l,i)*psif[l];
988 
989  //Loop over derivative directions
990  for(unsigned j=0;j<2;j++)
991  {
992  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
993  }
994  }
995  }
996 
997  const double Alpha = alpha();
998  //Can no longer loop over strain components
999  strainrate(0,0)=interpolated_dudx(0,0);
1000  strainrate(1,1)=(1./(Alpha*interpolated_x[0]))*interpolated_dudx(1,1)+(interpolated_u[0]/interpolated_x[0]);
1001  strainrate(1,0)=0.5*(interpolated_dudx(1,0)-(interpolated_u[1]/interpolated_x[0])+(1./(Alpha*interpolated_x[0]))*interpolated_dudx(0,1));
1002  strainrate(0,1)=strainrate(1,0);
1003 
1004  strainrate(0,0)*=interpolated_x[0];
1005  strainrate(1,1)*=interpolated_x[0];
1006  strainrate(1,0)*=interpolated_x[0];
1007  strainrate(0,1)*=interpolated_x[0];
1008  /*
1009  strainrate(0,0)*=interpolated_x[0];
1010  strainrate(1,1)*=interpolated_x[0];
1011  strainrate(1,0)*=interpolated_x[0];
1012  strainrate(0,1)*=interpolated_x[0];
1013  */
1014  }
1015 
1016 //==============================================================
1017 /// \short Get integral of kinetic energy over element:
1018 //==============================================================
1020 {
1021 
1022  // Initialise
1023  double kin_en=0.0;
1024 
1025  //Set the value of n_intpt
1026  unsigned n_intpt = integral_pt()->nweight();
1027 
1028  //Set the Vector to hold local coordinates
1029  Vector<double> s(2);
1030 
1031  //Loop over the integration points
1032  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1033  {
1034  //Assign values of s
1035  for(unsigned i=0;i<2;i++)
1036  {
1037  s[i] = integral_pt()->knot(ipt,i);
1038  }
1039 
1040  //Get the integral weight
1041  double w = integral_pt()->weight(ipt);
1042 
1043  //Get Jacobian of mapping
1044  double J = J_eulerian(s);
1045 
1046  // Loop over directions
1047  double veloc_squared=0.0;
1048  for(unsigned i=0;i<2;i++)
1049  {
1050  veloc_squared+=interpolated_u_pnst(s,i)*interpolated_u_pnst(s,i);
1051  }
1052 
1053  kin_en+=0.5*veloc_squared*w*J;
1054  }
1055 
1056  return kin_en;
1057 
1058 }
1059 
1060 //==============================================================
1061 /// Return pressure integrated over the element
1062 //==============================================================
1064 {
1065 
1066  // Initialise
1067  double press_int=0;
1068 
1069  //Set the value of n_intpt
1070  unsigned n_intpt = integral_pt()->nweight();
1071 
1072  //Set the Vector to hold local coordinates
1073  Vector<double> s(2);
1074 
1075  //Loop over the integration points
1076  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1077  {
1078 
1079  //Assign values of s
1080  for(unsigned i=0;i<2;i++)
1081  {
1082  s[i] = integral_pt()->knot(ipt,i);
1083  }
1084 
1085  //Get the integral weight
1086  double w = integral_pt()->weight(ipt);
1087 
1088  //Get Jacobian of mapping
1089  double J = J_eulerian(s);
1090 
1091  //Premultiply the weights and the Jacobian
1092  double W = w*J;
1093 
1094  // Get pressure
1095  double press=interpolated_p_pnst(s);
1096 
1097  // Add
1098  press_int+=press*W;
1099 
1100  }
1101 
1102  return press_int;
1103 
1104 }
1105 
1106 //////////////////////////////////////////////////////////////////////
1107 //// The finished version of the new equations /////
1108 //////////////////////////////////////////////////////////////////////
1109 
1110 //==============================================================
1111 /// Compute the residuals for the Navier--Stokes
1112 /// equations; flag=1(or 0): do (or don't) compute the
1113 /// Jacobian as well.
1114 /// flag=2 for Residuals, Jacobian and mass_matrix
1115 ///
1116 /// This is now my new version with Jacobian and
1117 /// dimensionless phi
1118 //==============================================================
1121  DenseMatrix<double> &jacobian,
1122  DenseMatrix<double> &mass_matrix,
1123  unsigned flag)
1124 {
1125 
1126  //Find out how many nodes there are
1127  unsigned n_node = nnode();
1128 
1129  //Find out how many pressure dofs there are
1130  unsigned n_pres = npres_pnst();
1131 
1132  //Find the indices at which the local velocities are stored
1133  unsigned u_nodal_index[2];
1134  for(unsigned i=0;i<2;i++) {u_nodal_index[i] = u_index_pnst(i);}
1135 
1136  //Set up memory for the shape and test functions
1137  Shape psif(n_node), testf(n_node);
1138  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
1139 
1140  //Set up memory for pressure shape and test functions
1141  Shape psip(n_pres), testp(n_pres);
1142 
1143  //Number of integration points
1144  unsigned n_intpt = integral_pt()->nweight();
1145 
1146  //Set the Vector to hold local coordinates
1147  Vector<double> s(2);
1148 
1149  //Get the reynolds number and Alpha
1150  const double Re = re();
1151  const double Alpha = alpha();
1152  const double Re_St = re_st();
1153 
1154  //Integers to store the local equations and unknowns
1155  int local_eqn=0, local_unknown=0;
1156 
1157  //Loop over the integration points
1158  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1159  {
1160  //Assign values of s
1161  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
1162  //Get the integral weight
1163  double w = integral_pt()->weight(ipt);
1164 
1165  //Call the derivatives of the shape and test functions
1166  double J =
1167  dshape_and_dtest_eulerian_at_knot_pnst(ipt,psif,dpsifdx,testf,dtestfdx);
1168 
1169  //Call the pressure shape and test functions
1170  pshape_pnst(s,psip,testp);
1171 
1172  //Premultiply the weights and the Jacobian
1173  double W = w*J;
1174 
1175  //Calculate local values of the pressure and velocity components
1176  //Allocate
1177  double interpolated_p=0.0;
1178  Vector<double> interpolated_u(2);
1180  //Vector<double> mesh_velocity(2);
1181  //Vector<double> dudt(2);
1182  DenseMatrix<double> interpolated_dudx(2,2);
1183 
1184  //Initialise to zero
1185  for(unsigned i=0;i<2;i++)
1186  {
1187  //dudt[i] = 0.0;
1188  //mesh_velocity[i] = 0.0;
1189  interpolated_u[i] = 0.0;
1190  interpolated_x[i] = 0.0;
1191  for(unsigned j=0;j<2;j++)
1192  {
1193  interpolated_dudx(i,j) = 0.0;
1194  }
1195  }
1196 
1197  //Calculate pressure
1198  for(unsigned l=0;l<n_pres;l++) interpolated_p += p_pnst(l)*psip[l];
1199 
1200  //Calculate velocities and derivatives:
1201 
1202  // Loop over nodes
1203  for(unsigned l=0;l<n_node;l++)
1204  {
1205  //Loop over directions
1206  for(unsigned i=0;i<2;i++)
1207  {
1208  //Get the nodal value
1209  double u_value = this->nodal_value(l,u_nodal_index[i]);
1210  interpolated_u[i] += u_value*psif[l];
1211  interpolated_x[i] += this->nodal_position(l,i)*psif[l];
1212  //dudt[i]+=du_dt_pnst(l,i)*psif[l];
1213  //mesh_velocity[i] +=dx_dt(l,i)*psif[l];
1214 
1215  //Loop over derivative directions
1216  for(unsigned j=0;j<2;j++)
1217  {
1218  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
1219  }
1220  }
1221  }
1222 
1223  //MOMENTUM EQUATIONS
1224  //------------------
1225 
1226  //Loop over the test functions
1227  for(unsigned l=0;l<n_node;l++)
1228  {
1229  // Can't loop over velocity components as don't have identical contributions
1230  // Do seperately for i = {0,1} instead
1231  unsigned i=0;
1232  {
1233  /*IF it's not a boundary condition*/
1234  local_eqn = nodal_local_eqn(l,u_nodal_index[i]);
1235  if(local_eqn >= 0)
1236  {
1237  //Add the testf[l] term of the stress tensor
1238  residuals[local_eqn] += ((interpolated_p/interpolated_x[0])
1239  - ((1.+Gamma[i])/pow(interpolated_x[0],2.))*((1./Alpha)*interpolated_dudx(1,1)+interpolated_u[0]))
1240  *testf[l]*interpolated_x[0]*Alpha*W;
1241 
1242  //Add the dtestfdx(l,0) term of the stress tensor
1243  residuals[local_eqn] += (interpolated_p - (1.+Gamma[i])*interpolated_dudx(0,0))
1244  *dtestfdx(l,0)*interpolated_x[0]*Alpha*W;
1245 
1246  //Add the dtestfdx(l,1) term of the stress tensor
1247  residuals[local_eqn] -= ((1./(interpolated_x[0]*Alpha))*interpolated_dudx(0,1)
1248  - (interpolated_u[1]/interpolated_x[0])
1249  + Gamma[i]*interpolated_dudx(1,0))
1250  *(1./(interpolated_x[0]*Alpha))*dtestfdx(l,1)*interpolated_x[0]*Alpha*W;
1251 
1252  //Convective terms
1253  residuals[local_eqn] -= Re*(interpolated_u[0]*interpolated_dudx(0,0)
1254  +(interpolated_u[1]/(interpolated_x[0]*Alpha))*interpolated_dudx(0,1)
1255  -(pow(interpolated_u[1],2.)/interpolated_x[0]))
1256  *testf[l]*interpolated_x[0]*Alpha*W;
1257 
1258 
1259  //CALCULATE THE JACOBIAN
1260  if(flag)
1261  {
1262  //Loop over the velocity shape functions again
1263  for(unsigned l2=0;l2<n_node;l2++)
1264  {
1265  // Again can't loop over velocity components due to loss of symmetry
1266  unsigned i2=0;
1267  {
1268  //If at a non-zero degree of freedom add in the entry
1269  local_unknown = nodal_local_eqn(l2,u_nodal_index[i2]);
1270  if(local_unknown >= 0)
1271  {
1272  //Add contribution to Elemental Matrix
1273  jacobian(local_eqn,local_unknown)
1274  -= (1.+Gamma[i])*(psif[l2]/pow(interpolated_x[0],2.))*testf[l]*interpolated_x[0]*Alpha*W;
1275 
1276  jacobian(local_eqn,local_unknown)
1277  -= (1.+Gamma[i])*dpsifdx(l2,0)*dtestfdx(l,0)*interpolated_x[0]*Alpha*W;
1278 
1279  jacobian(local_eqn,local_unknown)
1280  -= (1./(interpolated_x[0]*Alpha))*dpsifdx(l2,1)*(1./(interpolated_x[0]*Alpha))*dtestfdx(l,1)*interpolated_x[0]*Alpha*W;
1281 
1282  //Now add in the inertial terms
1283  jacobian(local_eqn,local_unknown)
1284  -= Re*(psif[l2]*interpolated_dudx(0,0)+interpolated_u[0]*dpsifdx(l2,0)
1285  +(interpolated_u[1]/(interpolated_x[0]*Alpha))*dpsifdx(l2,1))*testf[l]*interpolated_x[0]*Alpha*W;
1286 
1287  //extra bit for mass matrix
1288  if(flag==2)
1289  {
1290  mass_matrix(local_eqn, local_unknown) +=
1291  Re_St*psif[l2]*testf[l]*interpolated_x[0]*Alpha*W;
1292  }
1293 
1294  } //End of (Jacobian's) if not boundary condition statement
1295  } //End of i2=0 section
1296 
1297  i2=1;
1298  {
1299 
1300  //If at a non-zero degree of freedom add in the entry
1301  local_unknown = nodal_local_eqn(l2,u_nodal_index[i2]);
1302  if(local_unknown >= 0)
1303  {
1304  //Add contribution to Elemental Matrix
1305  jacobian(local_eqn,local_unknown)
1306  -= ((1.+Gamma[i])/(pow(interpolated_x[0],2.)*Alpha))*dpsifdx(l2,1)*testf[l]*interpolated_x[0]*Alpha*W;
1307 
1308  jacobian(local_eqn,local_unknown)
1309  -= (-(psif[l2]/interpolated_x[0])+Gamma[i]*dpsifdx(l2,0))
1310  *(1./(interpolated_x[0]*Alpha))*dtestfdx(l,1)*interpolated_x[0]*Alpha*W;
1311 
1312  //Now add in the inertial terms
1313  jacobian(local_eqn,local_unknown)
1314  -= Re*((psif[l2]/(interpolated_x[0]*Alpha))*interpolated_dudx(0,1)-2*interpolated_u[1]*(psif[l2]/interpolated_x[0]))
1315  *testf[l]*interpolated_x[0]*Alpha*W;
1316 
1317  } //End of (Jacobian's) if not boundary condition statement
1318  } //End of i2=1 section
1319 
1320  } //End of l2 loop
1321 
1322  /*Now loop over pressure shape functions*/
1323  /*This is the contribution from pressure gradient*/
1324  for(unsigned l2=0;l2<n_pres;l2++)
1325  {
1326  /*If we are at a non-zero degree of freedom in the entry*/
1327  local_unknown = p_local_eqn(l2);
1328  if(local_unknown >= 0)
1329  {
1330  jacobian(local_eqn,local_unknown)
1331  += (psip[l2]/interpolated_x[0])*testf[l]*interpolated_x[0]*Alpha*W;
1332 
1333  jacobian(local_eqn,local_unknown)
1334  += psip[l2]*dtestfdx(l,0)*interpolated_x[0]*Alpha*W;
1335  }
1336  }
1337  } /*End of Jacobian calculation*/
1338 
1339  } //End of if not boundary condition statement
1340  } //End of i=0 section
1341 
1342  i=1;
1343  {
1344  /*IF it's not a boundary condition*/
1345  local_eqn = nodal_local_eqn(l,u_nodal_index[i]);
1346  if(local_eqn >= 0)
1347  {
1348  //Add the testf[l] term of the stress tensor
1349  residuals[local_eqn] += ((1./(pow(interpolated_x[0],2.)*Alpha))*interpolated_dudx(0,1)
1350  -(interpolated_u[1]/pow(interpolated_x[0],2.))
1351  +Gamma[i]*(1./interpolated_x[0])*interpolated_dudx(1,0))
1352  *testf[l]*interpolated_x[0]*Alpha*W;
1353 
1354  //Add the dtestfdx(l,0) term of the stress tensor
1355  residuals[local_eqn] -= (interpolated_dudx(1,0)
1356  +Gamma[i]*((1./(interpolated_x[0]*Alpha))*interpolated_dudx(0,1)-(interpolated_u[1]/interpolated_x[0])))
1357  *dtestfdx(l,0)*interpolated_x[0]*Alpha*W;
1358 
1359  //Add the dtestfdx(l,1) term of the stress tensor
1360  residuals[local_eqn] += (interpolated_p
1361  -(1.+Gamma[i])*((1./(interpolated_x[0]*Alpha))*interpolated_dudx(1,1)+(interpolated_u[0]/interpolated_x[0])))
1362  *(1./(interpolated_x[0]*Alpha))*dtestfdx(l,1)*interpolated_x[0]*Alpha*W;
1363 
1364  //Convective terms
1365  residuals[local_eqn] -= Re*(interpolated_u[0]*interpolated_dudx(1,0)
1366  +(interpolated_u[1]/(interpolated_x[0]*Alpha))*interpolated_dudx(1,1)
1367  +((interpolated_u[0]*interpolated_u[1])/interpolated_x[0]))
1368  *testf[l]*interpolated_x[0]*Alpha*W;
1369 
1370  //CALCULATE THE JACOBIAN
1371  if(flag)
1372  {
1373  //Loop over the velocity shape functions again
1374  for(unsigned l2=0;l2<n_node;l2++)
1375  {
1376  // Again can't loop over velocity components due to loss of symmetry
1377  unsigned i2=0;
1378  {
1379  //If at a non-zero degree of freedom add in the entry
1380  local_unknown = nodal_local_eqn(l2,u_nodal_index[i2]);
1381  if(local_unknown >= 0)
1382  {
1383  //Add contribution to Elemental Matrix
1384  jacobian(local_eqn,local_unknown)
1385  += (1./(pow(interpolated_x[0],2.)*Alpha))*dpsifdx(l2,1)
1386  *testf[l]*interpolated_x[0]*Alpha*W;
1387 
1388  jacobian(local_eqn,local_unknown)
1389  -= Gamma[i]*(1./(interpolated_x[0]*Alpha))*dpsifdx(l2,1)
1390  *dtestfdx(l,0)*interpolated_x[0]*Alpha*W;
1391 
1392  jacobian(local_eqn,local_unknown)
1393  -= (1+Gamma[i])*(psif[l2]/interpolated_x[0])*(1./(interpolated_x[0]*Alpha))
1394  *dtestfdx(l,1)*interpolated_x[0]*Alpha*W;
1395 
1396  //Now add in the inertial terms
1397  jacobian(local_eqn,local_unknown)
1398  -= Re*(psif[l2]*interpolated_dudx(1,0)+(psif[l2]*interpolated_u[1]/interpolated_x[0]))
1399  *testf[l]*interpolated_x[0]*Alpha*W;
1400 
1401  } //End of (Jacobian's) if not boundary condition statement
1402  } //End of i2=0 section
1403 
1404  i2=1;
1405  {
1406  //If at a non-zero degree of freedom add in the entry
1407  local_unknown = nodal_local_eqn(l2,u_nodal_index[i2]);
1408  if(local_unknown >= 0)
1409  {
1410  //Add contribution to Elemental Matrix
1411  jacobian(local_eqn,local_unknown)
1412  += (-(psif[l2]/pow(interpolated_x[0],2.))+Gamma[i]*(1./interpolated_x[0])*dpsifdx(l2,0))
1413  *testf[l]*interpolated_x[0]*Alpha*W;
1414 
1415  jacobian(local_eqn,local_unknown)
1416  -= (dpsifdx(l2,0)-Gamma[i]*(psif[l2]/interpolated_x[0]))
1417  *dtestfdx(l,0)*interpolated_x[0]*Alpha*W;
1418 
1419  jacobian(local_eqn,local_unknown)
1420  -= (1.+Gamma[i])*(1./(interpolated_x[0]*Alpha))*dpsifdx(l2,1)
1421  *(1./(interpolated_x[0]*Alpha))*dtestfdx(l,1)*interpolated_x[0]*Alpha*W;
1422 
1423  //Now add in the inertial terms
1424  jacobian(local_eqn,local_unknown)
1425  -= Re*(interpolated_u[0]*dpsifdx(l2,0)+(psif[l2]/(interpolated_x[0]*Alpha))*interpolated_dudx(1,1)
1426  +(interpolated_u[1]/(interpolated_x[0]*Alpha))*dpsifdx(l2,1)+(interpolated_u[0]*psif[l2]/interpolated_x[0]))
1427  *testf[l]*interpolated_x[0]*Alpha*W;
1428 
1429  //extra bit for mass matrix
1430  if(flag==2)
1431  {
1432  mass_matrix(local_eqn, local_unknown)
1433  += Re_St*psif[l2]*testf[l]*interpolated_x[0]*Alpha*W;
1434  }
1435 
1436  } //End of (Jacobian's) if not boundary condition statement
1437  } //End of i2=1 section
1438 
1439  } //End of l2 loop
1440 
1441  /*Now loop over pressure shape functions*/
1442  /*This is the contribution from pressure gradient*/
1443  for(unsigned l2=0;l2<n_pres;l2++)
1444  {
1445  /*If we are at a non-zero degree of freedom in the entry*/
1446  local_unknown = p_local_eqn(l2);
1447  if(local_unknown >= 0)
1448  {
1449  jacobian(local_eqn,local_unknown)
1450  += (psip[l2]/interpolated_x[0])*(1./Alpha)*dtestfdx(l,1)*interpolated_x[0]*Alpha*W;
1451 
1452  }
1453  }
1454  } /*End of Jacobian calculation*/
1455 
1456  } //End of if not boundary condition statement
1457 
1458  } //End of i=1 section
1459 
1460  } //End of loop over shape functions
1461 
1462  //CONTINUITY EQUATION
1463  //-------------------
1464 
1465  //Loop over the shape functions
1466  for(unsigned l=0;l<n_pres;l++)
1467  {
1468  local_eqn = p_local_eqn(l);
1469  //If not a boundary conditions
1470  if(local_eqn >= 0)
1471  {
1472  residuals[local_eqn] += (interpolated_dudx(0,0)+(interpolated_u[0]/interpolated_x[0])+(1./(interpolated_x[0]*Alpha))*interpolated_dudx(1,1))
1473  *testp[l]*interpolated_x[0]*Alpha*W;
1474 
1475 
1476  /*CALCULATE THE JACOBIAN*/
1477  if(flag)
1478  {
1479  /*Loop over the velocity shape functions*/
1480  for(unsigned l2=0;l2<n_node;l2++)
1481  {
1482  unsigned i2=0;
1483  {
1484  /*If we're at a non-zero degree of freedom add it in*/
1485  local_unknown = nodal_local_eqn(l2,u_nodal_index[i2]);
1486  if(local_unknown >= 0)
1487  {
1488  jacobian(local_eqn,local_unknown)
1489  += (dpsifdx(l2,0)+(psif[l2]/interpolated_x[0]))*testp[l]*interpolated_x[0]*Alpha*W;
1490  }
1491  } //End of i2=0 section
1492 
1493  i2=1;
1494  {
1495  /*If we're at a non-zero degree of freedom add it in*/
1496  local_unknown = nodal_local_eqn(l2,u_nodal_index[i2]);
1497  if(local_unknown >= 0)
1498  {
1499  jacobian(local_eqn,local_unknown)
1500  += (1./(interpolated_x[0]*Alpha))*dpsifdx(l2,1)*testp[l]*interpolated_x[0]*Alpha*W;
1501  }
1502  } //End of i2=1 section
1503 
1504  } /*End of loop over l2*/
1505  } /*End of Jacobian calculation*/
1506 
1507  } //End of if not boundary condition
1508  } //End of loop over l
1509 
1510 
1511  } //End of loop over integration points
1512 
1513 } //End of add_generic_residual_contribution
1514 
1515 //////////////////////////////////////////////////////////////////////
1516 //////////////////////////////////////////////////////////////////////
1517 //////////////////////////////////////////////////////////////////////
1518 
1519 ///2D Crouzeix-Raviart elements
1520 //Set the data for the number of Variables at each node
1522 ={2,2,2,2,2,2,2,2,2};
1523 
1524 //========================================================================
1525 /// Number of values (pinned or dofs) required at node n.
1526 //========================================================================
1528  int &n) const
1529  {return Initial_Nvalue[n];}
1530 
1531 //=========================================================================
1532 /// Add pointers to Data and indices of the values
1533 /// that affect the potential load (traction) applied
1534 /// to the SolidElements to the set paired_load_data
1535 //=========================================================================
1537 get_load_data(std::set<std::pair<Data*,unsigned> > &paired_load_data)
1538 {
1539  //Loop over the nodes
1540  unsigned n_node = this->nnode();
1541  for(unsigned n=0;n<n_node;n++)
1542  {
1543  //Loop over the velocity components and add pointer to their data
1544  //and indices to the vectors
1545  for(unsigned i=0;i<2;i++)
1546  {
1547  paired_load_data.insert(std::make_pair(this->node_pt(n),i));
1548  }
1549  }
1550 
1551  //Loop over the internal data
1552  unsigned n_internal = this->ninternal_data();
1553  for(unsigned l=0;l<n_internal;l++)
1554  {
1555  unsigned nval=this->internal_data_pt(l)->nvalue();
1556  //Add internal data
1557  for (unsigned j=0;j<nval;j++)
1558  {
1559  paired_load_data.insert(std::make_pair(this->internal_data_pt(l),j));
1560  }
1561  }
1562 }
1563 
1564 ///////////////////////////////////////////////////////////////////////////
1565 ///////////////////////////////////////////////////////////////////////////
1566 ///////////////////////////////////////////////////////////////////////////
1567 
1568 //2D Taylor--Hood
1569 //Set the data for the number of Variables at each node
1570 const unsigned PolarTaylorHoodElement::Initial_Nvalue[9]={3,2,3,2,2,2,3,2,3};
1571 
1572 //Set the data for the pressure conversion array
1573 const unsigned PolarTaylorHoodElement::Pconv[4]={0,2,6,8};
1574 
1575 //=========================================================================
1576 /// Add pointers to Data and the indices of the values
1577 /// that affect the potential load (traction) applied
1578 /// to the SolidElements to the set paired_load_data.
1579 //=========================================================================
1581 get_load_data(std::set<std::pair<Data*,unsigned> > &paired_load_data)
1582 {
1583  //Loop over the nodes
1584  unsigned n_node = this->nnode();
1585  for(unsigned n=0;n<n_node;n++)
1586  {
1587  //Loop over the velocity components and add pointer to their data
1588  //and indices to the vectors
1589  for(unsigned i=0;i<2;i++)
1590  {
1591  paired_load_data.insert(std::make_pair(this->node_pt(n),i));
1592  }
1593  }
1594 
1595  //Loop over the pressure data
1596  unsigned n_pres= npres_pnst();
1597  for(unsigned l=0;l<n_pres;l++)
1598  {
1599  //The 2th entry in each nodal data is the pressure, which
1600  //affects the traction
1601  paired_load_data.insert(std::make_pair(this->node_pt(Pconv[l]),2));
1602  }
1603 }
1604 
1605 //========================================================================
1606 /// Compute traction at local coordinate s for outer unit normal N
1607 //========================================================================
1608 //template<unsigned 2>
1609 //void PolarTaylorHoodElement::get_traction(const Vector<double>& s,
1610 // const Vector<double>& N,
1611 // Vector<double>& traction)
1612 //{
1613 // PolarNavierStokesEquations::traction(s,N,traction);
1614 //}
1615 
1616 }
void get_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs of pointers to data objects and unsigned integers that index th...
void full_output(std::ostream &outfile)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format...
const double & re() const
Reynolds number.
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.
virtual unsigned npres_pnst() const =0
Function to return number of pressure degrees of freedom.
unsigned npres_pnst() const
Return number of pressure values.
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
cstr elem_len * i
Definition: cfortran.h:607
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
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
void strain_rate_by_r(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Function to return polar strain multiplied by r.
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
double interpolated_p_pnst(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
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...
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition: elements.h:2226
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
char t
Definition: cfortran.h:572
double kin_energy() const
Get integral of kinetic energy over element.
virtual void pshape_pnst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
virtual double p_pnst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes...
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
double du_dt_pnst(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes...
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
void get_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs of pointers to data objects and unsigned integers that index th...
void interpolated_u_pnst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:504
virtual void fill_in_generic_residual_contribution(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)
Compute traction (on the viscous scale) at local coordinate s for outer unit normal N...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: Now returns polar strain.
static char t char * s
Definition: cfortran.h:572
double dissipation() const
Return integral of dissipation over element.
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:507
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
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
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
virtual double dshape_and_dtest_eulerian_at_knot_pnst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and derivatives w.r.t. global coords at ipt-th integration point Return J...
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:828
void output(std::ostream &outfile)
Output functionget_vels(const Vector<double>& x_to_get, Vector<double>& vels): x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
virtual int p_local_eqn(const unsigned &n)=0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
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.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual unsigned u_index_pnst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.
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
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
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.
virtual double u_pnst(const unsigned &n, const unsigned &i) const =0
Velocity i at local node n. Uses suitably interpolated value for hanging nodes.
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
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...
double pressure_integral() const
Integral of pressure over element.