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