refineable_spacetime_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 space-time Navier-Stokes elements
28
29/// //////////////////////////////////////////////////////////////////////
30/// //////////////////////////////////////////////////////////////////////
31/// //////////////////////////////////////////////////////////////////////
32
33namespace oomph
34{
35 //===================================================================
36 /// Compute the diagonal of the velocity/pressure mass matrices.
37 /// If which one=0, both are computed, otherwise only the pressure
38 /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
39 /// LSC version of the preconditioner only needs that one)
40 //===================================================================
41 template<unsigned DIM>
44 Vector<double>& press_mass_diag,
45 Vector<double>& veloc_mass_diag,
46 const unsigned& which_one)
47 {
48 // Resize and initialise
49 unsigned n_dof = ndof();
50
51 // If the pressure mass matrix needs to be constructed
52 if ((which_one == 0) || (which_one == 1))
53 {
54 // Assign the appropriate amount of space
55 press_mass_diag.assign(n_dof, 0.0);
56 }
57
58 // If the velocity mass matrix needs to be constructed
59 if ((which_one == 0) || (which_one == 2))
60 {
61 // Assign the appropriate amount of space
62 veloc_mass_diag.assign(n_dof, 0.0);
63 }
64
65 // Pointer to hang info object
66 HangInfo* hang_info_pt = 0;
67
68 // Number of master nodes
69 unsigned n_master = 1;
70
71 // Hang weight for shape functions
72 double hang_weight = 1.0;
73
74 // Find out how many nodes there are in the element
75 unsigned n_node = nnode();
76
77 // Local coordinates
78 Vector<double> s(DIM + 1, 0.0);
79
80 // Storage for the local velocity indices
81 Vector<unsigned> u_nodal_index(DIM, 0.0);
82
83 // Find the indices at which the local velocities are stored
84 for (unsigned i = 0; i < DIM; i++)
85 {
86 // Calculate the i-th local velocity component
87 u_nodal_index[i] = this->u_index_nst(i);
88 }
89
90 // Set up memory for velocity shape functions
91 Shape psi(n_node);
92
93 // Find number of pressure dofs
94 unsigned n_pres = this->npres_nst();
95
96 // Pressure shape function
97 Shape psi_p(n_pres);
98
99 // Which nodal value represents the pressure? (Negative if pressure
100 // is not based on nodal interpolation).
101 int p_index = this->p_nodal_index_nst();
102
103 // Local array of booleans that are true if the l-th pressure value is
104 // hanging (avoid repeated virtual function calls)
105 bool pressure_dof_is_hanging[n_pres];
106
107 // If the pressure is stored at a node
108 if (p_index >= 0)
109 {
110 // Loop over the pressure nodes
111 for (unsigned l = 0; l < n_pres; l++)
112 {
113 // Read out the hang status of the pressure node
114 pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
115 }
116 }
117 // Otherwise the pressure is not stored at a node and so cannot hang
118 else
119 {
120 // Loop over the pressures nodes
121 for (unsigned l = 0; l < n_pres; l++)
122 {
123 // Indicate that the pressure node is not hanging
124 pressure_dof_is_hanging[l] = false;
125 }
126 } // if (p_index>=0)
127
128 // Number of integration points
129 unsigned n_intpt = integral_pt()->nweight();
130
131 // Integer to store the local equations no
132 int local_eqn = 0;
133
134 // Loop over the integration points
135 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
136 {
137 // Get the integral weight
138 double w = integral_pt()->weight(ipt);
139
140 // Get determinant of Jacobian of the mapping
141 double J = J_eulerian_at_knot(ipt);
142
143 // Assign values of s
144 for (unsigned i = 0; i < DIM + 1; i++)
145 {
146 // Calculate the i-th local coordinate at the ipt-th integration point
147 s[i] = integral_pt()->knot(ipt, i);
148 }
149
150 // Premultiply weights and Jacobian
151 double W = w * J;
152
153 // Do we want the velocity one?
154 if ((which_one == 0) || (which_one == 2))
155 {
156 // Get the velocity shape functions
157 shape_at_knot(ipt, psi);
158
159 // Loop over the velocity shape functions
160 for (unsigned l = 0; l < n_node; l++)
161 {
162 // Local boolean to indicate whether or not the node is hanging
163 bool is_node_hanging = node_pt(l)->is_hanging();
164
165 // If the node is hanging
166 if (is_node_hanging)
167 {
168 // Get the HangInfo pointer from the node
169 hang_info_pt = node_pt(l)->hanging_pt();
170
171 // Read out number of master nodes from hanging data
172 n_master = hang_info_pt->nmaster();
173 }
174 // Otherwise the node is its own master
175 else
176 {
177 // There is only one node to consider; the node itself
178 n_master = 1;
179 }
180
181 // Loop over the master nodes
182 for (unsigned m = 0; m < n_master; m++)
183 {
184 // Loop over velocity components for equations
185 for (unsigned i = 0; i < DIM; i++)
186 {
187 // If the node is hanging
188 if (is_node_hanging)
189 {
190 // Get the equation number from the master node
191 local_eqn = this->local_hang_eqn(
192 hang_info_pt->master_node_pt(m), u_nodal_index[i]);
193
194 // Get the hang weight from the master node
195 hang_weight = hang_info_pt->master_weight(m);
196 }
197 // If the node is not hanging
198 else
199 {
200 // Local equation number
201 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
202
203 // Node contributes with full weight
204 hang_weight = 1.0;
205 }
206
207 // If it's not a boundary condition
208 if (local_eqn >= 0)
209 {
210 // Add the contribution
211 veloc_mass_diag[local_eqn] += pow(psi[l] * hang_weight, 2) * W;
212 }
213 } // for (unsigned i=0;i<DIM;i++)
214 } // for (unsigned m=0;m<n_master;m++)
215 } // if ((which_one==0)||(which_one==2))
216 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
217
218 // Do we want the pressure one?
219 if ((which_one == 0) || (which_one == 1))
220 {
221 // Get the pressure shape functions
222 this->pshape_nst(s, psi_p);
223
224 // Loop over the pressure shape functions
225 for (unsigned l = 0; l < n_pres; l++)
226 {
227 // If the pressure dof is hanging
228 if (pressure_dof_is_hanging[l])
229 {
230 // Get the HangInfo object
231 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
232
233 // Get the number of master nodes from the pressure node
234 n_master = hang_info_pt->nmaster();
235 }
236 // Otherwise the node is its own master
237 else
238 {
239 // There is only one node to consider; the node itself
240 n_master = 1;
241 }
242
243 // Loop over the master nodes
244 for (unsigned m = 0; m < n_master; m++)
245 {
246 // If the pressure dof is hanging
247 if (pressure_dof_is_hanging[l])
248 {
249 // Get the local equation from the master node
250 local_eqn =
251 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
252
253 // Get the weight for the node
254 hang_weight = hang_info_pt->master_weight(m);
255 }
256 // If the pressure dof is not hanging
257 else
258 {
259 // Get the local equation number of this pressure dof
260 local_eqn = this->p_local_eqn(l);
261
262 // Get the hang weight associated with this pressure node
263 hang_weight = 1.0;
264 }
265
266 // If the equation is not pinned
267 if (local_eqn >= 0)
268 {
269 // Add the contribution
270 press_mass_diag[local_eqn] += pow(psi_p[l] * hang_weight, 2) * W;
271 }
272 } // for (unsigned m=0;m<n_master;m++)
273 } // for (unsigned l=0;l<n_pres;l++)
274 } // if ((which_one==0)||(which_one==1))
275 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
276 } // End of get_pressure_and_velocity_mass_matrix_diagonal
277
278
279 //==============================================================
280 /// Compute the residuals for the associated pressure advection
281 /// diffusion problem. Used by the Fp preconditioner.
282 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
283 //==============================================================
284 template<unsigned DIM>
287 Vector<double>& residuals,
288 DenseMatrix<double>& jacobian,
289 const unsigned& flag)
290 {
291 OomphLibWarning("I'm not sure this is correct yet...",
292 OOMPH_CURRENT_FUNCTION,
293 OOMPH_EXCEPTION_LOCATION);
294
295 // Return immediately if there are no dofs
296 if (ndof() == 0) return;
297
298 // Find out how many nodes there are
299 unsigned n_node = nnode();
300
301 // Find out how many pressure dofs there are
302 unsigned n_pres = this->npres_nst();
303
304 // Find the indices at which the local velocities are stored
305 unsigned u_nodal_index[DIM];
306 for (unsigned i = 0; i < DIM; i++)
307 {
308 u_nodal_index[i] = this->u_index_nst(i);
309 }
310
311
312 // Which nodal value represents the pressure? (Negative if pressure
313 // is not based on nodal interpolation).
314 int p_index = this->p_nodal_index_nst();
315
316 // Local array of booleans that are true if the l-th pressure value is
317 // hanging (avoid repeated virtual function calls)
318 bool pressure_dof_is_hanging[n_pres];
319 // If the pressure is stored at a node
320 if (p_index >= 0)
321 {
322 // Read out whether the pressure is hanging
323 for (unsigned l = 0; l < n_pres; ++l)
324 {
325 pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
326 }
327 }
328 // Otherwise the pressure is not stored at a node and so cannot hang
329 else
330 {
331 // pressure advection diffusion doesn't work for this one!
332 throw OomphLibError(
333 "Pressure advection diffusion does not work in this case\n",
334 OOMPH_CURRENT_FUNCTION,
335 OOMPH_EXCEPTION_LOCATION);
336
337 for (unsigned l = 0; l < n_pres; ++l)
338 {
339 pressure_dof_is_hanging[l] = false;
340 }
341 }
342
343 // Set up memory for the velocity shape fcts
344 Shape psif(n_node);
345 DShape dpsidx(n_node, DIM);
346
347 // Set up memory for pressure shape and test functions
348 Shape psip(n_pres), testp(n_pres);
349 DShape dpsip(n_pres, DIM);
350 DShape dtestp(n_pres, DIM);
351
352 // Number of integration points
353 unsigned n_intpt = integral_pt()->nweight();
354
355 // Set the Vector to hold local coordinates
356 Vector<double> s(DIM);
357
358 // Get Physical Variables from Element
359 // Reynolds number must be multiplied by the density ratio
360 double scaled_re = this->re() * this->density_ratio();
361
362 // Integers to store the local equations and unknowns
363 int local_eqn = 0, local_unknown = 0;
364
365 // Pointers to hang info objects
366 HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
367
368 // Loop over the integration points
369 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
370 {
371 // Assign values of s
372 for (unsigned i = 0; i < DIM; i++) s[i] = integral_pt()->knot(ipt, i);
373
374 // Get the integral weight
375 double w = integral_pt()->weight(ipt);
376
377 // Call the derivatives of the veloc shape functions
378 // (Derivs not needed but they are free)
379 double J = this->dshape_eulerian_at_knot(ipt, psif, dpsidx);
380
381 // Call the pressure shape and test functions
382 this->dpshape_and_dptest_eulerian_nst(s, psip, dpsip, testp, dtestp);
383
384 // Premultiply the weights and the Jacobian
385 double W = w * J;
386
387 // Calculate local values of the pressure and velocity components
388 // Allocate
389 Vector<double> interpolated_u(DIM, 0.0);
390 Vector<double> interpolated_x(DIM, 0.0);
391 Vector<double> interpolated_dpdx(DIM, 0.0);
392
393 // Calculate pressure gradient
394 for (unsigned l = 0; l < n_pres; l++)
395 {
396 for (unsigned i = 0; i < DIM; i++)
397 {
398 interpolated_dpdx[i] += this->p_nst(l) * dpsip(l, i);
399 }
400 }
401
402 // Calculate velocities
403
404 // Loop over nodes
405 for (unsigned l = 0; l < n_node; l++)
406 {
407 // Loop over directions
408 for (unsigned i = 0; i < DIM; i++)
409 {
410 // Get the nodal value
411 double u_value = nodal_value(l, u_nodal_index[i]);
412 interpolated_u[i] += u_value * psif[l];
413 interpolated_x[i] += nodal_position(l, i) * psif[l];
414 }
415 }
416
417 // Source function (for validaton only)
418 double source = 0.0;
419 if (this->Press_adv_diff_source_fct_pt != 0)
420 {
421 source = this->Press_adv_diff_source_fct_pt(interpolated_x);
422 }
423
424
425 // Number of master nodes and storage for the weight of the shape function
426 unsigned n_master = 1;
427 double hang_weight = 1.0;
428
429
430 // Loop over the pressure shape functions
431 for (unsigned l = 0; l < n_pres; l++)
432 {
433 // If the pressure dof is hanging
434 if (pressure_dof_is_hanging[l])
435 {
436 // Pressure dof is hanging so it must be nodal-based
437 // Get the hang info object
438 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
439
440 // Get the number of master nodes from the pressure node
441 n_master = hang_info_pt->nmaster();
442 }
443 // Otherwise the node is its own master
444 else
445 {
446 n_master = 1;
447 }
448
449 // Loop over the master nodes
450 for (unsigned m = 0; m < n_master; m++)
451 {
452 // Get the number of the unknown
453 // If the pressure dof is hanging
454 if (pressure_dof_is_hanging[l])
455 {
456 // Get the local equation from the master node
457 local_eqn =
458 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
459 // Get the weight for the node
460 hang_weight = hang_info_pt->master_weight(m);
461 }
462 else
463 {
464 local_eqn = this->p_local_eqn(l);
465 hang_weight = 1.0;
466 }
467
468 // If the equation is not pinned
469 if (local_eqn >= 0)
470 {
471 residuals[local_eqn] -= source * testp[l] * W * hang_weight;
472 for (unsigned k = 0; k < DIM; k++)
473 {
474 residuals[local_eqn] +=
475 interpolated_dpdx[k] *
476 (scaled_re * interpolated_u[k] * testp[l] + dtestp(l, k)) * W *
477 hang_weight;
478 }
479
480 // Jacobian too?
481 if (flag)
482 {
483 // Number of master nodes and weights
484 unsigned n_master2 = 1;
485 double hang_weight2 = 1.0;
486
487 // Loop over the pressure shape functions
488 for (unsigned l2 = 0; l2 < n_pres; l2++)
489 {
490 // If the pressure dof is hanging
491 if (pressure_dof_is_hanging[l2])
492 {
493 hang_info2_pt =
494 this->pressure_node_pt(l2)->hanging_pt(p_index);
495 // Pressure dof is hanging so it must be nodal-based
496 // Get the number of master nodes from the pressure node
497 n_master2 = hang_info2_pt->nmaster();
498 }
499 // Otherwise the node is its own master
500 else
501 {
502 n_master2 = 1;
503 }
504
505 // Loop over the master nodes
506 for (unsigned m2 = 0; m2 < n_master2; m2++)
507 {
508 // Get the number of the unknown
509 // If the pressure dof is hanging
510 if (pressure_dof_is_hanging[l2])
511 {
512 // Get the unknown from the master node
513 local_unknown = this->local_hang_eqn(
514 hang_info2_pt->master_node_pt(m2), p_index);
515 // Get the weight from the hanging object
516 hang_weight2 = hang_info2_pt->master_weight(m2);
517 }
518 else
519 {
520 local_unknown = this->p_local_eqn(l2);
521 hang_weight2 = 1.0;
522 }
523
524 // If the unknown is not pinned
525 if (local_unknown >= 0)
526 {
527 if ((int(eqn_number(local_eqn)) !=
528 this->Pinned_fp_pressure_eqn) &&
529 (int(eqn_number(local_unknown)) !=
530 this->Pinned_fp_pressure_eqn))
531 {
532 for (unsigned k = 0; k < DIM; k++)
533 {
534 jacobian(local_eqn, local_unknown) +=
535 dtestp(l2, k) *
536 (scaled_re * interpolated_u[k] * testp[l] +
537 dtestp(l, k)) *
538 W * hang_weight * hang_weight2;
539 }
540 }
541 else
542 {
543 if ((int(eqn_number(local_eqn)) ==
544 this->Pinned_fp_pressure_eqn) &&
545 (int(eqn_number(local_unknown)) ==
546 this->Pinned_fp_pressure_eqn))
547 {
548 jacobian(local_eqn, local_unknown) = 1.0;
549 }
550 }
551 }
552 }
553 }
554 } /*End of Jacobian calculation*/
555 } // End of if not boundary condition
556 } // End of loop over master nodes
557 } // End of loop over l
558 } // end of integration loop
559
560 // Now add boundary contributions from Robin BCs
561 unsigned nrobin =
562 this->Pressure_advection_diffusion_robin_element_pt.size();
563 for (unsigned e = 0; e < nrobin; e++)
564 {
565 this->Pressure_advection_diffusion_robin_element_pt[e]
566 ->fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(
567 residuals, jacobian, flag);
568 }
569 }
570
571
572 //========================================================================
573 /// Add element's contribution to the elemental
574 /// residual vector and/or Jacobian matrix.
575 /// flag=1: compute both
576 /// flag=0: compute only residual vector
577 //========================================================================
578 template<unsigned DIM>
580 fill_in_generic_residual_contribution_nst(Vector<double>& residuals,
581 DenseMatrix<double>& jacobian,
582 DenseMatrix<double>& mass_matrix,
583 const unsigned& flag)
584 {
585 // Return immediately if there are no dofs
586 if (ndof() == 0) return;
587
588 // Find out how many nodes there are
589 unsigned n_node = nnode();
590
591 // Find out how many pressure dofs there are
592 unsigned n_pres = this->npres_nst();
593
594 // Allocate storage for the indices of the velocity components
595 unsigned u_nodal_index[DIM];
596
597 // Loop over the velocity components
598 for (unsigned i = 0; i < DIM; i++)
599 {
600 // Find the index at which the i-th local velocity is stored
601 u_nodal_index[i] = this->u_index_nst(i);
602 }
603
604 // Which nodal value represents the pressure? (Negative if pressure
605 // is not based on nodal interpolation).
606 int p_index = this->p_nodal_index_nst();
607
608 // Local array of booleans that are true if the l-th pressure value is
609 // hanging (avoid repeated virtual function calls)
610 bool pressure_dof_is_hanging[n_pres];
611
612 // If the pressure is stored at a node
613 if (p_index >= 0)
614 {
615 // Loop over the pressure nodes
616 for (unsigned l = 0; l < n_pres; ++l)
617 {
618 // Get the hang status of the l-th pressure node
619 pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
620 }
621 }
622 // Otherwise the pressure is not stored at a node and so cannot hang
623 else
624 {
625 // Loop over the pressure nodes
626 for (unsigned l = 0; l < n_pres; ++l)
627 {
628 // Indicate that the l-th pressure node is not hangings
629 pressure_dof_is_hanging[l] = false;
630 }
631 } // if (p_index>=0)
632
633 // Set up memory for the shape functions
634 Shape psif(n_node);
635
636 // Set up memory for the test functions
637 Shape testf(n_node);
638
639 // Allocate space for the derivatives of the shape functions
640 DShape dpsifdx(n_node, DIM + 1);
641
642 // Allocate space for the derivatives of the test functions
643 DShape dtestfdx(n_node, DIM + 1);
644
645 // Set up memory for pressure shape functions
646 Shape psip(n_pres);
647
648 // Set up memory for pressure test functions
649 Shape testp(n_pres);
650
651 // Number of integration points
652 unsigned n_intpt = integral_pt()->nweight();
653
654 // Set the Vector to hold local coordinates
655 Vector<double> s(DIM + 1, 0.0);
656
657 //-------------------------------------
658 // Get physical variables from element:
659 //-------------------------------------
660 // Reynolds number must be multiplied by the density ratio
661 double scaled_re = this->re() * this->density_ratio();
662
663 // Get the scaled Womersley value
664 double scaled_re_st = this->re() * this->st() * this->density_ratio();
665
666 // Get the scaled Womersley value differentiated w.r.t. the Strouhal number
667 double scaled_dre_st_dst = this->re() * this->density_ratio();
668
669 // Get the scaled Reynolds / Froude number
670 double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
671
672 // Get the viscosity ratio
673 double visc_ratio = this->viscosity_ratio();
674
675 // Get the gravity vector
676 Vector<double> G = this->g();
677
678 // Integer to store the local equation number
679 int local_eqn = 0;
680
681 // Integer to store the local unknown number
682 int local_unknown = 0;
683
684 // Pointers to first hang info object
685 HangInfo* hang_info_pt = 0;
686
687 // Pointers to second hang info object
688 HangInfo* hang_info2_pt = 0;
689
690 // Loop over the integration points
691 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
692 {
693 // Assign values of s
694 for (unsigned i = 0; i < DIM + 1; i++)
695 {
696 // Calculate the i-th local coordinate
697 s[i] = integral_pt()->knot(ipt, i);
698 }
699
700 // Get the integral weight
701 double w = integral_pt()->weight(ipt);
702
703 // Call the derivatives of the shape and test functions
704 double J = this->dshape_and_dtest_eulerian_at_knot_nst(
705 ipt, psif, dpsifdx, testf, dtestfdx);
706
707 // Call the pressure shape and test functions
708 this->pshape_nst(s, psip, testp);
709
710 // Pre-multiply the weights and the Jacobian
711 double W = w * J;
712
713 // Storage for the interpolated time value
714 double interpolated_t = 0.0;
715
716 // Storage for the interpolated pressure value
717 double interpolated_p = 0.0;
718
719 // Storage for the spatial coordinates
720 Vector<double> interpolated_x(DIM, 0.0);
721
722 // Storage for the interpolated velocity components
723 Vector<double> interpolated_u(DIM, 0.0);
724
725 // Storage for the interpolated time-derivative of the velocities
726 Vector<double> interpolated_dudt(DIM, 0.0);
727
728 // Storage for the mesh velocity
729 Vector<double> mesh_velocity(DIM, 0.0);
730
731 // Storage for the spatial derivatives of the velocity components
732 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
733
734 // Calculate pressure
735 for (unsigned l = 0; l < n_pres; l++)
736 {
737 // Update the current approximation to the interpolated pressure
738 interpolated_p += this->p_nst(l) * psip[l];
739 }
740
741 //-------------------------------------------------------------------
742 // Calculate velocities, derivatives, source function and body force:
743 //-------------------------------------------------------------------
744 // Loop over nodes
745 for (unsigned l = 0; l < n_node; l++)
746 {
747 // Update the interpolated time value
748 interpolated_t += this->nodal_position(l, DIM) * psif(l);
749
750 // Loop over coordinate directions
751 for (unsigned i = 0; i < DIM; i++)
752 {
753 // Get the nodal value
754 double u_value = this->nodal_value(l, u_nodal_index[i]);
755
756 // Update the i-th interpolated velocity component
757 interpolated_u[i] += u_value * psif[l];
758
759 // Update the i-th interpolated coordinate value
760 interpolated_x[i] += this->nodal_position(l, i) * psif[l];
761
762 // Update the interpolated du_i/dt value
763 interpolated_dudt[i] += u_value * dpsifdx(l, DIM);
764
765 // Loop over derivative directions
766 for (unsigned j = 0; j < DIM; j++)
767 {
768 // Update the interpolated du_i/dx_j value
769 interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
770 }
771 } // for (unsigned i=0;i<DIM;i++)
772 } // for (unsigned l=0;l<n_node;l++)
773
774 // If ALE is enabled
775 if (!(this->ALE_is_disabled))
776 {
777 // Loop over nodes
778 for (unsigned l = 0; l < n_node; l++)
779 {
780 // Loop over directions
781 for (unsigned i = 0; i < DIM; i++)
782 {
783 // Update the i-th mesh velocity component
784 mesh_velocity[i] += this->nodal_position(l, i) * dpsifdx(l, DIM);
785 }
786 } // for (unsigned l=0;l<n_node;l++)
787 } // if (!(this->ALE_is_disabled))
788
789 // Allocate space for the body force
790 Vector<double> body_force(DIM, 0.0);
791
792 // Get the user-defined body force term
793 this->get_body_force_nst(
794 interpolated_t, ipt, s, interpolated_x, body_force);
795
796 // Get the user-defined source function
797 double source = this->get_source_nst(interpolated_t, ipt, interpolated_x);
798
799 //---------------------------------
800 // Assemble residuals and Jacobian:
801 //---------------------------------
802 //--------------------
803 // Momentum equations:
804 //--------------------
805 // Number of master nodes
806 unsigned n_master = 1;
807
808 // Storage for the weight of the shape function
809 double hang_weight = 1.0;
810
811 // Loop over the test functions
812 for (unsigned l = 0; l < n_node; l++)
813 {
814 // Local boolean to indicate whether the node is hanging
815 bool is_node_hanging = node_pt(l)->is_hanging();
816
817 // If the node is hanging
818 if (is_node_hanging)
819 {
820 // Get the HangInfo object associated with the l-th node
821 hang_info_pt = node_pt(l)->hanging_pt();
822
823 // Read out number of master nodes from hanging data
824 n_master = hang_info_pt->nmaster();
825 }
826 // Otherwise the node is its own master
827 else
828 {
829 // There is only one node to consider; the node itself
830 n_master = 1;
831 }
832
833 // Loop over the master nodes
834 for (unsigned m = 0; m < n_master; m++)
835 {
836 // Loop over velocity components
837 for (unsigned i = 0; i < DIM; i++)
838 {
839 // Check if the node is hanging
840 if (is_node_hanging)
841 {
842 // Get the equation number from the master node
843 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
844 u_nodal_index[i]);
845
846 // Get the hang weight from the master node
847 hang_weight = hang_info_pt->master_weight(m);
848 }
849 // If the node is not hanging
850 else
851 {
852 // Local equation number
853 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
854
855 // Node contributes with full weight
856 hang_weight = 1.0;
857 }
858
859 // If it's not a boundary condition
860 if (local_eqn >= 0)
861 {
862 // Add the user-defined body force terms
863 residuals[local_eqn] +=
864 body_force[i] * testf[l] * W * hang_weight;
865
866 // Add the gravitational body force term
867 residuals[local_eqn] +=
868 scaled_re_inv_fr * testf[l] * G[i] * W * hang_weight;
869
870 // Add the pressure gradient term
871 residuals[local_eqn] +=
872 interpolated_p * dtestfdx(l, i) * W * hang_weight;
873
874 // Add in the contribution from the time derivative
875 residuals[local_eqn] -= (scaled_re_st * interpolated_dudt[i] *
876 testf[l] * W * hang_weight);
877
878 // If ALE is enabled
879 if (!(this->ALE_is_disabled))
880 {
881 // Loop over the coordinate directions
882 for (unsigned k = 0; k < DIM; k++)
883 {
884 // Add in the mesh velocity contribution
885 residuals[local_eqn] +=
886 (scaled_re_st * mesh_velocity[k] * interpolated_dudx(i, k) *
887 testf[l] * W * hang_weight);
888 }
889 } // if (!(this->ALE_is_disabled))
890
891 // Loop over the coordinate directions
892 for (unsigned k = 0; k < DIM; k++)
893 {
894 // Add in the convective term contribution
895 residuals[local_eqn] -=
896 (scaled_re * interpolated_u[k] * interpolated_dudx(i, k) *
897 testf[l] * W * hang_weight);
898 }
899
900 // Loop over the coordinate directions
901 for (unsigned k = 0; k < DIM; k++)
902 {
903 // Add in the stress tensor terms:
904 // NOTE: The viscosity ratio needs to go in here to ensure
905 // continuity of normal stress is satisfied even in flows
906 // with zero pressure gradient!
907 residuals[local_eqn] -=
908 ((interpolated_dudx(i, k) +
909 this->Gamma[i] * interpolated_dudx(k, i)) *
910 visc_ratio * dtestfdx(l, k) * W * hang_weight);
911 }
912
913 //------------------------
914 // Calculate the Jacobian:
915 //------------------------
916 // If we also need to construct the Jacobian
917 if (flag)
918 {
919 // Number of master nodes
920 unsigned n_master2 = 1;
921
922 // Storage for the weight of the shape function
923 double hang_weight2 = 1.0;
924
925 // Loop over the velocity nodes for columns
926 for (unsigned l2 = 0; l2 < n_node; l2++)
927 {
928 // Local boolean to indicate whether the node is hanging
929 bool is_node2_hanging = node_pt(l2)->is_hanging();
930
931 // If the node is hanging
932 if (is_node2_hanging)
933 {
934 // Get the HangInfo object associated with the l-th node
935 hang_info2_pt = node_pt(l2)->hanging_pt();
936
937 // Read out number of master nodes from hanging data
938 n_master2 = hang_info2_pt->nmaster();
939 }
940 // Otherwise the node is its own master
941 else
942 {
943 // There is only one node to consider; the node itself
944 n_master2 = 1;
945 }
946
947 // Loop over the master nodes
948 for (unsigned m2 = 0; m2 < n_master2; m2++)
949 {
950 // Loop over the velocity components
951 for (unsigned i2 = 0; i2 < DIM; i2++)
952 {
953 // Check if the node is hanging
954 if (is_node2_hanging)
955 {
956 // Get the equation number from the master node
957 local_unknown = this->local_hang_eqn(
958 hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
959
960 // Get the hang weight from the master node
961 hang_weight2 = hang_info2_pt->master_weight(m2);
962 }
963 // If the node is not hanging
964 else
965 {
966 // Get the local equation number
967 local_unknown =
968 this->nodal_local_eqn(l2, u_nodal_index[i2]);
969
970 // Node contributes with full weight
971 hang_weight2 = 1.0;
972 }
973
974 // If the unknown is non-zero
975 if (local_unknown >= 0)
976 {
977 // Add contribution to elemental matrix
978 jacobian(local_eqn, local_unknown) -=
979 (visc_ratio * this->Gamma[i] * dpsifdx(l2, i) *
980 dtestfdx(l, i2) * W * hang_weight * hang_weight2);
981
982 // Now add in the inertial terms
983 jacobian(local_eqn, local_unknown) -=
984 (scaled_re * psif[l2] * interpolated_dudx(i, i2) *
985 testf[l] * W * hang_weight * hang_weight2);
986
987 // Extra component if i2=i
988 if (i2 == i)
989 {
990 // If we also need to construct the mass matrix (only
991 // diagonal entries)
992 if (flag == 2)
993 {
994 // NOTE: This is positive because the mass matrix is
995 // taken to the other side of the equation when
996 // formulating the generalised eigenproblem.
997 mass_matrix(local_eqn, local_unknown) +=
998 (scaled_re_st * psif[l2] * testf[l] * W *
999 hang_weight * hang_weight2);
1000 }
1001
1002 // Add in the time-derivative contribution
1003 jacobian(local_eqn, local_unknown) -=
1004 (scaled_re_st * dpsifdx(l2, DIM) * testf[l] * W *
1005 hang_weight * hang_weight2);
1006
1007 // Loop over the velocity components
1008 for (unsigned k = 0; k < DIM; k++)
1009 {
1010 // Add in the convective term contribution
1011 jacobian(local_eqn, local_unknown) -=
1012 (scaled_re * interpolated_u[k] * dpsifdx(l2, k) *
1013 testf[l] * W * hang_weight * hang_weight2);
1014 }
1015
1016 // If ALE is enabled
1017 if (!(this->ALE_is_disabled))
1018 {
1019 // Loop over the velocity components
1020 for (unsigned k = 0; k < DIM; k++)
1021 {
1022 // Add in the mesh velocity contribution
1023 jacobian(local_eqn, local_unknown) +=
1024 (scaled_re_st * mesh_velocity[k] *
1025 dpsifdx(l2, k) * testf[l] * W * hang_weight *
1026 hang_weight2);
1027 }
1028 } // if (!(this->ALE_is_disabled))
1029
1030 // Loop over velocity components
1031 for (unsigned k = 0; k < DIM; k++)
1032 {
1033 // Add in the velocity gradient terms
1034 jacobian(local_eqn, local_unknown) -=
1035 (visc_ratio * dpsifdx(l2, k) * dtestfdx(l, k) *
1036 W * hang_weight * hang_weight2);
1037 }
1038 } // if (i2==i)
1039 } // if (local_unknown>=0)
1040 } // for (unsigned i2=0;i2<DIM;i2++)
1041 } // for (unsigned m2=0;m2<n_master2;m2++)
1042 } // for (unsigned l2=0;l2<n_node;l2++)
1043
1044 // Loop over the pressure shape functions
1045 for (unsigned l2 = 0; l2 < n_pres; l2++)
1046 {
1047 // If the pressure dof is hanging
1048 if (pressure_dof_is_hanging[l2])
1049 {
1050 // Get the HangInfo object associated with the l2-th
1051 // pressure node
1052 hang_info2_pt =
1053 this->pressure_node_pt(l2)->hanging_pt(p_index);
1054
1055 // Get the number of master nodes from the pressure node
1056 n_master2 = hang_info2_pt->nmaster();
1057 }
1058 // Otherwise the node is its own master
1059 else
1060 {
1061 // There is only one node to consider; the node itself
1062 n_master2 = 1;
1063 }
1064
1065 // Loop over the master nodes
1066 for (unsigned m2 = 0; m2 < n_master2; m2++)
1067 {
1068 // If the pressure dof is hanging
1069 if (pressure_dof_is_hanging[l2])
1070 {
1071 // Get the unknown from the master node
1072 local_unknown = this->local_hang_eqn(
1073 hang_info2_pt->master_node_pt(m2), p_index);
1074
1075 // Get the weight from the HangInfo object
1076 hang_weight2 = hang_info2_pt->master_weight(m2);
1077 }
1078 // If the node is not hanging
1079 else
1080 {
1081 // Get the local equation number
1082 local_unknown = this->p_local_eqn(l2);
1083
1084 // Node contributes with full weight
1085 hang_weight2 = 1.0;
1086 }
1087
1088 // If the unknown is not pinned
1089 if (local_unknown >= 0)
1090 {
1091 jacobian(local_eqn, local_unknown) +=
1092 (psip[l2] * dtestfdx(l, i) * W * hang_weight *
1093 hang_weight2);
1094 }
1095 } // for (unsigned m2=0;m2<n_master2;m2++)
1096 } // for (unsigned l2=0;l2<n_pres;l2++)
1097
1098 //------------------------------------
1099 // Calculate external data information
1100 //------------------------------------
1101 // If we're storing the Strouhal number as external data then we
1102 // need to calculate dr/d(St) in the elemental Jacobian
1103 if (this->Strouhal_is_stored_as_external_data)
1104 {
1105 // The index of the external data (there's only one!)
1106 unsigned data_index = 0;
1107
1108 // The index of the unknown value stored in the external data
1109 unsigned value_index = 0;
1110
1111 // Get the local equation number associated with the extra
1112 // unknown
1113 local_unknown =
1114 this->external_local_eqn(data_index, value_index);
1115
1116 // If we're at a non-zero degree of freedom add in the entry
1117 if (local_unknown >= 0)
1118 {
1119 // Add in the contribution from the time derivative
1120 jacobian(local_eqn, local_unknown) -=
1121 (scaled_dre_st_dst * interpolated_dudt[i] * testf[l] * W *
1122 hang_weight);
1123
1124 // If ALE is enabled
1125 if (!this->ALE_is_disabled)
1126 {
1127 // Loop over the coordinate directions
1128 for (unsigned k = 0; k < DIM; k++)
1129 {
1130 // Add in the mesh velocity contribution
1131 jacobian(local_eqn, local_unknown) +=
1132 (scaled_dre_st_dst * mesh_velocity[k] *
1133 interpolated_dudx(i, k) * testf[l] * W *
1134 hang_weight);
1135 }
1136 } // if (!ALE_is_disabled)
1137 } // if (local_unknown>=0)
1138 } // if (Strouhal_is_stored_as_external_data)
1139 } // if (flag)
1140 } // if (local_eqn>=0)
1141 } // for (unsigned i=0;i<DIM;i++)
1142 } // for (unsigned m=0;m<n_master;m++)
1143 } // for (unsigned l=0;l<n_node;l++)
1144
1145 //---------------------
1146 // Continuity equation:
1147 //---------------------
1148 // Loop over the pressure shape functions
1149 for (unsigned l = 0; l < n_pres; l++)
1150 {
1151 // If the pressure dof is hanging
1152 if (pressure_dof_is_hanging[l])
1153 {
1154 // Get the HangInfo object associated with the l2-th pressure node
1155 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1156
1157 // Get the number of master nodes from the pressure node
1158 n_master = hang_info_pt->nmaster();
1159 }
1160 // Otherwise the node is its own master
1161 else
1162 {
1163 // There is only one node to consider; the node itself
1164 n_master = 1;
1165 }
1166
1167 // Loop over the master nodes
1168 for (unsigned m = 0; m < n_master; m++)
1169 {
1170 // If the pressure dof is hanging
1171 if (pressure_dof_is_hanging[l])
1172 {
1173 // Get the unknown from the master node
1174 local_eqn =
1175 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
1176
1177 // Get the weight from the HangInfo object
1178 hang_weight = hang_info_pt->master_weight(m);
1179 }
1180 // If the node is not hanging
1181 else
1182 {
1183 // Get the local equation number
1184 local_eqn = this->p_local_eqn(l);
1185
1186 // Node contributes with full weight
1187 hang_weight = 1.0;
1188 }
1189
1190 // If the equation is not pinned
1191 if (local_eqn >= 0)
1192 {
1193 // Add in the source term contribution
1194 residuals[local_eqn] -= source * testp[l] * W * hang_weight;
1195
1196 // Loop over velocity components
1197 for (unsigned k = 0; k < DIM; k++)
1198 {
1199 // Add in the velocity gradient terms
1200 residuals[local_eqn] +=
1201 interpolated_dudx(k, k) * testp[l] * W * hang_weight;
1202 }
1203
1204 //------------------------
1205 // Calculate the Jacobian:
1206 //------------------------
1207 // If we also need to construct the Jacobian
1208 if (flag)
1209 {
1210 // Number of master nodes
1211 unsigned n_master2 = 1;
1212
1213 // Storage for the weight of the shape function
1214 double hang_weight2 = 1.0;
1215
1216 // Loop over the velocity nodes for columns
1217 for (unsigned l2 = 0; l2 < n_node; l2++)
1218 {
1219 // Local boolean to indicate whether the node is hanging
1220 bool is_node2_hanging = node_pt(l2)->is_hanging();
1221
1222 // If the node is hanging
1223 if (is_node2_hanging)
1224 {
1225 // Get the HangInfo object associated with the l2-th node
1226 hang_info2_pt = node_pt(l2)->hanging_pt();
1227
1228 // Read out number of master nodes from hanging data
1229 n_master2 = hang_info2_pt->nmaster();
1230 }
1231 // Otherwise the node is its own master
1232 else
1233 {
1234 // There is only one node to consider; the node itself
1235 n_master2 = 1;
1236 }
1237
1238 // Loop over the master nodes
1239 for (unsigned m2 = 0; m2 < n_master2; m2++)
1240 {
1241 // Loop over the velocity components
1242 for (unsigned i2 = 0; i2 < DIM; i2++)
1243 {
1244 // If the node is hanging
1245 if (is_node2_hanging)
1246 {
1247 // Get the equation number from the master node
1248 local_unknown = this->local_hang_eqn(
1249 hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
1250
1251 // Get the weight from the HangInfo object
1252 hang_weight2 = hang_info2_pt->master_weight(m2);
1253 }
1254 // If the node is not hanging
1255 else
1256 {
1257 // Get the local equation number
1258 local_unknown =
1259 this->nodal_local_eqn(l2, u_nodal_index[i2]);
1260
1261 // Node contributes with full weight
1262 hang_weight2 = 1.0;
1263 }
1264
1265 // If the unknown is not pinned
1266 if (local_unknown >= 0)
1267 {
1268 // Add in the velocity gradient contribution
1269 jacobian(local_eqn, local_unknown) +=
1270 (dpsifdx(l2, i2) * testp[l] * W * hang_weight *
1271 hang_weight2);
1272 }
1273 } // for (unsigned i2=0;i2<DIM;i2++)s
1274 } // for (unsigned m2=0;m2<n_master2;m2++)
1275 } // for (unsigned l2=0;l2<n_node;l2++)
1276 } // if (flag)
1277 } // if (local_eqn>=0)
1278 } // for (unsigned m=0;m<n_master;m++)
1279 } // for (unsigned l=0;l<n_pres;l++)
1280 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
1281 } // End of fill_in_generic_residual_contribution_nst
1282
1283
1284 //======================================================================
1285 /// Compute derivatives of elemental residual vector with respect
1286 /// to nodal coordinates.
1287 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
1288 /// Overloads the FD-based version in the FE base class.
1289 //======================================================================
1290 template<unsigned DIM>
1291 void RefineableSpaceTimeNavierStokesEquations<
1292 DIM>::get_dresidual_dnodal_coordinates(RankThreeTensor<double>&
1293 dresidual_dnodal_coordinates)
1294 {
1295 // Throw a warning
1296 throw OomphLibError("Space-time update needs to be checked!",
1297 OOMPH_CURRENT_FUNCTION,
1298 OOMPH_EXCEPTION_LOCATION);
1299
1300 // Return immediately if there are no dofs
1301 if (ndof() == 0)
1302 {
1303 return;
1304 }
1305
1306 // Determine number of nodes in element
1307 const unsigned n_node = nnode();
1308
1309 // Get continuous time from timestepper of first node
1310 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
1311
1312 // Determine number of pressure dofs in element
1313 const unsigned n_pres = this->npres_nst();
1314
1315 // Find the indices at which the local velocities are stored
1316 unsigned u_nodal_index[DIM];
1317 for (unsigned i = 0; i < DIM; i++)
1318 {
1319 u_nodal_index[i] = this->u_index_nst(i);
1320 }
1321
1322 // Which nodal value represents the pressure? (Negative if pressure
1323 // is not based on nodal interpolation).
1324 const int p_index = this->p_nodal_index_nst();
1325
1326 // Local array of booleans that are true if the l-th pressure value is
1327 // hanging (avoid repeated virtual function calls)
1328 bool pressure_dof_is_hanging[n_pres];
1329
1330 // If the pressure is stored at a node
1331 if (p_index >= 0)
1332 {
1333 // Read out whether the pressure is hanging
1334 for (unsigned l = 0; l < n_pres; ++l)
1335 {
1336 pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
1337 }
1338 }
1339 // Otherwise the pressure is not stored at a node and so cannot hang
1340 else
1341 {
1342 for (unsigned l = 0; l < n_pres; ++l)
1343 {
1344 pressure_dof_is_hanging[l] = false;
1345 }
1346 }
1347
1348 // Set up memory for the shape and test functions
1349 Shape psif(n_node), testf(n_node);
1350 DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
1351
1352 // Set up memory for pressure shape and test functions
1353 Shape psip(n_pres), testp(n_pres);
1354
1355 // Determine number of shape controlling nodes
1356 const unsigned n_shape_controlling_node = nshape_controlling_nodes();
1357
1358 // Deriatives of shape fct derivatives w.r.t. nodal coords
1359 RankFourTensor<double> d_dpsifdx_dX(
1360 DIM, n_shape_controlling_node, n_node, DIM);
1361 RankFourTensor<double> d_dtestfdx_dX(
1362 DIM, n_shape_controlling_node, n_node, DIM);
1363
1364 // Derivative of Jacobian of mapping w.r.t. to nodal coords
1365 DenseMatrix<double> dJ_dX(DIM, n_shape_controlling_node);
1366
1367 // Derivatives of derivative of u w.r.t. nodal coords
1368 RankFourTensor<double> d_dudx_dX(DIM, n_shape_controlling_node, DIM, DIM);
1369
1370 // Derivatives of nodal velocities w.r.t. nodal coords:
1371 // Assumption: Interaction only local via no-slip so
1372 // X_ij only affects U_ij.
1373 DenseMatrix<double> d_U_dX(DIM, n_shape_controlling_node, 0.0);
1374
1375 // Determine the number of integration points
1376 const unsigned n_intpt = integral_pt()->nweight();
1377
1378 // Vector to hold local coordinates
1379 Vector<double> s(DIM);
1380
1381 // Get physical variables from element
1382 // (Reynolds number must be multiplied by the density ratio)
1383 double scaled_re = this->re() * this->density_ratio();
1384 double scaled_re_st = this->re_st() * this->density_ratio();
1385 double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
1386 double visc_ratio = this->viscosity_ratio();
1387 Vector<double> G = this->g();
1388
1389 // FD step
1391
1392 // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1393 // Assumption: Interaction only local via no-slip so
1394 // X_ij only affects U_ij.
1395 bool element_has_node_with_aux_node_update_fct = false;
1396
1397 std::map<Node*, unsigned> local_shape_controlling_node_lookup =
1398 shape_controlling_node_lookup();
1399
1400 // FD loop over shape-controlling nodes
1401 for (std::map<Node*, unsigned>::iterator it =
1402 local_shape_controlling_node_lookup.begin();
1403 it != local_shape_controlling_node_lookup.end();
1404 it++)
1405 {
1406 // Get node
1407 Node* nod_pt = it->first;
1408
1409 // Get its number
1410 unsigned q = it->second;
1411
1412 // Only compute if there's a node-update fct involved
1413 if (nod_pt->has_auxiliary_node_update_fct_pt())
1414 {
1415 element_has_node_with_aux_node_update_fct = true;
1416
1417 // Current nodal velocity
1418 Vector<double> u_ref(DIM);
1419 for (unsigned i = 0; i < DIM; i++)
1420 {
1421 u_ref[i] = *(nod_pt->value_pt(u_nodal_index[i]));
1422 }
1423
1424 // FD
1425 for (unsigned p = 0; p < DIM; p++)
1426 {
1427 // Make backup
1428 double backup = nod_pt->x(p);
1429
1430 // Do FD step. No node update required as we're
1431 // attacking the coordinate directly...
1432 nod_pt->x(p) += eps_fd;
1433
1434 // Do auxiliary node update (to apply no slip)
1435 nod_pt->perform_auxiliary_node_update_fct();
1436
1437 // Evaluate
1438 d_U_dX(p, q) =
1439 (*(nod_pt->value_pt(u_nodal_index[p])) - u_ref[p]) / eps_fd;
1440
1441 // Reset
1442 nod_pt->x(p) = backup;
1443
1444 // Do auxiliary node update (to apply no slip)
1445 nod_pt->perform_auxiliary_node_update_fct();
1446 }
1447 }
1448 }
1449
1450 // Integer to store the local equation number
1451 int local_eqn = 0;
1452
1453 // Pointers to hang info object
1454 HangInfo* hang_info_pt = 0;
1455
1456 // Loop over the integration points
1457 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1458 {
1459 // Assign values of s
1460 for (unsigned i = 0; i < DIM; i++)
1461 {
1462 s[i] = integral_pt()->knot(ipt, i);
1463 }
1464
1465 // Get the integral weight
1466 const double w = integral_pt()->weight(ipt);
1467
1468 // Call the derivatives of the shape and test functions
1469 const double J =
1470 this->dshape_and_dtest_eulerian_at_knot_nst(ipt,
1471 psif,
1472 dpsifdx,
1473 d_dpsifdx_dX,
1474 testf,
1475 dtestfdx,
1476 d_dtestfdx_dX,
1477 dJ_dX);
1478
1479 // Call the pressure shape and test functions
1480 this->pshape_nst(s, psip, testp);
1481
1482 // Calculate local values of the pressure and velocity components
1483 // Allocate
1484 double interpolated_p = 0.0;
1485 Vector<double> interpolated_u(DIM, 0.0);
1486 Vector<double> interpolated_x(DIM, 0.0);
1487 Vector<double> mesh_velocity(DIM, 0.0);
1488 Vector<double> dudt(DIM, 0.0);
1489 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1490
1491 // Calculate pressure
1492 for (unsigned l = 0; l < n_pres; l++)
1493 {
1494 interpolated_p += this->p_nst(l) * psip[l];
1495 }
1496
1497 // Calculate velocities and derivatives:
1498
1499 // Loop over nodes
1500 for (unsigned l = 0; l < n_node; l++)
1501 {
1502 // Loop over directions
1503 for (unsigned i = 0; i < DIM; i++)
1504 {
1505 // Get the nodal value
1506 const double u_value = nodal_value(l, u_nodal_index[i]);
1507 interpolated_u[i] += u_value * psif[l];
1508 interpolated_x[i] += nodal_position(l, i) * psif[l];
1509 dudt[i] += this->du_dt_nst(l, i) * psif[l];
1510
1511 // Loop over derivative directions
1512 for (unsigned j = 0; j < DIM; j++)
1513 {
1514 interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1515 }
1516 }
1517 }
1518
1519 if (!this->ALE_is_disabled)
1520 {
1521 // Loop over nodes
1522 for (unsigned l = 0; l < n_node; l++)
1523 {
1524 // Loop over directions
1525 for (unsigned i = 0; i < DIM; i++)
1526 {
1527 mesh_velocity[i] += this->dnodal_position_dt(l, i) * psif[l];
1528 }
1529 }
1530 }
1531
1532 // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1533
1534 // Loop over shape-controlling nodes
1535 for (unsigned q = 0; q < n_shape_controlling_node; q++)
1536 {
1537 // Loop over coordinate directions
1538 for (unsigned p = 0; p < DIM; p++)
1539 {
1540 for (unsigned i = 0; i < DIM; i++)
1541 {
1542 for (unsigned k = 0; k < DIM; k++)
1543 {
1544 double aux = 0.0;
1545 for (unsigned j = 0; j < n_node; j++)
1546 {
1547 aux +=
1548 nodal_value(j, u_nodal_index[i]) * d_dpsifdx_dX(p, q, j, k);
1549 }
1550 d_dudx_dX(p, q, i, k) = aux;
1551 }
1552 }
1553 }
1554 }
1555
1556 // Get weight of actual nodal position/value in computation of mesh
1557 // velocity from positional/value time stepper
1558 const double pos_time_weight =
1559 node_pt(0)->position_time_stepper_pt()->weight(1, 0);
1560 const double val_time_weight =
1561 node_pt(0)->time_stepper_pt()->weight(1, 0);
1562
1563 // Get the user-defined body force terms
1564 Vector<double> body_force(DIM);
1565 this->get_body_force_nst(time, ipt, s, interpolated_x, body_force);
1566
1567 // Get the user-defined source function
1568 const double source = this->get_source_nst(time, ipt, interpolated_x);
1569
1570 // Get gradient of body force function
1571 DenseMatrix<double> d_body_force_dx(DIM, DIM, 0.0);
1572 this->get_body_force_gradient_nst(
1573 time, ipt, s, interpolated_x, d_body_force_dx);
1574
1575 // Get gradient of source function
1576 Vector<double> source_gradient(DIM, 0.0);
1577 this->get_source_gradient_nst(time, ipt, interpolated_x, source_gradient);
1578
1579
1580 // Assemble shape derivatives
1581 //---------------------------
1582
1583 // MOMENTUM EQUATIONS
1584 // ------------------
1585
1586 // Number of master nodes and storage for the weight of the shape function
1587 unsigned n_master = 1;
1588 double hang_weight = 1.0;
1589
1590 // Loop over the test functions
1591 for (unsigned l = 0; l < n_node; l++)
1592 {
1593 // Local boolean to indicate whether the node is hanging
1594 bool is_node_hanging = node_pt(l)->is_hanging();
1595
1596 // If the node is hanging
1597 if (is_node_hanging)
1598 {
1599 hang_info_pt = node_pt(l)->hanging_pt();
1600
1601 // Read out number of master nodes from hanging data
1602 n_master = hang_info_pt->nmaster();
1603 }
1604 // Otherwise the node is its own master
1605 else
1606 {
1607 n_master = 1;
1608 }
1609
1610 // Loop over the master nodes
1611 for (unsigned m = 0; m < n_master; m++)
1612 {
1613 // Loop over coordinate directions
1614 for (unsigned i = 0; i < DIM; i++)
1615 {
1616 // Get the equation number
1617 // If the node is hanging
1618 if (is_node_hanging)
1619 {
1620 // Get the equation number from the master node
1621 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1622 u_nodal_index[i]);
1623 // Get the hang weight from the master node
1624 hang_weight = hang_info_pt->master_weight(m);
1625 }
1626 // If the node is not hanging
1627 else
1628 {
1629 // Local equation number
1630 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
1631
1632 // Node contributes with full weight
1633 hang_weight = 1.0;
1634 }
1635
1636 // IF it's not a boundary condition
1637 if (local_eqn >= 0)
1638 {
1639 // Loop over coordinate directions
1640 for (unsigned p = 0; p < DIM; p++)
1641 {
1642 // Loop over shape controlling nodes
1643 for (unsigned q = 0; q < n_shape_controlling_node; q++)
1644 {
1645 // Residual x deriv of Jacobian
1646 // ----------------------------
1647
1648 // Add the user-defined body force terms
1649 double sum = body_force[i] * testf[l];
1650
1651 // Add the gravitational body force term
1652 sum += scaled_re_inv_fr * testf[l] * G[i];
1653
1654 // Add the pressure gradient term
1655 sum += interpolated_p * dtestfdx(l, i);
1656
1657 // Add in the stress tensor terms
1658 // The viscosity ratio needs to go in here to ensure
1659 // continuity of normal stress is satisfied even in flows
1660 // with zero pressure gradient!
1661 for (unsigned k = 0; k < DIM; k++)
1662 {
1663 sum -= visc_ratio *
1664 (interpolated_dudx(i, k) +
1665 this->Gamma[i] * interpolated_dudx(k, i)) *
1666 dtestfdx(l, k);
1667 }
1668
1669 // Add in the inertial terms
1670
1671 // du/dt term
1672 sum -= scaled_re_st * dudt[i] * testf[l];
1673
1674 // Convective terms, including mesh velocity
1675 for (unsigned k = 0; k < DIM; k++)
1676 {
1677 double tmp = scaled_re * interpolated_u[k];
1678 if (!this->ALE_is_disabled)
1679 {
1680 tmp -= scaled_re_st * mesh_velocity[k];
1681 }
1682 sum -= tmp * interpolated_dudx(i, k) * testf[l];
1683 }
1684
1685 // Multiply throsugh by deriv of Jacobian and integration
1686 // weight
1687 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1688 sum * dJ_dX(p, q) * w * hang_weight;
1689
1690 // Derivative of residual x Jacobian
1691 // ---------------------------------
1692
1693 // Body force
1694 sum = d_body_force_dx(i, p) * psif(q) * testf(l);
1695
1696 // Pressure gradient term
1697 sum += interpolated_p * d_dtestfdx_dX(p, q, l, i);
1698
1699 // Viscous term
1700 for (unsigned k = 0; k < DIM; k++)
1701 {
1702 sum -=
1703 visc_ratio * ((interpolated_dudx(i, k) +
1704 this->Gamma[i] * interpolated_dudx(k, i)) *
1705 d_dtestfdx_dX(p, q, l, k) +
1706 (d_dudx_dX(p, q, i, k) +
1707 this->Gamma[i] * d_dudx_dX(p, q, k, i)) *
1708 dtestfdx(l, k));
1709 }
1710
1711 // Convective terms, including mesh velocity
1712 for (unsigned k = 0; k < DIM; k++)
1713 {
1714 double tmp = scaled_re * interpolated_u[k];
1715 if (!this->ALE_is_disabled)
1716 {
1717 tmp -= scaled_re_st * mesh_velocity[k];
1718 }
1719 sum -= tmp * d_dudx_dX(p, q, i, k) * testf(l);
1720 }
1721 if (!this->ALE_is_disabled)
1722 {
1723 sum += scaled_re_st * pos_time_weight * psif(q) *
1724 interpolated_dudx(i, p) * testf(l);
1725 }
1726
1727 // Multiply through by Jacobian and integration weight
1728 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1729 sum * J * w * hang_weight;
1730
1731 } // End of loop over shape controlling nodes q
1732 } // End of loop over coordinate directions p
1733
1734
1735 // Derivs w.r.t. to nodal velocities
1736 // ---------------------------------
1737 if (element_has_node_with_aux_node_update_fct)
1738 {
1739 // Loop over local nodes
1740 for (unsigned q_local = 0; q_local < n_node; q_local++)
1741 {
1742 // Number of master nodes and storage for the weight of
1743 // the shape function
1744 unsigned n_master2 = 1;
1745 double hang_weight2 = 1.0;
1746 HangInfo* hang_info2_pt = 0;
1747
1748 // Local boolean to indicate whether the node is hanging
1749 bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1750
1751 Node* actual_shape_controlling_node_pt = node_pt(q_local);
1752
1753 // If the node is hanging
1754 if (is_node_hanging2)
1755 {
1756 hang_info2_pt = node_pt(q_local)->hanging_pt();
1757
1758 // Read out number of master nodes from hanging data
1759 n_master2 = hang_info2_pt->nmaster();
1760 }
1761 // Otherwise the node is its own master
1762 else
1763 {
1764 n_master2 = 1;
1765 }
1766
1767 // Loop over the master nodes
1768 for (unsigned mm = 0; mm < n_master2; mm++)
1769 {
1770 if (is_node_hanging2)
1771 {
1772 actual_shape_controlling_node_pt =
1773 hang_info2_pt->master_node_pt(mm);
1774 hang_weight2 = hang_info2_pt->master_weight(mm);
1775 }
1776
1777 // Find the corresponding number
1778 unsigned q = local_shape_controlling_node_lookup
1779 [actual_shape_controlling_node_pt];
1780
1781 // Loop over coordinate directions
1782 for (unsigned p = 0; p < DIM; p++)
1783 {
1784 double sum = -visc_ratio * this->Gamma[i] *
1785 dpsifdx(q_local, i) * dtestfdx(l, p) -
1786 scaled_re * psif(q_local) *
1787 interpolated_dudx(i, p) * testf(l);
1788 if (i == p)
1789 {
1790 sum -= scaled_re_st * val_time_weight * psif(q_local) *
1791 testf(l);
1792 for (unsigned k = 0; k < DIM; k++)
1793 {
1794 sum -=
1795 visc_ratio * dpsifdx(q_local, k) * dtestfdx(l, k);
1796 double tmp = scaled_re * interpolated_u[k];
1797 if (!this->ALE_is_disabled)
1798 {
1799 tmp -= scaled_re_st * mesh_velocity[k];
1800 }
1801 sum -= tmp * dpsifdx(q_local, k) * testf(l);
1802 }
1803 }
1804
1805 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1806 sum * d_U_dX(p, q) * J * w * hang_weight * hang_weight2;
1807 }
1808 } // End of loop over master nodes
1809 } // End of loop over local nodes
1810 } // End of if(element_has_node_with_aux_node_update_fct)
1811
1812
1813 } // local_eqn>=0
1814 }
1815 }
1816 } // End of loop over test functions
1817
1818
1819 // CONTINUITY EQUATION
1820 // -------------------
1821
1822 // Loop over the shape functions
1823 for (unsigned l = 0; l < n_pres; l++)
1824 {
1825 // If the pressure dof is hanging
1826 if (pressure_dof_is_hanging[l])
1827 {
1828 // Pressure dof is hanging so it must be nodal-based
1829 // Get the hang info object
1830 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1831
1832 // Get the number of master nodes from the pressure node
1833 n_master = hang_info_pt->nmaster();
1834 }
1835 // Otherwise the node is its own master
1836 else
1837 {
1838 n_master = 1;
1839 }
1840
1841 // Loop over the master nodes
1842 for (unsigned m = 0; m < n_master; m++)
1843 {
1844 // Get the number of the unknown
1845 // If the pressure dof is hanging
1846 if (pressure_dof_is_hanging[l])
1847 {
1848 // Get the local equation from the master node
1849 local_eqn =
1850 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
1851 // Get the weight for the node
1852 hang_weight = hang_info_pt->master_weight(m);
1853 }
1854 else
1855 {
1856 local_eqn = this->p_local_eqn(l);
1857 hang_weight = 1.0;
1858 }
1859
1860 // If not a boundary conditions
1861 if (local_eqn >= 0)
1862 {
1863 // Loop over coordinate directions
1864 for (unsigned p = 0; p < DIM; p++)
1865 {
1866 // Loop over nodes
1867 for (unsigned q = 0; q < n_shape_controlling_node; q++)
1868 {
1869 // Residual x deriv of Jacobian
1870 //-----------------------------
1871
1872 // Source term
1873 double aux = -source;
1874
1875 // Loop over velocity components
1876 for (unsigned k = 0; k < DIM; k++)
1877 {
1878 aux += interpolated_dudx(k, k);
1879 }
1880
1881 // Multiply through by deriv of Jacobian and integration weight
1882 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1883 aux * dJ_dX(p, q) * testp[l] * w * hang_weight;
1884
1885
1886 // Derivative of residual x Jacobian
1887 // ---------------------------------
1888
1889 // Loop over velocity components
1890 aux = -source_gradient[p] * psif(q);
1891 for (unsigned k = 0; k < DIM; k++)
1892 {
1893 aux += d_dudx_dX(p, q, k, k);
1894 }
1895 // Multiply through by Jacobian and integration weight
1896 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1897 aux * testp[l] * J * w * hang_weight;
1898 }
1899 }
1900
1901
1902 // Derivs w.r.t. to nodal velocities
1903 // ---------------------------------
1904 if (element_has_node_with_aux_node_update_fct)
1905 {
1906 // Loop over local nodes
1907 for (unsigned q_local = 0; q_local < n_node; q_local++)
1908 {
1909 // Number of master nodes and storage for the weight of
1910 // the shape function
1911 unsigned n_master2 = 1;
1912 double hang_weight2 = 1.0;
1913 HangInfo* hang_info2_pt = 0;
1914
1915 // Local boolean to indicate whether the node is hanging
1916 bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1917
1918 Node* actual_shape_controlling_node_pt = node_pt(q_local);
1919
1920 // If the node is hanging
1921 if (is_node_hanging2)
1922 {
1923 hang_info2_pt = node_pt(q_local)->hanging_pt();
1924
1925 // Read out number of master nodes from hanging data
1926 n_master2 = hang_info2_pt->nmaster();
1927 }
1928 // Otherwise the node is its own master
1929 else
1930 {
1931 n_master2 = 1;
1932 }
1933
1934 // Loop over the master nodes
1935 for (unsigned mm = 0; mm < n_master2; mm++)
1936 {
1937 if (is_node_hanging2)
1938 {
1939 actual_shape_controlling_node_pt =
1940 hang_info2_pt->master_node_pt(mm);
1941 hang_weight2 = hang_info2_pt->master_weight(mm);
1942 }
1943
1944 // Find the corresponding number
1945 unsigned q = local_shape_controlling_node_lookup
1946 [actual_shape_controlling_node_pt];
1947
1948 // Loop over coordinate directions
1949 for (unsigned p = 0; p < DIM; p++)
1950 {
1951 double aux = dpsifdx(q_local, p) * testp(l);
1952 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1953 aux * d_U_dX(p, q) * J * w * hang_weight * hang_weight2;
1954 }
1955 } // End of loop over (mm) master nodes
1956 } // End of loop over local nodes q_local
1957 } // End of if(element_has_node_with_aux_node_update_fct)
1958 } // End of if it's not a boundary condition
1959 } // End of loop over (m) master nodes
1960 } // End of loop over shape functions for continuity eqn
1961
1962 } // End of loop over integration points
1963 }
1964
1965
1966 //====================================================================
1967 /// Force build of templates
1968 //====================================================================
1969 template class RefineableSpaceTimeNavierStokesEquations<2>;
1970 template class RefineableQTaylorHoodSpaceTimeElement<2>;
1971} // namespace oomph
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1198
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 fill_in_generic_pressure_advection_diffusion_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &compute_jacobian_flag)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
void fill_in_generic_residual_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, const unsigned &compute_jacobian_flag)
Add the elements contribution to elemental residual vector and/or Jacobian matrix....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...