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