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