womersley_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-2022 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
27// Non-inline functions for Womersley elements
28#include "womersley_elements.h"
29
30
31namespace oomph
32{
33 /// Default Womersley number
34 template<unsigned DIM>
36
37
38 //========================================================================
39 /// Instantiation of static bool to suppress warning; initialise to false.
40 //========================================================================
42 false;
43
44
45 //========================================================================
46 /// Zero!
47 //========================================================================
49
50
51 //======================================================================
52 // Set the data for the number of Variables at each node
53 //======================================================================
54 template<unsigned DIM, unsigned NNODE_1D>
56
57
58 //======================================================================
59 /// Compute element residual Vector and/or element Jacobian matrix
60 ///
61 /// flag=1: compute both
62 /// flag=0: compute only residual Vector
63 ///
64 /// Pure version without hanging nodes
65 //======================================================================
66 template<unsigned DIM>
68 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
69 {
70 // Find out how many nodes there are
71 unsigned n_node = nnode();
72
73 // Find the index at which the variable is stored
74 unsigned u_nodal_index = u_index_womersley();
75
76 // Set up memory for the shape and test functions
77 Shape psi(n_node), test(n_node);
78 DShape dpsidx(n_node, DIM), dtestdx(n_node, DIM);
79
80 // Set the value of n_intpt
81 unsigned n_intpt = integral_pt()->nweight();
82
83 // Set the Vector to hold local coordinates
84 Vector<double> s(DIM);
85
86 // Integers to hold the local equation and unknowns
87 int local_eqn = 0, local_unknown = 0;
88
89 // Loop over the integration points
90 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
91 {
92 // Assign values of s
93 for (unsigned i = 0; i < DIM; i++) s[i] = integral_pt()->knot(ipt, i);
94
95 // Get the integral weight
96 double w = integral_pt()->weight(ipt);
97
98 // Call the derivatives of the shape and test functions
99 double J = dshape_and_dtest_eulerian_at_knot_womersley(
100 ipt, psi, dpsidx, test, dtestdx);
101
102 // Premultiply the weights and the Jacobian
103 double W = w * J;
104
105 // Allocate memory for local quantities and initialise to zero
106 double interpolated_u = 0.0;
107 double dudt = 0.0;
108 Vector<double> interpolated_x(DIM, 0.0);
109 Vector<double> interpolated_dudx(DIM, 0.0);
110
111 // Calculate function value and derivatives:
112 // Loop over nodes
113 for (unsigned l = 0; l < n_node; l++)
114 {
115 // Calculate the value at the nodes
116 double u_value = raw_nodal_value(l, u_nodal_index);
117 interpolated_u += u_value * psi(l);
118 dudt += du_dt_womersley(l) * psi(l);
119 // Loop over directions
120 for (unsigned j = 0; j < DIM; j++)
121 {
122 interpolated_x[j] += raw_nodal_position(l, j) * psi(l);
123 interpolated_dudx[j] += u_value * dpsidx(l, j);
124 }
125 }
126
127 // Get pressure gradient
128 //---------------------
129 double dpdz;
130
131 // If no pressure gradient has been set, use zero
132 if (Pressure_gradient_data_pt == 0)
133 {
134 dpdz = 0.0;
135 }
136 else
137 {
138 // Get gressure gradient
139 dpdz = Pressure_gradient_data_pt->value(0);
140 }
141
142
143 // Assemble residuals and Jacobian
144 //--------------------------------
145
146 // Loop over the test functions
147 for (unsigned l = 0; l < n_node; l++)
148 {
149 local_eqn = nodal_local_eqn(l, u_nodal_index);
150 /*IF it's not a boundary condition*/
151 if (local_eqn >= 0)
152 {
153 // Add dpdz term and time derivative
154 residuals[local_eqn] += (re_st() * dudt + dpdz) * test(l) * W;
155
156 // Laplace operator
157 for (unsigned k = 0; k < DIM; k++)
158 {
159 residuals[local_eqn] += interpolated_dudx[k] * dtestdx(l, k) * W;
160 }
161
162
163 // Calculate the jacobian
164 //-----------------------
165 if (flag)
166 {
167 // Loop over the velocity shape functions again
168 for (unsigned l2 = 0; l2 < n_node; l2++)
169 {
170 local_unknown = nodal_local_eqn(l2, u_nodal_index);
171 // If at a non-zero degree of freedom add in the entry
172 if (local_unknown >= 0)
173 {
174 // Mass matrix
175 jacobian(local_eqn, local_unknown) +=
176 re_st() * test(l) * psi(l2) *
177 node_pt(l2)->time_stepper_pt()->weight(1, 0) * W;
178
179 // Laplace operator & mesh velocity bit
180 for (unsigned i = 0; i < DIM; i++)
181 {
182 double tmp = dtestdx(l, i);
183 jacobian(local_eqn, local_unknown) += dpsidx(l2, i) * tmp * W;
184 }
185 }
186 }
187
188 // Derivative w.r.t. pressure gradient data (if it's
189 // an unknown)
190 if ((Pressure_gradient_data_pt != 0) &&
191 (!Pressure_gradient_data_pt->is_pinned(0)))
192 {
193 local_unknown = external_local_eqn(0, 0);
194 if (local_unknown >= 0)
195 {
196 jacobian(local_eqn, local_unknown) += test(l) * W;
197 }
198
199 // Derivatives of the final eqn (volume flux constraint
200 // w.r.t. to this unknown)
201 unsigned final_local_eqn = external_local_eqn(0, 0);
202 unsigned local_unknown = local_eqn; // [from above just renaming
203 // // for clarity(!?)]
204 jacobian(final_local_eqn, local_unknown) += psi(l) * W;
205 }
206 }
207 }
208 }
209
210
211 } // End of loop over integration points
212 }
213
214
215 //======================================================================
216 /// Compute volume flux through element
217 //======================================================================
218 template<unsigned DIM>
220 {
221 // Find out how many nodes there are
222 unsigned n_node = nnode();
223
224 // Find the index at which the variable is stored
225 unsigned u_nodal_index = u_index_womersley();
226
227 // Set up memory for the shape fcs
228 Shape psi(n_node);
229 DShape dpsidx(n_node, DIM);
230
231 // Set the value of n_intpt
232 unsigned n_intpt = integral_pt()->nweight();
233
234 // Set the Vector to hold local coordinates
235 Vector<double> s(DIM);
236
237 // Initialise flux
238 double flux = 0.0;
239
240 // Loop over the integration points
241 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
242 {
243 // Assign values of s
244 for (unsigned i = 0; i < DIM; i++) s[i] = integral_pt()->knot(ipt, i);
245
246 // Get the integral weight
247 double w = integral_pt()->weight(ipt);
248
249 // Call the derivatives of the shape and test functions
250 double J = dshape_eulerian_at_knot(ipt, psi, dpsidx);
251
252 // Premultiply the weights and the Jacobian
253 double W = w * J;
254
255 // Allocate memory for local quantities and initialise to zero
256 double interpolated_u = 0.0;
257
258 // Calculate function value
259
260 // Loop over nodes
261 for (unsigned l = 0; l < n_node; l++)
262 {
263 // Calculate the value at the nodes (takes hanging node status
264 // into account
265 interpolated_u += nodal_value(l, u_nodal_index) * psi(l);
266 }
267
268 // Add to flux
269 flux += interpolated_u * W;
270
271 } // End of loop over integration points
272
273 return flux;
274 }
275
276
277 //======================================================================
278 /// Self-test: Return 0 for OK
279 //======================================================================
280 template<unsigned DIM>
282 {
283 bool passed = true;
284
285 // Check lower-level stuff
286 if (FiniteElement::self_test() != 0)
287 {
288 passed = false;
289 }
290
291 // Return verdict
292 if (passed)
293 {
294 return 0;
295 }
296 else
297 {
298 return 1;
299 }
300 }
301
302
303 //======================================================================
304 /// Output function: x,y,z_out,0,0,u,0 to allow comparison
305 /// against full Navier Stokes at n_nplot x n_plot points (2D)
306 //======================================================================
307 template<unsigned DIM>
308 void WomersleyEquations<DIM>::output_3d(std::ostream& outfile,
309 const unsigned& nplot,
310 const double& z_out)
311 {
312 // Vector of local coordinates
313 Vector<double> s(DIM);
314
315 // Tecplot header info
316 outfile << tecplot_zone_string(nplot);
317
318 // Loop over plot points
319 unsigned num_plot_points = nplot_points(nplot);
320 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
321 {
322 // Get local coordinates of plot point
323 get_s_plot(iplot, nplot, s);
324
325 for (unsigned i = 0; i < DIM; i++)
326 {
327 outfile << interpolated_x(s, i) << " ";
328 }
329 outfile << z_out << " 0.0 0.0 ";
330 outfile << interpolated_u_womersley(s);
331 outfile << " 0.0 " << std::endl;
332 }
333
334 // Write tecplot footer (e.g. FE connectivity lists)
335 write_tecplot_zone_footer(outfile, nplot);
336 }
337
338 //======================================================================
339 /// Output function:
340 ///
341 /// x,y,u or x,y,z,u
342 ///
343 /// nplot points in each coordinate direction
344 //======================================================================
345 template<unsigned DIM>
346 void WomersleyEquations<DIM>::output(std::ostream& outfile,
347 const unsigned& nplot)
348 {
349 // Vector of local coordinates
350 Vector<double> s(DIM);
351
352 // Tecplot header info
353 outfile << tecplot_zone_string(nplot);
354
355 // Loop over plot points
356 unsigned num_plot_points = nplot_points(nplot);
357 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
358 {
359 // Get local coordinates of plot point
360 get_s_plot(iplot, nplot, s);
361
362 for (unsigned i = 0; i < DIM; i++)
363 {
364 outfile << interpolated_x(s, i) << " ";
365 }
366 outfile << interpolated_u_womersley(s) << std::endl;
367 }
368
369 // Write tecplot footer (e.g. FE connectivity lists)
370 write_tecplot_zone_footer(outfile, nplot);
371 }
372
373
374 //======================================================================
375 /// C-style output function:
376 ///
377 /// x,y,u or x,y,z,u
378 ///
379 /// nplot points in each coordinate direction
380 //======================================================================
381 template<unsigned DIM>
382 void WomersleyEquations<DIM>::output(FILE* file_pt, const unsigned& nplot)
383 {
384 // Vector of local coordinates
385 Vector<double> s(DIM);
386
387 // Tecplot header info
388 fprintf(file_pt, "%s", tecplot_zone_string(nplot).c_str());
389
390 // Loop over plot points
391 unsigned num_plot_points = nplot_points(nplot);
392 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
393 {
394 // Get local coordinates of plot point
395 get_s_plot(iplot, nplot, s);
396
397 for (unsigned i = 0; i < DIM; i++)
398 {
399 fprintf(file_pt, "%g ", interpolated_x(s, i));
400 }
401 fprintf(file_pt, "%g \n", interpolated_u_womersley(s));
402 }
403
404 // Write tecplot footer (e.g. FE connectivity lists)
405 write_tecplot_zone_footer(file_pt, nplot);
406 }
407
408
409 //======================================================================
410 /// Output exact solution
411 ///
412 /// Solution is provided via function pointer.
413 /// Plot at a given number of plot points.
414 ///
415 /// x,y,u_exact or x,y,z,u_exact
416 //======================================================================
417 template<unsigned DIM>
419 std::ostream& outfile,
420 const unsigned& nplot,
422 {
423 // Vector of local coordinates
424 Vector<double> s(DIM);
425
426 // Vector for coordintes
427 Vector<double> x(DIM);
428
429 // Tecplot header info
430 outfile << tecplot_zone_string(nplot);
431
432 // Exact solution Vector (here a scalar)
433 Vector<double> exact_soln(1);
434
435 // Loop over plot points
436 unsigned num_plot_points = nplot_points(nplot);
437 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
438 {
439 // Get local coordinates of plot point
440 get_s_plot(iplot, nplot, s);
441
442 // Get x position as Vector
443 interpolated_x(s, x);
444
445 // Get exact solution at this point
446 (*exact_soln_pt)(x, exact_soln);
447
448 // Output x,y,...,u_exact
449 for (unsigned i = 0; i < DIM; i++)
450 {
451 outfile << x[i] << " ";
452 }
453 outfile << exact_soln[0] << std::endl;
454 }
455
456 // Write tecplot footer (e.g. FE connectivity lists)
457 write_tecplot_zone_footer(outfile, nplot);
458 }
459
460
461 //======================================================================
462 /// Output exact solution at time t
463 ///
464 /// Solution is provided via function pointer.
465 /// Plot at a given number of plot points.
466 ///
467 /// x,y,u_exact or x,y,z,u_exact
468 //======================================================================
469 template<unsigned DIM>
471 std::ostream& outfile,
472 const unsigned& nplot,
473 const double& time,
475
476 {
477 // Vector of local coordinates
478 Vector<double> s(DIM);
479
480 // Vector for coordintes
481 Vector<double> x(DIM);
482
483
484 // Tecplot header info
485 outfile << tecplot_zone_string(nplot);
486
487 // Exact solution Vector (here a scalar)
488 Vector<double> exact_soln(1);
489
490 // Loop over plot points
491 unsigned num_plot_points = nplot_points(nplot);
492 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
493 {
494 // Get local coordinates of plot point
495 get_s_plot(iplot, nplot, s);
496
497 // Get x position as Vector
498 interpolated_x(s, x);
499
500 // Get exact solution at this point
501 (*exact_soln_pt)(time, x, exact_soln);
502
503 // Output x,y,...,u_exact
504 for (unsigned i = 0; i < DIM; i++)
505 {
506 outfile << x[i] << " ";
507 }
508 outfile << exact_soln[0] << std::endl;
509 }
510
511 // Write tecplot footer (e.g. FE connectivity lists)
512 write_tecplot_zone_footer(outfile, nplot);
513 }
514
515
516 //======================================================================
517 /// Validate against exact solution
518 ///
519 /// Solution is provided via function pointer.
520 /// Plot error at a given number of plot points.
521 ///
522 //======================================================================
523 template<unsigned DIM>
525 std::ostream& outfile,
527 double& error,
528 double& norm)
529 {
530 // Initialise
531 error = 0.0;
532 norm = 0.0;
533
534 // Vector of local coordinates
535 Vector<double> s(DIM);
536
537 // Vector for coordintes
538 Vector<double> x(DIM);
539
540 // Find out how many nodes there are in the element
541 unsigned n_node = nnode();
542
543 Shape psi(n_node);
544
545 // Set the value of n_intpt
546 unsigned n_intpt = integral_pt()->nweight();
547
548 // Tecplot header info
549 outfile << "ZONE" << std::endl;
550
551 // Exact solution Vector (here a scalar)
552 Vector<double> exact_soln(1);
553
554 // Loop over the integration points
555 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
556 {
557 // Assign values of s
558 for (unsigned i = 0; i < DIM; i++)
559 {
560 s[i] = integral_pt()->knot(ipt, i);
561 }
562
563 // Get the integral weight
564 double w = integral_pt()->weight(ipt);
565
566 // Get jacobian of mapping
567 double J = J_eulerian(s);
568
569 // Premultiply the weights and the Jacobian
570 double W = w * J;
571
572 // Get x position as Vector
573 interpolated_x(s, x);
574
575 // Get FE function value
576 double u_fe = interpolated_u_womersley(s);
577
578 // Get exact solution at this point
579 (*exact_soln_pt)(x, exact_soln);
580
581 // Output x,y,...,error
582 for (unsigned i = 0; i < DIM; i++)
583 {
584 outfile << x[i] << " ";
585 }
586 outfile << exact_soln[0] << " " << exact_soln[0] - u_fe << std::endl;
587
588 // Add to error and norm
589 norm += exact_soln[0] * exact_soln[0] * W;
590 error += (exact_soln[0] - u_fe) * (exact_soln[0] - u_fe) * W;
591 }
592 }
593
594
595 //======================================================================
596 /// Validate against exact solution at time t.
597 ///
598 /// Solution is provided via function pointer.
599 /// Plot error at a given number of plot points.
600 ///
601 //======================================================================
602 template<unsigned DIM>
604 std::ostream& outfile,
606 const double& time,
607 double& error,
608 double& norm)
609
610 {
611 // Initialise
612 error = 0.0;
613 norm = 0.0;
614
615 // Vector of local coordinates
616 Vector<double> s(DIM);
617
618 // Vector for coordintes
619 Vector<double> x(DIM);
620
621
622 // Find out how many nodes there are in the element
623 unsigned n_node = nnode();
624
625 Shape psi(n_node);
626
627 // Set the value of n_intpt
628 unsigned n_intpt = integral_pt()->nweight();
629
630 // Tecplot
631 outfile << "ZONE" << std::endl;
632
633 // Exact solution Vector (here a scalar)
634 Vector<double> exact_soln(1);
635
636 // Loop over the integration points
637 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
638 {
639 // Assign values of s
640 for (unsigned i = 0; i < DIM; i++)
641 {
642 s[i] = integral_pt()->knot(ipt, i);
643 }
644
645 // Get the integral weight
646 double w = integral_pt()->weight(ipt);
647
648 // Get jacobian of mapping
649 double J = J_eulerian(s);
650
651 // Premultiply the weights and the Jacobian
652 double W = w * J;
653
654 // Get x position as Vector
655 interpolated_x(s, x);
656
657 // Get FE function value
658 double u_fe = interpolated_u_womersley(s);
659
660 // Get exact solution at this point
661 (*exact_soln_pt)(time, x, exact_soln);
662
663 // Output x,y,...,error
664 for (unsigned i = 0; i < DIM; i++)
665 {
666 outfile << x[i] << " ";
667 }
668 outfile << exact_soln[0] << " " << exact_soln[0] - u_fe << std::endl;
669
670 // Add to error and norm
671 norm += exact_soln[0] * exact_soln[0] * W;
672 error += (exact_soln[0] - u_fe) * (exact_soln[0] - u_fe) * W;
673 }
674 }
675
676
677 //====================================================================
678 // Force build of templates
679 //====================================================================
680
681 template class WomersleyEquations<1>;
682
683 template class QWomersleyElement<1, 2>;
684 template class QWomersleyElement<1, 3>;
685 template class QWomersleyElement<1, 4>;
686
687 template class WomersleyEquations<2>;
688
689 template class QWomersleyElement<2, 2>;
690 template class QWomersleyElement<2, 3>;
691 template class QWomersleyElement<2, 4>;
692
693 // template class QWomersleyElement<3,2>;
694 // template class QWomersleyElement<3,3>;
695 // template class QWomersleyElement<3,4>;
696
697} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1765
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK.
Definition: elements.cc:4440
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
static const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
static bool Suppress_warning_about_unpinned_nst_dofs
Static bool to suppress warning.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
virtual void fill_in_generic_residual_contribution_womersley(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
void output(std::ostream &outfile)
Output with default number of plot points.
static double Default_ReSt_value
Static default value for the Womersley number.
double get_volume_flux()
Compute total volume flux through element.
unsigned self_test()
Self-test: Return 0 for OK.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points.
void output_3d(std::ostream &outfile, const unsigned &n_plot, const double &z_out)
Output function: x,y,z_out,0,0,u,0 to allow comparison against full Navier Stokes at n_nplot x n_plot...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...