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