linearised_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 linearised axisymmetric Navier-Stokes elements
27
28// oomph-lib includes
30
31namespace oomph
32{
33 //=======================================================================
34 /// Linearised axisymmetric Navier--Stokes equations static data
35 //=======================================================================
36
37 // Use the stress-divergence form by default (Gamma=1)
38 Vector<double> LinearisedNavierStokesEquations::Gamma(2, 1.0);
39
40 // "Magic" number to indicate pressure is not stored at node
42
43 // Physical constants default to zero
45
46 // Density/viscosity ratios default to one
48
49
50 //=======================================================================
51 /// Output function in tecplot format: Velocities only
52 /// r, z, U^C, U^S, V^C, V^S, W^C, W^S
53 /// at specified previous timestep (t=0 present; t>0 previous timestep).
54 /// Specified number of plot points in each coordinate direction.
55 //=======================================================================
57 const unsigned& nplot,
58 const unsigned& t)
59 {
60 // Determine number of nodes in element
61 const unsigned n_node = nnode();
62
63 // Provide storage for local shape functions
64 Shape psi(n_node);
65
66 // Provide storage for vectors of local coordinates and
67 // global coordinates and velocities
68 Vector<double> s(DIM);
70 Vector<double> interpolated_u(2 * DIM);
71
72 // Tecplot header info
73 outfile << tecplot_zone_string(nplot);
74
75 // Determine number of plot points
76 const unsigned n_plot_points = nplot_points(nplot);
77
78 // Loop over plot points
79 for (unsigned iplot = 0; iplot < n_plot_points; iplot++)
80 {
81 // Get local coordinates of plot point
82 get_s_plot(iplot, nplot, s);
83
84 // Get shape functions
85 shape(s, psi);
86
87 // Loop over coordinate directions
88 for (unsigned i = 0; i < DIM; i++)
89 {
90 // Initialise global coordinate
91 interpolated_x[i] = 0.0;
92
93 // Loop over the local nodes and sum
94 for (unsigned l = 0; l < n_node; l++)
95 {
96 interpolated_x[i] += nodal_position(t, l, i) * psi[l];
97 }
98 }
99
100 // Loop over the velocity components
101 for (unsigned i = 0; i < (2 * DIM); i++)
102 {
103 // Get the index at which the velocity is stored
104 const unsigned u_nodal_index = u_index_linearised_nst(i);
105
106 // Initialise velocity
107 interpolated_u[i] = 0.0;
108
109 // Loop over the local nodes and sum
110 for (unsigned l = 0; l < n_node; l++)
111 {
112 interpolated_u[i] += nodal_value(t, l, u_nodal_index) * psi[l];
113 }
114 }
115
116 // Output global coordinates to file
117 for (unsigned i = 0; i < DIM; i++)
118 {
119 outfile << interpolated_x[i] << " ";
120 }
121
122 // Output velocities to file
123 for (unsigned i = 0; i < (2 * DIM); i++)
124 {
125 outfile << interpolated_u[i] << " ";
126 }
127
128 outfile << std::endl;
129 }
130
131 // Write tecplot footer (e.g. FE connectivity lists)
132 write_tecplot_zone_footer(outfile, nplot);
133
134 } // End of output_veloc
135
136
137 //=======================================================================
138 /// Output function in tecplot format:
139 /// r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
140 /// in tecplot format. Specified number of plot points in each
141 /// coordinate direction.
142 //=======================================================================
143 void LinearisedNavierStokesEquations::output(std::ostream& outfile,
144 const unsigned& nplot)
145 {
146 // Provide storage for vector of local coordinates
147 Vector<double> s(DIM);
148
149 // Tecplot header info
150 outfile << tecplot_zone_string(nplot);
151
152 // Determine number of plot points
153 const unsigned n_plot_points = nplot_points(nplot);
154
155 // Loop over plot points
156 for (unsigned iplot = 0; iplot < n_plot_points; iplot++)
157 {
158 // Get local coordinates of plot point
159 get_s_plot(iplot, nplot, s);
160
161 // Output global coordinates to file
162 for (unsigned i = 0; i < DIM; i++)
163 {
164 outfile << interpolated_x(s, i) << " ";
165 }
166
167 // Output velocities (and normalisation) to file
168 for (unsigned i = 0; i < (4 * DIM); i++)
169 {
170 outfile << interpolated_u_linearised_nst(s, i) << " ";
171 }
172
173 // Output pressure to file
174 for (unsigned i = 0; i < 2; i++)
175 {
176 outfile << interpolated_p_linearised_nst(s, i) << " ";
177 }
178
179 outfile << std::endl;
180 }
181 outfile << std::endl;
182
183 // Write tecplot footer (e.g. FE connectivity lists)
184 write_tecplot_zone_footer(outfile, nplot);
185
186 } // End of output
187
188
189 //=======================================================================
190 /// Output function in tecplot format:
191 /// r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
192 /// Specified number of plot points in each coordinate direction.
193 //=======================================================================
195 const unsigned& nplot)
196 {
197 // Provide storage for vector of local coordinates
198 Vector<double> s(2);
199
200 // Tecplot header info
201 fprintf(file_pt, "%s ", tecplot_zone_string(nplot).c_str());
202
203 // Determine number of plot points
204 const unsigned n_plot_points = nplot_points(nplot);
205
206 // Loop over plot points
207 for (unsigned iplot = 0; iplot < n_plot_points; iplot++)
208 {
209 // Get local coordinates of plot point
210 get_s_plot(iplot, nplot, s);
211
212 // Output global coordinates to file
213 for (unsigned i = 0; i < 2; i++)
214 {
215 fprintf(file_pt, "%g ", interpolated_x(s, i));
216 }
217
218 // Output velocities to file
219 for (unsigned i = 0; i < (2 * DIM); i++)
220 {
221 fprintf(file_pt, "%g ", interpolated_u_linearised_nst(s, i));
222 }
223
224 // Output pressure to file
225 for (unsigned i = 0; i < 2; i++)
226 {
227 fprintf(file_pt, "%g ", interpolated_p_linearised_nst(s, i));
228 }
229
230 fprintf(file_pt, "\n");
231 }
232
233 fprintf(file_pt, "\n");
234
235 // Write tecplot footer (e.g. FE connectivity lists)
236 write_tecplot_zone_footer(file_pt, nplot);
237
238 } // End of output
239
240
241 //=======================================================================
242 /// Get strain-rate tensor: \f$ e_{ij} \f$ where
243 /// \f$ i,j = r,z,\theta \f$ (in that order). We evaluate this tensor
244 /// at a value of theta such that the product of theta and the azimuthal
245 /// mode number (k) gives \f$ \pi/4 \f$. Therefore
246 /// \f$ \cos(k \theta) = \sin(k \theta) = 1/\sqrt{2} \f$.
247 //=======================================================================
249 const Vector<double>& s,
250 DenseMatrix<double>& strainrate,
251 const unsigned& real) const
252 {
253#ifdef PARANOID
254 if ((strainrate.ncol() != DIM) || (strainrate.nrow() != DIM))
255 {
256 std::ostringstream error_message;
257 error_message << "The strain rate has incorrect dimensions "
258 << strainrate.ncol() << " x " << strainrate.nrow()
259 << " Not " << DIM << std::endl;
260
261 throw OomphLibError(
262 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
263 }
264#endif
265
266 // The question is what to do about the real and imaginary parts
267 // The answer is that we know that the "real" velocity is
268 // U_r cos(omega t) - U_i sin(omega t) and we can choose omega t
269 // to be 7pi/4 so that cos = 1/sqrt(2) and sin = -1/sqrt(2)
270 // to get equal contributions
271 double cosomegat = 1.0 / sqrt(2.0);
272 double sinomegat = cosomegat;
273
274 // Find out how many nodes there are in the element
275 unsigned n_node = nnode();
276
277 // Set up memory for the shape and test functions
278 Shape psif(n_node);
279 DShape dpsifdx(n_node, DIM);
280
281 // Call the derivatives of the shape functions
282 dshape_eulerian(s, psif, dpsifdx);
283
284 // Velocity gradient matrix
286 // Initialise to zero
287 for (unsigned i = 0; i < DIM; ++i)
288 {
289 dudx[i].resize(DIM);
290 for (unsigned j = 0; j < DIM; ++j)
291 {
292 dudx[i][j].real(0.0);
293 dudx[i][j].imag(0.0);
294 }
295 }
296
297 // Get the nodal indices at which the velocity is stored
298 unsigned n_veloc = 2 * DIM;
299 unsigned u_nodal_index[n_veloc];
300 for (unsigned i = 0; i < n_veloc; ++i)
301 {
302 u_nodal_index[i] = u_index_linearised_nst(i);
303 }
304
305 // Loop over the element's nodes
306 for (unsigned l = 0; l < n_node; l++)
307 {
308 // Loop over the DIM complex velocity components
309 for (unsigned i = 0; i < DIM; i++)
310 {
311 // Get the value
312 const std::complex<double> u_value(
313 this->raw_nodal_value(l, u_nodal_index[2 * i + 0]),
314 this->raw_nodal_value(l, u_nodal_index[2 * i + 1]));
315
316 // Loop over two coordinate directions (for derivatives)
317 for (unsigned j = 0; j < DIM; j++)
318 {
319 dudx[i][j] += u_value * dpsifdx(l, j);
320 }
321 }
322 }
323
324 // Take the real part of the velocity gradient matrix
325 DenseMatrix<double> real_dudx(DIM);
326 if (real == 0)
327 {
328 cosomegat = 1.0;
329 sinomegat = 0.0;
330 }
331 else
332 {
333 cosomegat = 0.0;
334 sinomegat = 1.0;
335 }
336
337 for (unsigned i = 0; i < DIM; ++i)
338 {
339 for (unsigned j = 0; j < DIM; ++j)
340 {
341 real_dudx(i, j) =
342 dudx[i][j].real() * cosomegat + dudx[i][j].imag() * sinomegat;
343 }
344 }
345
346
347 // Now set the strain rate
348 // Loop over veclocity components
349 for (unsigned i = 0; i < DIM; i++)
350 {
351 // Loop over derivative directions
352 for (unsigned j = 0; j < DIM; j++)
353 {
354 strainrate(i, j) = 0.5 * (real_dudx(i, j) + real_dudx(j, i));
355 }
356 }
357
358 } // End of strain_rate
359
360
361 //=======================================================================
362 /// Compute the residuals for the linearised axisymmetric Navier--Stokes
363 /// equations; flag=1(or 0): do (or don't) compute the Jacobian as well.
364 //=======================================================================
367 Vector<double>& residuals,
368 DenseMatrix<double>& jacobian,
369 DenseMatrix<double>& mass_matrix,
370 unsigned flag)
371 {
372 // Get the time from the first node in the element
373 const double time = this->node_pt(0)->time_stepper_pt()->time();
374
375 // Determine number of nodes in the element
376 const unsigned n_node = nnode();
377
378 // Determine how many pressure values there are associated with
379 // a single pressure component
380 const unsigned n_pres = npres_linearised_nst();
381
382 const unsigned n_veloc = 4 * DIM;
383
384 // Get the nodal indices at which the velocity is stored
385 unsigned u_nodal_index[n_veloc];
386 for (unsigned i = 0; i < n_veloc; ++i)
387 {
388 u_nodal_index[i] = u_index_linearised_nst(i);
389 }
390
391 // Set up memory for the fluid shape and test functions
392 Shape psif(n_node), testf(n_node);
393 DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
394
395 // Set up memory for the pressure shape and test functions
396 Shape psip(n_pres), testp(n_pres);
397
398 // Determine number of integration points
399 const unsigned n_intpt = integral_pt()->nweight();
400
401 // Set up memory for the vector to hold local coordinates (two dimensions)
402 Vector<double> s(DIM);
403
404 // Get physical variables from the element
405 // (Reynolds number must be multiplied by the density ratio)
406 const double scaled_re = re() * density_ratio();
407 const double scaled_re_st = re_st() * density_ratio();
408 const double visc_ratio = viscosity_ratio();
409
410 const double eval_real = lambda();
411 const double eval_imag = omega();
412
413 const std::complex<double> eigenvalue(eval_real, eval_imag);
414
415 // Integers used to store the local equation and unknown numbers
416 int local_eqn = 0;
417
418 // Loop over the integration points
419 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
420 {
421 // Assign values of the local coordinates s
422 for (unsigned i = 0; i < DIM; i++)
423 {
424 s[i] = integral_pt()->knot(ipt, i);
425 }
426
427 // Get the integral weight
428 const double w = integral_pt()->weight(ipt);
429
430 // Calculate the fluid shape and test functions, and their derivatives
431 // w.r.t. the global coordinates
433 ipt, psif, dpsifdx, testf, dtestfdx);
434
435 // Calculate the pressure shape and test functions
436 pshape_linearised_nst(s, psip, testp);
437
438 // Premultiply the weights and the Jacobian of the mapping between
439 // local and global coordinates
440 const double W = w * J;
441
442 // Allocate storage for the position and the derivative of the
443 // mesh positions w.r.t. time
445 // Vector<double> mesh_velocity(2,0.0);
446
447 // Allocate storage for the velocity components (six of these)
448 // and their derivatives w.r.t. time
449 Vector<std::complex<double>> interpolated_u(DIM);
450 // Vector<double> dudt(6,0.0);
451 // Allocate storage for the eigen function normalisation
452 Vector<std::complex<double>> interpolated_u_normalisation(DIM);
453 for (unsigned i = 0; i < DIM; ++i)
454 {
455 interpolated_u[i].real(0.0);
456 interpolated_u[i].imag(0.0);
457 interpolated_u_normalisation[i].real(0.0);
458 interpolated_u_normalisation[i].imag(0.0);
459 }
460
461 // Allocate storage for the pressure components (two of these
462 std::complex<double> interpolated_p(0.0, 0.0);
463 std::complex<double> interpolated_p_normalisation(0.0, 0.0);
464
465 // Allocate storage for the derivatives of the velocity components
466 // w.r.t. global coordinates (r and z)
467 Vector<Vector<std::complex<double>>> interpolated_dudx(DIM);
468 for (unsigned i = 0; i < DIM; ++i)
469 {
470 interpolated_dudx[i].resize(DIM);
471 for (unsigned j = 0; j < DIM; ++j)
472 {
473 interpolated_dudx[i][j].real(0.0);
474 interpolated_dudx[i][j].imag(0.0);
475 }
476 }
477
478 // Calculate pressure at the integration point
479 // -------------------------------------------
480
481 // Loop over pressure degrees of freedom (associated with a single
482 // pressure component) in the element
483 for (unsigned l = 0; l < n_pres; l++)
484 {
485 // Cache the shape function
486 const double psip_ = psip(l);
487
488 // Get the complex nodal pressure values
489 const std::complex<double> p_value(this->p_linearised_nst(l, 0),
490 this->p_linearised_nst(l, 1));
491
492 // Add contribution
493 interpolated_p += p_value * psip_;
494
495 // Repeat for normalisation
496 const std::complex<double> p_norm_value(this->p_linearised_nst(l, 2),
497 this->p_linearised_nst(l, 3));
498 interpolated_p_normalisation += p_norm_value * psip_;
499 }
500 // End of loop over the pressure degrees of freedom in the element
501
502 // Calculate velocities and their derivatives at the integration point
503 // -------------------------------------------------------------------
504
505 // Loop over the element's nodes
506 for (unsigned l = 0; l < n_node; l++)
507 {
508 // Cache the shape function
509 const double psif_ = psif(l);
510
511 // Loop over the DIM coordinate directions
512 for (unsigned i = 0; i < DIM; i++)
513 {
514 interpolated_x[i] += this->raw_nodal_position(l, i) * psif_;
515 }
516
517 // Loop over the DIM complex velocity components
518 for (unsigned i = 0; i < DIM; i++)
519 {
520 // Get the value
521 const std::complex<double> u_value(
522 this->raw_nodal_value(l, u_nodal_index[2 * i + 0]),
523 this->raw_nodal_value(l, u_nodal_index[2 * i + 1]));
524
525 // Add contribution
526 interpolated_u[i] += u_value * psif_;
527
528 // Add contribution to dudt
529 // dudt[i] += du_dt_linearised_nst(l,i)*psif_;
530
531 // Loop over two coordinate directions (for derivatives)
532 for (unsigned j = 0; j < DIM; j++)
533 {
534 interpolated_dudx[i][j] += u_value * dpsifdx(l, j);
535 }
536
537 // Interpolate the normalisation function
538 const std::complex<double> normalisation_value(
539 this->raw_nodal_value(l, u_nodal_index[2 * (DIM + i)]),
540 this->raw_nodal_value(l, u_nodal_index[2 * (DIM + i) + 1]));
541 interpolated_u_normalisation[i] += normalisation_value * psif_;
542 }
543 } // End of loop over the element's nodes
544
545 // Get the mesh velocity if ALE is enabled
546 /*if(!ALE_is_disabled)
547 {
548 // Loop over the element's nodes
549 for(unsigned l=0;l<n_node;l++)
550 {
551 // Loop over the two coordinate directions
552 for(unsigned i=0;i<2;i++)
553 {
554 mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif(l);
555 }
556 }
557 }*/
558
559 // Get velocities and their derivatives from base flow problem
560 // -----------------------------------------------------------
561
562 // Allocate storage for the velocity components of the base state
563 // solution (initialise to zero)
564 Vector<double> base_flow_u(DIM, 0.0);
565
566 // Get the user-defined base state solution velocity components
567 get_base_flow_u(time, ipt, interpolated_x, base_flow_u);
568
569 // Allocate storage for the derivatives of the base state solution's
570 // velocity components w.r.t. global coordinate (r and z)
571 // N.B. the derivatives of the base flow components w.r.t. the
572 // azimuthal coordinate direction (theta) are always zero since the
573 // base flow is axisymmetric
574 DenseMatrix<double> base_flow_dudx(DIM, DIM, 0.0);
575
576 // Get the derivatives of the user-defined base state solution
577 // velocity components w.r.t. global coordinates
578 get_base_flow_dudx(time, ipt, interpolated_x, base_flow_dudx);
579
580
581 // MOMENTUM EQUATIONS
582 //------------------
583
584 // Loop over the test functions
585 for (unsigned l = 0; l < n_node; l++)
586 {
587 // Loop over the velocity components
588 for (unsigned i = 0; i < DIM; i++)
589 {
590 // Assemble the residuals
591 // Time dependent term
592 std::complex<double> residual_contribution =
593 -scaled_re_st * eigenvalue * interpolated_u[i] * testf[l] * W;
594 // Pressure term
595 residual_contribution += interpolated_p * dtestfdx(l, i) * W;
596 // Viscous terms
597 for (unsigned k = 0; k < DIM; ++k)
598 {
599 residual_contribution -=
600 visc_ratio *
601 (interpolated_dudx[i][k] + Gamma[i] * interpolated_dudx[k][i]) *
602 dtestfdx(l, k) * W;
603 }
604
605 // Advective terms
606 for (unsigned k = 0; k < DIM; ++k)
607 {
608 residual_contribution -=
609 scaled_re *
610 (base_flow_u[k] * interpolated_dudx[i][k] +
611 interpolated_u[k] * base_flow_dudx(i, k)) *
612 testf[l] * W;
613 }
614
615 // Now separate real and imaginary parts
616
617 /*IF it's not a boundary condition*/
618 // Here assume that we're only going to pin entire complex
619 // number or not
620 local_eqn = nodal_local_eqn(l, u_nodal_index[2 * i]);
621 if (local_eqn >= 0)
622 {
623 residuals[local_eqn] += residual_contribution.real();
624 }
625
626 local_eqn = nodal_local_eqn(l, u_nodal_index[2 * i + 1]);
627 if (local_eqn >= 0)
628 {
629 residuals[local_eqn] += residual_contribution.imag();
630 }
631
632
633 // CALCULATE THE JACOBIAN
634 /*if(flag)
635 {
636 //Loop over the velocity shape functions again
637 for(unsigned l2=0;l2<n_node;l2++)
638 {
639 //Loop over the velocity components again
640 for(unsigned i2=0;i2<DIM;i2++)
641 {
642 //If at a non-zero degree of freedom add in the entry
643 local_unknown = nodal_local_eqn(l2,u_nodal_index[i2]);
644 if(local_unknown >= 0)
645 {
646 //Add contribution to Elemental Matrix
647 jacobian(local_eqn,local_unknown)
648 -= visc_ratio*Gamma[i]*dpsifdx(l2,i)*dtestfdx(l,i2)*W;
649
650 //Extra component if i2 = i
651 if(i2 == i)
652 {
653 //Loop over velocity components
654 for(unsigned k=0;k<DIM;k++)
655 {
656 jacobian(local_eqn,local_unknown)
657 -= visc_ratio*dpsifdx(l2,k)*dtestfdx(l,k)*W;
658 }
659 }
660
661 //Now add in the inertial terms
662 jacobian(local_eqn,local_unknown)
663 -= scaled_re*psif[l2]*interpolated_dudx(i,i2)*testf[l]*W;
664
665 //Extra component if i2=i
666 if(i2 == i)
667 {
668 //Add the mass matrix term (only diagonal entries)
669 //Note that this is positive because the mass matrix
670 //is taken to the other side of the equation when
671 //formulating the generalised eigenproblem.
672 if(flag==2)
673 {
674 mass_matrix(local_eqn,local_unknown) +=
675 scaled_re_st*psif[l2]*testf[l]*W;
676 }
677
678 //du/dt term
679 jacobian(local_eqn,local_unknown)
680 -= scaled_re_st*
681 node_pt(l2)->time_stepper_pt()->weight(1,0)*
682 psif[l2]*testf[l]*W;
683
684 //Loop over the velocity components
685 for(unsigned k=0;k<DIM;k++)
686 {
687 double tmp=scaled_re*interpolated_u[k];
688 if (!ALE_is_disabled) tmp-=scaled_re_st*mesh_velocity[k];
689 jacobian(local_eqn,local_unknown) -=
690 tmp*dpsifdx(l2,k)*testf[l]*W;
691 }
692 }
693
694 }
695 }
696 }
697
698 //Now loop over pressure shape functions
699 //This is the contribution from pressure gradient
700 for(unsigned l2=0;l2<n_pres;l2++)
701 {
702 //If we are at a non-zero degree of freedom in the entry
703 local_unknown = p_local_eqn(l2);
704 if(local_unknown >= 0)
705 {
706 jacobian(local_eqn,local_unknown)
707 += psip[l2]*dtestfdx(l,i)*W;
708 }
709 }
710 } //End of Jacobian calculation
711
712 }*/ //End of if not boundary condition statement
713
714 } // End of loop over dimension
715 } // End of loop over shape functions
716
717
718 // CONTINUITY EQUATION
719 //-------------------
720
721 // Loop over the shape functions
722 for (unsigned l = 0; l < n_pres; l++)
723 {
724 // Assemble the residuals
725 std::complex<double> residual_contribution = interpolated_dudx[0][0];
726 for (unsigned k = 1; k < DIM; ++k)
727 {
728 residual_contribution += interpolated_dudx[k][k];
729 }
730
731 local_eqn = p_local_eqn(l, 0);
732 // If not a boundary conditions
733 if (local_eqn >= 0)
734 {
735 residuals[local_eqn] += residual_contribution.real() * testp[l] * W;
736 }
737
738 local_eqn = p_local_eqn(l, 1);
739 // If not a boundary conditions
740 if (local_eqn >= 0)
741 {
742 residuals[local_eqn] += residual_contribution.imag() * testp[l] * W;
743 }
744
745 } // End of loop over l
746
747 // Normalisation condition
748 std::complex<double> residual_contribution =
749 interpolated_p_normalisation * interpolated_p;
750 for (unsigned k = 0; k < DIM; ++k)
751 {
752 residual_contribution +=
753 interpolated_u_normalisation[k] * interpolated_u[k];
754 }
755
756 local_eqn = this->eigenvalue_local_eqn(0);
757 if (local_eqn >= 0)
758 {
759 residuals[local_eqn] += residual_contribution.real() * W;
760 }
761
762 local_eqn = this->eigenvalue_local_eqn(1);
763 if (local_eqn >= 0)
764 {
765 residuals[local_eqn] += residual_contribution.imag() * W;
766 }
767
768 } // End of loop over the integration points
769
770 } // End of fill_in_generic_residual_contribution_linearised_nst
771
772
773 /// ////////////////////////////////////////////////////////////////////////
774 /// ////////////////////////////////////////////////////////////////////////
775 /// ////////////////////////////////////////////////////////////////////////
776
777
778 /// Linearised axisymmetric Crouzeix-Raviart elements
779 /// -------------------------------------------------
780
781
782 //=======================================================================
783 /// Set the data for the number of variables at each node
784 //=======================================================================
786 8, 8, 8, 8, 8, 8, 8, 8, 8};
787
788
789 //========================================================================
790 /// Number of values (pinned or dofs) required at node n
791 //========================================================================
793 const unsigned& n) const
794 {
795 return Initial_Nvalue[n];
796 }
797
798
799 /// ////////////////////////////////////////////////////////////////////////
800 /// ////////////////////////////////////////////////////////////////////////
801 /// ////////////////////////////////////////////////////////////////////////
802
803
804 /// Linearised axisymmetric Taylor-Hood elements
805 /// --------------------------------------------
806
807
808 //=======================================================================
809 /// Set the data for the number of variables at each node
810 //=======================================================================
812 12, 8, 12, 8, 8, 8, 12, 8, 12};
813
814
815 //=======================================================================
816 /// Set the data for the pressure conversion array
817 //=======================================================================
818 const unsigned LinearisedQTaylorHoodElement::Pconv[4] = {0, 2, 6, 8};
819
820
821} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
unsigned 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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3298
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2576
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation....
Definition: elements.cc:1686
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
const double & viscosity_ratio() const
Viscosity ratio for element: element's viscosity relative to the viscosity used in the definition of ...
virtual unsigned u_index_linearised_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
virtual unsigned npres_linearised_nst() const =0
Return the number of pressure degrees of freedom associated with a single pressure component in the e...
const double & density_ratio() const
Density ratio for element: element's density relative to the viscosity used in the definition of the ...
virtual int p_local_eqn(const unsigned &n, const unsigned &i)=0
Access function for the local equation number information for the i-th component of the pressure....
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate, const unsigned &real) const
Strain-rate tensor: where (in that order)
void output_veloc(std::ostream &outfile, const unsigned &nplot, const unsigned &t)
Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, in tecplot format. nplot points in each coordina...
virtual void pshape_linearised_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
virtual void get_base_flow_u(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result) const
Calculate the velocity components of the base flow solution at a given time and Eulerian position.
virtual double p_linearised_nst(const unsigned &n_p, const unsigned &i) const =0
Return the i-th pressure value at local pressure "node" n_p. Uses suitably interpolated value for han...
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
virtual double dshape_and_dtest_eulerian_at_knot_linearised_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and their derivatives w.r.t. global coordinates at the ipt-th integration...
void output(std::ostream &outfile)
Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S in tecplot format. Default number of pl...
virtual void fill_in_generic_residual_contribution_linearised_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 Jacobi...
double interpolated_p_linearised_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated pressure p[i] at local coordinate s.
virtual void get_base_flow_dudx(const double &time, const unsigned &ipt, const Vector< double > &x, DenseMatrix< double > &result) const
Calculate the derivatives of the velocity components of the base flow solution w.r....
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
double interpolated_u_linearised_nst(const Vector< double > &s, const unsigned &i) const
Compute the element's residual Vector and the jacobian matrix. Virtual function can be overloaded by ...
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all initialised to one)
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
virtual unsigned required_nvalue(const unsigned &n) const
Return number of values (pinned or dofs) required at local node n.
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
An OomphLibError object which should be thrown when an run-time error is encountered....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
double & time()
Return current value of continous time.
Definition: timesteppers.h:332
//////////////////////////////////////////////////////////////////// ////////////////////////////////...