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