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