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