womersley_elements.h
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
28// Header file for Womersley elements
29#ifndef OOMPH_WOMERSLEY_ELEMENTS_HEADER
30#define OOMPH_WOMERSLEY_ELEMENTS_HEADER
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37
38// OOMPH-LIB headers
39#include "../generic/nodes.h"
40#include "../generic/Qelements.h"
41#include "../generic/mesh.h"
42#include "../generic/problem.h"
43#include "../generic/oomph_utilities.h"
44#include "../navier_stokes/navier_stokes_flux_control_elements.h"
45
46
47namespace oomph
48{
49 /// //////////////////////////////////////////////////////////////////////
50 /// //////////////////////////////////////////////////////////////////////
51 /// //////////////////////////////////////////////////////////////////////
52
53
54 //===============================================================
55 /// Template-free base class for Impedance Tube -- to faciliate
56 /// interactions between the Womersley elements and the
57 /// Navier Stokes impedance traction elements.
58 //===============================================================
60 {
61 public:
62 /// Empty constructor
64
65 /// Empty virtual destructor
67
68 /// Empty virtual dummy member function -- every base class
69 /// needs at least one virtual member function if it's
70 /// to be used as a base class for a polymorphic object.
71 // virtual void dummy(){};
72
73 /// Pure virtual function to compute inlet pressure, p_in,
74 /// required to achieve the currently imposed, instantaneous
75 /// volume flux q prescribed by total_volume_flux_into_impedance_tube(),
76 /// and its derivative, dp_in/dq.
77 virtual void get_response(double& p_in, double& dp_in_dq) = 0;
78
79 /// Zero!
80 static double Zero;
81 };
82
83
84 /// //////////////////////////////////////////////////////////////////////
85 /// //////////////////////////////////////////////////////////////////////
86 /// //////////////////////////////////////////////////////////////////////
87
88
89 //======================================================================
90 /// A base class for elements that allow the imposition of an impedance
91 /// type boundary condition to the Navier--Stokes equations. Establishes
92 /// the template-free common functionality, that they must have to be able
93 /// to compute the volume flux that passes through them, etc.
94 //======================================================================
96 {
97 public:
98 // Empty constructor
100
101 /// Empty vitual destructor
103
104 /// Pure virtual function that must be implemented to compute
105 /// the volume flux that passes through this element
106 virtual double get_volume_flux() = 0;
107
108 /// Add the element's contribution to the auxiliary integral
109 /// used in the element's Jacobian. The aux_integral contains
110 /// the derivative of the total volume flux through the
111 /// outflow boundary of the (higher-dimensional) Navier-Stokes mesh w.r.t.
112 /// to the discrete (global) (velocity) degrees of freedom.
114 std::map<unsigned, double>* aux_integral_pt) = 0;
115
116 /// Pass the pointer to the pre-computed auxiliary integral
117 /// to the element so it can be accessed when computing the
118 /// elemental Jacobian.
120 std::map<unsigned, double>* aux_integral_pt) = 0;
121
122 /// Pass the pointer to the "impedance tube" that computes
123 /// the flow resistance via the solution of Womersley's equations
124 /// to the element.
126 TemplateFreeWomersleyImpedanceTubeBase* impedance_tube_pt) = 0;
127
128
129 /// Pass the pointer to the mesh containing all
130 /// NavierStokesImpedanceTractionElements that contribute
131 /// to the volume flux into the downstream "impedance tube"
132 /// to the element and classify all nodes in that mesh
133 /// as external Data for this element (unless the nodes
134 /// are also the element's own nodes, of course).
136 Mesh* navier_stokes_outflow_mesh_pt_mesh_pt) = 0;
137
138 protected:
139 /// Pointer to mesh containing the
140 /// NavierStokesImpedanceTractionElements
141 /// that contribute to the volume flux into the "downstream tube" that
142 /// provides the flow resistance
144 };
145
146
147 /// //////////////////////////////////////////////////////////////////////
148 /// //////////////////////////////////////////////////////////////////////
149 /// //////////////////////////////////////////////////////////////////////
150
151
152 //=============================================================
153 /// A class for all isoparametric elements that solve the
154 /// Womersley (parallel flow) equations.
155 /// \f[ Re St \frac{\partial u}{\partial t} = - g + \frac{\partial^2 u}{\partial x_i^2} \f]
156 /// which may be derived from the full Navier-Stokes equations
157 /// (with a viscous scaling of the pressure) under the assumption of
158 /// parallel flow in the z direction. u then represents the axial
159 /// velocity and g is the (spatially constant) axial component of
160 /// the pressure gradient.
161 ///
162 /// This class contains the generic maths. Shape functions, geometric
163 /// mapping etc. must get implemented in derived class.
164 /// Note that this class assumes an isoparametric formulation, i.e. that
165 /// the scalar unknown is interpolated using the same shape functions
166 /// as the position.
167 ///
168 /// Generally, the instantaneous value of the pressure gradient, g,
169 /// is prescribed (and specified via a pointer to a single-valued
170 /// Data object whose current (pinned) value contains the pressure.
171 ///
172 /// It is also possible to prescribe the flow rate through
173 /// a mesh of Womersley elements and to determine the pressure
174 /// gradient required to achieve this flow rate as an unknown.
175 /// In that case the external pressure is treated as
176 /// an external Data object that an associated
177 /// ImposeFluxForWomersleyElement is in charge of. Note that only
178 /// the ImposeFluxForWomersleyElement can set the pressure gradient
179 /// Data object as external Data. This is because (counter to
180 /// general practice) the WomersleyEquations make contributions
181 /// to the residuals of the ImposeFluxForWomersleyElements in order
182 /// to keep the elemental Jacobians as small as possible.
183 //=============================================================
184 template<unsigned DIM>
185 class WomersleyEquations : public virtual FiniteElement
186 {
187 public:
188 /// Constructor: Initialises the Pressure_gradient_data_pt to null
190 {
192 }
193
194
195 /// Broken copy constructor
197
198 /// Broken assignment operator
199 void operator=(const WomersleyEquations&) = delete;
200
201 /// Set pointer to pressure gradient (single-valued Data)
202 void set_pressure_gradient_pt(Data*& pressure_gradient_data_pt)
203 {
204#ifdef PARANOID
205 if (pressure_gradient_data_pt->nvalue() != 1)
206 {
207 throw OomphLibError(
208 "Pressure gradient Data must only contain a single value!\n",
209 OOMPH_CURRENT_FUNCTION,
210 OOMPH_EXCEPTION_LOCATION);
211 }
212#endif
213 Pressure_gradient_data_pt = pressure_gradient_data_pt;
214 }
215
216
217 /// Read-only access to pointer to pressure gradient
219 {
221 }
222
223
224 /// Product of Reynolds and Strouhal number (=Womersley number)
225 const double& re_st() const
226 {
227 return *ReSt_pt;
228 }
229
230
231 /// Pointer to product of Reynolds and Strouhal number (=Womersley number)
232 double*& re_st_pt()
233 {
234 return ReSt_pt;
235 }
236
237
238 /// Return the index at which the unknown value
239 /// is stored. The default value, 0, is appropriate for single-physics
240 /// problems, when there is only one variable, the value that satisfies the
241 /// Womersley equation.
242 /// In derived multi-physics elements, this function should be overloaded
243 /// to reflect the chosen storage scheme. Note that these equations require
244 /// that the unknown is always stored at the same index at each node.
245 virtual inline unsigned u_index_womersley() const
246 {
247 return 0;
248 }
249
250
251 /// du/dt at local node n.
252 /// Uses suitably interpolated value for hanging nodes.
253 double du_dt_womersley(const unsigned& n) const
254 {
255 // Get the data's timestepper
257
258 // Initialise dudt
259 double dudt = 0.0;
260
261 // Loop over the timesteps, if there is a non Steady timestepper
263 {
264 // Find the index at which the variable is stored
265 const unsigned u_nodal_index = u_index_womersley();
266
267 // Number of timsteps (past & present)
268 const unsigned n_time = time_stepper_pt->ntstorage();
269
270 // Add the contributions to the time derivative
271 for (unsigned t = 0; t < n_time; t++)
272 {
273 dudt +=
274 time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
275 }
276 }
277 return dudt;
278 }
279
280
281 /// Output with default number of plot points
282 void output(std::ostream& outfile)
283 {
284 unsigned nplot = 5;
285 output(outfile, nplot);
286 }
287
288 /// Output function: x,y,z_out,0,0,u,0 to allow comparison
289 /// against full Navier Stokes at n_nplot x n_plot points (2D)
290 void output_3d(std::ostream& outfile,
291 const unsigned& n_plot,
292 const double& z_out);
293
294 /// Output FE representation of soln: x,y,u or x,y,z,u at
295 /// n_plot^DIM plot points
296 void output(std::ostream& outfile, const unsigned& nplot);
297
298
299 /// C_style output with default number of plot points
300 void output(FILE* file_pt)
301 {
302 unsigned n_plot = 5;
303 output(file_pt, n_plot);
304 }
305
306
307 /// C-style output FE representation of soln: x,y,u or x,y,z,u at
308 /// n_plot^DIM plot points
309 void output(FILE* file_pt, const unsigned& n_plot);
310
311
312 /// Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points
313 void output_fct(std::ostream& outfile,
314 const unsigned& nplot,
316
317
318 /// Output exact soln: x,y,u_exact or x,y,z,u_exact at
319 /// nplot^DIM plot points (time-dependent version)
320 virtual void output_fct(
321 std::ostream& outfile,
322 const unsigned& nplot,
323 const double& time,
325
326
327 /// Get error against and norm of exact solution
328 void compute_error(std::ostream& outfile,
330 double& error,
331 double& norm);
332
333
334 /// Get error against and norm of exact solution
335 void compute_error(std::ostream& outfile,
337 const double& time,
338 double& error,
339 double& norm);
340
341 /// Get flux: flux[i] = du/dx_i
342 void get_flux(const Vector<double>& s, Vector<double>& flux) const
343 {
344 // Find out how many nodes there are in the element
345 unsigned n_node = nnode();
346
347 // Find the index at which the variable is stored
348 unsigned u_nodal_index = u_index_womersley();
349
350 // Set up memory for the shape and test functions
351 Shape psi(n_node);
352 DShape dpsidx(n_node, DIM);
353
354 // Call the derivatives of the shape and test functions
355 dshape_eulerian(s, psi, dpsidx);
356
357 // Initialise to zero
358 for (unsigned j = 0; j < DIM; j++)
359 {
360 flux[j] = 0.0;
361 }
362
363 // Loop over nodes
364 for (unsigned l = 0; l < n_node; l++)
365 {
366 // Loop over derivative directions
367 for (unsigned j = 0; j < DIM; j++)
368 {
369 flux[j] += nodal_value(l, u_nodal_index) * dpsidx(l, j);
370 }
371 }
372 }
373
374
375 /// Compute element residual Vector (wrapper)
377 {
378 // Call the generic residuals function with flag set to 0
379 // using a dummy matrix argument
382 }
383
384
385 /// Compute element residual Vector and element Jacobian matrix (wrapper)
387 DenseMatrix<double>& jacobian)
388 {
389 // Call the generic routine with the flag set to 1
391 }
392
393
394 /// Return FE representation of function value u(s) at local coordinate s
395 inline double interpolated_u_womersley(const Vector<double>& s) const
396 {
397 // Find number of nodes
398 unsigned n_node = nnode();
399
400 // Find the index at which the variable is stored
401 unsigned u_nodal_index = u_index_womersley();
402
403 // Local shape function
404 Shape psi(n_node);
405
406 // Find values of shape function
407 shape(s, psi);
408
409 // Initialise value of u
410 double interpolated_u = 0.0;
411
412 // Loop over the local nodes and sum
413 for (unsigned l = 0; l < n_node; l++)
414 {
415 interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
416 }
417
418 return (interpolated_u);
419 }
420
421 /// Self-test: Return 0 for OK
422 unsigned self_test();
423
424 /// Compute total volume flux through element
425 double get_volume_flux();
426
427 protected:
428 /// Shape/test functions and derivs w.r.t. to global coords at
429 /// local coord. s; return Jacobian of mapping
431 const Vector<double>& s,
432 Shape& psi,
433 DShape& dpsidx,
434 Shape& test,
435 DShape& dtestdx) const = 0;
436
437
438 /// Shape/test functions and derivs w.r.t. to global coords at
439 /// integration point ipt; return Jacobian of mapping
441 const unsigned& ipt,
442 Shape& psi,
443 DShape& dpsidx,
444 Shape& test,
445 DShape& dtestdx) const = 0;
446
447 /// Compute element residual Vector only (if flag=and/or element
448 /// Jacobian matrix
450 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
451
452 /// Pointer to pressure gradient Data (single value Data item)
454
455 /// Pointer to global Reynolds number x Strouhal number (=Womersley)
456 double* ReSt_pt;
457
458 /// Static default value for the Womersley number
459 static double Default_ReSt_value;
460
461 private:
462 template<unsigned DIMM>
464
465 /// Set pointer to pressure gradient (single-valued Data) and
466 /// treat it as external data -- this can only be called
467 /// by the friend class ImposeFluxForWomersleyElement which
468 /// imposes a volume flux constraint and trades it
469 /// for the now unknown pressure gradient Data that is
470 /// treated as external Data for this element.
471 /// This slightly convoluted (private/friend) construction
472 /// is necessary because, counter to practice, the current
473 /// element adds contributions to the equation that determines
474 /// the external data. This obviously requires that the
475 /// ImposeFluxForWomersleyElement doesn't do the same. We know
476 /// that it doesn't and therefore we make it a friend that's allowed
477 /// to collaborate with this element...
479 Data* pressure_gradient_data_pt)
480 {
481 Pressure_gradient_data_pt = pressure_gradient_data_pt;
482 add_external_data(pressure_gradient_data_pt);
483 }
484 };
485
486
487 /// ////////////////////////////////////////////////////////////////////////
488 /// ////////////////////////////////////////////////////////////////////////
489 /// ////////////////////////////////////////////////////////////////////////
490
491
492 //========================================================================
493 /// Element to impose volume flux through collection of Womersley
494 /// elements, in exchange for treating the pressure gradient
495 /// as an unknown. The pressure gradient is created (as a single-valued
496 /// Data item) in the constructor for this element which also
497 /// takes a pointer to the Mesh containing the Womersley elements whose
498 /// total flux is being controlled. While doing this we tell them
499 /// that their pressure gradient is now an unknown and must be
500 /// treated as external Data.
501 //========================================================================
502 template<unsigned DIM>
504 {
505 public:
506 /// Constructor: Pass pointer to mesh that contains the
507 /// Womersley elements whose volume flux is controlled, and pointer to
508 /// double that contains the instantaneous value of the
509 /// prescribed flux
511 double* prescribed_flux_pt)
512 : Prescribed_flux_pt(prescribed_flux_pt)
513 {
514 // Store the mesh of the flux-controlled Womerersley elements
515 Womersley_mesh_pt = womersley_mesh_pt;
516
517 // Create the pressure gradient Data
520
521 // Pressure gradient is internal data of this element
523
524 // Find number of elements in the mesh of Womersley elements
525 // whose total flux is controlled
526 unsigned n_element = womersley_mesh_pt->nelement();
527
528 // Loop over the elements to tell them that the pressure
529 // gradient is given by the newly created Data object
530 // which is to be treated as external Data.
531 for (unsigned e = 0; e < n_element; e++)
532 {
533 // Upcast from FiniteElement to the present element
534 WomersleyEquations<DIM>* el_pt = dynamic_cast<WomersleyEquations<DIM>*>(
535 womersley_mesh_pt->element_pt(e));
536
537 // Set the pressure gradient function pointer
540 }
541 }
542
543 /// Read-only access to the single-valued Data item that
544 /// stores the pressure gradient (to be determined via the
545 /// flux control)
547 {
549 }
550
551
552 /// Get volume flux through all Womersley elements
554 {
555 // Initialise
556 double flux = 0.0;
557
558 // Assemble contributions from elements
559 unsigned nelem = Womersley_mesh_pt->nelement();
560 for (unsigned e = 0; e < nelem; e++)
561 {
562 WomersleyEquations<DIM>* el_pt = dynamic_cast<WomersleyEquations<DIM>*>(
564 if (el_pt != 0) flux += el_pt->get_volume_flux();
565 }
566
567 // Return total volume flux
568 return flux;
569 }
570
571
572 /// Compute residual vector: the volume flux constraint
573 /// determines this element's one-and-only internal Data which represents
574 /// the pressure gradient
576 {
577 // Local equation number of volume flux constraint -- associated
578 // with the internal data (the unknown pressure gradient)
579 int local_eqn = internal_local_eqn(0, 0);
580 if (local_eqn >= 0)
581 {
582 residuals[local_eqn] += total_volume_flux() - (*Prescribed_flux_pt);
583 }
584 }
585
586
587 /// Compute element residual Vector and element Jacobian matrix
588 /// Note: Jacobian is zero because the derivatives w.r.t. to
589 /// velocity dofs are added by the Womersley elements; the current
590 /// element's internal Data (the pressure gradient) does not feature
591 /// in the volume constraint.
593 {
594 // Initialise Jacobian
595 unsigned n_dof = ndof();
596 for (unsigned i = 0; i < n_dof; i++)
597 {
598 for (unsigned j = 0; j < n_dof; j++)
599 {
600 jacobian(i, j) = 0.0;
601 }
602 }
603 // Get residuals
604 get_residuals(residuals);
605 }
606
607 private:
608 /// Pointer to mesh that contains the Womersley elements
610
611 /// Data item whose one and only value contains the pressure
612 /// gradient
614
615 /// Pointer to current value of prescribed flux
617 };
618
619
620 /// ////////////////////////////////////////////////////////////////////////
621 /// ////////////////////////////////////////////////////////////////////////
622 /// ////////////////////////////////////////////////////////////////////////
623
624
625 //======================================================================
626 /// QWomersleyElement elements are linear/quadrilateral/brick-shaped
627 /// Womersley elements with isoparametric interpolation for the function.
628 //======================================================================
629 template<unsigned DIM, unsigned NNODE_1D>
630 class QWomersleyElement : public virtual QElement<DIM, NNODE_1D>,
631 public virtual WomersleyEquations<DIM>
632 {
633 private:
634 /// Static array of ints to hold number of variables at
635 /// nodes: Initial_Nvalue[n]
636 static const unsigned Initial_Nvalue;
637
638 public:
639 /// Constructor: Call constructors for QElement and
640 /// Womersley equations
641 QWomersleyElement() : QElement<DIM, NNODE_1D>(), WomersleyEquations<DIM>()
642 {
643 }
644
645 /// Broken copy constructor
647
648 /// Broken assignment operator
650
651 /// Required # of `values' (pinned or dofs)
652 /// at node n
653 inline unsigned required_nvalue(const unsigned& n) const
654 {
655 return Initial_Nvalue;
656 }
657
658 /// Output function:
659 /// x,y,u or x,y,z,u
660 void output(std::ostream& outfile)
661 {
663 }
664
665
666 /// Output function:
667 /// x,y,u or x,y,z,u at n_plot^DIM plot points
668 void output(std::ostream& outfile, const unsigned& n_plot)
669 {
670 WomersleyEquations<DIM>::output(outfile, n_plot);
671 }
672
673
674 /// C-style output function:
675 /// x,y,u or x,y,z,u
676 void output(FILE* file_pt)
677 {
679 }
680
681
682 /// C-style output function:
683 /// x,y,u or x,y,z,u at n_plot^DIM plot points
684 void output(FILE* file_pt, const unsigned& n_plot)
685 {
686 WomersleyEquations<DIM>::output(file_pt, n_plot);
687 }
688
689
690 /// Output function for an exact solution:
691 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
692 void output_fct(std::ostream& outfile,
693 const unsigned& n_plot,
695 {
696 WomersleyEquations<DIM>::output_fct(outfile, n_plot, exact_soln_pt);
697 }
698
699
700 /// Output function for a time-dependent exact solution.
701 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
702 /// (Calls the steady version)
703 void output_fct(std::ostream& outfile,
704 const unsigned& n_plot,
705 const double& time,
707 {
708 WomersleyEquations<DIM>::output_fct(outfile, n_plot, time, exact_soln_pt);
709 }
710
711
712 protected:
713 /// Shape, test functions & derivs. w.r.t. to global coords. Return
714 /// Jacobian.
716 Shape& psi,
717 DShape& dpsidx,
718 Shape& test,
719 DShape& dtestdx) const;
720
721
722 /// Shape/test functions and derivs w.r.t. to global coords at
723 /// integration point ipt; return Jacobian of mapping
725 const unsigned& ipt,
726 Shape& psi,
727 DShape& dpsidx,
728 Shape& test,
729 DShape& dtestdx) const;
730 };
731
732
733 /// /////////////////////////////////////////////////////////////////////
734 /// /////////////////////////////////////////////////////////////////////
735 /// /////////////////////////////////////////////////////////////////////
736
737
738 //=====start_of_problem_class=========================================
739 /// Womersley problem
740 //====================================================================
741 template<class ELEMENT, unsigned DIM>
743 {
744 public:
745 /// Function pointer to fct that prescribes pressure gradient
746 /// g=fct(t)
747 typedef double (*PrescribedPressureGradientFctPt)(const double& time);
748
749
750 /// Constructor: Pass pointer to Womersley number, pointer to the
751 /// double that stores the currently imposed flow rate, the pointer
752 /// to the mesh of WomersleyElements (of the type specified by the template
753 /// argument) and the TimeStepper used in these
754 /// elements. (Note that the mesh must be created, and boundary
755 /// conditions applied BEFORE creating the Womersley problem.
756 /// This is to facilitate the re-use of this class
757 /// for different geometries.
758 WomersleyProblem(double* re_st_pt,
759 double* prescribed_volume_flux_pt,
761 Mesh* womersley_mesh_pt);
762
763
764 /// Constructor: Pass pointer to Womersley number, pointer to the
765 /// function that returns the imposed pressure gradient, the pointer
766 /// to the mesh of WomersleyElements and the TimeStepper used in these
767 /// elements. (Note that the mesh must be created, and boundary
768 /// conditions applied BEFORE creating the Womersley problem.
769 /// This is to allow the facilitate the re-use of this class
770 /// for different geometries.)
771 WomersleyProblem(double* re_st_pt,
772 PrescribedPressureGradientFctPt pressure_gradient_fct_pt,
774 Mesh* womersley_mesh_pt);
775
776 /// Destructor to clean up memory
778
779 /// Update the problem specs after solve (empty)
781
782 /// Update the problem specs before solve (empty)
784
785
786 /// Update the problem specs before next timestep:
787 /// Update time-varying pressure gradient (if prescribed)
789 {
790 /// Assign current prescribed pressure gradient to Data
792 {
795 }
796 }
797
798 /// Doc the solution incl. trace file for various quantities of
799 /// interest (to some...)
800 void doc_solution(DocInfo& doc_info,
801 std::ofstream& trace_file,
802 const double& z_out = 0.0);
803
804 /// Doc the solution
805 void doc_solution(DocInfo& doc_info, const double& z_out = 0.0)
806 {
807 std::ofstream trace_file;
808 doc_solution(doc_info, trace_file, z_out);
809 }
810
811 /// Access function to the single-valued Data object that
812 /// contains the unknown pressure gradient (used if flux is prescribed)
814 {
816 }
817
818
819 private:
820 /// Pointer to currently prescribed volume flux
822
823 /// Pointer to element that imposes the flux through the collection
824 /// of Womersley elements
826
827 /// Fct pointer to fct that prescribes pressure gradient
829
830 /// Pointer to single-valued Data item that stores pressure gradient
832
833 }; // end of problem class
834
835
836 //========start_of_constructor============================================
837 /// Constructor: Pass pointer to Womersley number, fct pointer to the
838 /// function that returns the prescribed pressure gradient, the pointer
839 /// to the mesh of WomersleyElements (of the type specified by the
840 /// template argument), and the TimeStepper used in these
841 /// elements. (Note that the mesh must be created, and boundary
842 /// conditions applied BEFORE creating the Womersley problem.
843 /// This is to facilitate the re-use of this class
844 /// for different geometries.)
845 //========================================================================
846 template<class ELEMENT, unsigned DIM>
848 double* re_st_pt,
849 PrescribedPressureGradientFctPt pressure_gradient_fct_pt,
850 TimeStepper* time_stepper_pt,
851 Mesh* womersley_mesh_pt)
852 : Prescribed_volume_flux_pt(0),
853 Flux_el_pt(0),
854 Prescribed_pressure_gradient_fct_pt(pressure_gradient_fct_pt)
855 {
856 // Problem is linear: Skip convergence check in Newton solver
857 Problem_is_nonlinear = false;
858
859 // Set the timestepper
861
862 // Set the mesh (bcs have already been allocated!)
863 mesh_pt() = womersley_mesh_pt;
864
865 // Complete the build of all elements so they are fully functional
866 //----------------------------------------------------------------
867
868 // Find number of elements in mesh
869 unsigned n_element = mesh_pt()->nelement();
870
871 // Loop over the elements to set up element-specific
872 // things that cannot be handled by constructor
873 for (unsigned i = 0; i < n_element; i++)
874 {
875 // Upcast from FiniteElement to the present element
876 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
877
878 // Set pointer to Womersley number
879 el_pt->re_st_pt() = re_st_pt;
880 }
881
882 // Create pressure gradient as pinned, single-valued Data item
885
886 // Pass pointer to pressure gradient Data to elements
887 unsigned nelem = mesh_pt()->nelement();
888 for (unsigned e = 0; e < nelem; e++)
889 {
890 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
891 if (el_pt != 0)
892 {
893 el_pt->set_pressure_gradient_pt(Pressure_gradient_data_pt);
894 }
895 }
896
897
898 // Do equation numbering
899 oomph_info << "Number of equations in WomersleyProblem: "
900 << assign_eqn_numbers() << std::endl;
901
902 } // end of constructor
903
904
905 //========start_of_constructor============================================
906 /// Constructor: Pass pointer to Womersley number, pointer to the
907 /// double that stores the currently imposed flow rate, the pointer
908 /// to the mesh of WomersleyElements and the TimeStepper used in these
909 /// elements. (Note that the mesh must be created, and boundary
910 /// conditions applied BEFORE creating the Womersley problem.
911 /// This is to facilitate the re-use of this class
912 /// for different geometries.)
913 //========================================================================
914 template<class ELEMENT, unsigned DIM>
916 double* re_st_pt,
917 double* prescribed_volume_flux_pt,
918 TimeStepper* time_stepper_pt,
919 Mesh* womersley_mesh_pt)
920 : Prescribed_volume_flux_pt(prescribed_volume_flux_pt),
921 Flux_el_pt(0),
922 Prescribed_pressure_gradient_fct_pt(0)
923 {
924 // Problem is linear: Skip convergence check in Newton solver
925 Problem_is_nonlinear = false;
926
927 // Set the timestepper
929
930 // Set the mesh (bcs have already been allocated!)
931 mesh_pt() = womersley_mesh_pt;
932
933
934 // Complete the build of all elements so they are fully functional
935 //----------------------------------------------------------------
936
937 // Find number of elements in mesh
938 unsigned n_element = mesh_pt()->nelement();
939
940 // Loop over the elements to set up element-specific
941 // things that cannot be handled by constructor
942 for (unsigned i = 0; i < n_element; i++)
943 {
944 // Upcast from FiniteElement to the present element
945 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
946
947 // Set pointer to Womersley number
948 el_pt->re_st_pt() = re_st_pt;
949 }
950
951 // Create element that imposes the flux -- this element creates
952 // the single-valued Data object that represents the (unknown)
953 // pressure gradient internally. It also passes the pointer to
954 // this Data object to the Womersley elements contained in the
955 // mesh. The Womersley elements treat this Data as external Data.
958
959 // Add the ImposeFluxForWomersleyElement to the mesh
961
962 // Store pressure gradient data that was
963 // created in the ImposeFluxForWomersleyElement
964 Pressure_gradient_data_pt = Flux_el_pt->pressure_gradient_data_pt();
965
966 // Do equation numbering
967 oomph_info << "Number of equations in WomersleyProblem: "
968 << assign_eqn_numbers() << std::endl;
969
970 } // end of constructor
971
972
973 //======start_of_destructor===============================================
974 /// Destructor for Womersley problem
975 //========================================================================
976 template<class ELEMENT, unsigned DIM>
978 {
979 // Timestepper gets killed in general problem destructor
980
981 // Mesh gets killed in general problem destructor
982
983 } // end of destructor
984
985
986 //=======start_of_doc_solution============================================
987 /// Doc the solution
988 //========================================================================
989 template<class ELEMENT, unsigned DIM>
991 std::ofstream& trace_file,
992 const double& z_out)
993 {
994 std::ofstream some_file;
995 std::ofstream some_file1;
996 std::ostringstream filename;
997
998 // Number of plot points
999 unsigned npts;
1000 npts = 5;
1001
1002
1003 // Compute total volume flux directly
1004 double flux = 0.0;
1005
1006 // Output solution
1007 //-----------------
1008 filename << doc_info.directory() << "/womersley_soln" << doc_info.number()
1009 << ".dat";
1010 some_file1.open(filename.str().c_str());
1011
1012 filename.str("");
1013 filename << doc_info.directory() << "/womersley_soln_3d_"
1014 << doc_info.number() << ".dat";
1015 some_file.open(filename.str().c_str());
1016
1017 // Assemble contributions from elements and output 3D solution
1018 unsigned nelem = mesh_pt()->nelement();
1019 for (unsigned e = 0; e < nelem; e++)
1020 {
1021 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
1022 if (el_pt != 0)
1023 {
1024 flux += el_pt->get_volume_flux();
1025 el_pt->output_3d(some_file, npts, z_out);
1026 el_pt->output(some_file1, npts);
1027 }
1028 }
1029 some_file.close();
1030 some_file1.close();
1031
1032 double prescribed_g = 0.0;
1033 if (Prescribed_pressure_gradient_fct_pt != 0)
1034 {
1035 prescribed_g = Prescribed_pressure_gradient_fct_pt(time_pt()->time());
1036 }
1037
1038
1039 double prescribed_q = 0.0;
1040 if (Prescribed_volume_flux_pt != 0)
1041 {
1042 prescribed_q = *Prescribed_volume_flux_pt;
1043 }
1044
1045 if (trace_file.is_open())
1046 {
1047 trace_file << time_pt()->time() << " "
1048 << pressure_gradient_data_pt()->value(0) << " " << flux << " "
1049 << prescribed_g << " " << prescribed_q << " " << std::endl;
1050 }
1051
1052 } // end of doc_solution
1053
1054
1055 /// /////////////////////////////////////////////////////////////////////
1056 /// /////////////////////////////////////////////////////////////////////
1057 /// /////////////////////////////////////////////////////////////////////
1058
1059
1060 //====================================================================
1061 /// Base class for Womersley impedance tube. Allows the computation
1062 /// of the inlet pressure p_in into a uniform tube of specified length
1063 /// that is assumed to convey fully-developed, but time-dependent flow with a
1064 /// presribed instantaneous flow rate, q. Also computes the derivative
1065 /// dp_in/dq required when this is used to determine impedance-type
1066 /// outlet boundary conditions in a Navier-Stokes computation.
1067 //====================================================================
1068 template<class ELEMENT, unsigned DIM>
1071 {
1072 public:
1073 /// Function pointer to fct that prescribes volume flux
1074 /// q=fct(t) -- mainly used for validation purposes.
1075 typedef double (*PrescribedVolumeFluxFctPt)(const double& time);
1076
1077
1078 /// Constructor: Specify length of tube and pointer to function that
1079 /// specifies the prescribed volume flux. Outlet pressure is set to zero.
1081 const double& length,
1082 PrescribedVolumeFluxFctPt prescribed_volume_flux_fct_pt)
1083 : Length(length),
1084 P_out(0.0),
1085 Prescribed_volume_flux_fct_pt(prescribed_volume_flux_fct_pt),
1087 {
1088 // Initialise currently prescribed flux
1089 Current_volume_flux_pt = new double(0.0);
1090
1091 // Auxiliary integral isn't used if flux isn't prescribed
1092 // via outflow through NavierStokesImpedanceTractionElements
1093 Aux_integral_pt = 0;
1094 }
1095
1096
1097 /// Constructor: Specify length of tube and the pointer to the
1098 /// mesh of either NavierStokesImpedanceTractionElements or
1099 /// NavierStokesFluxControlElements that are attached
1100 /// to the outflow cross-section of a (higher-dimensional)
1101 /// Navier Stokes mesh and provide the inflow into the ImpedanceTube.
1102 /// Outlet pressure is set to zero.
1103 WomersleyImpedanceTubeBase(const double& length,
1104 Mesh* navier_stokes_outflow_mesh_pt)
1105 : Length(length),
1106 P_out(0.0),
1108 Navier_stokes_outflow_mesh_pt(navier_stokes_outflow_mesh_pt)
1109 {
1110 // Initialise currently prescribed flux
1111 Current_volume_flux_pt = new double(0.0);
1112
1113 // Initialise flag to record if NavierStokesFluxControlElement
1114 // or NavierStokesImpedanceTractionElement elements are being used
1116
1117 // Attempt to cast 1st element to NavierStokesImpedanceTractionElementBase
1119 navier_stokes_outflow_mesh_pt->element_pt(0)))
1120 {
1122
1123 // Create map used to store the non-zero entries of the
1124 // auxiliary integral, containing the derivative of the total
1125 // volume flux through the outflow boundary of the (higher-dimensional)
1126 // Navier-Stokes mesh w.r.t. to the discrete (global) (velocity)
1127 // degrees of freedom.
1128 Aux_integral_pt = new std::map<unsigned, double>;
1129
1130 // Pass pointer to Navier_stokes_outflow_mesh_pt to the Navier
1131 // Stokes traction elements
1132 unsigned nelem = navier_stokes_outflow_mesh_pt->nelement();
1133 for (unsigned e = 0; e < nelem; e++)
1134 {
1137 navier_stokes_outflow_mesh_pt->element_pt(e));
1138
1139 // Pass the mesh of all NavierStokesImpedanceTractionElements to
1140 // each NavierStokesImpedanceTractionElements in that mesh
1141 // and treat nodes in that mesh that are not part of the element
1142 // itself as external data (since they affect the total volume
1143 // flux and therefore the traction onto the element).
1145 navier_stokes_outflow_mesh_pt);
1146 }
1147 }
1148#ifdef PARANOID
1149 // Test to make sure the elements in the mesh are valid
1150 else
1151 {
1153 navier_stokes_outflow_mesh_pt->element_pt(0)))
1154 {
1155 std::ostringstream error_message;
1156 error_message
1157 << "WomersleyImpedanceTubeBase requires a Navier-Stokes\n"
1158 << "outflow mesh of elements which inherit from either\n"
1159 << "TemplateFreeNavierStokesFluxControlElementBase or\n"
1160 << "NavierStokesImpedanceTractionElementBase.\n";
1161 throw OomphLibError(error_message.str(),
1162 OOMPH_CURRENT_FUNCTION,
1163 OOMPH_EXCEPTION_LOCATION);
1164 }
1165 }
1166#endif
1167 }
1168
1169 /// Access fct to outlet pressure
1170 double& p_out()
1171 {
1172 return P_out;
1173 }
1174
1175 /// Pure virtual function in which the user of a derived class
1176 /// must create the mesh of WomersleyElements (of the type specified
1177 /// by the class's template argument) and apply the boundary conditions.
1178 /// The Womersley elements use the timestepper specified as the
1179 /// input argument.
1181 TimeStepper* time_stepper_pt) = 0;
1182
1183
1184 /// Set up the Womersley tubes so that a subsequent call
1185 /// to get_response(...) computes the inlet pressure for the currently
1186 /// prescribed instantaneous flow rate. Steady version!
1187 void setup()
1188 {
1189 // Dummy parameters
1190 double* re_st_pt = &Zero;
1191 double dt = 0.0;
1192 double q_initial = 0;
1193 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper;
1194 setup(re_st_pt, dt, q_initial, time_stepper_pt);
1195 }
1196
1197
1198 /// Set up the Womersley tubes so that a subsequent call
1199 /// to get_response(...) computes the inlet pressure for the currently
1200 /// prescribed instantaneous flow rate, assuming that at all previous times
1201 /// the tube conveyed steady, fully-developed flow with flowrate q_initial.
1202 /// dt specifies the timestep for the subsequent time integration.
1203 /// Specify: Womersley number, (constant) timestep, the initial volume
1204 /// flux (from which the subsequent impulsive start is performed) and,
1205 /// optionally the pointer to the timestepper to be used in the Womersley
1206 /// elements (defaults to BDF<2>).
1207 void setup(double* re_st_pt,
1208 const double& dt,
1209 const double& q_initial,
1210 TimeStepper* time_stepper_pt = 0)
1211 {
1212 // Create timestepper if none specified so far
1213 if (time_stepper_pt == 0)
1214 {
1215 time_stepper_pt = new BDF<2>;
1216 }
1217
1218 // Build mesh and apply bcs
1219 Mesh* my_mesh_pt =
1221
1222 // Build problem
1224 re_st_pt, Current_volume_flux_pt, time_stepper_pt, my_mesh_pt);
1225
1226 /// By default, we do want to suppress the output from the
1227 /// Newton solver
1228 Womersley_problem_pt->disable_info_in_newton_solve();
1229 oomph_info << "NOTE: We're suppressing timings etc from \n"
1230 << " Newton solver in WomersleyImpedanceTubeBase. "
1231 << std::endl;
1232
1233 // Precompute the auxiliary integrals for the Navier-Stokes
1234 // impedance traction elements (if they're used to specify the inflow
1237 {
1239 }
1240
1241 // Initialise timestep -- also sets the weights for all timesteppers
1242 // in the problem.
1243 Womersley_problem_pt->initialise_dt(dt);
1244
1245 // Set currently imposed flux
1246 *Current_volume_flux_pt = q_initial;
1247
1248 // Assign steady initial solution for this flux
1249 Womersley_problem_pt->steady_newton_solve();
1250
1251 // Allow for resolve
1252 Womersley_problem_pt->linear_solver_pt()->enable_resolve();
1253
1254 // Re-use Jacobian
1255 Womersley_problem_pt->enable_jacobian_reuse();
1256
1257 // Shut up
1258 Womersley_problem_pt->linear_solver_pt()->disable_doc_time();
1259
1260 // Do a dummy solve with time-dependent terms switched on
1261 // to generate (and store) the Jacobian. (We're not using
1262 // a Newton solve because the initial residual may be zero
1263 // in which case the Jacobian would never be computed!)
1264 unsigned n_dof = Womersley_problem_pt->ndof();
1265
1266 // Local scope to make sure dx goes out of scope
1267 {
1268 DoubleVector dx;
1269 Womersley_problem_pt->linear_solver_pt()->solve(Womersley_problem_pt,
1270 dx);
1271 }
1272
1273
1274 // Pre-compute derivative of p_in w.r.t. q
1275
1276 // Setup vector of derivatives of residuals & unknowns w.r.t. Q
1278 Womersley_problem_pt->communicator_pt(), n_dof, false);
1279 DoubleVector drdq(&dist, 0.0);
1280 DoubleVector dxdq(&dist, 0.0);
1281
1282 // What's the global equation number of the equation that
1283 // determines the pressure gradient
1284 unsigned g_eqn =
1285 Womersley_problem_pt->pressure_gradient_data_pt()->eqn_number(0);
1286
1287 // Derivative of volume constraint residual w.r.t. prescribed
1288 // instantaenous volume flux (in ImposeFluxForWomersleyElement)
1289 drdq[g_eqn] = -1.0;
1290
1291 // Solve for derivatives of unknowns in Womersley problem, w.r.t.
1292 // instantaenous volume flux (in ImposeFluxForWomersleyElement)
1293 Womersley_problem_pt->linear_solver_pt()->resolve(drdq, dxdq);
1294
1295 // Rate of change of inflow pressure w.r.t to instantaneous
1296 // volume flux
1297 Dp_in_dq = dxdq[g_eqn] * Length;
1298 }
1299
1300
1301 /// Access to underlying Womersley problem
1303 {
1304 return Womersley_problem_pt;
1305 }
1306
1307
1308 /// Shift history values to allow coputation of next timestep.
1309 /// Note: When used with a full Navier-Stokes problem this function
1310 /// must be called in actions_before_implicit_timestep()
1311 void shift_time_values(const double& dt)
1312 {
1313 // Shift the history values in the Womersley problem
1314 Womersley_problem_pt->shift_time_values();
1315
1316 // Advance global time and set current value of dt
1317 Womersley_problem_pt->time_pt()->time() += dt;
1318 Womersley_problem_pt->time_pt()->dt() = dt;
1319
1320 // Find out how many timesteppers there are
1321 unsigned n_time_steppers = Womersley_problem_pt->ntime_stepper();
1322
1323 // Loop over them all and set the weights
1324 for (unsigned i = 0; i < n_time_steppers; i++)
1325 {
1326 Womersley_problem_pt->time_stepper_pt(i)->set_weights();
1327 }
1328 }
1329
1330
1331 /// Compute total current volume flux into the "impedance tube" that
1332 /// provides the flow resistance (flux is either obtained from
1333 /// the function that specifies it externally or by by adding up the flux
1334 /// through all NavierStokesImpedanceTractionElements in
1335 /// the mesh pointed to by the Navier_stokes_outflow_mesh_pt.
1337 {
1339 {
1341 Womersley_problem_pt->time_pt()->time());
1342 }
1343 else
1344 {
1345 unsigned nelem = Navier_stokes_outflow_mesh_pt->nelement();
1346 double flux = 0.0;
1348 {
1349 for (unsigned e = 0; e < nelem; e++)
1350 {
1351 flux +=
1354 ->get_volume_flux();
1355 }
1356 }
1357 else
1358 {
1359 for (unsigned e = 0; e < nelem; e++)
1360 {
1361 flux += dynamic_cast<NavierStokesImpedanceTractionElementBase*>(
1363 ->get_volume_flux();
1364 }
1365 }
1366 return flux;
1367 }
1368 }
1369
1370
1371 /// Compute inlet pressure, p_in, required to achieve the currently
1372 /// imposed, instantaneous volume flux q prescribed by
1373 /// total_volume_flux_into_impedance_tube(), and its
1374 /// derivative, dp_in/dq.
1375 void get_response(double& p_in, double& dp_in_dq)
1376 {
1377 // Set currently imposed flux
1379
1380 // Do a Newton solve to compute the pressure gradient
1381 // required to achieve the imposed instantaneous flow rate
1382 Womersley_problem_pt->newton_solve();
1383
1384 // Compute inflow pressure based on computed pressure gradient,
1385 // the length of tube, and the outlet pressure
1386 p_in =
1387 -Womersley_problem_pt->pressure_gradient_data_pt()->value(0) * Length +
1388 P_out;
1389
1390 // Return pre-computed value for dp_in/dq
1391 dp_in_dq = Dp_in_dq;
1392 }
1393
1394
1395 protected:
1396 /// Precompute auxiliary integrals required for the computation of
1397 /// the Jacobian in the NavierStokesImpedanceTractionElement. Also pass the
1398 /// pointer to the pre-computed integrals to the elements in the
1399 /// Navier_stokes_outflow_mesh_pt so they can refer to it.
1401 {
1402 // Loop over all elements
1403 unsigned nelem = Navier_stokes_outflow_mesh_pt->nelement();
1404 for (unsigned e = 0; e < nelem; e++)
1405 {
1409
1410 // Add the element's contribution
1412
1413 // Tell the elements who's setting their flow resistance
1414 el_pt->set_impedance_tube_pt(this);
1415 }
1416
1417 // Pass pointer to Aux_integral to the elements so they can
1418 // use it in the computation of the Jacobian
1419 for (unsigned e = 0; e < nelem; e++)
1420 {
1424
1425 // Pass pointer to elements
1427 }
1428 }
1429
1430 /// Length of the tube
1431 double Length;
1432
1433 /// Derivative of inflow pressure w.r.t. instantaenous volume flux
1434 /// (Note: Can be pre-computed)
1435 double Dp_in_dq;
1436
1437 /// Pointer to double that specifies the currently imposed
1438 /// instantaneous volume flux into the impedance tube. This is
1439 /// used to communicate with the Womersley elements which require
1440 /// access to the flux via a pointer to a double.
1442
1443 /// Pointer to Womersley problem that determines the
1444 /// pressure gradient along the tube
1446
1447 /// Outlet pressure
1448 double P_out;
1449
1450 /// Pointer to function that specifies the prescribed volume flux
1452
1453 /// Pointer to the mesh of NavierStokesImpedanceTractionElements
1454 /// that are attached to the outflow cross-section of the higher-dimensional
1455 /// Navier Stokes mesh and provide the inflow into the Impedance tube.
1457
1458 /// Pointer to auxiliary integral, containing
1459 /// the derivative of the total volume flux through the
1460 /// outflow boundary of the (higher-dimensional) Navier-Stokes mesh w.r.t.
1461 /// to the discrete (global) (velocity) degrees of freedom.
1462 std::map<unsigned, double>* Aux_integral_pt;
1463
1464 private:
1465 // Flag to record if NavierStokesFluxControlElement
1466 // or NavierStokesImpedanceTractionElement elements are being used
1468 };
1469
1470 /// /////////////////////////////////////////////////////////////////////
1471 /// /////////////////////////////////////////////////////////////////////
1472 /// /////////////////////////////////////////////////////////////////////
1473
1474 // Inline functions:
1475
1476
1477 //======================================================================
1478 /// Define the shape functions and test functions and derivatives
1479 /// w.r.t. global coordinates and return Jacobian of mapping.
1480 ///
1481 /// Galerkin: Test functions = shape functions
1482 //======================================================================
1483 template<unsigned DIM, unsigned NNODE_1D>
1485 const Vector<double>& s,
1486 Shape& psi,
1487 DShape& dpsidx,
1488 Shape& test,
1489 DShape& dtestdx) const
1490 {
1491 // Call the geometrical shape functions and derivatives
1492 double J = this->dshape_eulerian(s, psi, dpsidx);
1493
1494 // Loop over the test functions and derivatives and set them equal to the
1495 // shape functions
1496 for (unsigned i = 0; i < NNODE_1D; i++)
1497 {
1498 test[i] = psi[i];
1499 for (unsigned j = 0; j < DIM; j++)
1500 {
1501 dtestdx(i, j) = dpsidx(i, j);
1502 }
1503 }
1504
1505 // Return the jacobian
1506 return J;
1507 }
1508
1509
1510 //======================================================================
1511 /// Define the shape functions and test functions and derivatives
1512 /// w.r.t. global coordinates and return Jacobian of mapping.
1513 ///
1514 /// Galerkin: Test functions = shape functions
1515 //======================================================================
1516 template<unsigned DIM, unsigned NNODE_1D>
1519 Shape& psi,
1520 DShape& dpsidx,
1521 Shape& test,
1522 DShape& dtestdx) const
1523 {
1524 // Call the geometrical shape functions and derivatives
1525 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
1526
1527 // Set the test functions equal to the shape functions
1528 //(sets internal pointers)
1529 test = psi;
1530 dtestdx = dpsidx;
1531
1532 // Return the jacobian
1533 return J;
1534 }
1535
1536
1537 /// /////////////////////////////////////////////////////////////////////
1538 /// /////////////////////////////////////////////////////////////////////
1539
1540
1541 //=======================================================================
1542 /// Face geometry for the QWomersleyElement elements: The spatial
1543 /// dimension of the face elements is one lower than that of the
1544 /// bulk element but they have the same number of points
1545 /// along their 1D edges.
1546 //=======================================================================
1547 template<unsigned DIM, unsigned NNODE_1D>
1548 class FaceGeometry<QWomersleyElement<DIM, NNODE_1D>>
1549 : public virtual QElement<DIM - 1, NNODE_1D>
1550 {
1551 public:
1552 /// Constructor: Call the constructor for the
1553 /// appropriate lower-dimensional QElement
1554 FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
1555 };
1556
1557
1558 /// /////////////////////////////////////////////////////////////////////
1559 /// /////////////////////////////////////////////////////////////////////
1560 /// /////////////////////////////////////////////////////////////////////
1561
1562
1563 //=======================================================================
1564 /// Face geometry for the 1D QWomersleyElement elements: Point elements
1565 //=======================================================================
1566 template<unsigned NNODE_1D>
1568 : public virtual PointElement
1569 {
1570 public:
1571 /// Constructor: Call the constructor for the
1572 /// appropriate lower-dimensional QElement
1574 };
1575
1576
1577 /// /////////////////////////////////////////////////////////////////////
1578 /// /////////////////////////////////////////////////////////////////////
1579 /// /////////////////////////////////////////////////////////////////////
1580
1581
1582 //====================================================================
1583 /// Template-free base class
1584 //====================================================================
1586 {
1587 public:
1588 /// Static bool to suppress warning
1590 };
1591
1592
1593 /// /////////////////////////////////////////////////////////////////////
1594 /// /////////////////////////////////////////////////////////////////////
1595 /// /////////////////////////////////////////////////////////////////////
1596
1597
1598 //====================================================================
1599 /// Mesh of Womersley elements whose topology, nodal position etc.
1600 /// matches that of a given mesh of face elements in the outflow
1601 /// cross-section of a full Navier-Stokes mesh.
1602 //====================================================================
1603 template<class WOMERSLEY_ELEMENT>
1604 class WomersleyMesh : public virtual Mesh,
1605 public virtual TemplateFreeWomersleyMeshBase
1606 {
1607 public:
1608 /// Constructor: Pass pointer to mesh of face elements in the
1609 /// outflow cross-section of a full Navier-Stokes mesh, the timestepper to
1610 /// be used for the Womersley elements, the coordinate (in the
1611 /// higher-dimensional Navier-Stokes mesh) that is constant
1612 /// in the outflow cross-section and the velocity component
1613 /// (in terms of the nodal index) that represents the outflow
1614 /// component -- the latter is used to automatically apply
1615 /// the boundary conditions for the Womersley problem.
1616 WomersleyMesh(Mesh* n_st_outflow_mesh_pt,
1617 TimeStepper* time_stepper_pt,
1618 const unsigned& fixed_coordinate,
1619 const unsigned& w_index)
1620 {
1621 /// How many elements and nodes are there in the original mesh?
1622 unsigned nelem = n_st_outflow_mesh_pt->nelement();
1623
1624 // Navier-Stokes outflow mesh may not have any nodes stored (it usually
1625 // just acts as a container for traction elements) -->
1626 // Count number of distinct Navier-Stokes nodes by adding
1627 // the elements' nodes to a set
1628 std::set<Node*> n_st_nodes;
1629 for (unsigned e = 0; e < nelem; e++)
1630 {
1631 FiniteElement* el_pt = n_st_outflow_mesh_pt->finite_element_pt(e);
1632 unsigned nnod_el = el_pt->nnode();
1633 for (unsigned j = 0; j < nnod_el; j++)
1634 {
1635 n_st_nodes.insert(el_pt->node_pt(j));
1636
1637 // Careful: It there are hanging nodes this won't work!
1638 if (el_pt->node_pt(j)->is_hanging())
1639 {
1640 throw OomphLibError(
1641 "Cannot build WomersleyMesh from mesh with hanging nodes!",
1642 OOMPH_CURRENT_FUNCTION,
1643 OOMPH_EXCEPTION_LOCATION);
1644 }
1645 }
1646 }
1647
1648 // Extract size then wipe
1649 unsigned nnode_n_st = n_st_nodes.size();
1650 n_st_nodes.clear();
1651
1652 // Create enough storage
1653 Node_pt.resize(nnode_n_st);
1654
1655 /// Create new elements
1656 for (unsigned e = 0; e < nelem; e++)
1657 {
1658 add_element_pt(new WOMERSLEY_ELEMENT);
1659#ifdef PARANOID
1660 if (finite_element_pt(e)->nnode() !=
1661 n_st_outflow_mesh_pt->finite_element_pt(e)->nnode())
1662 {
1663 throw OomphLibError(
1664 "Number of nodes in existing and new elements don't match",
1665 OOMPH_CURRENT_FUNCTION,
1666 OOMPH_EXCEPTION_LOCATION);
1667 }
1668#endif
1669 }
1670
1671 // Map to record which Navier-Stokes nodes have been visited (default
1672 // return is false)
1673 std::map<Node*, bool> n_st_node_done;
1674
1675 // Map to store the Womersley node that corresponds to a
1676 // Navier Stokes node
1677 std::map<Node*, Node*> equivalent_womersley_node_pt;
1678
1679 // Initialise count of newly created nodes
1680 unsigned node_count = 0;
1681
1682
1683 // This is awkward do diagnose: We're assuming that
1684 // the boundary conditions have been applied for the
1685 // underlying Navier-Stokes problem before calling
1686 // this function (otherwise it's really tricky to
1687 // apply the right boundary conditions here), but it's
1688 // hard to police. Issue definite (but suppressable)
1689 // warning if nothing has been pinned at all.
1690 unsigned n_pinned_nodes = 0;
1691
1692 // Loop over nst and womersley elements in tandem to sort out
1693 // which new nodes are required
1694 for (unsigned e = 0; e < nelem; e++)
1695 {
1696 FiniteElement* n_st_el_pt = n_st_outflow_mesh_pt->finite_element_pt(e);
1697 unsigned nnod_el = n_st_el_pt->nnode();
1698 for (unsigned j = 0; j < nnod_el; j++)
1699 {
1700 // Has the Navier Stokes node been done yet?
1701 Node* n_st_node_pt = n_st_el_pt->node_pt(j);
1702
1703 // Hasn't been done: Create new node in Womersley element
1704 if (!n_st_node_done[n_st_node_pt])
1705 {
1706 // Create a new node in the Womersley element
1707 Node* new_node_pt =
1708 finite_element_pt(e)->construct_node(j, time_stepper_pt);
1709
1710 // Add newly created node
1711 Node_pt[node_count] = new_node_pt;
1712 node_count++;
1713
1714
1715 // Set coordinates
1716 unsigned dim = n_st_node_pt->ndim();
1717 unsigned icount = 0;
1718 for (unsigned i = 0; i < dim; i++)
1719 {
1720 if (i != fixed_coordinate)
1721 {
1722 new_node_pt->x(icount) = n_st_node_pt->x(i);
1723 icount++;
1724 }
1725 }
1726
1727 // Set pin status
1728 if (n_st_node_pt->is_pinned(w_index))
1729 {
1730 new_node_pt->pin(0);
1731 n_pinned_nodes++;
1732 }
1733 else
1734 {
1735 // shouldn't need this, but...
1736 new_node_pt->unpin(0);
1737 }
1738
1739 // Record which Womersley node the
1740 // Navier Stokes node is associated with
1741 equivalent_womersley_node_pt[n_st_node_pt] = new_node_pt;
1742
1743 // Now the Navier-Stokes node has been done
1744 n_st_node_done[n_st_node_pt] = true;
1745 }
1746 // The node has already been done -- set pointer to existing
1747 // Womersley node
1748 else
1749 {
1751 equivalent_womersley_node_pt[n_st_node_pt];
1752 }
1753 }
1754
1755 bool passed = true;
1757 if (!passed)
1758 {
1759 // Reverse the nodes
1760 unsigned nnod = finite_element_pt(e)->nnode();
1761 Vector<Node*> orig_nod_pt(nnod);
1762 for (unsigned j = 0; j < nnod; j++)
1763 {
1764 orig_nod_pt[j] = finite_element_pt(e)->node_pt(j);
1765 }
1766 for (unsigned j = 0; j < nnod; j++)
1767 {
1768 finite_element_pt(e)->node_pt(j) = orig_nod_pt[nnod - j - 1];
1769 }
1770 bool passed = true;
1772 if (!passed)
1773 {
1774 throw OomphLibError("Element remains inverted even after reversing "
1775 "the local node numbers",
1776 OOMPH_CURRENT_FUNCTION,
1777 OOMPH_EXCEPTION_LOCATION);
1778 }
1779 }
1780 }
1781
1782
1783#ifdef PARANOID
1785 {
1786 if (n_pinned_nodes == 0)
1787 {
1788 std::stringstream bla;
1789 bla << "Boundary conditions must be applied in Navier-Stokes\n"
1790 << "problem before attaching impedance elements.\n"
1791 << "Note: This warning can be suppressed by setting the\n"
1792 << "global static boolean\n\n"
1793 << " "
1794 "TemplateFreeWomersleyMeshBase::Suppress_warning_about_"
1795 "unpinned_nst_dofs\n\n"
1796 << "to true\n";
1798 bla.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1799 }
1800 }
1801#endif
1802
1803#ifdef PARANOID
1804 if (nnode() != nnode_n_st)
1805 {
1806 throw OomphLibError(
1807 "Number of nodes in the new mesh don't match that in the old one",
1808 OOMPH_CURRENT_FUNCTION,
1809 OOMPH_EXCEPTION_LOCATION);
1810 }
1811#endif
1812 }
1813 };
1814
1815
1816 /// /////////////////////////////////////////////////////////////////////
1817 /// //////////////////////////////////////////////////////////////////////
1818 /// //////////////////////////////////////////////////////////////////////
1819
1820
1821 //====================================================================
1822 /// WomersleyImpedanceTube that attaches itself to the outflow
1823 /// of a Navier-Stokes mesh.
1824 //====================================================================
1825 template<class ELEMENT, unsigned DIM>
1827 : public WomersleyImpedanceTubeBase<ELEMENT, DIM>
1828 {
1829 public:
1830 /// Constructor: Pass length and mesh of face elements that
1831 /// are attached to the outflow cross-section of the Navier Stokes mesh
1832 /// to constructor of underlying base class. Also specify
1833 /// the coordinate (in the higher-dimensional
1834 /// Navier-Stokes mesh) that is constant
1835 /// in the outflow cross-section and the velocity component
1836 /// (in terms of the nodal index) that represents the outflow
1837 /// component -- the latter is used to automatically apply
1838 /// the boundary conditions for the Womersley problem.
1839 WomersleyOutflowImpedanceTube(const double& length,
1840 Mesh* navier_stokes_outflow_mesh_pt,
1841 const unsigned& fixed_coordinate,
1842 const unsigned& w_index)
1843 : WomersleyImpedanceTubeBase<ELEMENT, DIM>(length,
1844 navier_stokes_outflow_mesh_pt),
1845 Fixed_coordinate(fixed_coordinate),
1846 W_index(w_index)
1847 {
1848 }
1849
1850 /// Implement pure virtual fct (defined in the base class
1851 /// WomersleyImpedanceTubeBase) that builds the mesh of Womersley elements
1852 /// (of the type specified by the template argument), using the
1853 /// specified timestepper. Also applies the boundary condition.
1855 {
1856 // Build mesh and automatically apply the same boundary
1857 // conditions as those that are applied to the W_index-th
1858 // value of the nodes in the Navier-Stokes mesh.
1859 WomersleyMesh<ELEMENT>* womersley_mesh_pt =
1861 time_stepper_pt,
1863 W_index);
1864
1865 return womersley_mesh_pt;
1866 }
1867
1868 private:
1869 /// The coordinate (in the higher-dimensional Navier-Stokes
1870 /// mesh) that is constant in the outflow cross-section
1872
1873 /// The velocity component
1874 /// (in terms of the nodal index) that represents the outflow
1875 /// component -- the latter is used to automatically apply
1876 /// the boundary conditions for the Womersley problem.
1877 unsigned W_index;
1878 };
1879
1880 /// /////////////////////////////////////////////////////////////////////
1881 /// //////////////////////////////////////////////////////////////////////
1882 /// //////////////////////////////////////////////////////////////////////
1883
1884
1885 /* //==================================================================== */
1886 /* /// WomersleyImpedanceTube that attaches itself to the outflow */
1887 /* /// of a Navier-Stokes mesh for use with . */
1888 /* //==================================================================== */
1889 /* template<class ELEMENT, unsigned DIM> */
1890 /* class WomersleyOutflowImpedanceTubeForNavierStokesBlockPreconditioner : */
1891 /* public WomersleyImpedanceTubeBaseForNavierStokesBlockPreconditioner */
1892 /* <ELEMENT,DIM> */
1893 /* { */
1894
1895 /* public: */
1896
1897 /* /// Constructor: Pass length and mesh of face elements that */
1898 /* /// are attached to the outflow cross-section of the Navier Stokes mesh */
1899 /* /// to constructor of underlying base class. Also specify */
1900 /* /// the coordinate (in the higher-dimensional */
1901 /* /// Navier-Stokes mesh) that is constant */
1902 /* /// in the outflow cross-section and the velocity component */
1903 /* /// (in terms of the nodal index) that represents the outflow */
1904 /* /// component -- the latter is used to automatically apply */
1905 /* /// the boundary conditions for the Womersley problem. */
1906 /* WomersleyOutflowImpedanceTubeForNavierStokesBlockPreconditioner( */
1907 /* const double& length, */
1908 /* Mesh* navier_stokes_outflow_mesh_pt, */
1909 /* const unsigned& fixed_coordinate, */
1910 /* const unsigned& w_index) : */
1911 /* WomersleyImpedanceTubeBaseForNavierStokesBlockPreconditioner<ELEMENT,DIM>
1912 */
1913 /* (length, */
1914 /* navier_stokes_outflow_mesh_pt), */
1915 /* Fixed_coordinate(fixed_coordinate), W_index(w_index) */
1916 /* {} */
1917
1918 /* /// Implement pure virtual fct (defined in the base class */
1919 /* /// WomersleyImpedanceTubeBase) that builds the mesh of Womersley elements
1920 */
1921 /* /// (of the type specified by the template argument), using the */
1922 /* /// specified timestepper. Also applies the boundary condition. */
1923 /* Mesh* build_mesh_and_apply_boundary_conditions(TimeStepper*
1924 * time_stepper_pt) */
1925 /* { */
1926 /* // Build mesh and automatically apply the same boundary */
1927 /* // conditions as those that are applied to the W_index-th */
1928 /* // value of the nodes in the Navier-Stokes mesh. */
1929 /* WomersleyMesh<ELEMENT>* womersley_mesh_pt= // CHANGED TO USE
1930 * ELEMENT */
1931 /* new WomersleyMesh<ELEMENT>( */
1932 /* this->Navier_stokes_outflow_mesh_pt, */
1933 /* time_stepper_pt, */
1934 /* Fixed_coordinate, */
1935 /* W_index); */
1936
1937 /* return womersley_mesh_pt; */
1938
1939 /* } */
1940
1941 /* private: */
1942
1943 /* /// The coordinate (in the higher-dimensional Navier-Stokes */
1944 /* /// mesh) that is constant in the outflow cross-section */
1945 /* unsigned Fixed_coordinate; */
1946
1947 /* /// The velocity component */
1948 /* /// (in terms of the nodal index) that represents the outflow */
1949 /* /// component -- the latter is used to automatically apply */
1950 /* /// the boundary conditions for the Womersley problem. */
1951 /* unsigned W_index; */
1952
1953
1954 /* }; */
1955
1956
1957 /// //////////////////////////////////////////////////////////////////////
1958 /// //////////////////////////////////////////////////////////////////////
1959 /// //////////////////////////////////////////////////////////////////////
1960
1961
1962 //======================================================================
1963 /// A class for elements that allow the imposition of an impedance type
1964 /// traction boundary condition to the Navier--Stokes equations
1965 /// The geometrical information can be read from the FaceGeometery<ELEMENT>
1966 /// class and and thus, we can be generic enough without the need to have
1967 /// a separate equations class. Template arguments specify the
1968 /// type of the bulk Navier Stokes elements that the elements are
1969 /// attached to, and the type of the Womersley element used
1970 /// to compute the flow resistance in the downstream "impedance tube".
1971 //======================================================================
1972 template<class BULK_NAVIER_STOKES_ELEMENT,
1973 class WOMERSLEY_ELEMENT,
1974 unsigned DIM>
1976 : public virtual FaceGeometry<BULK_NAVIER_STOKES_ELEMENT>,
1977 public virtual FaceElement,
1979 {
1980 private:
1981 /// Pointer to auxiliary integral, containing
1982 /// the derivative of the total volume flux through the
1983 /// outflow boundary of the (higher-dimensional) Navier-Stokes mesh w.r.t.
1984 /// to the discrete (global) (velocity) degrees of freedom.
1985 std::map<unsigned, double>* Aux_integral_pt;
1986
1987 /// Pointer to ImpedanceTubeProblem that computes the flow resistance
1989
1990
1991 protected:
1992 /// Access function that returns the local equation numbers
1993 /// for velocity components.
1994 /// u_local_eqn(n,i) = local equation number or < 0 if pinned.
1995 /// The default is to asssume that n is the local node number
1996 /// and the i-th velocity component is the i-th unknown stored at the node.
1997 virtual inline int u_local_eqn(const unsigned& n, const unsigned& i)
1998 {
1999 return nodal_local_eqn(n, i);
2000 }
2001
2002 /// Function to compute the shape and test functions and to return
2003 /// the Jacobian of mapping
2004 inline double shape_and_test_at_knot(const unsigned& ipt,
2005 Shape& psi,
2006 Shape& test) const
2007 {
2008 // Find number of nodes
2009 unsigned n_node = nnode();
2010
2011 // Calculate the shape functions
2012 shape_at_knot(ipt, psi);
2013
2014 // Set the test functions to be the same as the shape functions
2015 for (unsigned i = 0; i < n_node; i++)
2016 {
2017 test[i] = psi[i];
2018 }
2019
2020 // Return the value of the jacobian
2021 return J_eulerian_at_knot(ipt);
2022 }
2023
2024
2025 /// This function returns the residuals for the
2026 /// traction function.
2027 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
2029 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
2030
2031
2032 public:
2033 /// Constructor, which takes a "bulk" element and the value of the index
2034 /// and its limit
2036 const int& face_index)
2037 : FaceGeometry<BULK_NAVIER_STOKES_ELEMENT>(), FaceElement()
2038 {
2039 // Attach the geometrical information to the element. N.B. This function
2040 // also assigns nbulk_value from the required_nvalue of the bulk element
2041 element_pt->build_face_element(face_index, this);
2042
2043 // Initialise pointer to mesh containing the
2044 // NavierStokesImpedanceTractionElements
2045 // that contribute to the volume flux into the "downstream tube" that
2046 // provides the impedance
2048
2049 // Initialise pointer to impedance tube
2051
2052 // Initialise pointer to auxiliary integral
2053 Aux_integral_pt = 0;
2054
2055 // Set the dimension from the dimension of the first node
2056 // Dim = this->node_pt(0)->ndim();
2057
2058#ifdef PARANOID
2059 {
2060 // Check that the element is not a refineable 3d element
2061 BULK_NAVIER_STOKES_ELEMENT* elem_pt =
2062 dynamic_cast<BULK_NAVIER_STOKES_ELEMENT*>(element_pt);
2063 // If it's three-d
2064 if (elem_pt->dim() == 3)
2065 {
2066 // Is it refineable
2067 RefineableElement* ref_el_pt =
2068 dynamic_cast<RefineableElement*>(elem_pt);
2069 if (ref_el_pt != 0)
2070 {
2071 if (this->has_hanging_nodes())
2072 {
2073 throw OomphLibError("This flux element will not work correctly "
2074 "if nodes are hanging\n",
2075 OOMPH_CURRENT_FUNCTION,
2076 OOMPH_EXCEPTION_LOCATION);
2077 }
2078 }
2079 }
2080 }
2081#endif
2082 }
2083
2084
2085 /// Destructor should not delete anything
2087
2088 /// Access to mesh containing all
2089 /// NavierStokesImpedanceTractionElements that contribute to the volume flux
2090 /// into the "downstream tube" that provides the impedance
2092 {
2094 }
2095
2096 /// Get integral of volume flux through element
2098 {
2099 // Initialise
2100 double volume_flux_integral = 0.0;
2101
2102 // Vector of local coordinates in face element
2103 Vector<double> s(DIM);
2104
2105 // Vector for global Eulerian coordinates
2106 Vector<double> x(DIM + 1);
2107
2108 // Vector for local coordinates in bulk element
2109 Vector<double> s_bulk(DIM + 1);
2110
2111 // Set the value of n_intpt
2112 unsigned n_intpt = integral_pt()->nweight();
2113
2114 // Get pointer to assocated bulk element
2115 BULK_NAVIER_STOKES_ELEMENT* bulk_el_pt =
2116 dynamic_cast<BULK_NAVIER_STOKES_ELEMENT*>(bulk_element_pt());
2117
2118 // Loop over the integration points
2119 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2120 {
2121 // Assign values of s in FaceElement and local coordinates in bulk
2122 // element
2123 for (unsigned i = 0; i < DIM; i++)
2124 {
2125 s[i] = integral_pt()->knot(ipt, i);
2126 }
2127
2128 // Get the bulk coordinates
2129 this->get_local_coordinate_in_bulk(s, s_bulk);
2130
2131 // Get the integral weight
2132 double w = integral_pt()->weight(ipt);
2133
2134 // Get jacobian of mapping
2135 double J = J_eulerian(s);
2136
2137 // Premultiply the weights and the Jacobian
2138 double W = w * J;
2139
2140
2141#ifdef PARANOID
2142
2143 // Get x position as Vector
2144 interpolated_x(s, x);
2145
2146 // Get x position as Vector from bulk element
2147 Vector<double> x_bulk(DIM + 1);
2148 bulk_el_pt->interpolated_x(s_bulk, x_bulk);
2149
2150 double max_legal_error = 1.0e-12;
2151 double error = 0.0;
2152 for (unsigned i = 0; i < DIM + 1; i++)
2153 {
2154 error += std::fabs(x[i] - x_bulk[i]);
2155 }
2156 if (error > max_legal_error)
2157 {
2158 std::ostringstream error_stream;
2159 error_stream << "difference in Eulerian posn from bulk and face: "
2160 << error << " exceeds threshold " << max_legal_error
2161 << std::endl;
2162 throw OomphLibError(error_stream.str(),
2163 OOMPH_CURRENT_FUNCTION,
2164 OOMPH_EXCEPTION_LOCATION);
2165 }
2166#endif
2167
2168 // Outer unit normal
2169 Vector<double> normal(DIM + 1);
2170 outer_unit_normal(s, normal);
2171
2172 // Get velocity from bulk element
2173 Vector<double> veloc(DIM + 1);
2174 bulk_el_pt->interpolated_u_nst(s_bulk, veloc);
2175
2176 // Volume flux
2177 double volume_flux = 0.0;
2178 for (unsigned i = 0; i < DIM + 1; i++)
2179 {
2180 volume_flux += normal[i] * veloc[i];
2181 }
2182
2183 // Add to integral
2184 volume_flux_integral += volume_flux * W;
2185 }
2186
2187 return volume_flux_integral;
2188 }
2189
2190
2191 /// NavierStokesImpedanceTractionElements that contribute
2192 /// to the volume flux into the downstream "impedance tube"
2193 /// to the element and classify all nodes in that mesh
2194 /// as external Data for this element (unless the nodes
2195 /// are also the element's own nodes, of course).
2198 {
2199 // Store pointer to mesh of NavierStokesImpedanceTractionElement
2200 // that contribute to the volume flux into the "impedance tube" that
2201 // provides the flow resistance
2203
2204 // Create a set the contains all nodal Data in the flux mesh
2205 std::set<Data*> external_data_set;
2206 unsigned nelem = Navier_stokes_outflow_mesh_pt->nelement();
2207 for (unsigned e = 0; e < nelem; e++)
2208 {
2209 FiniteElement* el_pt =
2211 unsigned nnod = el_pt->nnode();
2212 for (unsigned j = 0; j < nnod; j++)
2213 {
2214 external_data_set.insert(el_pt->node_pt(j));
2215 }
2216 }
2217
2218 // Remove the element's own nodes
2219 unsigned nnod = nnode();
2220 for (unsigned j = 0; j < nnod; j++)
2221 {
2222 external_data_set.erase(node_pt(j));
2223 }
2224
2225 // Copy across
2226 for (std::set<Data*>::iterator it = external_data_set.begin();
2227 it != external_data_set.end();
2228 it++)
2229 {
2230 add_external_data(*it);
2231 }
2232 }
2233
2234
2235 /// Set pointer to the precomputed auxiliary integral that contains
2236 /// the derivative of the total volume flux through the
2237 /// outflow boundary of the (higher-dimensional) Navier-Stokes mesh w.r.t.
2238 /// to the discrete (global) (velocity) degrees of freedom.
2239 void set_aux_integral_pt(std::map<unsigned, double>* aux_integral_pt)
2240 {
2241 Aux_integral_pt = aux_integral_pt;
2242 }
2243
2244
2245 /// Compute total volume flux into the "downstream tube" that
2246 /// provides the impedance (computed by adding up the flux
2247 /// through all NavierStokesImpedanceTractionElements in
2248 /// the mesh specified by volume_flux_mesh_pt().
2250 {
2251#ifdef PARANOID
2253 {
2254 throw OomphLibError(
2255 "Navier_stokes_outflow_mesh_pt==0 -- set it with \n "
2256 "set_external_data_from_navier_stokes_outflow_mesh() before calling "
2257 "this function!\n",
2258 OOMPH_CURRENT_FUNCTION,
2259 OOMPH_EXCEPTION_LOCATION);
2260 }
2261#endif
2262
2263
2264 double total_flux = 0.0;
2265 unsigned nelem = Navier_stokes_outflow_mesh_pt->nelement();
2266 for (unsigned e = 0; e < nelem; e++)
2267 {
2268 NavierStokesImpedanceTractionElement<BULK_NAVIER_STOKES_ELEMENT,
2269 WOMERSLEY_ELEMENT,
2270 DIM>* el_pt =
2271 dynamic_cast<
2272 NavierStokesImpedanceTractionElement<BULK_NAVIER_STOKES_ELEMENT,
2273 WOMERSLEY_ELEMENT,
2274 DIM>*>(
2276 total_flux += el_pt->get_volume_flux();
2277 }
2278 return total_flux;
2279 }
2280
2281
2282 /// Set pointer to "impedance tube" that provides the flow
2283 /// resistance
2285 TemplateFreeWomersleyImpedanceTubeBase* impedance_tube_pt)
2286 {
2289 impedance_tube_pt);
2290 }
2291
2292
2293 /// Add the element's contribution to the auxiliary integral
2294 /// that contains the derivative of the total volume flux through the
2295 /// outflow boundary of the (higher-dimensional) Navier-Stokes mesh w.r.t.
2296 /// to the discrete (global) (velocity) degrees of freedom.
2298 std::map<unsigned, double>* aux_integral_pt)
2299 {
2300 // Spatial dimension of element
2301 // unsigned ndim=dim();
2302
2303 // Vector of local coordinates in face element
2304 Vector<double> s(DIM);
2305
2306 // Create storage for shape functions
2307 unsigned nnod = nnode();
2308 Shape psi(nnod);
2309
2310 // Set the value of n_intpt
2311 unsigned n_intpt = integral_pt()->nweight();
2312
2313 // Loop over the integration points
2314 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2315 {
2316 // Assign values of s in FaceElement and local coordinates in bulk
2317 // element
2318 for (unsigned i = 0; i < DIM; i++)
2319 {
2320 s[i] = integral_pt()->knot(ipt, i);
2321 }
2322
2323 // Get the integral weight
2324 double w = integral_pt()->weight(ipt);
2325
2326 // Get jacobian of mapping
2327 double J = J_eulerian(s);
2328
2329 // Get shape functions
2330 shape(s, psi);
2331
2332 // Premultiply the weights and the Jacobian
2333 double W = w * J;
2334
2335 // Outer unit normal
2336 Vector<double> normal(DIM + 1);
2337 outer_unit_normal(s, normal);
2338
2339 // Loop over nodes
2340 for (unsigned j = 0; j < nnod; j++)
2341 {
2342 // Get pointer to Node
2343 Node* nod_pt = node_pt(j);
2344
2345 // Loop over directions
2346 for (unsigned i = 0; i < (DIM + 1); i++)
2347 {
2348 // Get global equation number
2349 int i_global = nod_pt->eqn_number(i);
2350
2351 // Real dof or bc?
2352 if (i_global >= 0)
2353 {
2354 (*aux_integral_pt)[i_global] += psi[j] * normal[i] * W;
2355 }
2356 }
2357 }
2358 }
2359 }
2360
2361
2362 /// Fill in the element's contribution to the element's residual vector
2364 {
2365 // Call the generic residuals function with flag set to 0
2366 // using a dummy matrix argument
2368 residuals, GeneralisedElement::Dummy_matrix, 0);
2369 }
2370
2371
2372 /// Fill in the element's contribution to the element's residual
2373 /// vector and Jacobian matrix
2375 DenseMatrix<double>& jacobian)
2376 {
2377 // Call the generic routine with the flag set to 1
2379 residuals, jacobian, 1);
2380 }
2381
2382
2383 /// Specify the value of nodal zeta from the face geometry
2384 /// The "global" intrinsic coordinate of the element when
2385 /// viewed as part of a geometric object should be given by
2386 /// the FaceElement representation, by default (needed to break
2387 /// indeterminacy if bulk element is SolidElement)
2388 double zeta_nodal(const unsigned& n,
2389 const unsigned& k,
2390 const unsigned& i) const
2391 {
2392 return FaceElement::zeta_nodal(n, k, i);
2393 }
2394
2395
2396 /// Overload the output function
2397 void output(std::ostream& outfile)
2398 {
2399 FiniteElement::output(outfile);
2400 }
2401
2402 /// Output function: x,y,[z],u,v,[w],p in tecplot format
2403 void output(std::ostream& outfile, const unsigned& nplot)
2404 {
2405 FiniteElement::output(outfile, nplot);
2406 }
2407 };
2408
2409
2410 /// ////////////////////////////////////////////////////////////////////
2411 /// ////////////////////////////////////////////////////////////////////
2412 /// ////////////////////////////////////////////////////////////////////
2413
2414
2415 //============================================================================
2416 /// Function that returns the residuals for the imposed traction Navier_Stokes
2417 /// equations
2418 //============================================================================
2419 template<class BULK_NAVIER_STOKES_ELEMENT,
2420 class WOMERSLEY_ELEMENT,
2421 unsigned DIM>
2422 void NavierStokesImpedanceTractionElement<BULK_NAVIER_STOKES_ELEMENT,
2423 WOMERSLEY_ELEMENT,
2424 DIM>::
2425 fill_in_generic_residual_contribution_fluid_traction(
2426 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
2427 {
2428 // Find out how many nodes there are
2429 unsigned n_node = nnode();
2430
2431 // Set up memory for the shape and test functions
2432 Shape psif(n_node), testf(n_node);
2433
2434 // Set the value of n_intpt
2435 unsigned n_intpt = integral_pt()->nweight();
2436
2437 // Integers to store local equation numbers
2438 int local_eqn = 0;
2439 int local_unknown = 0;
2440
2441 // Loop over the integration points
2442 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2443 {
2444 // Get the integral weight
2445 double w = integral_pt()->weight(ipt);
2446
2447 // Find the shape and test functions and return the Jacobian
2448 // of the mapping
2449 double J = shape_and_test_at_knot(ipt, psif, testf);
2450
2451 // Premultiply the weights and the Jacobian
2452 double W = w * J;
2453
2454 // Traction vector
2455 Vector<double> traction(DIM + 1);
2456
2457 // Initialise response
2458 double p_in = 0.0;
2459 double dp_in_dq = 0.0;
2460
2461 // Traction= outer unit normal x pressure at upstream end of
2462 // impedance tube
2463 if (Navier_stokes_outflow_mesh_pt != 0)
2464 {
2465 // Get response of the impedance tube:
2466 Impedance_tube_pt->get_response(p_in, dp_in_dq);
2467 }
2468
2469 // Get outer unit normal at current integration point
2470 Vector<double> unit_normal(DIM + 1);
2471 outer_unit_normal(ipt, unit_normal);
2472
2473 // Loop over the directions
2474 for (unsigned i = 0; i < DIM + 1; i++)
2475 {
2476 traction[i] = -unit_normal[i] * p_in;
2477 }
2478
2479
2480 // Loop over the test functions
2481 for (unsigned l = 0; l < n_node; l++)
2482 {
2483 // Loop over the velocity components
2484 for (unsigned i = 0; i < DIM + 1; i++)
2485 {
2486 local_eqn = u_local_eqn(l, i);
2487 /*IF it's not a boundary condition*/
2488 if (local_eqn >= 0)
2489 {
2490 // Add the user-defined traction terms
2491 residuals[local_eqn] += traction[i] * testf[l] * W;
2492
2493 // Compute the Jacobian too?
2494 if (flag && (Navier_stokes_outflow_mesh_pt != 0))
2495 {
2496 // Loop over the nodes
2497 for (unsigned j = 0; j < n_node; j++)
2498 {
2499 // Get pointer to Node
2500 Node* nod_pt = node_pt(j);
2501
2502 // Loop over the velocity components
2503 for (unsigned ii = 0; ii < DIM + 1; ii++)
2504 {
2505 local_unknown = u_local_eqn(j, ii);
2506
2507 /*IF it's not a boundary condition*/
2508 if (local_unknown >= 0)
2509 {
2510 // Get corresponding global unknown number
2511 unsigned global_unknown = nod_pt->eqn_number(ii);
2512
2513 // Add contribution
2514 jacobian(local_eqn, local_unknown) -=
2515 (*Aux_integral_pt)[global_unknown] * psif[l] *
2516 unit_normal[i] * dp_in_dq * W;
2517 }
2518 }
2519 }
2520
2521
2522 // Loop over external dofs for unknowns
2523 unsigned n_ext = nexternal_data();
2524 for (unsigned j = 0; j < n_ext; j++)
2525 {
2526 // Get pointer to external Data (=other nodes)
2527 Data* ext_data_pt = external_data_pt(j);
2528
2529 // Loop over directions for equation
2530 for (unsigned ii = 0; ii < DIM + 1; ii++)
2531 {
2532 // Get local unknown number
2533 int local_unknown = external_local_eqn(j, ii);
2534
2535 // Real dof or bc?
2536 if (local_unknown >= 0)
2537 {
2538 // Get corresponding global unknown number
2539 unsigned global_unknown = ext_data_pt->eqn_number(ii);
2540
2541 // Add contribution
2542 jacobian(local_eqn, local_unknown) -=
2543 (*Aux_integral_pt)[global_unknown] * psif[l] *
2544 unit_normal[i] * dp_in_dq * W;
2545 }
2546 }
2547 }
2548 } // end of computation of Jacobian terms
2549 }
2550 } // End of loop over dimension
2551 } // End of loop over shape functions
2552 }
2553 }
2554
2555 /// ///////////////////////////////////////////////////////////////////////
2556 /// ///////////////////////////////////////////////////////////////////////
2557 /// ///////////////////////////////////////////////////////////////////////
2558
2559
2560 //======================================================================
2561 /// An element to impose a fluid pressure obtained from a Womersley
2562 /// impedance tube at a boundary. This element is used in conjunction with a
2563 /// NetFluxControlElementForWomersleyPressureControl element, and is
2564 /// passed to the NetFluxControlElementForWomersleyPressureControl element's
2565 /// constructor. The volume flux across the boundary is then an
2566 /// unknown of the problem. The constructor argument for this element
2567 /// is a suitable Womersley impedance tube to give the pressure via
2568 /// its get_response(...) function.
2569 ///
2570 /// Note: the NavierStokesWomersleyPressureControlElement element calculates
2571 /// Jacobian entries BOTH for itself AND for the
2572 /// NetFluxControlElementForWomersleyPressureControl with respect to
2573 /// the unknowns in this (NavierStokesWomersleyPressureControlElement)
2574 /// element.
2575 //======================================================================
2577 : public virtual GeneralisedElement
2578 {
2579 public:
2580 /// Constructor takes a pointer to a suitable Womersley
2581 /// impedance tube which defines the pressure via get_response(...)
2583 TemplateFreeWomersleyImpedanceTubeBase* womersley_tube_pt)
2584 : Womersley_tube_pt(womersley_tube_pt)
2585 {
2586 // Create the new Data which contains the volume flux.
2587 Volume_flux_data_pt = new Data(1);
2588
2589 // Add new Data to internal data
2591 }
2592
2593 /// Destructor should not delete anything
2595
2596 /// This function returns the residuals
2598 {
2599 // Call the generic residuals function using a dummy matrix argument
2601 residuals, GeneralisedElement::Dummy_matrix, 0);
2602 }
2603
2604 /// This function returns the residuals and the Jacobian,
2605 /// plus the Jacobian contribution for the
2606 /// NetFluxControlElementForWomersleyPressureControl
2607 /// with respect to unknowns in this element
2609 DenseMatrix<double>& jacobian)
2610 {
2611 // Call the generic routine
2613 residuals, jacobian, 1);
2614 }
2615
2616 /// Function to return a pointer to the Data object whose
2617 /// single value is the flux degree of freedom
2619 {
2620 return Volume_flux_data_pt;
2621 }
2622
2623 /// Function to add to external data the Data object whose
2624 /// single value is the pressure applied at the boundary
2625 void add_pressure_data(Data* pressure_data_pt)
2626 {
2627 Pressure_data_id = add_external_data(pressure_data_pt);
2628 }
2629
2630 /// The number of "DOF types" that degrees of freedom in this element
2631 /// are sub-divided into - set to 1
2632 unsigned ndof_types() const
2633 {
2634 return 1;
2635 }
2636
2637 /// Create a list of pairs for all unknowns in this element,
2638 /// so that the first entry in each pair contains the global equation
2639 /// number of the unknown, while the second one contains the number
2640 /// of the "DOF type" that this unknown is associated with.
2641 /// (Function can obviously only be called if the equation numbering
2642 /// scheme has been set up.)
2644 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
2645 {
2646 // pair to store dof lookup prior to being added to list
2647 std::pair<unsigned, unsigned> dof_lookup;
2648
2649 dof_lookup.first = this->eqn_number(0);
2650 dof_lookup.second = 0;
2651
2652 // add to list
2653 dof_lookup_list.push_front(dof_lookup);
2654 }
2655
2656 protected:
2657 /// This function returns the residuals.
2658 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
2659 /// Note that this function also calculates the Jacobian contribution
2660 /// for the NetFluxControlElementForWomersleyPressureControl
2662 Vector<double>& residuals,
2663 DenseMatrix<double>& jacobian,
2664 const unsigned& flag)
2665 {
2666 // Get Womersley pressure and derivative with respect to the flux
2667 double womersley_pressure = 0.0;
2668 double d_womersley_pressure_d_q = 0.0;
2669
2670 // Get response of impedance tube
2671 Womersley_tube_pt->get_response(womersley_pressure,
2672 d_womersley_pressure_d_q);
2673
2674 // Get the current pressure
2675 double pressure = external_data_pt(Pressure_data_id)->value(0);
2676
2677 // Get equation number of the volume flux unknown
2678 int local_eq = internal_local_eqn(Volume_flux_data_id, 0);
2679
2680 // Calculate residuals
2681 residuals[local_eq] += pressure - womersley_pressure;
2682
2683 // Calculate Jacobian contributions if required
2684 if (flag)
2685 {
2686 // Get equation number of the pressure data unknown
2687 int local_unknown = external_local_eqn(Pressure_data_id, 0);
2688
2689 // Add the Jacobian contriburions
2690 jacobian(local_eq, local_eq) -= d_womersley_pressure_d_q;
2691 jacobian(local_eq, local_unknown) += 1.0;
2692 jacobian(local_unknown, local_eq) += 1.0;
2693 }
2694 }
2695
2696 private:
2697 /// Data object whose single value is the volume flux
2698 /// applied by the elements in the Flux_control_mesh_pt
2700
2701 /// Pointer to the Womersley impedance tube
2703
2704 /// Id of external Data object whose single value is the pressure
2706
2707 /// Id of internal Data object whose single value is the volume
2708 /// flux
2710 };
2711
2712
2713 /// ///////////////////////////////////////////////////////////////////////
2714 /// ///////////////////////////////////////////////////////////////////////
2715 /// ///////////////////////////////////////////////////////////////////////
2716
2717
2718 //======================================================================
2719 /// A class for an element to control net fluid flux across a boundary
2720 /// by imposing an applied pressure to the Navier-Stokes equations.
2721 /// This element is used with a mesh of NavierStokesFluxControlElements
2722 /// attached to the boundary. The flux imposed by this element is given
2723 /// by a NavierStokesWomersleyPressureControlElement.
2724 /// Note: fill_in_contribution_to_jacobian() does not calculate any
2725 /// Jacobian contributions for this element as they are calculated by
2726 /// NavierStokesFluxControlElement::fill_in_contribution_to_jacobian(...)
2727 /// and
2728 /// NavierStokesWomersleyPressureControlElement::
2729 /// fill_in_contribution_to_jacobian(...)
2730 //======================================================================
2732 : public virtual NetFluxControlElement
2733 {
2734 public:
2735 /// Constructor takes the mesh of
2736 /// TemplateFreeNavierStokesFluxControlElementBase which impose
2737 /// the pressure to controls the flux, plus a pointer to
2738 /// the PressureControlElement whoes internal data is the prescribed
2739 /// flux.
2741 Mesh* flux_control_mesh_pt,
2742 NavierStokesWomersleyPressureControlElement* pressure_control_element_pt)
2744 flux_control_mesh_pt,
2745 pressure_control_element_pt->volume_flux_data_pt()->value_pt(0))
2746 {
2747 // There's no need to add external data to this element since
2748 // this element's Jacobian contributions are calculated by the
2749 // NavierStokesFluxControlElements and the P
2750 // NavierStokesWomersleyPressureControlElement
2751
2752 // Add this elements Data to the external data of the
2753 // PressureControlElement
2754 pressure_control_element_pt->add_pressure_data(pressure_data_pt());
2755 }
2756
2757 /// Empty Destructor - Data gets deleted automatically
2759
2760 /// Broken copy constructor
2763
2764
2765 /// Broken assignment operator
2767 delete;
2768
2769
2770 /// The number of "DOF types" that degrees of freedom in this element
2771 /// are sub-divided into - set to 1
2772 unsigned ndof_types() const
2773 {
2774 return 1;
2775 }
2776
2777 /// Create a list of pairs for all unknowns in this element,
2778 /// so that the first entry in each pair contains the global equation
2779 /// number of the unknown, while the second one contains the number
2780 /// of the "DOF type" that this unknown is associated with.
2781 /// (Function can obviously only be called if the equation numbering
2782 /// scheme has been set up.)
2784 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
2785 {
2786 // pair to store dof lookup prior to being added to list
2787 std::pair<unsigned, unsigned> dof_lookup;
2788
2789 dof_lookup.first = this->eqn_number(0);
2790 dof_lookup.second = 0;
2791
2792 // add to list
2793 dof_lookup_list.push_front(dof_lookup);
2794 }
2795 };
2796
2797 /// //////////////////////////////////////////////////////////////////
2798 /// //////////////////////////////////////////////////////////////////
2799 /// //////////////////////////////////////////////////////////////////
2800
2801} // namespace oomph
2802
2803#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
////////////////////////////////////////////////////////////////////// //////////////////////////////...
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
void unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:391
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:271
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:293
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
bool is_pinned(const unsigned &i) const
Test whether the i-th variable is pinned (1: true; 0: false).
Definition: nodes.h:417
Information for documentation of results: Directory and file number to enable output in the form RESL...
std::string directory() const
Output directory.
unsigned & number()
Number used (e.g.) for labeling output files.
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:6006
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary,...
Definition: elements.h:4497
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4528
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4735
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s....
Definition: elements.cc:5242
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5328
void get_local_coordinate_in_bulk(const Vector< double > &s, Vector< double > &s_bulk) const
Calculate the vector of local coordinate in the bulk element given the local coordinates in this Face...
Definition: elements.cc:6384
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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:2593
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
Definition: elements.h:3050
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object.
Definition: elements.h:2509
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
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:1432
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3298
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
void check_J_eulerian_at_knots(bool &passed) const
Check that Jacobian of mapping between local and Eulerian coordinates at all integration points is po...
Definition: elements.cc:4237
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement....
Definition: elements.cc:5132
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 void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3220
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2470
A Generalised Element class.
Definition: elements.h:73
Data *& external_data_pt(const unsigned &i)
Return a pointer to i-th external data object.
Definition: elements.h:659
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
Definition: elements.cc:62
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Definition: elements.h:704
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:227
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data.
Definition: elements.h:267
int external_local_eqn(const unsigned &i, const unsigned &j)
Return the local equation number corresponding to the j-th value stored at the i-th external data.
Definition: elements.h:311
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) external data object to the element and return its index (i....
Definition: elements.cc:307
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:192
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
double * Prescribed_flux_pt
Pointer to current value of prescribed flux.
Data * Pressure_gradient_data_pt
Data item whose one and only value contains the pressure gradient.
void get_residuals(Vector< double > &residuals)
Compute residual vector: the volume flux constraint determines this element's one-and-only internal D...
void get_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute element residual Vector and element Jacobian matrix Note: Jacobian is zero because the deriva...
double total_volume_flux()
Get volume flux through all Womersley elements.
Mesh * Womersley_mesh_pt
Pointer to mesh that contains the Womersley elements.
ImposeFluxForWomersleyElement(Mesh *womersley_mesh_pt, double *prescribed_flux_pt)
Constructor: Pass pointer to mesh that contains the Womersley elements whose volume flux is controlle...
Data * pressure_gradient_data_pt()
Read-only access to the single-valued Data item that stores the pressure gradient (to be determined v...
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
A general mesh class.
Definition: mesh.h:67
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors.
Definition: mesh.h:75
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
void add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
Definition: mesh.h:617
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
////////////////////////////////////////////////////////////////////// //////////////////////////////...
virtual ~NavierStokesImpedanceTractionElementBase()
Empty vitual destructor.
Mesh * Navier_stokes_outflow_mesh_pt
Pointer to mesh containing the NavierStokesImpedanceTractionElements that contribute to the volume fl...
virtual void add_element_contribution_to_aux_integral(std::map< unsigned, double > *aux_integral_pt)=0
Add the element's contribution to the auxiliary integral used in the element's Jacobian....
virtual void set_external_data_from_navier_stokes_outflow_mesh(Mesh *navier_stokes_outflow_mesh_pt_mesh_pt)=0
Pass the pointer to the mesh containing all NavierStokesImpedanceTractionElements that contribute to ...
virtual double get_volume_flux()=0
Pure virtual function that must be implemented to compute the volume flux that passes through this el...
virtual void set_aux_integral_pt(std::map< unsigned, double > *aux_integral_pt)=0
Pass the pointer to the pre-computed auxiliary integral to the element so it can be accessed when com...
virtual void set_impedance_tube_pt(TemplateFreeWomersleyImpedanceTubeBase *impedance_tube_pt)=0
Pass the pointer to the "impedance tube" that computes the flow resistance via the solution of Womers...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
void fill_in_generic_residual_contribution_fluid_traction(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
This function returns the residuals for the traction function. flag=1(or 0): do (or don't) compute th...
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
Mesh *& navier_stokes_outflow_mesh_pt()
Access to mesh containing all NavierStokesImpedanceTractionElements that contribute to the volume flu...
~NavierStokesImpedanceTractionElement()
Destructor should not delete anything.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
void output(std::ostream &outfile)
Overload the output function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the element's contribution to the element's residual vector.
WomersleyImpedanceTubeBase< WOMERSLEY_ELEMENT, DIM > * Impedance_tube_pt
Pointer to ImpedanceTubeProblem that computes the flow resistance.
NavierStokesImpedanceTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in the element's contribution to the element's residual vector and Jacobian matrix.
std::map< unsigned, double > * Aux_integral_pt
Pointer to auxiliary integral, containing the derivative of the total volume flux through the outflow...
double get_volume_flux()
Get integral of volume flux through element.
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping.
void set_impedance_tube_pt(TemplateFreeWomersleyImpedanceTubeBase *impedance_tube_pt)
Set pointer to "impedance tube" that provides the flow resistance.
void set_aux_integral_pt(std::map< unsigned, double > *aux_integral_pt)
Set pointer to the precomputed auxiliary integral that contains the derivative of the total volume fl...
virtual int u_local_eqn(const unsigned &n, const unsigned &i)
Access function that returns the local equation numbers for velocity components. u_local_eqn(n,...
void add_element_contribution_to_aux_integral(std::map< unsigned, double > *aux_integral_pt)
Add the element's contribution to the auxiliary integral that contains the derivative of the total vo...
double total_volume_flux_into_downstream_tube()
Compute total volume flux into the "downstream tube" that provides the impedance (computed by adding ...
void set_external_data_from_navier_stokes_outflow_mesh(Mesh *navier_stokes_outflow_mesh_pt)
NavierStokesImpedanceTractionElements that contribute to the volume flux into the downstream "impedan...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns the residuals.
void fill_in_generic_residual_contribution_pressure_control(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
This function returns the residuals. flag=1(or 0): do (or don't) compute the Jacobian as well....
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into - set to 1.
Data * volume_flux_data_pt() const
Function to return a pointer to the Data object whose single value is the flux degree of freedom.
void add_pressure_data(Data *pressure_data_pt)
Function to add to external data the Data object whose single value is the pressure applied at the bo...
TemplateFreeWomersleyImpedanceTubeBase * Womersley_tube_pt
Pointer to the Womersley impedance tube.
unsigned Volume_flux_data_id
Id of internal Data object whose single value is the volume flux.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the Jacobian, plus the Jacobian contribution for the NetFluxC...
~NavierStokesWomersleyPressureControlElement()
Destructor should not delete anything.
Data * Volume_flux_data_pt
Data object whose single value is the volume flux applied by the elements in the Flux_control_mesh_pt...
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
unsigned Pressure_data_id
Id of external Data object whose single value is the pressure.
NavierStokesWomersleyPressureControlElement(TemplateFreeWomersleyImpedanceTubeBase *womersley_tube_pt)
Constructor takes a pointer to a suitable Womersley impedance tube which defines the pressure via get...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
void operator=(const NetFluxControlElementForWomersleyPressureControl &)=delete
Broken assignment operator.
NetFluxControlElementForWomersleyPressureControl(const NetFluxControlElementForWomersleyPressureControl &dummy)=delete
Broken copy constructor.
NetFluxControlElementForWomersleyPressureControl(Mesh *flux_control_mesh_pt, NavierStokesWomersleyPressureControlElement *pressure_control_element_pt)
Constructor takes the mesh of TemplateFreeNavierStokesFluxControlElementBase which impose the pressur...
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into - set to 1.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
~NetFluxControlElementForWomersleyPressureControl()
Empty Destructor - Data gets deleted automatically.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Data * pressure_data_pt() const
Function to return a pointer to the Data object whose single value is the pressure applied by the Nav...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: elements.h:3439
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
Definition: problem.cc:1989
void add_time_stepper_pt(TimeStepper *const &time_stepper_pt)
Add a timestepper to the problem. The function will automatically create or resize the Time object so...
Definition: problem.cc:1545
TimeStepper *& time_stepper_pt()
Access function for the pointer to the first (presumably only) timestepper.
Definition: problem.h:1524
double & time()
Return the current value of continuous time.
Definition: problem.cc:11531
bool Problem_is_nonlinear
Boolean flag indicating if we're dealing with a linear or nonlinear Problem – if set to false the New...
Definition: problem.h:628
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1504
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
void operator=(const QWomersleyElement< DIM, NNODE_1D > &)=delete
Broken assignment operator.
void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot ...
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
unsigned required_nvalue(const unsigned &n) const
Required # of ‘values’ (pinned or dofs) at node n.
QWomersleyElement()
Constructor: Call constructors for QElement and Womersley equations.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
double dshape_and_dtest_eulerian_at_knot_womersley(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
QWomersleyElement(const QWomersleyElement< DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points.
static const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
double dshape_and_dtest_eulerian_womersley(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
virtual ~TemplateFreeWomersleyImpedanceTubeBase()
Empty virtual destructor.
virtual void get_response(double &p_in, double &dp_in_dq)=0
Empty virtual dummy member function – every base class needs at least one virtual member function if ...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
static bool Suppress_warning_about_unpinned_nst_dofs
Static bool to suppress warning.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:389
////////////////////////////////////////////////////////////////////// //////////////////////////////...
void operator=(const WomersleyEquations &)=delete
Broken assignment operator.
WomersleyEquations()
Constructor: Initialises the Pressure_gradient_data_pt to null.
virtual double dshape_and_dtest_eulerian_womersley(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at local coord. s; return Jacobian of mapping...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute element residual Vector and element Jacobian matrix (wrapper)
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 get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: flux[i] = du/dx_i.
virtual double dshape_and_dtest_eulerian_at_knot_womersley(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
double du_dt_womersley(const unsigned &n) const
du/dt at local node n. Uses suitably interpolated value for hanging nodes.
double interpolated_u_womersley(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
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 *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
double get_volume_flux()
Compute total volume flux through element.
void output(FILE *file_pt)
C_style output with default number of plot points.
unsigned self_test()
Self-test: Return 0 for OK.
Data * Pressure_gradient_data_pt
Pointer to pressure gradient Data (single value Data item)
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute element residual Vector (wrapper)
void set_pressure_gradient_pt(Data *&pressure_gradient_data_pt)
Set pointer to pressure gradient (single-valued Data)
void set_pressure_gradient_and_add_as_external_data(Data *pressure_gradient_data_pt)
Set pointer to pressure gradient (single-valued Data) and treat it as external data – this can only b...
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.
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
WomersleyEquations(const WomersleyEquations &dummy)=delete
Broken copy constructor.
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...
Data * set_pressure_gradient_pt() const
Read-only access to pointer to pressure gradient.
virtual unsigned u_index_womersley() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
double(* PrescribedVolumeFluxFctPt)(const double &time)
Function pointer to fct that prescribes volume flux q=fct(t) – mainly used for validation purposes.
double Length
Length of the tube.
double Dp_in_dq
Derivative of inflow pressure w.r.t. instantaenous volume flux (Note: Can be pre-computed)
WomersleyImpedanceTubeBase(const double &length, PrescribedVolumeFluxFctPt prescribed_volume_flux_fct_pt)
Constructor: Specify length of tube and pointer to function that specifies the prescribed volume flux...
Mesh * Navier_stokes_outflow_mesh_pt
Pointer to the mesh of NavierStokesImpedanceTractionElements that are attached to the outflow cross-s...
void get_response(double &p_in, double &dp_in_dq)
Compute inlet pressure, p_in, required to achieve the currently imposed, instantaneous volume flux q ...
double & p_out()
Access fct to outlet pressure.
void precompute_aux_integrals()
Precompute auxiliary integrals required for the computation of the Jacobian in the NavierStokesImpeda...
double * Current_volume_flux_pt
Pointer to double that specifies the currently imposed instantaneous volume flux into the impedance t...
virtual Mesh * build_mesh_and_apply_boundary_conditions(TimeStepper *time_stepper_pt)=0
Pure virtual function in which the user of a derived class must create the mesh of WomersleyElements ...
WomersleyImpedanceTubeBase(const double &length, Mesh *navier_stokes_outflow_mesh_pt)
Constructor: Specify length of tube and the pointer to the mesh of either NavierStokesImpedanceTracti...
void setup()
Set up the Womersley tubes so that a subsequent call to get_response(...) computes the inlet pressure...
void setup(double *re_st_pt, const double &dt, const double &q_initial, TimeStepper *time_stepper_pt=0)
Set up the Womersley tubes so that a subsequent call to get_response(...) computes the inlet pressure...
WomersleyProblem< ELEMENT, DIM > * womersley_problem_pt()
Access to underlying Womersley problem.
double total_volume_flux_into_impedance_tube()
Compute total current volume flux into the "impedance tube" that provides the flow resistance (flux i...
WomersleyProblem< ELEMENT, DIM > * Womersley_problem_pt
Pointer to Womersley problem that determines the pressure gradient along the tube.
std::map< unsigned, double > * Aux_integral_pt
Pointer to auxiliary integral, containing the derivative of the total volume flux through the outflow...
void shift_time_values(const double &dt)
Shift history values to allow coputation of next timestep. Note: When used with a full Navier-Stokes ...
PrescribedVolumeFluxFctPt Prescribed_volume_flux_fct_pt
Pointer to function that specifies the prescribed volume flux.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
WomersleyMesh(Mesh *n_st_outflow_mesh_pt, TimeStepper *time_stepper_pt, const unsigned &fixed_coordinate, const unsigned &w_index)
Constructor: Pass pointer to mesh of face elements in the outflow cross-section of a full Navier-Stok...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Mesh * build_mesh_and_apply_boundary_conditions(TimeStepper *time_stepper_pt)
Implement pure virtual fct (defined in the base class WomersleyImpedanceTubeBase) that builds the mes...
WomersleyOutflowImpedanceTube(const double &length, Mesh *navier_stokes_outflow_mesh_pt, const unsigned &fixed_coordinate, const unsigned &w_index)
Constructor: Pass length and mesh of face elements that are attached to the outflow cross-section of ...
unsigned Fixed_coordinate
The coordinate (in the higher-dimensional Navier-Stokes mesh) that is constant in the outflow cross-s...
unsigned W_index
The velocity component (in terms of the nodal index) that represents the outflow component – the latt...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void doc_solution(DocInfo &doc_info, const double &z_out=0.0)
Doc the solution.
ImposeFluxForWomersleyElement< DIM > * Flux_el_pt
Pointer to element that imposes the flux through the collection of Womersley elements.
PrescribedPressureGradientFctPt Prescribed_pressure_gradient_fct_pt
Fct pointer to fct that prescribes pressure gradient.
void doc_solution(DocInfo &doc_info, std::ofstream &trace_file, const double &z_out=0.0)
Doc the solution incl. trace file for various quantities of interest (to some...)
void actions_before_newton_solve()
Update the problem specs before solve (empty)
double * Prescribed_volume_flux_pt
Pointer to currently prescribed volume flux.
~WomersleyProblem()
Destructor to clean up memory.
WomersleyProblem(double *re_st_pt, double *prescribed_volume_flux_pt, TimeStepper *time_stepper_pt, Mesh *womersley_mesh_pt)
Constructor: Pass pointer to Womersley number, pointer to the double that stores the currently impose...
double(* PrescribedPressureGradientFctPt)(const double &time)
Function pointer to fct that prescribes pressure gradient g=fct(t)
Data * Pressure_gradient_data_pt
Pointer to single-valued Data item that stores pressure gradient.
Data * pressure_gradient_data_pt()
Access function to the single-valued Data object that contains the unknown pressure gradient (used if...
void actions_before_implicit_timestep()
Update the problem specs before next timestep: Update time-varying pressure gradient (if prescribed)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...