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