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