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