refineable_spherical_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-2024 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_spherical_navier_stokes_elements.h"
27 
28 
29 namespace oomph
30 {
31  //========================================================================
32  /// Add element's contribution to the elemental
33  /// residual vector and/or Jacobian matrix.
34  /// flag=1: compute both
35  /// flag=0: compute only residual vector
36  //=======================================================================
37  void RefineableSphericalNavierStokesEquations::
38  fill_in_generic_residual_contribution_spherical_nst(
39  Vector<double>& residuals,
40  DenseMatrix<double>& jacobian,
41  DenseMatrix<double>& mass_matrix,
42  unsigned flag)
43  {
44  // Find out how many nodes there are
45  const unsigned n_node = nnode();
46 
47  // Get continuous time from timestepper of first node
48  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
49 
50  // Find out how many pressure dofs there are
51  const unsigned n_pres = npres_spherical_nst();
52 
53  // Get the local indices of the nodal coordinates
54  unsigned u_nodal_index[3];
55  for (unsigned i = 0; i < 3; ++i)
56  {
57  u_nodal_index[i] = this->u_index_spherical_nst(i);
58  }
59 
60  // Which nodal value represents the pressure? (Negative if pressure
61  // is not based on nodal interpolation).
62  const int p_index = this->p_nodal_index_spherical_nst();
63 
64  // Local array of booleans that are true if the l-th pressure value is
65  // hanging (avoid repeated virtual function calls)
66  bool pressure_dof_is_hanging[n_pres];
67  // If the pressure is stored at a node
68  if (p_index >= 0)
69  {
70  // Read out whether the pressure is hanging
71  for (unsigned l = 0; l < n_pres; ++l)
72  {
73  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
74  }
75  }
76  // Otherwise the pressure is not stored at a node and so cannot hang
77  else
78  {
79  for (unsigned l = 0; l < n_pres; ++l)
80  {
81  pressure_dof_is_hanging[l] = false;
82  }
83  }
84 
85 
86  // Set up memory for the shape and test functions
87  Shape psif(n_node), testf(n_node);
88  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
89 
90  // Set up memory for pressure shape and test functions
91  Shape psip(n_pres), testp(n_pres);
92 
93  // Set the value of n_intpt
94  const unsigned n_intpt = integral_pt()->nweight();
95 
96  // Set the Vector to hold local coordinates
97  Vector<double> s(2);
98 
99  // Get Physical Variables from Element
100  // Reynolds number must be multiplied by the density ratio
101  const double dens_ratio = density_ratio();
102  const double scaled_re = re() * dens_ratio;
103  const double scaled_re_st = re_st() * dens_ratio;
104  const double scaled_re_inv_ro = re_invro() * dens_ratio;
105  // const double scaled_re_inv_fr = re_invfr()*dens_ratio;
106  // const double visc_ratio = viscosity_ratio();
107  Vector<double> G = g();
108 
109  // Integers to store the local equation and unknown numbers
110  int local_eqn = 0, local_unknown = 0;
111 
112  // Local storage for 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++)
120  {
121  s[i] = integral_pt()->knot(ipt, i);
122  }
123 
124  // Get the integral weight
125  double w = integral_pt()->weight(ipt);
126 
127  // Call the derivatives of the shape and test functions
128  double J = dshape_and_dtest_eulerian_at_knot_spherical_nst(
129  ipt, psif, dpsifdx, testf, dtestfdx);
130 
131  // Call the pressure shape and test functions
132  pshape_spherical_nst(s, psip, testp);
133 
134  // Premultiply the weights and the Jacobian
135  double W = w * J;
136 
137  // Calculate local values of the pressure and velocity components
138  //--------------------------------------------------------------
139  double interpolated_p = 0.0;
140  Vector<double> interpolated_u(3, 0.0);
141  Vector<double> interpolated_x(2, 0.0);
142  Vector<double> mesh_velocity(2, 0.0);
143  Vector<double> dudt(3, 0.0);
144  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
145 
146  // Calculate pressure
147  for (unsigned l = 0; l < n_pres; l++)
148  {
149  interpolated_p += this->p_spherical_nst(l) * psip(l);
150  }
151 
152  // Calculate velocities and derivatives
153 
154  // Loop over nodes
155  for (unsigned l = 0; l < n_node; l++)
156  {
157  // Cache the shape function
158  double psi_ = psif(l);
159  // Loop over positions separately
160  for (unsigned i = 0; i < 2; i++)
161  {
162  interpolated_x[i] += nodal_position(l, i) * psi_;
163  }
164 
165  // Loop over velocity directions (three of these)
166  for (unsigned i = 0; i < 3; i++)
167  {
168  const double u_value = this->nodal_value(l, u_nodal_index[i]);
169  interpolated_u[i] += u_value * psi_;
170  dudt[i] += du_dt_spherical_nst(l, i) * psi_;
171 
172  // Loop over derivative directions for gradients
173  for (unsigned j = 0; j < 2; j++)
174  {
175  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
176  }
177  }
178 
179  // Only bother to calculate the mesh velocity if we are using an ALE
180  // method
181  if (!ALE_is_disabled)
182  {
183  throw OomphLibError(
184  "ALE is not properly implemented for Refineable Spherical NS yet",
185  OOMPH_CURRENT_FUNCTION,
186  OOMPH_EXCEPTION_LOCATION);
187 
188  // Loop over directions (only DIM (2) of these)
189  for (unsigned i = 0; i < 2; i++)
190  {
191  mesh_velocity[i] += this->dnodal_position_dt(l, i) * psi_;
192  }
193  }
194  } // End of loop over the nodes
195 
196  // Get the user-defined body force terms
197  Vector<double> body_force(3);
198  this->get_body_force_spherical_nst(
199  time, ipt, s, interpolated_x, body_force);
200 
201  // Get the user-defined source function
202  // double
203  // source=this->get_source_spherical_nst(time(),ipt,interpolated_x);
204 
205  // r is the first postition component
206  const double r = interpolated_x[0];
207  const double r2 = r * r;
208  // const double theta = interpolated_x[1];
209  const double sin_theta = sin(interpolated_x[1]);
210  const double cos_theta = cos(interpolated_x[1]);
211  const double cot_theta = cos_theta / sin_theta;
212 
213  const double u_r = interpolated_u[0];
214  const double u_theta = interpolated_u[1];
215  const double u_phi = interpolated_u[2];
216 
217  // Pre-calculate the scaling factor of the jacobian
218  // double W_pure = W;
219 
220  // W *= r*r*sin(theta);
221 
222  // MOMENTUM EQUATIONS
223  //==================
224  // Number of master nodes and storage for the weight of the shape function
225  unsigned n_master = 1;
226  double hang_weight = 1.0;
227 
228  // Loop over the nodes for the test functions/equations
229  //----------------------------------------------------
230  for (unsigned l = 0; l < n_node; l++)
231  {
232  // Local boolean that indicates the hanging status of the node
233  const bool is_node_hanging = node_pt(l)->is_hanging();
234 
235  // If the node is hanging
236  if (is_node_hanging)
237  {
238  // Get the hanging pointer
239  hang_info_pt = node_pt(l)->hanging_pt();
240  // Read out number of master nodes from hanging data
241  n_master = hang_info_pt->nmaster();
242  }
243  // Otherwise the node is its own master
244  else
245  {
246  n_master = 1;
247  }
248 
249  // Loop over the master nodes
250  for (unsigned m = 0; m < n_master; m++)
251  {
252  // Loop over velocity components for equations
253  for (unsigned i = 0; i < 2 + 1; i++)
254  {
255  // Get the equation number
256  // If the node is hanging
257  if (is_node_hanging)
258  {
259  // Get the equation number from the master node
260  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
261  u_nodal_index[i]);
262  // Get the hang weight from the master node
263  hang_weight = hang_info_pt->master_weight(m);
264  }
265  // If the node is not hanging
266  else
267  {
268  // Local equation number
269  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
270  // Node contributes with full weight
271  hang_weight = 1.0;
272  }
273 
274  // If it's not a boundary condition...
275  if (local_eqn >= 0)
276  {
277  // initialise for residual calculation
278  double sum = 0.0;
279 
280  switch (i)
281  {
282  // RADIAL MOMENTUM EQUATION
283  case 0:
284  {
285  // Convective r-terms
286  double conv = r * u_r * interpolated_dudx(0, 0);
287 
288  // Convective theta-terms
289  conv += u_theta * interpolated_dudx(0, 1);
290 
291  // Remaining convective terms
292  conv -= (u_theta * u_theta + u_phi * u_phi);
293 
294  // Add the time-derivative and convective terms
295  sum += (scaled_re_st * dudt[0] * r2 + scaled_re * r * conv) *
296  sin_theta * testf(l);
297 
298  // Add the user-defined body force term
299  sum -= r2 * sin_theta * body_force[0] * testf(l);
300 
301  // Add the Coriolis term
302  sum -= 2.0 * scaled_re_inv_ro * sin_theta * u_phi * r2 *
303  sin_theta * testf(l);
304 
305  // r-derivative test function component of stress tensor
306  sum += (-interpolated_p + 2 * interpolated_dudx(0, 0)) * r2 *
307  sin_theta * dtestfdx(l, 0);
308 
309  // theta-derivative test function component of stress tensor
310  sum += (r * interpolated_dudx(1, 0) - u_theta +
311  interpolated_dudx(0, 1)) *
312  sin_theta * dtestfdx(l, 1);
313 
314  // Undifferentiated test function component of stress tensor
315  sum += 2.0 *
316  ((-r * interpolated_p + +interpolated_dudx(1, 1) +
317  2.0 * u_r) *
318  sin_theta +
319  u_theta * cos_theta) *
320  testf(l);
321  }
322  break;
323 
324  // THETA-COMPONENT MOMENTUM EQUATION
325  case 1:
326  {
327  // All convective terms
328  double conv =
329  (u_r * interpolated_dudx(1, 0) * r +
330  u_theta * interpolated_dudx(1, 1) + u_r * u_theta) *
331  sin_theta -
332  u_phi * u_phi * cos_theta;
333 
334  // Add the time-derivative and convective terms to the
335  // residual
336  sum += (scaled_re_st * r2 * sin_theta * dudt[1] +
337  scaled_re * r * conv) *
338  testf(l);
339 
340  // Add the user-defined body force term
341  sum -= r2 * sin_theta * body_force[1] * testf(l);
342 
343  // Add the Coriolis term
344  sum -= 2.0 * scaled_re_inv_ro * cos_theta * u_phi * r2 *
345  sin_theta * testf(l);
346 
347  // r-derivative test function component of stress tensor
348  sum += (r * interpolated_dudx(1, 0) - u_theta +
349  interpolated_dudx(0, 1)) *
350  r * sin_theta * dtestfdx(l, 0);
351 
352  // theta-derivative test function component of stress tensor
353  sum += (-r * interpolated_p + 2.0 * interpolated_dudx(1, 1) +
354  2 * u_r) *
355  sin_theta * dtestfdx(l, 1);
356 
357  // Undifferentiated test function component of stress tensor
358  sum -= ((r * interpolated_dudx(1, 0) - u_theta +
359  interpolated_dudx(0, 1)) *
360  sin_theta -
361  (-r * interpolated_p + 2.0 * u_r +
362  2.0 * u_theta * cot_theta) *
363  cos_theta) *
364  testf(l);
365  }
366  break;
367 
368  // PHI-COMPONENT MOMENTUM EQUATION
369  case 2:
370 
371  {
372  // Convective r-terms
373  double conv = u_r * interpolated_dudx(2, 0) * r * sin_theta;
374 
375  // Convective theta-terms
376  conv += u_theta * interpolated_dudx(2, 1) * sin_theta;
377 
378  // Remaining convective terms
379  conv += u_phi * (u_r * sin_theta + u_theta * cos_theta);
380 
381  // Add the time-derivative and convective terms
382  sum += (scaled_re_st * r2 * dudt[2] * sin_theta +
383  scaled_re * conv * r) *
384  testf(l);
385 
386  sum -= r2 * sin_theta * body_force[2] * testf(l);
387 
388  // Add the Coriolis term
389  sum += 2.0 * scaled_re_inv_ro *
390  (cos_theta * u_theta + sin_theta * u_r) * r2 *
391  sin_theta * testf(l);
392 
393 
394  // r-derivative test function component of stress tensor
395  sum += (r2 * interpolated_dudx(2, 0) - r * u_phi) *
396  dtestfdx(l, 0) * sin_theta;
397 
398  // theta-derivative test function component of stress tensor
399  sum +=
400  (interpolated_dudx(2, 1) * sin_theta - u_phi * cos_theta) *
401  dtestfdx(l, 1);
402 
403  // Undifferentiated test function component of stress tensor
404  sum -= ((r * interpolated_dudx(2, 0) - u_phi) * sin_theta +
405  (interpolated_dudx(2, 1) - u_phi * cot_theta) *
406  cos_theta) *
407  testf(l);
408  }
409  break;
410  }
411 
412  // Add contribution
413  // Sign changed to be consistent with other NS implementations
414  residuals[local_eqn] -= sum * W * hang_weight;
415 
416  // CALCULATE THE JACOBIAN
417  if (flag)
418  {
419  // Number of master nodes and weights
420  unsigned n_master2 = 1;
421  double hang_weight2 = 1.0;
422  // Loop over the velocity nodes for columns
423  for (unsigned l2 = 0; l2 < n_node; l2++)
424  {
425  // Local boolean for hanging status
426  const bool is_node2_hanging = node_pt(l2)->is_hanging();
427 
428  // If the node is hanging
429  if (is_node2_hanging)
430  {
431  hang_info2_pt = node_pt(l2)->hanging_pt();
432  // Read out number of master nodes from hanging data
433  n_master2 = hang_info2_pt->nmaster();
434  }
435  // Otherwise the node is its own master
436  else
437  {
438  n_master2 = 1;
439  }
440 
441  // Loop over the master nodes
442  for (unsigned m2 = 0; m2 < n_master2; m2++)
443  {
444  // Loop over the velocity components
445  for (unsigned i2 = 0; i2 < 2 + 1; i2++)
446  {
447  // Get the number of the unknown
448  // If the node is hanging
449  if (is_node2_hanging)
450  {
451  // Get the equation number from the master node
452  local_unknown = this->local_hang_eqn(
453  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
454  hang_weight2 = hang_info2_pt->master_weight(m2);
455  }
456  else
457  {
458  local_unknown =
459  this->nodal_local_eqn(l2, u_nodal_index[i2]);
460  hang_weight2 = 1.0;
461  }
462 
463  // If the unknown is non-zero
464  if (local_unknown >= 0)
465  {
466  // Different results for i and i2
467  switch (i)
468  {
469  // RADIAL MOMENTUM EQUATION
470  case 0:
471  switch (i2)
472  {
473  // radial component
474  case 0:
475 
476  // Add the mass matrix entries
477  if (flag == 2)
478  {
479  mass_matrix(local_eqn, local_unknown) +=
480  scaled_re_st * psif(l2) * testf(l) * r2 *
481  sin_theta * W * hang_weight * hang_weight2;
482  }
483 
484  {
485  // Convective r-term contribution
486  double jac_conv =
487  r * (psif(l2) * interpolated_dudx(0, 0) +
488  u_r * dpsifdx(l2, 0));
489 
490  // Convective theta-term contribution
491  jac_conv += u_theta * dpsifdx(l2, 1);
492 
493  // Add the time-derivative contribution and
494  // the convective
495  // contribution to the sum
496  double jac_sum =
497  (scaled_re_st *
498  node_pt(l2)->time_stepper_pt()->weight(
499  1, 0) *
500  psif(l2) * r2 +
501  scaled_re * jac_conv * r) *
502  testf(l);
503 
504 
505  // Contribution from the r-derivative test
506  // function
507  // part of stress tensor
508  jac_sum +=
509  2.0 * dpsifdx(l2, 0) * dtestfdx(l, 0) * r2;
510 
511  // Contribution from the theta-derivative
512  // test function part of stress tensor
513  jac_sum += dpsifdx(l2, 1) * dtestfdx(l, 1);
514 
515 
516  // Contribution from the undifferentiated
517  // test function part
518  // of stress tensor
519  jac_sum += 4.0 * psif[l2] * testf(l);
520 
521  // Add the total contribution to the
522  // jacobian multiplied
523  // by the jacobian weight
524  jacobian(local_eqn, local_unknown) -=
525  jac_sum * sin_theta * W * hang_weight *
526  hang_weight2;
527  }
528 
529  break;
530 
531  // axial component
532  case 1:
533  {
534  // No time derivative contribution
535 
536  // Convective contribution
537  double jac_sum =
538  scaled_re *
539  (interpolated_dudx(0, 1) - 2.0 * u_theta) *
540  psif(l2) * r * sin_theta * testf(l);
541 
542  // Contribution from the theta-derivative
543  // test function
544  // part of stress tensor
545  jac_sum += (r * dpsifdx(l2, 0) - psif(l2)) *
546  dtestfdx(l, 1) * sin_theta;
547 
548  // Contribution from the undifferentiated
549  // test function
550  // part of stress tensor
551  jac_sum += 2.0 *
552  (dpsifdx(l2, 1) * sin_theta +
553  psif(l2) * cos_theta) *
554  testf(l);
555 
556  // Add the full contribution to the jacobian
557  // matrix
558  jacobian(local_eqn, local_unknown) -=
559  jac_sum * W * hang_weight * hang_weight2;
560 
561  } // End of i2 = 1 section
562 
563  break;
564 
565  // azimuthal component
566  case 2:
567  {
568  // Single convective-term contribution
569  jacobian(local_eqn, local_unknown) +=
570  2.0 * scaled_re * u_phi * psif[l2] * r *
571  sin_theta * testf[l] * W * hang_weight *
572  hang_weight2;
573 
574  // Add the Coriolis term
575  jacobian(local_eqn, local_unknown) +=
576  2.0 * scaled_re_inv_ro * sin_theta *
577  psif(l2) * r2 * sin_theta * testf[l] * W *
578  hang_weight * hang_weight2;
579  }
580 
581  break;
582  } // End of contribution radial momentum eqn
583  break;
584 
585  // AXIAL MOMENTUM EQUATION
586  case 1:
587  switch (i2)
588  {
589  // radial component
590  case 0:
591  {
592  // Convective terms
593  double jac_sum =
594  scaled_re *
595  (r2 * interpolated_dudx(1, 0) + r * u_theta) *
596  psif(l2) * sin_theta * testf(l);
597 
598  // Contribution from the r-derivative
599  // test function part of stress tensor
600  jac_sum += dpsifdx(l2, 1) * dtestfdx(l, 0) *
601  sin_theta * r;
602 
603  // Contribution from the theta-derivative
604  // test function
605  // part of stress tensor
606  jac_sum +=
607  2.0 * psif(l2) * dtestfdx(l, 1) * sin_theta;
608 
609  // Contribution from the undifferentiated
610  // test function
611  // part of stress tensor
612  jac_sum -= (dpsifdx(l2, 1) * sin_theta -
613  2.0 * psif(l2) * cos_theta) *
614  testf(l);
615 
616  // Add the sum to the jacobian
617  jacobian(local_eqn, local_unknown) -=
618  jac_sum * W * hang_weight * hang_weight2;
619  }
620 
621  break;
622 
623  // axial component
624  case 1:
625 
626  // Add the mass matrix terms
627  if (flag == 2)
628  {
629  mass_matrix(local_eqn, local_unknown) +=
630  scaled_re_st * psif[l2] * testf[l] * W *
631  r2 * sin_theta * hang_weight * hang_weight2;
632  }
633 
634 
635  {
636  // Contribution from the convective terms
637  double jac_conv =
638  r * u_r * dpsifdx(l2, 0) +
639  u_theta * dpsifdx(l2, 1) +
640  (interpolated_dudx(1, 1) + u_r) * psif(l2);
641 
642  // Add the time-derivative term and the
643  // convective terms
644  double jac_sum =
645  (scaled_re_st *
646  node_pt(l2)->time_stepper_pt()->weight(
647  1, 0) *
648  psif(l2) * r2 +
649  scaled_re * r * jac_conv) *
650  testf(l) * sin_theta;
651 
652 
653  // Contribution from the r-derivative test
654  // function
655  // part of stress tensor
656  jac_sum += (r * dpsifdx(l2, 0) - psif(l2)) *
657  r * dtestfdx(l, 0) * sin_theta;
658 
659  // Contribution from the theta-derivative
660  // test function
661  // part of stress tensor
662  jac_sum += 2.0 * dpsifdx(l2, 1) *
663  dtestfdx(l, 1) * sin_theta;
664 
665  // Contribution from the undifferentiated
666  // test function
667  // part of stress tensor
668  jac_sum -=
669  ((r * dpsifdx(l2, 0) - psif(l2)) *
670  sin_theta -
671  2.0 * psif(l2) * cot_theta * cos_theta) *
672  testf(l);
673 
674  // Add the contribution to the jacobian
675  jacobian(local_eqn, local_unknown) -=
676  jac_sum * W * hang_weight * hang_weight2;
677  }
678 
679  break;
680 
681  // azimuthal component
682  case 2:
683  {
684  // Only a convective term contribution
685  jacobian(local_eqn, local_unknown) +=
686  scaled_re * 2.0 * r * psif(l2) * u_phi *
687  cos_theta * testf(l) * W * hang_weight *
688  hang_weight2;
689 
690  // Add the Coriolis term
691  jacobian(local_eqn, local_unknown) +=
692  2.0 * scaled_re_inv_ro * cos_theta *
693  psif(l2) * r2 * sin_theta * testf[l] * W *
694  hang_weight * hang_weight2;
695  }
696 
697  break;
698  }
699  break;
700 
701  // AZIMUTHAL MOMENTUM EQUATION
702  case 2:
703  switch (i2)
704  {
705  // radial component
706  case 0:
707  {
708  // Contribution from convective terms
709  jacobian(local_eqn, local_unknown) -=
710  scaled_re *
711  (r * interpolated_dudx(2, 0) + u_phi) *
712  psif(l2) * testf(l) * r * sin_theta * W *
713  hang_weight * hang_weight2;
714 
715  // Coriolis term
716  jacobian(local_eqn, local_unknown) -=
717  2.0 * scaled_re_inv_ro * sin_theta *
718  psif(l2) * r2 * sin_theta * testf[l] * W *
719  hang_weight * hang_weight2;
720  }
721  break;
722 
723  // axial component
724  case 1:
725  {
726  // Contribution from convective terms
727  jacobian(local_eqn, local_unknown) -=
728  scaled_re *
729  (interpolated_dudx(2, 1) * sin_theta +
730  u_phi * cos_theta) *
731  r * psif(l2) * testf(l) * W * hang_weight *
732  hang_weight2;
733 
734  // Coriolis term
735  jacobian(local_eqn, local_unknown) -=
736  2.0 * scaled_re_inv_ro * cos_theta *
737  psif(l2) * r2 * sin_theta * testf[l] * W *
738  hang_weight * hang_weight2;
739  }
740 
741  break;
742 
743  // azimuthal component
744  case 2:
745 
746  // Add the mass matrix terms
747  if (flag == 2)
748  {
749  mass_matrix(local_eqn, local_unknown) +=
750  scaled_re_st * psif[l2] * testf[l] * W *
751  r2 * sin_theta * hang_weight * hang_weight2;
752  }
753 
754  {
755  // Convective terms
756  double jac_conv =
757  r * u_r * dpsifdx(l2, 0) * sin_theta;
758 
759  // Convective theta-term contribution
760  jac_conv +=
761  u_theta * dpsifdx(l2, 1) * sin_theta;
762 
763  // Contribution from the remaining convective
764  // terms
765  jac_conv +=
766  (u_r * sin_theta + u_theta * cos_theta) *
767  psif(l2);
768 
769  // Add the convective and time derivatives
770  double jac_sum =
771  (scaled_re_st *
772  node_pt(l2)->time_stepper_pt()->weight(
773  1, 0) *
774  psif(l2) * r2 * sin_theta +
775  scaled_re * r * jac_conv) *
776  testf(l);
777 
778 
779  // Contribution from the r-derivative test
780  // function
781  // part of stress tensor
782  jac_sum += (r * dpsifdx(l2, 0) - psif(l2)) *
783  dtestfdx(l, 0) * r * sin_theta;
784 
785  // Contribution from the theta-derivative
786  // test function
787  // part of stress tensor
788  jac_sum += (dpsifdx(l2, 1) * sin_theta -
789  psif(l2) * cos_theta) *
790  dtestfdx(l, 1);
791 
792  // Contribution from the undifferentiated
793  // test function
794  // part of stress tensor
795  jac_sum -=
796  ((r * dpsifdx(l2, 0) - psif(l2)) *
797  sin_theta +
798  (dpsifdx(l2, 1) - psif(l2) * cot_theta) *
799  cos_theta) *
800  testf(l);
801 
802  // Add to the jacobian
803  jacobian(local_eqn, local_unknown) -=
804  jac_sum * W * hang_weight * hang_weight2;
805  }
806 
807  break;
808  }
809  break;
810  }
811  }
812  }
813  }
814  } // End of loop over the nodes
815 
816 
817  // Loop over the pressure shape functions
818  for (unsigned l2 = 0; l2 < n_pres; l2++)
819  {
820  // If the pressure dof is hanging
821  if (pressure_dof_is_hanging[l2])
822  {
823  // Pressure dof is hanging so it must be nodal-based
824  hang_info2_pt = pressure_node_pt(l2)->hanging_pt(p_index);
825 
826  // Get the number of master nodes from the pressure node
827  n_master2 = hang_info2_pt->nmaster();
828  }
829  // Otherwise the node is its own master
830  else
831  {
832  n_master2 = 1;
833  }
834 
835  // Loop over the master nodes
836  for (unsigned m2 = 0; m2 < n_master2; m2++)
837  {
838  // Get the number of the unknown
839  // If the pressure dof is hanging
840  if (pressure_dof_is_hanging[l2])
841  {
842  // Get the unknown from the node
843  local_unknown = local_hang_eqn(
844  hang_info2_pt->master_node_pt(m2), p_index);
845  // Get the weight from the hanging object
846  hang_weight2 = hang_info2_pt->master_weight(m2);
847  }
848  else
849  {
850  local_unknown = p_local_eqn(l2);
851  hang_weight2 = 1.0;
852  }
853 
854  // If the unknown is not pinned
855  if (local_unknown >= 0)
856  {
857  // Add in contributions to different equations
858  switch (i)
859  {
860  // RADIAL MOMENTUM EQUATION
861  case 0:
862  jacobian(local_eqn, local_unknown) +=
863  psip(l2) *
864  (r2 * dtestfdx(l, 0) + 2.0 * r * testf[l]) * W *
865  sin_theta * hang_weight * hang_weight2;
866 
867 
868  break;
869 
870  // AXIAL MOMENTUM EQUATION
871  case 1:
872  jacobian(local_eqn, local_unknown) +=
873  psip(l2) * r *
874  (dtestfdx(l, 1) * sin_theta +
875  cos_theta * testf(l)) *
876  W * hang_weight * hang_weight2;
877 
878  break;
879 
880  // AZIMUTHAL MOMENTUM EQUATION
881  case 2:
882  break;
883  }
884  }
885  }
886  } // End of loop over pressure dofs
887  } // End of Jacobian calculation
888  } // End of if not boundary condition statement
889  } // End of loop over velocity components
890  } // End of loop over master nodes
891  } // End of loop over nodes
892 
893 
894  // CONTINUITY EQUATION
895  //===================
896 
897  // Loop over the pressure shape functions
898  for (unsigned l = 0; l < n_pres; l++)
899  {
900  // If the pressure dof is hanging
901  if (pressure_dof_is_hanging[l])
902  {
903  // Pressure dof is hanging so it must be nodal-based
904  hang_info_pt = pressure_node_pt(l)->hanging_pt(p_index);
905  // Get the number of master nodes from the pressure node
906  n_master = hang_info_pt->nmaster();
907  }
908  // Otherwise the node is its own master
909  else
910  {
911  n_master = 1;
912  }
913 
914  // Loop over the master nodes
915  for (unsigned m = 0; m < n_master; m++)
916  {
917  // Get the number of the unknown
918  // If the pressure dof is hanging
919  if (pressure_dof_is_hanging[l])
920  {
921  local_eqn =
922  local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
923  hang_weight = hang_info_pt->master_weight(m);
924  }
925  else
926  {
927  local_eqn = p_local_eqn(l);
928  hang_weight = 1.0;
929  }
930 
931  // If the equation is not pinned
932  if (local_eqn >= 0)
933  {
934  // The entire continuity equation
935  residuals[local_eqn] += ((2.0 * u_r + r * interpolated_dudx(0, 0) +
936  interpolated_dudx(1, 1)) *
937  sin_theta +
938  u_theta * cos_theta) *
939  r * testp(l) * W * hang_weight;
940 
941  // CALCULATE THE JACOBIAN
942  //======================
943  if (flag)
944  {
945  unsigned n_master2 = 1;
946  double hang_weight2 = 1.0;
947  // Loop over the velocity nodes for columns
948  for (unsigned l2 = 0; l2 < n_node; l2++)
949  {
950  // Local boolean to indicate whether the node is hanging
951  bool is_node2_hanging = node_pt(l2)->is_hanging();
952 
953  // If the node is hanging
954  if (is_node2_hanging)
955  {
956  hang_info2_pt = node_pt(l2)->hanging_pt();
957  // Read out number of master nodes from hanging data
958  n_master2 = hang_info2_pt->nmaster();
959  }
960  // Otherwise the node is its own master
961  else
962  {
963  n_master2 = 1;
964  }
965 
966  // Loop over the master nodes
967  for (unsigned m2 = 0; m2 < n_master2; m2++)
968  {
969  // Loop over the velocity components
970  for (unsigned i2 = 0; i2 < 2 + 1; i2++)
971  {
972  // Get the number of the unknown
973  // If the node is hanging
974  if (is_node2_hanging)
975  {
976  // Get the equation number from the master node
977  local_unknown = local_hang_eqn(
978  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
979  hang_weight2 = hang_info2_pt->master_weight(m2);
980  }
981  else
982  {
983  local_unknown =
984  this->nodal_local_eqn(l2, u_nodal_index[i2]);
985  hang_weight2 = 1.0;
986  }
987 
988  // If the unknown is not pinned
989  if (local_unknown >= 0)
990  {
991  switch (i2)
992  {
993  // radial component
994  case 0:
995  jacobian(local_eqn, local_unknown) +=
996  (2.0 * psif(l2) + r * dpsifdx(l2, 0)) * r *
997  sin_theta * testp(l) * W * hang_weight *
998  hang_weight2;
999 
1000 
1001  break;
1002 
1003  // axial component
1004  case 1:
1005  jacobian(local_eqn, local_unknown) +=
1006  r *
1007  (dpsifdx(l2, 1) * sin_theta +
1008  psif(l2) * cos_theta) *
1009  testp(l) * W * hang_weight * hang_weight2;
1010 
1011 
1012  break;
1013 
1014  // azimuthal component
1015  case 2:
1016  break;
1017  }
1018  }
1019  }
1020  }
1021  } // End of loop over nodes
1022 
1023  // NO PRESSURE CONTRIBUTIONS TO CONTINUITY EQUATION
1024  } // End of Jacobian calculation
1025  }
1026  }
1027  } // End of loop over pressure nodes
1028 
1029  } // end of loop over integration points
1030  }
1031 
1032 } // 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::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2179
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:2597
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:3992
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:1436
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2214
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1967
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:2321
oomph::FiniteElement::dnodal_position_dt
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition: elements.h:2337
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::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition: oomph_definitions.h:222
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::RefineableSphericalNavierStokesEquations::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_spherical_navier_stokes_elements.h:53
oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst
void fill_in_generic_residual_contribution_spherical_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add element's contribution to the elemental residual vector and/or Jacobian matrix flag=1: compute bo...
Definition: refineable_spherical_navier_stokes_elements.cc:38
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::SphericalNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_spherical_nst
virtual double dshape_and_dtest_eulerian_at_knot_spherical_nst(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::SphericalNavierStokesEquations::pshape_spherical_nst
virtual void pshape_spherical_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
oomph::SphericalNavierStokesEquations::u_index_spherical_nst
virtual unsigned u_index_spherical_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
Definition: spherical_navier_stokes_elements.h:387
oomph::SphericalNavierStokesEquations::re
const double & re() const
Reynolds number.
Definition: spherical_navier_stokes_elements.h:249
oomph::SphericalNavierStokesEquations::re_invro
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
Definition: spherical_navier_stokes_elements.h:311
oomph::SphericalNavierStokesEquations::p_spherical_nst
virtual double p_spherical_nst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes.
oomph::SphericalNavierStokesEquations::density_ratio
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
Definition: spherical_navier_stokes_elements.h:287
oomph::SphericalNavierStokesEquations::get_body_force_spherical_nst
virtual void get_body_force_spherical_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Calculate the body force at a given time and local and/or Eulerian position. This function is virtual...
Definition: spherical_navier_stokes_elements.h:166
oomph::SphericalNavierStokesEquations::du_dt_spherical_nst
double du_dt_spherical_nst(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes.
Definition: spherical_navier_stokes_elements.h:395
oomph::SphericalNavierStokesEquations::p_nodal_index_spherical_nst
virtual int p_nodal_index_spherical_nst() const
Return the index at which the pressure is stored if it is stored at the nodes. If not stored at the n...
Definition: spherical_navier_stokes_elements.h:448
oomph::SphericalNavierStokesEquations::g
const Vector< double > & g() const
Vector of gravitational components.
Definition: spherical_navier_stokes_elements.h:323
oomph::SphericalNavierStokesEquations::p_local_eqn
virtual int p_local_eqn(const unsigned &n) const =0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
oomph::SphericalNavierStokesEquations::npres_spherical_nst
virtual unsigned npres_spherical_nst() const =0
Function to return number of pressure degrees of freedom.
oomph::SphericalNavierStokesEquations::ALE_is_disabled
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
Definition: spherical_navier_stokes_elements.h:124
oomph::SphericalNavierStokesEquations::re_st
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: spherical_navier_stokes_elements.h:255
oomph::TimeStepper::weight
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
oomph::TimeStepper::time_pt
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
oomph::Time::time
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
oomph
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: advection_diffusion_elements.cc:30
refineable_spherical_navier_stokes_elements.h