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