refineable_polar_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//====================================================================
26 #include "refineable_polar_navier_stokes_elements.h"
27 
28 namespace oomph
29 {
30  /// ///////////////////////////////////////////////////////////////////////
31  //======================================================================//
32  /// Start of what would've been refineable_navier_stokes_elements.cc //
33  //======================================================================//
34  /// ///////////////////////////////////////////////////////////////////////
35 
36  //==============================================================
37  /// Compute the residuals for the Navier--Stokes
38  /// equations; flag=1(or 0): do (or don't) compute the
39  /// Jacobian as well.
40  /// flag=2 for Residuals, Jacobian and mass_matrix
41  ///
42  /// This is now my new version with Jacobian and
43  /// dimensionless phi
44  ///
45  /// This version supports hanging nodes
46  //==============================================================
47  void RefineablePolarNavierStokesEquations::
48  fill_in_generic_residual_contribution(Vector<double>& residuals,
49  DenseMatrix<double>& jacobian,
50  DenseMatrix<double>& mass_matrix,
51  unsigned flag)
52  {
53  // Find out how many nodes there are
54  unsigned n_node = nnode();
55 
56  // Find out how many pressure dofs there are
57  unsigned n_pres = npres_pnst();
58 
59  // Find the indices at which the local velocities are stored
60  unsigned u_nodal_index[2];
61  for (unsigned i = 0; i < 2; i++)
62  {
63  u_nodal_index[i] = u_index_pnst(i);
64  }
65 
66  // Which nodal value represents the pressure? (Negative if pressure
67  // is not based on nodal interpolation).
68  int p_index = this->p_nodal_index_pnst();
69 
70  // Local array of booleans that are true if the l-th pressure value is
71  // hanging (avoid repeated virtual function calls)
72  bool pressure_dof_is_hanging[n_pres];
73  // If the pressure is stored at a node
74  if (p_index >= 0)
75  {
76  // Read out whether the pressure is hanging
77  for (unsigned l = 0; l < n_pres; ++l)
78  {
79  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
80  }
81  }
82  // Otherwise the pressure is not stored at a node and so cannot hang
83  else
84  {
85  for (unsigned l = 0; l < n_pres; ++l)
86  {
87  pressure_dof_is_hanging[l] = false;
88  }
89  }
90 
91  // Set up memory for the shape and test functions
92  Shape psif(n_node), testf(n_node);
93  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
94 
95  // Set up memory for pressure shape and test functions
96  Shape psip(n_pres), testp(n_pres);
97 
98  // Number of integration points
99  unsigned n_intpt = integral_pt()->nweight();
100 
101  // Set the Vector to hold local coordinates
102  Vector<double> s(2);
103 
104  // Get the reynolds number and Alpha
105  const double Re = re();
106  const double Alpha = alpha();
107  const double Re_St = re_st();
108 
109  // Integers to store the local equations and unknowns
110  int local_eqn = 0, local_unknown = 0;
111 
112  // Pointers to hang info objects
113  HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
114 
115  // Loop over the integration points
116  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
117  {
118  // Assign values of s
119  for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
120  // Get the integral weight
121  double w = integral_pt()->weight(ipt);
122 
123  // Call the derivatives of the shape and test functions
124  double J = dshape_and_dtest_eulerian_at_knot_pnst(
125  ipt, psif, dpsifdx, testf, dtestfdx);
126 
127  // Call the pressure shape and test functions
128  this->pshape_pnst(s, psip, testp);
129 
130  // Premultiply the weights and the Jacobian
131  double W = w * J;
132 
133  // Calculate local values of the pressure and velocity components
134  // Allocate storage initialised to zero
135  double interpolated_p = 0.0;
136  Vector<double> interpolated_u(2, 0.0);
137  Vector<double> interpolated_x(2, 0.0);
138  // Vector<double> dudt(2);
139  DenseMatrix<double> interpolated_dudx(2, 2, 0.0);
140 
141  // Initialise to zero
142  for (unsigned i = 0; i < 2; i++)
143  {
144  // dudt[i] = 0.0;
145  interpolated_u[i] = 0.0;
146  interpolated_x[i] = 0.0;
147  for (unsigned j = 0; j < 2; j++)
148  {
149  interpolated_dudx(i, j) = 0.0;
150  }
151  }
152 
153  // Calculate pressure
154  for (unsigned l = 0; l < n_pres; l++)
155  interpolated_p += this->p_pnst(l) * psip[l];
156 
157  // Calculate velocities and derivatives:
158 
159  // Loop over nodes
160  for (unsigned l = 0; l < n_node; l++)
161  {
162  // Loop over directions
163  for (unsigned i = 0; i < 2; i++)
164  {
165  // Get the nodal value
166  double u_value = this->nodal_value(l, u_nodal_index[i]);
167  interpolated_u[i] += u_value * psif[l];
168  interpolated_x[i] += this->nodal_position(l, i) * psif[l];
169  // dudt[i]+=du_dt_pnst(l,i)*psif[l];
170 
171  // Loop over derivative directions
172  for (unsigned j = 0; j < 2; j++)
173  {
174  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
175  }
176  }
177  }
178 
179  // MOMENTUM EQUATIONS
180  //------------------
181  // Number of master nodes and storage for the weight of the shape function
182  unsigned n_master = 1;
183  double hang_weight = 1.0;
184  // Loop over the nodes for the test functions
185  for (unsigned l = 0; l < n_node; l++)
186  {
187  // Local boolean to indicate whether the node is hanging
188  bool is_node_hanging = node_pt(l)->is_hanging();
189 
190  // If the node is hanging
191  if (is_node_hanging)
192  {
193  hang_info_pt = node_pt(l)->hanging_pt();
194  // Read out number of master nodes from hanging data
195  n_master = hang_info_pt->nmaster();
196  }
197  // Otherwise the node is its own master
198  else
199  {
200  n_master = 1;
201  }
202 
203  // Now add in a new loop over the master nodes
204  for (unsigned m = 0; m < n_master; m++)
205  {
206  // Can't loop over velocity components as don't have identical
207  // contributions Do seperately for i = {0,1} instead
208  unsigned i = 0;
209  {
210  // Get the equation number
211  // If the node is hanging
212  if (is_node_hanging)
213  {
214  // Get the equation number from the master node
215  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
216  u_nodal_index[i]);
217  // Get the hang weight from the master node
218  hang_weight = hang_info_pt->master_weight(m);
219  }
220  // If the node is not hanging
221  else
222  {
223  // Local equation number
224  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
225  // Node contributes with full weight
226  hang_weight = 1.0;
227  }
228 
229  /*IF it's not a boundary condition*/
230  if (local_eqn >= 0)
231  {
232  // Add the testf[l] term of the stress tensor
233  residuals[local_eqn] +=
234  ((interpolated_p / interpolated_x[0]) -
235  ((1. + Gamma[i]) / pow(interpolated_x[0], 2.)) *
236  ((1. / Alpha) * interpolated_dudx(1, 1) +
237  interpolated_u[0])) *
238  testf[l] * interpolated_x[0] * Alpha * W * hang_weight;
239 
240  // Add the dtestfdx(l,0) term of the stress tensor
241  residuals[local_eqn] +=
242  (interpolated_p - (1. + Gamma[i]) * interpolated_dudx(0, 0)) *
243  dtestfdx(l, 0) * interpolated_x[0] * Alpha * W * hang_weight;
244 
245  // Add the dtestfdx(l,1) term of the stress tensor
246  residuals[local_eqn] -=
247  ((1. / (interpolated_x[0] * Alpha)) * interpolated_dudx(0, 1) -
248  (interpolated_u[1] / interpolated_x[0]) +
249  Gamma[i] * interpolated_dudx(1, 0)) *
250  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
251  interpolated_x[0] * Alpha * W * hang_weight;
252 
253  // Convective terms
254  residuals[local_eqn] -=
255  Re *
256  (interpolated_u[0] * interpolated_dudx(0, 0) +
257  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
258  interpolated_dudx(0, 1) -
259  (pow(interpolated_u[1], 2.) / interpolated_x[0])) *
260  testf[l] * interpolated_x[0] * Alpha * W * hang_weight;
261 
262 
263  // CALCULATE THE JACOBIAN
264  if (flag)
265  {
266  // Number of master nodes and weights
267  unsigned n_master2 = 1;
268  double hang_weight2 = 1.0;
269  // Loop over the velocity shape functions again
270  for (unsigned l2 = 0; l2 < n_node; l2++)
271  {
272  // Local boolean to indicate whether the node is hanging
273  bool is_node2_hanging = node_pt(l2)->is_hanging();
274 
275  // If the node is hanging
276  if (is_node2_hanging)
277  {
278  hang_info2_pt = node_pt(l2)->hanging_pt();
279  // Read out number of master nodes from hanging data
280  n_master2 = hang_info2_pt->nmaster();
281  }
282  // Otherwise the node is its own master
283  else
284  {
285  n_master2 = 1;
286  }
287 
288  // Loop over the master nodes
289  for (unsigned m2 = 0; m2 < n_master2; m2++)
290  {
291  // Again can't loop over velocity components due to loss of
292  // symmetry
293  unsigned i2 = 0;
294  {
295  // Get the number of the unknown
296  // If the node is hanging
297  if (is_node2_hanging)
298  {
299  // Get the equation number from the master node
300  local_unknown = this->local_hang_eqn(
301  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
302  // Get the hang weights
303  hang_weight2 = hang_info2_pt->master_weight(m2);
304  }
305  else
306  {
307  local_unknown =
308  this->nodal_local_eqn(l2, u_nodal_index[i2]);
309  hang_weight2 = 1.0;
310  }
311 
312  // If at a non-zero degree of freedom add in the entry
313  if (local_unknown >= 0)
314  {
315  // Add contribution to Elemental Matrix
316  jacobian(local_eqn, local_unknown) -=
317  (1. + Gamma[i]) *
318  (psif[l2] / pow(interpolated_x[0], 2.)) * testf[l] *
319  interpolated_x[0] * Alpha * W * hang_weight *
320  hang_weight2;
321 
322  jacobian(local_eqn, local_unknown) -=
323  (1. + Gamma[i]) * dpsifdx(l2, 0) * dtestfdx(l, 0) *
324  interpolated_x[0] * Alpha * W * hang_weight *
325  hang_weight2;
326 
327  jacobian(local_eqn, local_unknown) -=
328  (1. / (interpolated_x[0] * Alpha)) * dpsifdx(l2, 1) *
329  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
330  interpolated_x[0] * Alpha * W * hang_weight *
331  hang_weight2;
332 
333  // Now add in the inertial terms
334  jacobian(local_eqn, local_unknown) -=
335  Re *
336  (psif[l2] * interpolated_dudx(0, 0) +
337  interpolated_u[0] * dpsifdx(l2, 0) +
338  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
339  dpsifdx(l2, 1)) *
340  testf[l] * interpolated_x[0] * Alpha * W *
341  hang_weight * hang_weight2;
342 
343  // extra bit for mass matrix
344  if (flag == 2)
345  {
346  mass_matrix(local_eqn, local_unknown) +=
347  Re_St * psif[l2] * testf[l] * interpolated_x[0] *
348  Alpha * W * hang_weight * hang_weight2;
349  }
350 
351  } // End of (Jacobian's) if not boundary condition
352  // statement
353  } // End of i2=0 section
354 
355  i2 = 1;
356  {
357  // Get the number of the unknown
358  // If the node is hanging
359  if (is_node2_hanging)
360  {
361  // Get the equation number from the master node
362  local_unknown = this->local_hang_eqn(
363  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
364  // Get the hang weights
365  hang_weight2 = hang_info2_pt->master_weight(m2);
366  }
367  else
368  {
369  local_unknown =
370  this->nodal_local_eqn(l2, u_nodal_index[i2]);
371  hang_weight2 = 1.0;
372  }
373 
374  // If at a non-zero degree of freedom add in the entry
375  if (local_unknown >= 0)
376  {
377  // Add contribution to Elemental Matrix
378  jacobian(local_eqn, local_unknown) -=
379  ((1. + Gamma[i]) /
380  (pow(interpolated_x[0], 2.) * Alpha)) *
381  dpsifdx(l2, 1) * testf[l] * interpolated_x[0] *
382  Alpha * W * hang_weight * hang_weight2;
383 
384  jacobian(local_eqn, local_unknown) -=
385  (-(psif[l2] / interpolated_x[0]) +
386  Gamma[i] * dpsifdx(l2, 0)) *
387  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
388  interpolated_x[0] * Alpha * W * hang_weight *
389  hang_weight2;
390 
391  // Now add in the inertial terms
392  jacobian(local_eqn, local_unknown) -=
393  Re *
394  ((psif[l2] / (interpolated_x[0] * Alpha)) *
395  interpolated_dudx(0, 1) -
396  2 * interpolated_u[1] *
397  (psif[l2] / interpolated_x[0])) *
398  testf[l] * interpolated_x[0] * Alpha * W *
399  hang_weight * hang_weight2;
400 
401  } // End of (Jacobian's) if not boundary condition
402  // statement
403  } // End of i2=1 section
404 
405  } // End of loop over master nodes m2
406  } // End of l2 loop
407 
408  /*Now loop over pressure shape functions*/
409  /*This is the contribution from pressure gradient*/
410  for (unsigned l2 = 0; l2 < n_pres; l2++)
411  {
412  // If the pressure dof is hanging
413  if (pressure_dof_is_hanging[l2])
414  {
415  hang_info2_pt =
416  this->pressure_node_pt(l2)->hanging_pt(p_index);
417  // Pressure dof is hanging so it must be nodal-based
418  // Get the number of master nodes from the pressure node
419  n_master2 = hang_info2_pt->nmaster();
420  }
421  // Otherwise the node is its own master
422  else
423  {
424  n_master2 = 1;
425  }
426 
427  // Loop over the master nodes
428  for (unsigned m2 = 0; m2 < n_master2; m2++)
429  {
430  // Get the number of the unknown
431  // If the pressure dof is hanging
432  if (pressure_dof_is_hanging[l2])
433  {
434  // Get the unknown from the master node
435  local_unknown = this->local_hang_eqn(
436  hang_info2_pt->master_node_pt(m2), p_index);
437  // Get the weight from the hanging object
438  hang_weight2 = hang_info2_pt->master_weight(m2);
439  }
440  else
441  {
442  local_unknown = this->p_local_eqn(l2);
443  hang_weight2 = 1.0;
444  }
445 
446  /*If we are at a non-zero degree of freedom in the entry*/
447  if (local_unknown >= 0)
448  {
449  jacobian(local_eqn, local_unknown) +=
450  (psip[l2] / interpolated_x[0]) * testf[l] *
451  interpolated_x[0] * Alpha * W * hang_weight *
452  hang_weight2;
453 
454  jacobian(local_eqn, local_unknown) +=
455  psip[l2] * dtestfdx(l, 0) * interpolated_x[0] * Alpha *
456  W * hang_weight * hang_weight2;
457 
458  } // End of Jacobian pressure if not a boundary condition
459  // statement
460 
461  } // End of loop over master nodes m2
462  } // End of loop over pressure shape functions l2
463 
464  } /*End of Jacobian calculation*/
465 
466  } // End of if not boundary condition statement
467  } // End of i=0 section
468 
469  i = 1;
470  {
471  // Get the equation number
472  // If the node is hanging
473  if (is_node_hanging)
474  {
475  // Get the equation number from the master node
476  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
477  u_nodal_index[i]);
478  // Get the hang weight from the master node
479  hang_weight = hang_info_pt->master_weight(m);
480  }
481  // If the node is not hanging
482  else
483  {
484  // Local equation number
485  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
486  // Node contributes with full weight
487  hang_weight = 1.0;
488  }
489 
490  /*IF it's not a boundary condition*/
491  if (local_eqn >= 0)
492  {
493  // Add the testf[l] term of the stress tensor
494  residuals[local_eqn] +=
495  ((1. / (pow(interpolated_x[0], 2.) * Alpha)) *
496  interpolated_dudx(0, 1) -
497  (interpolated_u[1] / pow(interpolated_x[0], 2.)) +
498  Gamma[i] * (1. / interpolated_x[0]) *
499  interpolated_dudx(1, 0)) *
500  testf[l] * interpolated_x[0] * Alpha * W * hang_weight;
501 
502  // Add the dtestfdx(l,0) term of the stress tensor
503  residuals[local_eqn] -=
504  (interpolated_dudx(1, 0) +
505  Gamma[i] * ((1. / (interpolated_x[0] * Alpha)) *
506  interpolated_dudx(0, 1) -
507  (interpolated_u[1] / interpolated_x[0]))) *
508  dtestfdx(l, 0) * interpolated_x[0] * Alpha * W * hang_weight;
509 
510  // Add the dtestfdx(l,1) term of the stress tensor
511  residuals[local_eqn] +=
512  (interpolated_p -
513  (1. + Gamma[i]) * ((1. / (interpolated_x[0] * Alpha)) *
514  interpolated_dudx(1, 1) +
515  (interpolated_u[0] / interpolated_x[0]))) *
516  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
517  interpolated_x[0] * Alpha * W * hang_weight;
518 
519  // Convective terms
520  residuals[local_eqn] -=
521  Re *
522  (interpolated_u[0] * interpolated_dudx(1, 0) +
523  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
524  interpolated_dudx(1, 1) +
525  ((interpolated_u[0] * interpolated_u[1]) /
526  interpolated_x[0])) *
527  testf[l] * interpolated_x[0] * Alpha * W * hang_weight;
528 
529  // CALCULATE THE JACOBIAN
530  if (flag)
531  {
532  // Number of master nodes and weights
533  unsigned n_master2 = 1;
534  double hang_weight2 = 1.0;
535 
536  // Loop over the velocity shape functions again
537  for (unsigned l2 = 0; l2 < n_node; l2++)
538  {
539  // Local boolean to indicate whether the node is hanging
540  bool is_node2_hanging = node_pt(l2)->is_hanging();
541 
542  // If the node is hanging
543  if (is_node2_hanging)
544  {
545  hang_info2_pt = node_pt(l2)->hanging_pt();
546  // Read out number of master nodes from hanging data
547  n_master2 = hang_info2_pt->nmaster();
548  }
549  // Otherwise the node is its own master
550  else
551  {
552  n_master2 = 1;
553  }
554 
555  // Loop over the master nodes
556  for (unsigned m2 = 0; m2 < n_master2; m2++)
557  {
558  // Again can't loop over velocity components due to loss of
559  // symmetry
560  unsigned i2 = 0;
561  {
562  // Get the number of the unknown
563  // If the node is hanging
564  if (is_node2_hanging)
565  {
566  // Get the equation number from the master node
567  local_unknown = this->local_hang_eqn(
568  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
569  // Get the hang weights
570  hang_weight2 = hang_info2_pt->master_weight(m2);
571  }
572  else
573  {
574  local_unknown =
575  this->nodal_local_eqn(l2, u_nodal_index[i2]);
576  hang_weight2 = 1.0;
577  }
578 
579  // If at a non-zero degree of freedom add in the entry
580  if (local_unknown >= 0)
581  {
582  // Add contribution to Elemental Matrix
583  jacobian(local_eqn, local_unknown) +=
584  (1. / (pow(interpolated_x[0], 2.) * Alpha)) *
585  dpsifdx(l2, 1) * testf[l] * interpolated_x[0] *
586  Alpha * W * hang_weight * hang_weight2;
587 
588  jacobian(local_eqn, local_unknown) -=
589  Gamma[i] * (1. / (interpolated_x[0] * Alpha)) *
590  dpsifdx(l2, 1) * dtestfdx(l, 0) * interpolated_x[0] *
591  Alpha * W * hang_weight * hang_weight2;
592 
593  jacobian(local_eqn, local_unknown) -=
594  (1 + Gamma[i]) * (psif[l2] / interpolated_x[0]) *
595  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
596  interpolated_x[0] * Alpha * W * hang_weight *
597  hang_weight2;
598 
599  // Now add in the inertial terms
600  jacobian(local_eqn, local_unknown) -=
601  Re *
602  (psif[l2] * interpolated_dudx(1, 0) +
603  (psif[l2] * interpolated_u[1] / interpolated_x[0])) *
604  testf[l] * interpolated_x[0] * Alpha * W *
605  hang_weight * hang_weight2;
606 
607  } // End of (Jacobian's) if not boundary condition
608  // statement
609  } // End of i2=0 section
610 
611  i2 = 1;
612  {
613  // Get the number of the unknown
614  // If the node is hanging
615  if (is_node2_hanging)
616  {
617  // Get the equation number from the master node
618  local_unknown = this->local_hang_eqn(
619  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
620  // Get the hang weights
621  hang_weight2 = hang_info2_pt->master_weight(m2);
622  }
623  else
624  {
625  local_unknown =
626  this->nodal_local_eqn(l2, u_nodal_index[i2]);
627  hang_weight2 = 1.0;
628  }
629 
630  // If at a non-zero degree of freedom add in the entry
631  if (local_unknown >= 0)
632  {
633  // Add contribution to Elemental Matrix
634  jacobian(local_eqn, local_unknown) +=
635  (-(psif[l2] / pow(interpolated_x[0], 2.)) +
636  Gamma[i] * (1. / interpolated_x[0]) *
637  dpsifdx(l2, 0)) *
638  testf[l] * interpolated_x[0] * Alpha * W *
639  hang_weight * hang_weight2;
640 
641  jacobian(local_eqn, local_unknown) -=
642  (dpsifdx(l2, 0) -
643  Gamma[i] * (psif[l2] / interpolated_x[0])) *
644  dtestfdx(l, 0) * interpolated_x[0] * Alpha * W *
645  hang_weight * hang_weight2;
646 
647  jacobian(local_eqn, local_unknown) -=
648  (1. + Gamma[i]) * (1. / (interpolated_x[0] * Alpha)) *
649  dpsifdx(l2, 1) * (1. / (interpolated_x[0] * Alpha)) *
650  dtestfdx(l, 1) * interpolated_x[0] * Alpha * W *
651  hang_weight * hang_weight2;
652 
653  // Now add in the inertial terms
654  jacobian(local_eqn, local_unknown) -=
655  Re *
656  (interpolated_u[0] * dpsifdx(l2, 0) +
657  (psif[l2] / (interpolated_x[0] * Alpha)) *
658  interpolated_dudx(1, 1) +
659  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
660  dpsifdx(l2, 1) +
661  (interpolated_u[0] * psif[l2] / interpolated_x[0])) *
662  testf[l] * interpolated_x[0] * Alpha * W *
663  hang_weight * hang_weight2;
664 
665  // extra bit for mass matrix
666  if (flag == 2)
667  {
668  mass_matrix(local_eqn, local_unknown) +=
669  Re_St * psif[l2] * testf[l] * interpolated_x[0] *
670  Alpha * W * hang_weight * hang_weight2;
671  }
672 
673  } // End of (Jacobian's) if not boundary condition
674  // statement
675  } // End of i2=1 section
676 
677  } // End of loop over master nodes m2
678  } // End of l2 loop
679 
680  /*Now loop over pressure shape functions*/
681  /*This is the contribution from pressure gradient*/
682  for (unsigned l2 = 0; l2 < n_pres; l2++)
683  {
684  // If the pressure dof is hanging
685  if (pressure_dof_is_hanging[l2])
686  {
687  hang_info2_pt =
688  this->pressure_node_pt(l2)->hanging_pt(p_index);
689  // Pressure dof is hanging so it must be nodal-based
690  // Get the number of master nodes from the pressure node
691  n_master2 = hang_info2_pt->nmaster();
692  }
693  // Otherwise the node is its own master
694  else
695  {
696  n_master2 = 1;
697  }
698 
699  // Loop over the master nodes
700  for (unsigned m2 = 0; m2 < n_master2; m2++)
701  {
702  // Get the number of the unknown
703  // If the pressure dof is hanging
704  if (pressure_dof_is_hanging[l2])
705  {
706  // Get the unknown from the master node
707  local_unknown = this->local_hang_eqn(
708  hang_info2_pt->master_node_pt(m2), p_index);
709  // Get the weight from the hanging object
710  hang_weight2 = hang_info2_pt->master_weight(m2);
711  }
712  else
713  {
714  local_unknown = this->p_local_eqn(l2);
715  hang_weight2 = 1.0;
716  }
717 
718 
719  /*If we are at a non-zero degree of freedom in the entry*/
720  if (local_unknown >= 0)
721  {
722  jacobian(local_eqn, local_unknown) +=
723  (psip[l2] / interpolated_x[0]) * (1. / Alpha) *
724  dtestfdx(l, 1) * interpolated_x[0] * Alpha * W *
725  hang_weight * hang_weight2;
726 
727  } // End of if not boundary condition for pressure in
728  // jacobian
729 
730  } // End of loop over master nodes m2
731  } // End of loop over pressure test functions l2
732 
733  } /*End of Jacobian calculation*/
734 
735  } // End of if not boundary condition statement
736  } // End of i=1 section
737 
738  } // End of loop over master nodes m
739  } // End of loop over shape functions
740 
741 
742  // CONTINUITY EQUATION
743  //-------------------
744 
745  // Loop over the shape functions
746  for (unsigned l = 0; l < n_pres; l++)
747  {
748  // If the pressure dof is hanging
749  if (pressure_dof_is_hanging[l])
750  {
751  // Pressure dof is hanging so it must be nodal-based
752  // Get the hang info object
753  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
754  // Get the number of master nodes from the pressure node
755  n_master = hang_info_pt->nmaster();
756  }
757  // Otherwise the node is its own master
758  else
759  {
760  n_master = 1;
761  }
762 
763  // Loop over the master nodes
764  for (unsigned m = 0; m < n_master; m++)
765  {
766  // Get the number of the unknown
767  // If the pressure dof is hanging
768  if (pressure_dof_is_hanging[l])
769  {
770  // Get the local equation from the master node
771  local_eqn =
772  this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
773  // Get the weight for the node
774  hang_weight = hang_info_pt->master_weight(m);
775  }
776  else
777  {
778  local_eqn = this->p_local_eqn(l);
779  hang_weight = 1.0;
780  }
781 
782  // If not a boundary conditions
783  if (local_eqn >= 0)
784  {
785  residuals[local_eqn] +=
786  (interpolated_dudx(0, 0) +
787  (interpolated_u[0] / interpolated_x[0]) +
788  (1. / (interpolated_x[0] * Alpha)) * interpolated_dudx(1, 1)) *
789  testp[l] * interpolated_x[0] * Alpha * W * hang_weight;
790 
791  /*CALCULATE THE JACOBIAN*/
792  if (flag)
793  {
794  unsigned n_master2 = 1;
795  double hang_weight2 = 1.0;
796  /*Loop over the velocity shape functions*/
797  for (unsigned l2 = 0; l2 < n_node; l2++)
798  {
799  // Local boolean to indicate whether the node is hanging
800  bool is_node2_hanging = node_pt(l2)->is_hanging();
801 
802  // If the node is hanging
803  if (is_node2_hanging)
804  {
805  hang_info2_pt = node_pt(l2)->hanging_pt();
806  // Read out number of master nodes from hanging data
807  n_master2 = hang_info2_pt->nmaster();
808  }
809  // Otherwise the node is its own master
810  else
811  {
812  n_master2 = 1;
813  }
814 
815  // Loop over the master nodes
816  for (unsigned m2 = 0; m2 < n_master2; m2++)
817  {
818  unsigned i2 = 0;
819  {
820  // Get the number of the unknown
821  // If the node is hanging
822  if (is_node2_hanging)
823  {
824  // Get the equation number from the master node
825  local_unknown = this->local_hang_eqn(
826  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
827  hang_weight2 = hang_info2_pt->master_weight(m2);
828  }
829  else
830  {
831  local_unknown =
832  this->nodal_local_eqn(l2, u_nodal_index[i2]);
833  hang_weight2 = 1.0;
834  }
835  /*If we're at a non-zero degree of freedom add it in*/
836  if (local_unknown >= 0)
837  {
838  jacobian(local_eqn, local_unknown) +=
839  (dpsifdx(l2, 0) + (psif[l2] / interpolated_x[0])) *
840  testp[l] * interpolated_x[0] * Alpha * W * hang_weight *
841  hang_weight2;
842  }
843  } // End of i2=0 section
844 
845  i2 = 1;
846  {
847  // Get the number of the unknown
848  // If the node is hanging
849  if (is_node2_hanging)
850  {
851  // Get the equation number from the master node
852  local_unknown = this->local_hang_eqn(
853  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
854  hang_weight2 = hang_info2_pt->master_weight(m2);
855  }
856  else
857  {
858  local_unknown =
859  this->nodal_local_eqn(l2, u_nodal_index[i2]);
860  hang_weight2 = 1.0;
861  }
862 
863  /*If we're at a non-zero degree of freedom add it in*/
864  if (local_unknown >= 0)
865  {
866  jacobian(local_eqn, local_unknown) +=
867  (1. / (interpolated_x[0] * Alpha)) * dpsifdx(l2, 1) *
868  testp[l] * interpolated_x[0] * Alpha * W * hang_weight *
869  hang_weight2;
870  }
871  } // End of i2=1 section
872 
873  } // End of loop over master nodes m2
874  } /*End of loop over l2*/
875  } /*End of Jacobian calculation*/
876 
877  } // End of if not boundary condition
878  } // End of loop over master nodes m
879  } // End of loop over pressure test functions l
880 
881  } // End of loop over integration points
882 
883  } // End of add_generic_residual_contribution
884 
885 
886 } // End of namespace oomph
s
static char t char * s
Definition: cfortran.h:568
i
cstr elem_len * i
Definition: cfortran.h:603
oomph::DShape
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::nodal_value
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
oomph::FiniteElement::interpolated_x
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
oomph::FiniteElement::nodal_local_eqn
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
oomph::FiniteElement::nodal_position
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
oomph::HangInfo
Class that contains data for hanging nodes.
Definition: nodes.h:742
oomph::HangInfo::master_weight
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
oomph::HangInfo::master_node_pt
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
oomph::HangInfo::nmaster
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
oomph::Integral::knot
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
oomph::Integral::nweight
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
oomph::Integral::weight
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
oomph::Node::hanging_pt
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
oomph::Node::is_hanging
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
oomph::PolarNavierStokesEquations::pshape_pnst
virtual void pshape_pnst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
oomph::PolarNavierStokesEquations::Gamma
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
Definition: polar_navier_stokes_elements.h:261
oomph::PolarNavierStokesEquations::p_local_eqn
virtual int p_local_eqn(const unsigned &n)=0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
oomph::PolarNavierStokesEquations::p_pnst
virtual double p_pnst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes.
oomph::PolarNavierStokesEquations::p_nodal_index_pnst
virtual int p_nodal_index_pnst()
Which nodal value represents the pressure? (Default: negative, indicating that pressure is not based ...
Definition: polar_navier_stokes_elements.h:432
oomph::PolarNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_pnst
virtual double dshape_and_dtest_eulerian_at_knot_pnst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and derivatives w.r.t. global coords at ipt-th integration point Return J...
oomph::PolarNavierStokesEquations::u_index_pnst
virtual unsigned u_index_pnst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
Definition: polar_navier_stokes_elements.h:394
oomph::PolarNavierStokesEquations::re
const double & re() const
Reynolds number.
Definition: polar_navier_stokes_elements.h:266
oomph::PolarNavierStokesEquations::re_st
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: polar_navier_stokes_elements.h:278
oomph::PolarNavierStokesEquations::npres_pnst
virtual unsigned npres_pnst() const =0
Function to return number of pressure degrees of freedom.
oomph::RefineableElement::local_hang_eqn
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
Definition: refineable_elements.h:278
oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution
virtual void fill_in_generic_residual_contribution(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...
Definition: refineable_polar_navier_stokes_elements.cc:48
oomph::RefineablePolarNavierStokesEquations::pressure_node_pt
virtual Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interp...
Definition: refineable_polar_navier_stokes_elements.h:62
oomph::Shape
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
oomph
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: advection_diffusion_elements.cc:30
refineable_polar_navier_stokes_elements.h