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 /// Compute the residuals for the associated pressure advection
286 /// diffusion problem. Used by the Fp preconditioner.
287 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
288 //==============================================================
289 template<unsigned DIM>
292 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
293 {
294 // Return immediately if there are no dofs
295 if (ndof() == 0) return;
296
297 // Find out how many nodes there are
298 unsigned n_node = nnode();
299
300 // Find out how many pressure dofs there are
301 unsigned n_pres = this->npres_nst();
302
303 // Find the indices at which the local velocities are stored
304 unsigned u_nodal_index[DIM];
305 for (unsigned i = 0; i < DIM; i++)
306 {
307 u_nodal_index[i] = this->u_index_nst(i);
308 }
309
310
311 // Which nodal value represents the pressure? (Negative if pressure
312 // is not based on nodal interpolation).
313 int p_index = this->p_nodal_index_nst();
314
315 // Local array of booleans that are true if the l-th pressure value is
316 // hanging (avoid repeated virtual function calls)
317 bool pressure_dof_is_hanging[n_pres];
318 // If the pressure is stored at a node
319 if (p_index >= 0)
320 {
321 // Read out whether the pressure is hanging
322 for (unsigned l = 0; l < n_pres; ++l)
323 {
324 pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
325 }
326 }
327 // Otherwise the pressure is not stored at a node and so cannot hang
328 else
329 {
330 // pressure advection diffusion doesn't work for this one!
331 throw OomphLibError(
332 "Pressure advection diffusion does not work in this case\n",
333 OOMPH_CURRENT_FUNCTION,
334 OOMPH_EXCEPTION_LOCATION);
335
336 for (unsigned l = 0; l < n_pres; ++l)
337 {
338 pressure_dof_is_hanging[l] = false;
339 }
340 }
341
342 // Set up memory for the velocity shape fcts
343 Shape psif(n_node);
344 DShape dpsidx(n_node, DIM);
345
346 // Set up memory for pressure shape and test functions
347 Shape psip(n_pres), testp(n_pres);
348 DShape dpsip(n_pres, DIM);
349 DShape dtestp(n_pres, DIM);
350
351 // Number of integration points
352 unsigned n_intpt = integral_pt()->nweight();
353
354 // Set the Vector to hold local coordinates
355 Vector<double> s(DIM);
356
357 // Get Physical Variables from Element
358 // Reynolds number must be multiplied by the density ratio
359 double scaled_re = this->re() * this->density_ratio();
360
361 // Integers to store the local equations and unknowns
362 int local_eqn = 0, local_unknown = 0;
363
364 // Pointers to hang info objects
365 HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
366
367 // Loop over the integration points
368 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
369 {
370 // Assign values of s
371 for (unsigned i = 0; i < DIM; i++) s[i] = integral_pt()->knot(ipt, i);
372
373 // Get the integral weight
374 double w = integral_pt()->weight(ipt);
375
376 // Call the derivatives of the veloc shape functions
377 // (Derivs not needed but they are free)
378 double J = this->dshape_eulerian_at_knot(ipt, psif, dpsidx);
379
380 // Call the pressure shape and test functions
381 this->dpshape_and_dptest_eulerian_nst(s, psip, dpsip, testp, dtestp);
382
383 // Premultiply the weights and the Jacobian
384 double W = w * J;
385
386 // Calculate local values of the pressure and velocity components
387 // Allocate
388 Vector<double> interpolated_u(DIM, 0.0);
389 Vector<double> interpolated_x(DIM, 0.0);
390 Vector<double> interpolated_dpdx(DIM, 0.0);
391
392 // Calculate pressure gradient
393 for (unsigned l = 0; l < n_pres; l++)
394 {
395 for (unsigned i = 0; i < DIM; i++)
396 {
397 interpolated_dpdx[i] += this->p_nst(l) * dpsip(l, i);
398 }
399 }
400
401 // Calculate velocities
402
403 // Loop over nodes
404 for (unsigned l = 0; l < n_node; l++)
405 {
406 // Loop over directions
407 for (unsigned i = 0; i < DIM; i++)
408 {
409 // Get the nodal value
410 double u_value = nodal_value(l, u_nodal_index[i]);
411 interpolated_u[i] += u_value * psif[l];
412 interpolated_x[i] += nodal_position(l, i) * psif[l];
413 }
414 }
415
416 // Source function (for validaton only)
417 double source = 0.0;
418 if (this->Press_adv_diff_source_fct_pt != 0)
419 {
420 source = this->Press_adv_diff_source_fct_pt(interpolated_x);
421 }
422
423
424 // Number of master nodes and storage for the weight of the shape function
425 unsigned n_master = 1;
426 double hang_weight = 1.0;
427
428
429 // Loop over the pressure shape functions
430 for (unsigned l = 0; l < n_pres; l++)
431 {
432 // If the pressure dof is hanging
433 if (pressure_dof_is_hanging[l])
434 {
435 // Pressure dof is hanging so it must be nodal-based
436 // Get the hang info object
437 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
438
439 // Get the number of master nodes from the pressure node
440 n_master = hang_info_pt->nmaster();
441 }
442 // Otherwise the node is its own master
443 else
444 {
445 n_master = 1;
446 }
447
448 // Loop over the master nodes
449 for (unsigned m = 0; m < n_master; m++)
450 {
451 // Get the number of the unknown
452 // If the pressure dof is hanging
453 if (pressure_dof_is_hanging[l])
454 {
455 // Get the local equation from the master node
456 local_eqn =
457 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
458 // Get the weight for the node
459 hang_weight = hang_info_pt->master_weight(m);
460 }
461 else
462 {
463 local_eqn = this->p_local_eqn(l);
464 hang_weight = 1.0;
465 }
466
467 // If the equation is not pinned
468 if (local_eqn >= 0)
469 {
470 residuals[local_eqn] -= source * testp[l] * W * hang_weight;
471 for (unsigned k = 0; k < DIM; k++)
472 {
473 residuals[local_eqn] +=
474 interpolated_dpdx[k] *
475 (scaled_re * interpolated_u[k] * testp[l] + dtestp(l, k)) * W *
476 hang_weight;
477 }
478
479 // Jacobian too?
480 if (flag)
481 {
482 // Number of master nodes and weights
483 unsigned n_master2 = 1;
484 double hang_weight2 = 1.0;
485
486 // Loop over the pressure shape functions
487 for (unsigned l2 = 0; l2 < n_pres; l2++)
488 {
489 // If the pressure dof is hanging
490 if (pressure_dof_is_hanging[l2])
491 {
492 hang_info2_pt =
493 this->pressure_node_pt(l2)->hanging_pt(p_index);
494 // Pressure dof is hanging so it must be nodal-based
495 // Get the number of master nodes from the pressure node
496 n_master2 = hang_info2_pt->nmaster();
497 }
498 // Otherwise the node is its own master
499 else
500 {
501 n_master2 = 1;
502 }
503
504 // Loop over the master nodes
505 for (unsigned m2 = 0; m2 < n_master2; m2++)
506 {
507 // Get the number of the unknown
508 // If the pressure dof is hanging
509 if (pressure_dof_is_hanging[l2])
510 {
511 // Get the unknown from the master node
512 local_unknown = this->local_hang_eqn(
513 hang_info2_pt->master_node_pt(m2), p_index);
514 // Get the weight from the hanging object
515 hang_weight2 = hang_info2_pt->master_weight(m2);
516 }
517 else
518 {
519 local_unknown = this->p_local_eqn(l2);
520 hang_weight2 = 1.0;
521 }
522
523 // If the unknown is not pinned
524 if (local_unknown >= 0)
525 {
526 if ((int(eqn_number(local_eqn)) !=
527 this->Pinned_fp_pressure_eqn) &&
528 (int(eqn_number(local_unknown)) !=
529 this->Pinned_fp_pressure_eqn))
530 {
531 for (unsigned k = 0; k < DIM; k++)
532 {
533 jacobian(local_eqn, local_unknown) +=
534 dtestp(l2, k) *
535 (scaled_re * interpolated_u[k] * testp[l] +
536 dtestp(l, k)) *
537 W * hang_weight * hang_weight2;
538 }
539 }
540 else
541 {
542 if ((int(eqn_number(local_eqn)) ==
543 this->Pinned_fp_pressure_eqn) &&
544 (int(eqn_number(local_unknown)) ==
545 this->Pinned_fp_pressure_eqn))
546 {
547 jacobian(local_eqn, local_unknown) = 1.0;
548 }
549 }
550 }
551 }
552 }
553 } /*End of Jacobian calculation*/
554 } // End of if not boundary condition
555 } // End of loop over master nodes
556 } // End of loop over l
557 } // end of integration loop
558
559 // Now add boundary contributions from Robin BCs
560 unsigned nrobin =
561 this->Pressure_advection_diffusion_robin_element_pt.size();
562 for (unsigned e = 0; e < nrobin; e++)
563 {
564 this->Pressure_advection_diffusion_robin_element_pt[e]
565 ->fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(
566 residuals, jacobian, flag);
567 }
568 }
569
570
571 //========================================================================
572 /// Add element's contribution to the elemental
573 /// residual vector and/or Jacobian matrix.
574 /// flag=1: compute both
575 /// flag=0: compute only residual vector
576 //========================================================================
577 template<unsigned DIM>
580 DenseMatrix<double>& jacobian,
581 DenseMatrix<double>& mass_matrix,
582 unsigned flag)
583 {
584 // Find out how many nodes there are
585 unsigned n_node = nnode();
586
587 // Get continuous time from timestepper of first node
588 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
589
590 // Find out how many pressure dofs there are
591 unsigned n_pres = this->npres_nst();
592
593 // Get the indices at which the velocity components are stored
594 unsigned u_nodal_index[DIM];
595 for (unsigned i = 0; i < DIM; i++)
596 {
597 u_nodal_index[i] = this->u_index_nst(i);
598 }
599
600 // Which nodal value represents the pressure? (Negative if pressure
601 // is not based on nodal interpolation).
602 int p_index = this->p_nodal_index_nst();
603
604 // Local array of booleans that are true if the l-th pressure value is
605 // hanging (avoid repeated virtual function calls)
606 bool pressure_dof_is_hanging[n_pres];
607 // If the pressure is stored at a node
608 if (p_index >= 0)
609 {
610 // Read out whether the pressure is hanging
611 for (unsigned l = 0; l < n_pres; ++l)
612 {
613 pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
614 }
615 }
616 // Otherwise the pressure is not stored at a node and so cannot hang
617 else
618 {
619 for (unsigned l = 0; l < n_pres; ++l)
620 {
621 pressure_dof_is_hanging[l] = false;
622 }
623 }
624
625 // Set up memory for the shape and test functions
626 Shape psif(n_node), testf(n_node);
627 DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
628
629
630 // Set up memory for pressure shape and test functions
631 Shape psip(n_pres), testp(n_pres);
632
633 // Set the value of n_intpt
634 unsigned n_intpt = integral_pt()->nweight();
635
636 // Set the Vector to hold local coordinates
637 Vector<double> s(DIM);
638
639 // Get Physical Variables from Element
640 // Reynolds number must be multiplied by the density ratio
641 double scaled_re = this->re() * this->density_ratio();
642 double scaled_re_st = this->re_st() * this->density_ratio();
643 double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
644 double visc_ratio = this->viscosity_ratio();
645 Vector<double> G = this->g();
646
647 // Integers that store the local equations and unknowns
648 int local_eqn = 0, local_unknown = 0;
649
650 // Pointers to hang info objects
651 HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
652
653 // Local boolean for ALE (or not)
654 bool ALE_is_disabled_flag = this->ALE_is_disabled;
655
656 // Loop over the integration points
657 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
658 {
659 // Assign values of s
660 for (unsigned i = 0; i < DIM; i++)
661 {
662 s[i] = integral_pt()->knot(ipt, i);
663 }
664
665 // Get the integral weight
666 double w = integral_pt()->weight(ipt);
667
668 // Call the derivatives of the shape and test functions
669 double J = this->dshape_and_dtest_eulerian_at_knot_nst(
670 ipt, psif, dpsifdx, testf, dtestfdx);
671
672 // Call the pressure shape and test functions
673 this->pshape_nst(s, psip, testp);
674
675 // Premultiply the weights and the Jacobian
676 double W = w * J;
677
678 // Calculate local values of the pressure and velocity components
679 //--------------------------------------------------------------
680 double interpolated_p = 0.0;
681 Vector<double> interpolated_u(DIM, 0.0);
682 Vector<double> interpolated_x(DIM, 0.0);
683 Vector<double> mesh_veloc(DIM, 0.0);
684 Vector<double> dudt(DIM, 0.0);
685 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
686
687 // Calculate pressure
688 for (unsigned l = 0; l < n_pres; l++)
689 {
690 interpolated_p += this->p_nst(l) * psip[l];
691 }
692
693
694 // Calculate velocities and derivatives
695
696 // Loop over nodes
697 for (unsigned l = 0; l < n_node; l++)
698 {
699 // Loop over directions
700 for (unsigned i = 0; i < DIM; i++)
701 {
702 // Get the nodal value
703 double u_value = this->nodal_value(l, u_nodal_index[i]);
704 interpolated_u[i] += u_value * psif[l];
705 interpolated_x[i] += this->nodal_position(l, i) * psif[l];
706 dudt[i] += this->du_dt_nst(l, i) * psif[l];
707
708 // Loop over derivative directions for velocity gradients
709 for (unsigned j = 0; j < DIM; j++)
710 {
711 interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
712 }
713 }
714 }
715
716 if (!ALE_is_disabled_flag)
717 {
718 // Loop over nodes
719 for (unsigned l = 0; l < n_node; l++)
720 {
721 // Loop over directions
722 for (unsigned i = 0; i < DIM; i++)
723 {
724 mesh_veloc[i] += this->dnodal_position_dt(l, i) * psif[l];
725 }
726 }
727 }
728
729 // Get the user-defined body force terms
730 Vector<double> body_force(DIM);
731 this->get_body_force_nst(time, ipt, s, interpolated_x, body_force);
732
733 // Get the user-defined source function
734 double source = this->get_source_nst(time, ipt, interpolated_x);
735
736 // MOMENTUM EQUATIONS
737 //==================
738
739 // Number of master nodes and storage for the weight of the shape function
740 unsigned n_master = 1;
741 double hang_weight = 1.0;
742
743 // Loop over the nodes for the test functions/equations
744 //----------------------------------------------------
745 for (unsigned l = 0; l < n_node; l++)
746 {
747 // Local boolean to indicate whether the node is hanging
748 bool is_node_hanging = node_pt(l)->is_hanging();
749
750 // If the node is hanging
751 if (is_node_hanging)
752 {
753 hang_info_pt = node_pt(l)->hanging_pt();
754
755 // Read out number of master nodes from hanging data
756 n_master = hang_info_pt->nmaster();
757 }
758 // Otherwise the node is its own master
759 else
760 {
761 n_master = 1;
762 }
763
764 // Loop over the master nodes
765 for (unsigned m = 0; m < n_master; m++)
766 {
767 // Loop over velocity components for equations
768 for (unsigned i = 0; i < DIM; i++)
769 {
770 // Get the equation number
771 // If the node is hanging
772 if (is_node_hanging)
773 {
774 // Get the equation number from the master node
775 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
776 u_nodal_index[i]);
777 // Get the hang weight from the master node
778 hang_weight = hang_info_pt->master_weight(m);
779 }
780 // If the node is not hanging
781 else
782 {
783 // Local equation number
784 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
785
786 // Node contributes with full weight
787 hang_weight = 1.0;
788 }
789
790 // If it's not a boundary condition...
791 if (local_eqn >= 0)
792 {
793 // Temporary variable to hold the residuals
794 double sum = 0.0;
795
796 // Add the user-defined body force terms
797 sum += body_force[i];
798
799 // Add the gravitational body force term
800 sum += scaled_re_inv_fr * G[i];
801
802 // Add in the inertial term
803 sum -= scaled_re_st * dudt[i];
804
805 // Convective terms, including mesh velocity
806 for (unsigned k = 0; k < DIM; k++)
807 {
808 double tmp = scaled_re * interpolated_u[k];
809 if (!ALE_is_disabled_flag)
810 {
811 tmp -= scaled_re_st * mesh_veloc[k];
812 }
813 sum -= tmp * interpolated_dudx(i, k);
814 }
815
816 // Add the pressure gradient term
817 sum = (sum * testf[l] + interpolated_p * dtestfdx(l, i)) * W *
818 hang_weight;
819
820 // Add in the stress tensor terms
821 // The viscosity ratio needs to go in here to ensure
822 // continuity of normal stress is satisfied even in flows
823 // with zero pressure gradient!
824 for (unsigned k = 0; k < DIM; k++)
825 {
826 sum -= visc_ratio *
827 (interpolated_dudx(i, k) +
828 this->Gamma[i] * interpolated_dudx(k, i)) *
829 dtestfdx(l, k) * W * hang_weight;
830 }
831
832 // Add contribution
833 residuals[local_eqn] += sum;
834
835 // CALCULATE THE JACOBIAN
836 if (flag)
837 {
838 // Number of master nodes and weights
839 unsigned n_master2 = 1;
840 double hang_weight2 = 1.0;
841 // Loop over the velocity nodes for columns
842 for (unsigned l2 = 0; l2 < n_node; l2++)
843 {
844 // Local boolean to indicate whether the node is hanging
845 bool is_node2_hanging = node_pt(l2)->is_hanging();
846
847 // If the node is hanging
848 if (is_node2_hanging)
849 {
850 hang_info2_pt = node_pt(l2)->hanging_pt();
851 // Read out number of master nodes from hanging data
852 n_master2 = hang_info2_pt->nmaster();
853 }
854 // Otherwise the node is its own master
855 else
856 {
857 n_master2 = 1;
858 }
859
860 // Loop over the master nodes
861 for (unsigned m2 = 0; m2 < n_master2; m2++)
862 {
863 // Loop over the velocity components
864 for (unsigned i2 = 0; i2 < DIM; i2++)
865 {
866 // Get the number of the unknown
867 // If the node is hanging
868 if (is_node2_hanging)
869 {
870 // Get the equation number from the master node
871 local_unknown = this->local_hang_eqn(
872 hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
873 // Get the hang weights
874 hang_weight2 = hang_info2_pt->master_weight(m2);
875 }
876 else
877 {
878 local_unknown =
879 this->nodal_local_eqn(l2, u_nodal_index[i2]);
880 hang_weight2 = 1.0;
881 }
882
883 // If the unknown is non-zero
884 if (local_unknown >= 0)
885 {
886 // Add contribution to Elemental Matrix
887 jacobian(local_eqn, local_unknown) -=
888 visc_ratio * this->Gamma[i] * dpsifdx(l2, i) *
889 dtestfdx(l, i2) * W * hang_weight * hang_weight2;
890
891 // Now add in the inertial terms
892 jacobian(local_eqn, local_unknown) -=
893 scaled_re * psif[l2] * interpolated_dudx(i, i2) *
894 testf[l] * W * hang_weight * hang_weight2;
895
896 // Extra diagonal components if i2=i
897 if (i2 == i)
898 {
899 // Mass matrix entries
900 // Again note the positive sign because the mass
901 // matrix is taken on the other side of the equation
902 if (flag == 2)
903 {
904 mass_matrix(local_eqn, local_unknown) +=
905 scaled_re_st * psif[l2] * testf[l] * W *
906 hang_weight * hang_weight2;
907 }
908
909 // du/dt term
910 jacobian(local_eqn, local_unknown) -=
911 scaled_re_st *
912 node_pt(l2)->time_stepper_pt()->weight(1, 0) *
913 psif[l2] * testf[l] * W * hang_weight *
914 hang_weight2;
915
916 // Extra advective terms
917 for (unsigned k = 0; k < DIM; k++)
918 {
919 double tmp = scaled_re * interpolated_u[k];
920 if (!ALE_is_disabled_flag)
921 {
922 tmp -= scaled_re_st * mesh_veloc[k];
923 }
924
925 jacobian(local_eqn, local_unknown) -=
926 tmp * dpsifdx(l2, k) * testf[l] * W *
927 hang_weight * hang_weight2;
928 }
929
930 // Extra viscous terms
931 for (unsigned k = 0; k < DIM; k++)
932 {
933 jacobian(local_eqn, local_unknown) -=
934 visc_ratio * dpsifdx(l2, k) * dtestfdx(l, k) * W *
935 hang_weight * hang_weight2;
936 }
937 }
938 }
939 }
940 }
941 }
942
943 // Loop over the pressure shape functions
944 for (unsigned l2 = 0; l2 < n_pres; l2++)
945 {
946 // If the pressure dof is hanging
947 if (pressure_dof_is_hanging[l2])
948 {
949 hang_info2_pt =
950 this->pressure_node_pt(l2)->hanging_pt(p_index);
951 // Pressure dof is hanging so it must be nodal-based
952 // Get the number of master nodes from the pressure node
953 n_master2 = hang_info2_pt->nmaster();
954 }
955 // Otherwise the node is its own master
956 else
957 {
958 n_master2 = 1;
959 }
960
961 // Loop over the master nodes
962 for (unsigned m2 = 0; m2 < n_master2; m2++)
963 {
964 // Get the number of the unknown
965 // If the pressure dof is hanging
966 if (pressure_dof_is_hanging[l2])
967 {
968 // Get the unknown from the master node
969 local_unknown = this->local_hang_eqn(
970 hang_info2_pt->master_node_pt(m2), p_index);
971 // Get the weight from the hanging object
972 hang_weight2 = hang_info2_pt->master_weight(m2);
973 }
974 else
975 {
976 local_unknown = this->p_local_eqn(l2);
977 hang_weight2 = 1.0;
978 }
979
980 // If the unknown is not pinned
981 if (local_unknown >= 0)
982 {
983 jacobian(local_eqn, local_unknown) +=
984 psip[l2] * dtestfdx(l, i) * W * hang_weight *
985 hang_weight2;
986 }
987 }
988 }
989
990 } // End of Jacobian calculation
991
992 } // End of if not boundary condition statement
993
994 } // End of loop over components of non-hanging node
995
996 } // End of loop over master nodes
997
998 } // End of loop over nodes for equations
999
1000
1001 // CONTINUITY EQUATION
1002 //===================
1003
1004 // Loop over the pressure shape functions
1005 for (unsigned l = 0; l < n_pres; l++)
1006 {
1007 // If the pressure dof is hanging
1008 if (pressure_dof_is_hanging[l])
1009 {
1010 // Pressure dof is hanging so it must be nodal-based
1011 // Get the hang info object
1012 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1013
1014 // Get the number of master nodes from the pressure node
1015 n_master = hang_info_pt->nmaster();
1016 }
1017 // Otherwise the node is its own master
1018 else
1019 {
1020 n_master = 1;
1021 }
1022
1023 // Loop over the master nodes
1024 for (unsigned m = 0; m < n_master; m++)
1025 {
1026 // Get the number of the unknown
1027 // If the pressure dof is hanging
1028 if (pressure_dof_is_hanging[l])
1029 {
1030 // Get the local equation from the master node
1031 local_eqn =
1032 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
1033 // Get the weight for the node
1034 hang_weight = hang_info_pt->master_weight(m);
1035 }
1036 else
1037 {
1038 local_eqn = this->p_local_eqn(l);
1039 hang_weight = 1.0;
1040 }
1041
1042 // If the equation is not pinned
1043 if (local_eqn >= 0)
1044 {
1045 // Source term
1046 residuals[local_eqn] -= source * testp[l] * W * hang_weight;
1047
1048 // Loop over velocity components
1049 for (unsigned k = 0; k < DIM; k++)
1050 {
1051 residuals[local_eqn] +=
1052 interpolated_dudx(k, k) * testp[l] * W * hang_weight;
1053 }
1054
1055 // CALCULATE THE JACOBIAN
1056 if (flag)
1057 {
1058 unsigned n_master2 = 1;
1059 double hang_weight2 = 1.0;
1060 // Loop over the velocity nodes for columns
1061 for (unsigned l2 = 0; l2 < n_node; l2++)
1062 {
1063 // Local boolean to indicate whether the node is hanging
1064 bool is_node2_hanging = node_pt(l2)->is_hanging();
1065
1066 // If the node is hanging
1067 if (is_node2_hanging)
1068 {
1069 hang_info2_pt = node_pt(l2)->hanging_pt();
1070 // Read out number of master nodes from hanging data
1071 n_master2 = hang_info2_pt->nmaster();
1072 }
1073 // Otherwise the node is its own master
1074 else
1075 {
1076 n_master2 = 1;
1077 }
1078
1079 // Loop over the master nodes
1080 for (unsigned m2 = 0; m2 < n_master2; m2++)
1081 {
1082 // Loop over the velocity components
1083 for (unsigned i2 = 0; i2 < DIM; i2++)
1084 {
1085 // Get the number of the unknown
1086 // If the node is hanging
1087 if (is_node2_hanging)
1088 {
1089 // Get the equation number from the master node
1090 local_unknown = this->local_hang_eqn(
1091 hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
1092 hang_weight2 = hang_info2_pt->master_weight(m2);
1093 }
1094 else
1095 {
1096 local_unknown =
1097 this->nodal_local_eqn(l2, u_nodal_index[i2]);
1098 hang_weight2 = 1.0;
1099 }
1100
1101 // If the unknown is not pinned
1102 if (local_unknown >= 0)
1103 {
1104 jacobian(local_eqn, local_unknown) +=
1105 dpsifdx(l2, i2) * testp[l] * W * hang_weight *
1106 hang_weight2;
1107 }
1108 }
1109 }
1110 }
1111
1112 // NO PRESSURE CONTRIBUTION TO THE JACOBIAN
1113
1114 } // End of jacobian calculation
1115 }
1116 }
1117 } // End of loop over pressure variables
1118
1119 } // End of loop over integration points
1120 }
1121
1122
1123 //======================================================================
1124 /// Compute derivatives of elemental residual vector with respect
1125 /// to nodal coordinates.
1126 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
1127 /// Overloads the FD-based version in the FE base class.
1128 //======================================================================
1129 template<unsigned DIM>
1131 RankThreeTensor<double>& dresidual_dnodal_coordinates)
1132 {
1133 // Return immediately if there are no dofs
1134 if (ndof() == 0)
1135 {
1136 return;
1137 }
1138
1139 // Determine number of nodes in element
1140 const unsigned n_node = nnode();
1141
1142 // Get continuous time from timestepper of first node
1143 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
1144
1145 // Determine number of pressure dofs in element
1146 const unsigned n_pres = this->npres_nst();
1147
1148 // Find the indices at which the local velocities are stored
1149 unsigned u_nodal_index[DIM];
1150 for (unsigned i = 0; i < DIM; i++)
1151 {
1152 u_nodal_index[i] = this->u_index_nst(i);
1153 }
1154
1155 // Which nodal value represents the pressure? (Negative if pressure
1156 // is not based on nodal interpolation).
1157 const int p_index = this->p_nodal_index_nst();
1158
1159 // Local array of booleans that are true if the l-th pressure value is
1160 // hanging (avoid repeated virtual function calls)
1161 bool pressure_dof_is_hanging[n_pres];
1162
1163 // If the pressure is stored at a node
1164 if (p_index >= 0)
1165 {
1166 // Read out whether the pressure is hanging
1167 for (unsigned l = 0; l < n_pres; ++l)
1168 {
1169 pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
1170 }
1171 }
1172 // Otherwise the pressure is not stored at a node and so cannot hang
1173 else
1174 {
1175 for (unsigned l = 0; l < n_pres; ++l)
1176 {
1177 pressure_dof_is_hanging[l] = false;
1178 }
1179 }
1180
1181 // Set up memory for the shape and test functions
1182 Shape psif(n_node), testf(n_node);
1183 DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
1184
1185 // Set up memory for pressure shape and test functions
1186 Shape psip(n_pres), testp(n_pres);
1187
1188 // Determine number of shape controlling nodes
1189 const unsigned n_shape_controlling_node = nshape_controlling_nodes();
1190
1191 // Deriatives of shape fct derivatives w.r.t. nodal coords
1192 RankFourTensor<double> d_dpsifdx_dX(
1193 DIM, n_shape_controlling_node, n_node, DIM);
1194 RankFourTensor<double> d_dtestfdx_dX(
1195 DIM, n_shape_controlling_node, n_node, DIM);
1196
1197 // Derivative of Jacobian of mapping w.r.t. to nodal coords
1198 DenseMatrix<double> dJ_dX(DIM, n_shape_controlling_node);
1199
1200 // Derivatives of derivative of u w.r.t. nodal coords
1201 RankFourTensor<double> d_dudx_dX(DIM, n_shape_controlling_node, DIM, DIM);
1202
1203 // Derivatives of nodal velocities w.r.t. nodal coords:
1204 // Assumption: Interaction only local via no-slip so
1205 // X_ij only affects U_ij.
1206 DenseMatrix<double> d_U_dX(DIM, n_shape_controlling_node, 0.0);
1207
1208 // Determine the number of integration points
1209 const unsigned n_intpt = integral_pt()->nweight();
1210
1211 // Vector to hold local coordinates
1212 Vector<double> s(DIM);
1213
1214 // Get physical variables from element
1215 // (Reynolds number must be multiplied by the density ratio)
1216 double scaled_re = this->re() * this->density_ratio();
1217 double scaled_re_st = this->re_st() * this->density_ratio();
1218 double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
1219 double visc_ratio = this->viscosity_ratio();
1220 Vector<double> G = this->g();
1221
1222 // FD step
1224
1225 // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1226 // Assumption: Interaction only local via no-slip so
1227 // X_ij only affects U_ij.
1228 bool element_has_node_with_aux_node_update_fct = false;
1229
1230 std::map<Node*, unsigned> local_shape_controlling_node_lookup =
1231 shape_controlling_node_lookup();
1232
1233 // FD loop over shape-controlling nodes
1234 for (std::map<Node*, unsigned>::iterator it =
1235 local_shape_controlling_node_lookup.begin();
1236 it != local_shape_controlling_node_lookup.end();
1237 it++)
1238 {
1239 // Get node
1240 Node* nod_pt = it->first;
1241
1242 // Get its number
1243 unsigned q = it->second;
1244
1245 // Only compute if there's a node-update fct involved
1247 {
1248 element_has_node_with_aux_node_update_fct = true;
1249
1250 // Current nodal velocity
1251 Vector<double> u_ref(DIM);
1252 for (unsigned i = 0; i < DIM; i++)
1253 {
1254 u_ref[i] = *(nod_pt->value_pt(u_nodal_index[i]));
1255 }
1256
1257 // FD
1258 for (unsigned p = 0; p < DIM; p++)
1259 {
1260 // Make backup
1261 double backup = nod_pt->x(p);
1262
1263 // Do FD step. No node update required as we're
1264 // attacking the coordinate directly...
1265 nod_pt->x(p) += eps_fd;
1266
1267 // Do auxiliary node update (to apply no slip)
1269
1270 // Evaluate
1271 d_U_dX(p, q) =
1272 (*(nod_pt->value_pt(u_nodal_index[p])) - u_ref[p]) / eps_fd;
1273
1274 // Reset
1275 nod_pt->x(p) = backup;
1276
1277 // Do auxiliary node update (to apply no slip)
1279 }
1280 }
1281 }
1282
1283 // Integer to store the local equation number
1284 int local_eqn = 0;
1285
1286 // Pointers to hang info object
1287 HangInfo* hang_info_pt = 0;
1288
1289 // Loop over the integration points
1290 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1291 {
1292 // Assign values of s
1293 for (unsigned i = 0; i < DIM; i++)
1294 {
1295 s[i] = integral_pt()->knot(ipt, i);
1296 }
1297
1298 // Get the integral weight
1299 const double w = integral_pt()->weight(ipt);
1300
1301 // Call the derivatives of the shape and test functions
1302 const double J =
1303 this->dshape_and_dtest_eulerian_at_knot_nst(ipt,
1304 psif,
1305 dpsifdx,
1306 d_dpsifdx_dX,
1307 testf,
1308 dtestfdx,
1309 d_dtestfdx_dX,
1310 dJ_dX);
1311
1312 // Call the pressure shape and test functions
1313 this->pshape_nst(s, psip, testp);
1314
1315 // Calculate local values of the pressure and velocity components
1316 // Allocate
1317 double interpolated_p = 0.0;
1318 Vector<double> interpolated_u(DIM, 0.0);
1319 Vector<double> interpolated_x(DIM, 0.0);
1320 Vector<double> mesh_velocity(DIM, 0.0);
1321 Vector<double> dudt(DIM, 0.0);
1322 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1323
1324 // Calculate pressure
1325 for (unsigned l = 0; l < n_pres; l++)
1326 {
1327 interpolated_p += this->p_nst(l) * psip[l];
1328 }
1329
1330 // Calculate velocities and derivatives:
1331
1332 // Loop over nodes
1333 for (unsigned l = 0; l < n_node; l++)
1334 {
1335 // Loop over directions
1336 for (unsigned i = 0; i < DIM; i++)
1337 {
1338 // Get the nodal value
1339 const double u_value = nodal_value(l, u_nodal_index[i]);
1340 interpolated_u[i] += u_value * psif[l];
1341 interpolated_x[i] += nodal_position(l, i) * psif[l];
1342 dudt[i] += this->du_dt_nst(l, i) * psif[l];
1343
1344 // Loop over derivative directions
1345 for (unsigned j = 0; j < DIM; j++)
1346 {
1347 interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1348 }
1349 }
1350 }
1351
1352 if (!this->ALE_is_disabled)
1353 {
1354 // Loop over nodes
1355 for (unsigned l = 0; l < n_node; l++)
1356 {
1357 // Loop over directions
1358 for (unsigned i = 0; i < DIM; i++)
1359 {
1360 mesh_velocity[i] += this->dnodal_position_dt(l, i) * psif[l];
1361 }
1362 }
1363 }
1364
1365 // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1366
1367 // Loop over shape-controlling nodes
1368 for (unsigned q = 0; q < n_shape_controlling_node; q++)
1369 {
1370 // Loop over coordinate directions
1371 for (unsigned p = 0; p < DIM; p++)
1372 {
1373 for (unsigned i = 0; i < DIM; i++)
1374 {
1375 for (unsigned k = 0; k < DIM; k++)
1376 {
1377 double aux = 0.0;
1378 for (unsigned j = 0; j < n_node; j++)
1379 {
1380 aux +=
1381 nodal_value(j, u_nodal_index[i]) * d_dpsifdx_dX(p, q, j, k);
1382 }
1383 d_dudx_dX(p, q, i, k) = aux;
1384 }
1385 }
1386 }
1387 }
1388
1389 // Get weight of actual nodal position/value in computation of mesh
1390 // velocity from positional/value time stepper
1391 const double pos_time_weight =
1392 node_pt(0)->position_time_stepper_pt()->weight(1, 0);
1393 const double val_time_weight =
1394 node_pt(0)->time_stepper_pt()->weight(1, 0);
1395
1396 // Get the user-defined body force terms
1397 Vector<double> body_force(DIM);
1398 this->get_body_force_nst(time, ipt, s, interpolated_x, body_force);
1399
1400 // Get the user-defined source function
1401 const double source = this->get_source_nst(time, ipt, interpolated_x);
1402
1403 // Get gradient of body force function
1404 DenseMatrix<double> d_body_force_dx(DIM, DIM, 0.0);
1405 this->get_body_force_gradient_nst(
1406 time, ipt, s, interpolated_x, d_body_force_dx);
1407
1408 // Get gradient of source function
1409 Vector<double> source_gradient(DIM, 0.0);
1410 this->get_source_gradient_nst(time, ipt, interpolated_x, source_gradient);
1411
1412
1413 // Assemble shape derivatives
1414 //---------------------------
1415
1416 // MOMENTUM EQUATIONS
1417 // ------------------
1418
1419 // Number of master nodes and storage for the weight of the shape function
1420 unsigned n_master = 1;
1421 double hang_weight = 1.0;
1422
1423 // Loop over the test functions
1424 for (unsigned l = 0; l < n_node; l++)
1425 {
1426 // Local boolean to indicate whether the node is hanging
1427 bool is_node_hanging = node_pt(l)->is_hanging();
1428
1429 // If the node is hanging
1430 if (is_node_hanging)
1431 {
1432 hang_info_pt = node_pt(l)->hanging_pt();
1433
1434 // Read out number of master nodes from hanging data
1435 n_master = hang_info_pt->nmaster();
1436 }
1437 // Otherwise the node is its own master
1438 else
1439 {
1440 n_master = 1;
1441 }
1442
1443 // Loop over the master nodes
1444 for (unsigned m = 0; m < n_master; m++)
1445 {
1446 // Loop over coordinate directions
1447 for (unsigned i = 0; i < DIM; i++)
1448 {
1449 // Get the equation number
1450 // If the node is hanging
1451 if (is_node_hanging)
1452 {
1453 // Get the equation number from the master node
1454 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1455 u_nodal_index[i]);
1456 // Get the hang weight from the master node
1457 hang_weight = hang_info_pt->master_weight(m);
1458 }
1459 // If the node is not hanging
1460 else
1461 {
1462 // Local equation number
1463 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
1464
1465 // Node contributes with full weight
1466 hang_weight = 1.0;
1467 }
1468
1469 // IF it's not a boundary condition
1470 if (local_eqn >= 0)
1471 {
1472 // Loop over coordinate directions
1473 for (unsigned p = 0; p < DIM; p++)
1474 {
1475 // Loop over shape controlling nodes
1476 for (unsigned q = 0; q < n_shape_controlling_node; q++)
1477 {
1478 // Residual x deriv of Jacobian
1479 // ----------------------------
1480
1481 // Add the user-defined body force terms
1482 double sum = body_force[i] * testf[l];
1483
1484 // Add the gravitational body force term
1485 sum += scaled_re_inv_fr * testf[l] * G[i];
1486
1487 // Add the pressure gradient term
1488 sum += interpolated_p * dtestfdx(l, i);
1489
1490 // Add in the stress tensor terms
1491 // The viscosity ratio needs to go in here to ensure
1492 // continuity of normal stress is satisfied even in flows
1493 // with zero pressure gradient!
1494 for (unsigned k = 0; k < DIM; k++)
1495 {
1496 sum -= visc_ratio *
1497 (interpolated_dudx(i, k) +
1498 this->Gamma[i] * interpolated_dudx(k, i)) *
1499 dtestfdx(l, k);
1500 }
1501
1502 // Add in the inertial terms
1503
1504 // du/dt term
1505 sum -= scaled_re_st * dudt[i] * testf[l];
1506
1507 // Convective terms, including mesh velocity
1508 for (unsigned k = 0; k < DIM; k++)
1509 {
1510 double tmp = scaled_re * interpolated_u[k];
1511 if (!this->ALE_is_disabled)
1512 {
1513 tmp -= scaled_re_st * mesh_velocity[k];
1514 }
1515 sum -= tmp * interpolated_dudx(i, k) * testf[l];
1516 }
1517
1518 // Multiply throsugh by deriv of Jacobian and integration
1519 // weight
1520 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1521 sum * dJ_dX(p, q) * w * hang_weight;
1522
1523 // Derivative of residual x Jacobian
1524 // ---------------------------------
1525
1526 // Body force
1527 sum = d_body_force_dx(i, p) * psif(q) * testf(l);
1528
1529 // Pressure gradient term
1530 sum += interpolated_p * d_dtestfdx_dX(p, q, l, i);
1531
1532 // Viscous term
1533 for (unsigned k = 0; k < DIM; k++)
1534 {
1535 sum -=
1536 visc_ratio * ((interpolated_dudx(i, k) +
1537 this->Gamma[i] * interpolated_dudx(k, i)) *
1538 d_dtestfdx_dX(p, q, l, k) +
1539 (d_dudx_dX(p, q, i, k) +
1540 this->Gamma[i] * d_dudx_dX(p, q, k, i)) *
1541 dtestfdx(l, k));
1542 }
1543
1544 // Convective terms, including mesh velocity
1545 for (unsigned k = 0; k < DIM; k++)
1546 {
1547 double tmp = scaled_re * interpolated_u[k];
1548 if (!this->ALE_is_disabled)
1549 {
1550 tmp -= scaled_re_st * mesh_velocity[k];
1551 }
1552 sum -= tmp * d_dudx_dX(p, q, i, k) * testf(l);
1553 }
1554 if (!this->ALE_is_disabled)
1555 {
1556 sum += scaled_re_st * pos_time_weight * psif(q) *
1557 interpolated_dudx(i, p) * testf(l);
1558 }
1559
1560 // Multiply through by Jacobian and integration weight
1561 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1562 sum * J * w * hang_weight;
1563
1564 } // End of loop over shape controlling nodes q
1565 } // End of loop over coordinate directions p
1566
1567
1568 // Derivs w.r.t. to nodal velocities
1569 // ---------------------------------
1570 if (element_has_node_with_aux_node_update_fct)
1571 {
1572 // Loop over local nodes
1573 for (unsigned q_local = 0; q_local < n_node; q_local++)
1574 {
1575 // Number of master nodes and storage for the weight of
1576 // the shape function
1577 unsigned n_master2 = 1;
1578 double hang_weight2 = 1.0;
1579 HangInfo* hang_info2_pt = 0;
1580
1581 // Local boolean to indicate whether the node is hanging
1582 bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1583
1584 Node* actual_shape_controlling_node_pt = node_pt(q_local);
1585
1586 // If the node is hanging
1587 if (is_node_hanging2)
1588 {
1589 hang_info2_pt = node_pt(q_local)->hanging_pt();
1590
1591 // Read out number of master nodes from hanging data
1592 n_master2 = hang_info2_pt->nmaster();
1593 }
1594 // Otherwise the node is its own master
1595 else
1596 {
1597 n_master2 = 1;
1598 }
1599
1600 // Loop over the master nodes
1601 for (unsigned mm = 0; mm < n_master2; mm++)
1602 {
1603 if (is_node_hanging2)
1604 {
1605 actual_shape_controlling_node_pt =
1606 hang_info2_pt->master_node_pt(mm);
1607 hang_weight2 = hang_info2_pt->master_weight(mm);
1608 }
1609
1610 // Find the corresponding number
1611 unsigned q = local_shape_controlling_node_lookup
1612 [actual_shape_controlling_node_pt];
1613
1614 // Loop over coordinate directions
1615 for (unsigned p = 0; p < DIM; p++)
1616 {
1617 double sum = -visc_ratio * this->Gamma[i] *
1618 dpsifdx(q_local, i) * dtestfdx(l, p) -
1619 scaled_re * psif(q_local) *
1620 interpolated_dudx(i, p) * testf(l);
1621 if (i == p)
1622 {
1623 sum -= scaled_re_st * val_time_weight * psif(q_local) *
1624 testf(l);
1625 for (unsigned k = 0; k < DIM; k++)
1626 {
1627 sum -=
1628 visc_ratio * dpsifdx(q_local, k) * dtestfdx(l, k);
1629 double tmp = scaled_re * interpolated_u[k];
1630 if (!this->ALE_is_disabled)
1631 {
1632 tmp -= scaled_re_st * mesh_velocity[k];
1633 }
1634 sum -= tmp * dpsifdx(q_local, k) * testf(l);
1635 }
1636 }
1637
1638 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1639 sum * d_U_dX(p, q) * J * w * hang_weight * hang_weight2;
1640 }
1641 } // End of loop over master nodes
1642 } // End of loop over local nodes
1643 } // End of if(element_has_node_with_aux_node_update_fct)
1644
1645
1646 } // local_eqn>=0
1647 }
1648 }
1649 } // End of loop over test functions
1650
1651
1652 // CONTINUITY EQUATION
1653 // -------------------
1654
1655 // Loop over the shape functions
1656 for (unsigned l = 0; l < n_pres; l++)
1657 {
1658 // If the pressure dof is hanging
1659 if (pressure_dof_is_hanging[l])
1660 {
1661 // Pressure dof is hanging so it must be nodal-based
1662 // Get the hang info object
1663 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1664
1665 // Get the number of master nodes from the pressure node
1666 n_master = hang_info_pt->nmaster();
1667 }
1668 // Otherwise the node is its own master
1669 else
1670 {
1671 n_master = 1;
1672 }
1673
1674 // Loop over the master nodes
1675 for (unsigned m = 0; m < n_master; m++)
1676 {
1677 // Get the number of the unknown
1678 // If the pressure dof is hanging
1679 if (pressure_dof_is_hanging[l])
1680 {
1681 // Get the local equation from the master node
1682 local_eqn =
1683 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
1684 // Get the weight for the node
1685 hang_weight = hang_info_pt->master_weight(m);
1686 }
1687 else
1688 {
1689 local_eqn = this->p_local_eqn(l);
1690 hang_weight = 1.0;
1691 }
1692
1693 // If not a boundary conditions
1694 if (local_eqn >= 0)
1695 {
1696 // Loop over coordinate directions
1697 for (unsigned p = 0; p < DIM; p++)
1698 {
1699 // Loop over nodes
1700 for (unsigned q = 0; q < n_shape_controlling_node; q++)
1701 {
1702 // Residual x deriv of Jacobian
1703 //-----------------------------
1704
1705 // Source term
1706 double aux = -source;
1707
1708 // Loop over velocity components
1709 for (unsigned k = 0; k < DIM; k++)
1710 {
1711 aux += interpolated_dudx(k, k);
1712 }
1713
1714 // Multiply through by deriv of Jacobian and integration weight
1715 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1716 aux * dJ_dX(p, q) * testp[l] * w * hang_weight;
1717
1718
1719 // Derivative of residual x Jacobian
1720 // ---------------------------------
1721
1722 // Loop over velocity components
1723 aux = -source_gradient[p] * psif(q);
1724 for (unsigned k = 0; k < DIM; k++)
1725 {
1726 aux += d_dudx_dX(p, q, k, k);
1727 }
1728 // Multiply through by Jacobian and integration weight
1729 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1730 aux * testp[l] * J * w * hang_weight;
1731 }
1732 }
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 aux = dpsifdx(q_local, p) * testp(l);
1785 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1786 aux * d_U_dX(p, q) * J * w * hang_weight * hang_weight2;
1787 }
1788 } // End of loop over (mm) master nodes
1789 } // End of loop over local nodes q_local
1790 } // End of if(element_has_node_with_aux_node_update_fct)
1791 } // End of if it's not a boundary condition
1792 } // End of loop over (m) master nodes
1793 } // End of loop over shape functions for continuity eqn
1794
1795 } // End of loop over integration points
1796 }
1797
1798 //======================================================================
1799 /// 2D Further build for Crouzeix_Raviart interpolates the internal
1800 /// pressure dofs from father element: Make sure pressure values and
1801 /// dp/ds agree between fathers and sons at the midpoints of the son
1802 /// elements.
1803 //======================================================================
1804 template<>
1806 {
1807 if (this->tree_pt()->father_pt() != 0)
1808 {
1809 // Call the generic further build (if there is a father)
1811 }
1812 // Now do the PRefineableQElement further_build()
1814
1815 // Resize internal pressure storage
1816 if (this->internal_data_pt(this->P_nst_internal_index)->nvalue() <=
1817 this->npres_nst())
1818 {
1819 this->internal_data_pt(this->P_nst_internal_index)
1820 ->resize(this->npres_nst());
1821 }
1822 else
1823 {
1824 Data* new_data_pt = new Data(this->npres_nst());
1825 delete internal_data_pt(this->P_nst_internal_index);
1826 internal_data_pt(this->P_nst_internal_index) = new_data_pt;
1827 }
1828
1829 if (this->tree_pt()->father_pt() != 0)
1830 {
1831 // Pointer to my father (in C-R element impersonation)
1832 PRefineableQCrouzeixRaviartElement<2>* father_element_pt =
1834 quadtree_pt()->father_pt()->object_pt());
1835
1836 // If element has same p-order as father then do the projection problem
1837 // (called after h-refinement)
1838 if (father_element_pt->p_order() == this->p_order())
1839 {
1840 using namespace QuadTreeNames;
1841
1842 // What type of son am I? Ask my quadtree representation...
1843 int son_type = quadtree_pt()->son_type();
1844
1845 Vector<double> s_father(2);
1846
1847 // Son midpoint is located at the following coordinates in father
1848 // element:
1849 switch (son_type)
1850 {
1851 case SW:
1852 // South west son
1853 s_father[0] = -0.5;
1854 s_father[1] = -0.5;
1855 break;
1856 case SE:
1857 // South east son
1858 s_father[0] = 0.5;
1859 s_father[1] = -0.5;
1860 break;
1861 case NE:
1862 // North east son
1863 s_father[0] = 0.5;
1864 s_father[1] = 0.5;
1865 break;
1866 case NW:
1867 // North west son
1868 s_father[0] = -0.5;
1869 s_father[1] = 0.5;
1870 break;
1871 default:
1872 throw OomphLibError("Invalid son type in",
1873 OOMPH_CURRENT_FUNCTION,
1874 OOMPH_EXCEPTION_LOCATION);
1875 break;
1876 }
1877
1878 // Get pressure value in father element
1879 double press = father_element_pt->interpolated_p_nst(s_father);
1880
1881 // Reset all pressures to zero
1882 for (unsigned p = 0; p < this->npres_nst(); p++)
1883 {
1884 internal_data_pt(this->P_nst_internal_index)->set_value(p, 0.0);
1885 }
1886
1887 // Set pressure values from father (BENFLAG: projection problem hack)
1888 if (this->npres_nst() == 1)
1889 {
1890 internal_data_pt(this->P_nst_internal_index)->set_value(0, press);
1891 }
1892 else if (this->npres_nst() == 3)
1893 {
1894 internal_data_pt(this->P_nst_internal_index)->set_value(0, press);
1895 internal_data_pt(this->P_nst_internal_index)->set_value(1, press);
1896 internal_data_pt(this->P_nst_internal_index)->set_value(2, press);
1897 }
1898 else
1899 {
1900 internal_data_pt(this->P_nst_internal_index)->set_value(0, press);
1901 internal_data_pt(this->P_nst_internal_index)->set_value(1, press);
1902 internal_data_pt(this->P_nst_internal_index)->set_value(2, press);
1903 internal_data_pt(this->P_nst_internal_index)->set_value(3, press);
1904 }
1905 } // Otherwise this is called after p-refinement
1906 }
1907 }
1908
1909 //======================================================================
1910 /// 2D Further build for Crouzeix_Raviart interpolates the internal
1911 /// pressure dofs from father element: Make sure pressure values and
1912 /// dp/ds agree between fathers and sons at the midpoints of the son
1913 /// elements.
1914 //======================================================================
1915 template<>
1917 {
1918 if (this->tree_pt()->father_pt() != 0)
1919 {
1920 // Call the generic further build (if there is a father)
1922 }
1923 // Now do the PRefineableQElement further_build()
1925
1926 // Resize internal pressure storage
1927 if (this->internal_data_pt(this->P_nst_internal_index)->nvalue() <=
1928 this->npres_nst())
1929 {
1930 this->internal_data_pt(this->P_nst_internal_index)
1931 ->resize(this->npres_nst());
1932 }
1933 else
1934 {
1935 Data* new_data_pt = new Data(this->npres_nst());
1936 delete internal_data_pt(this->P_nst_internal_index);
1937 internal_data_pt(this->P_nst_internal_index) = new_data_pt;
1938 }
1939
1940 if (this->tree_pt()->father_pt() != 0)
1941 {
1942 // Pointer to my father (in C-R element impersonation)
1943 PRefineableQCrouzeixRaviartElement<3>* father_element_pt =
1945 octree_pt()->father_pt()->object_pt());
1946
1947 // If element has same p-order as father then do the projection problem
1948 // (called after h-refinement)
1949 if (father_element_pt->p_order() == this->p_order())
1950 {
1951 using namespace OcTreeNames;
1952
1953 // What type of son am I? Ask my quadtree representation...
1954 int son_type = octree_pt()->son_type();
1955
1956 Vector<double> s_father(3);
1957
1958
1959 // Son midpoint is located at the following coordinates in father
1960 // element:
1961 for (unsigned i = 0; i < 3; i++)
1962 {
1963 s_father[i] = 0.5 * OcTree::Direction_to_vector[son_type][i];
1964 }
1965
1966 // Get pressure value in father element
1967 double press = father_element_pt->interpolated_p_nst(s_father);
1968
1969 // Reset all pressures to zero
1970 for (unsigned p = 0; p < this->npres_nst(); p++)
1971 {
1972 internal_data_pt(this->P_nst_internal_index)->set_value(p, 0.0);
1973 }
1974
1975 // Set pressure values from father (BENFLAG: projection problem hack)
1976 if (this->npres_nst() == 1)
1977 {
1978 internal_data_pt(this->P_nst_internal_index)->set_value(0, press);
1979 }
1980 else
1981 {
1982 internal_data_pt(this->P_nst_internal_index)->set_value(0, press);
1983 internal_data_pt(this->P_nst_internal_index)->set_value(1, press);
1984 internal_data_pt(this->P_nst_internal_index)->set_value(2, press);
1985 internal_data_pt(this->P_nst_internal_index)->set_value(3, press);
1986 internal_data_pt(this->P_nst_internal_index)->set_value(4, press);
1987 internal_data_pt(this->P_nst_internal_index)->set_value(5, press);
1988 internal_data_pt(this->P_nst_internal_index)->set_value(6, press);
1989 internal_data_pt(this->P_nst_internal_index)->set_value(7, press);
1990 }
1991 } // Otherwise this is called after p-refinement
1992 }
1993 }
1994
1995
1996 //====================================================================
1997 /// / Force build of templates
1998 //====================================================================
2001 template class RefineableQTaylorHoodElement<2>;
2002 template class RefineableQTaylorHoodElement<3>;
2007} // namespace oomph
e
Definition: cfortran.h:571
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
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
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
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
void further_build()
Further build for Crouzeix_Raviart interpolates the internal pressure dofs from father element: Make ...
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 further_build()
Further build, pass the pointers down to the sons.
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, unsigned flag)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...