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-2023 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
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
29 namespace 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
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
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
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;
419  {
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 =
562  for (unsigned e = 0; e < nrobin; e++)
563  {
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
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
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
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
1246  if (nod_pt->has_auxiliary_node_update_fct_pt())
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);
1406  time, ipt, s, interpolated_x, d_body_force_dx);
1407
1408  // Get gradient of source function
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
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
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 =
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  {
1841
1842  // What type of son am I? Ask my quadtree representation...
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  //====================================================================
1999  template class RefineableNavierStokesEquations<2>;
2000  template class RefineableNavierStokesEquations<3>;
2001  template class RefineableQTaylorHoodElement<2>;
2002  template class RefineableQTaylorHoodElement<3>;
2003  template class RefineableQCrouzeixRaviartElement<2>;
2004  template class RefineableQCrouzeixRaviartElement<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
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
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
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
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
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
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...