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