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