axisym_fvk_elements.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Header file for axisymmetric FoepplvonKarman elements
27 #ifndef OOMPH_AXISYM_FOEPPLVONKARMAN_ELEMENTS_HEADER
28 #define OOMPH_AXISYM_FOEPPLVONKARMAN_ELEMENTS_HEADER
29 
30 
31 // Config header generated by autoconfig
32 #ifdef HAVE_CONFIG_H
33 #include <oomph-lib-config.h>
34 #endif
35 
36 #include "../generic/nodes.h"
37 #include "../generic/Qelements.h"
38 #include "../generic/oomph_utilities.h"
39 #include "../generic/element_with_external_element.h"
40 
41 namespace oomph
42 {
43  //=============================================================
44  /// A class for all isoparametric elements that solve the
45  /// axisum Foeppl von Karman equations.
46  ///
47  /// This contains the generic maths. Shape functions, geometric
48  /// mapping etc. must get implemented in derived class.
49  //=============================================================
50  class AxisymFoepplvonKarmanEquations : public virtual FiniteElement
51  {
52  public:
53  /// Function pointer to pressure function fct(r,f(r)) --
54  /// r is a Vector!
55  typedef void (*AxisymFoepplvonKarmanPressureFctPt)(const double& r,
56  double& f);
57 
58  /// Constructor (must initialise the Pressure_fct_pt and
59  /// Airy_forcing_fct_pt to null). Also set physical parameters to their
60  /// default values.
63  {
64  // Set all the physical constants to the default value (zero)
66  Linear_bending_model = false;
67  }
68 
69  /// Broken copy constructor
71  const AxisymFoepplvonKarmanEquations& dummy) = delete;
72 
73  /// Broken assignment operator
75 
76  /// FvK parameter
77  const double& eta() const
78  {
79  return *Eta_pt;
80  }
81 
82  /// Pointer to FvK parameter
83  double*& eta_pt()
84  {
85  return Eta_pt;
86  }
87 
88  /// Return the index at which the i-th unknown value
89  /// is stored. The default value, i, is appropriate for single-physics
90  /// problems. By default, these are:
91  /// 0: w
92  /// 1: laplacian w
93  /// 2: phi
94  /// 3: laplacian phi
95  /// 4-5: smooth first derivatives
96  /// In derived multi-physics elements, this function should be overloaded
97  /// to reflect the chosen storage scheme. Note that these equations require
98  /// that the unknown is always stored at the same index at each node.
99  virtual inline unsigned nodal_index_fvk(const unsigned& i = 0) const
100  {
101  return i;
102  }
103 
104  /// Output with default number of plot points
105  void output(std::ostream& outfile)
106  {
107  const unsigned n_plot = 5;
108  output(outfile, n_plot);
109  }
110 
111  /// Output FE representation of soln: r,w,sigma_r_r,sigma_phi_phi
112  /// at n_plot plot points
113  void output(std::ostream& outfile, const unsigned& n_plot);
114 
115  /// C_style output with default number of plot points
116  void output(FILE* file_pt)
117  {
118  const unsigned n_plot = 5;
119  output(file_pt, n_plot);
120  }
121 
122  /// C-style output FE representation of soln: r,w at
123  /// n_plot plot points
124  void output(FILE* file_pt, const unsigned& n_plot);
125 
126  /// Output exact soln: r,w_exact at n_plot plot points
127  void output_fct(std::ostream& outfile,
128  const unsigned& n_plot,
130 
131  /// Output exact soln: r,w_exact at
132  /// n_plot plot points (dummy time-dependent version to
133  /// keep intel compiler happy)
134  virtual void output_fct(
135  std::ostream& outfile,
136  const unsigned& n_plot,
137  const double& time,
139  {
140  throw OomphLibError(
141  "There is no time-dependent output_fct() for Foeppl von Karman"
142  "elements ",
143  OOMPH_CURRENT_FUNCTION,
144  OOMPH_EXCEPTION_LOCATION);
145  }
146 
147  /// Get error against and norm of exact solution
148  void compute_error(std::ostream& outfile,
150  double& error,
151  double& norm);
152 
153 
154  /// Dummy, time dependent error checker
155  void compute_error(std::ostream& outfile,
157  const double& time,
158  double& error,
159  double& norm)
160  {
161  throw OomphLibError(
162  "There is no time-dependent compute_error() for Foeppl von Karman"
163  "elements",
164  OOMPH_CURRENT_FUNCTION,
165  OOMPH_EXCEPTION_LOCATION);
166  }
167 
168  /// Access function: Pointer to pressure function
170  {
171  return Pressure_fct_pt;
172  }
173 
174  /// Access function: Pointer to pressure function. Const version
176  {
177  return Pressure_fct_pt;
178  }
179 
180  /// Access function: Pointer to Airy forcing function
182  {
183  return Airy_forcing_fct_pt;
184  }
185 
186  /// Access function: Pointer to Airy forcing function. Const version
188  {
189  return Airy_forcing_fct_pt;
190  }
191 
192  /// Get pressure term at (Eulerian) position r. This function is
193  /// virtual to allow overloading in multi-physics problems where
194  /// the strength of the pressure function might be determined by
195  /// another system of equations.
196  inline virtual void get_pressure_fvk(const unsigned& ipt,
197  const double& r,
198  double& pressure) const
199  {
200  // If no pressure function has been set, return zero
201  if (Pressure_fct_pt == 0)
202  {
203  pressure = 0.0;
204  }
205  else
206  {
207  // Get pressure strength
208  (*Pressure_fct_pt)(r, pressure);
209  }
210  }
211 
212  /// Get Airy forcing term at (Eulerian) position r. This function is
213  /// virtual to allow overloading in multi-physics problems where
214  /// the strength of the pressure function might be determined by
215  /// another system of equations.
216  inline virtual void get_airy_forcing_fvk(const unsigned& ipt,
217  const double& r,
218  double& airy_forcing) const
219  {
220  // If no Airy forcing function has been set, return zero
221  if (Airy_forcing_fct_pt == 0)
222  {
223  airy_forcing = 0.0;
224  }
225  else
226  {
227  // Get Airy forcing strength
228  (*Airy_forcing_fct_pt)(r, airy_forcing);
229  }
230  }
231 
232  /// Get gradient of deflection: gradient[i] = dw/dr_i */
234  Vector<double>& gradient) const
235  {
236  // Find out how many nodes there are in the element
237  const unsigned n_node = nnode();
238 
239  // Get the index at which the unknown is stored
240  const unsigned w_nodal_index = nodal_index_fvk(0);
241 
242  // Set up memory for the shape and test functions
243  Shape psi(n_node);
244  DShape dpsidr(n_node, 1);
245 
246  // Call the derivatives of the shape and test functions
247  dshape_eulerian(s, psi, dpsidr);
248 
249  // Initialise to zero
250  gradient[0] = 0.0;
251 
252  // Loop over nodes
253  for (unsigned l = 0; l < n_node; l++)
254  {
255  gradient[0] += this->nodal_value(l, w_nodal_index) * dpsidr(l, 0);
256  }
257  }
258 
259  /// Fill in the residuals with this element's contribution
261 
262 
263  // hierher Jacobian not yet implemented
264  // void fill_in_contribution_to_jacobian(Vector<double> &residuals,
265  // DenseMatrix<double> &jacobian);
266 
267  /// Return FE representation of vertical displacement, w_fvk(s)
268  /// at local coordinate s
269  inline double interpolated_w_fvk(const Vector<double>& s) const
270  {
271  // Find number of nodes
272  const unsigned n_node = nnode();
273 
274  // Get the index at which the unknown is stored
275  const unsigned w_nodal_index = nodal_index_fvk(0);
276 
277  // Local shape function
278  Shape psi(n_node);
279 
280  // Find values of shape function
281  shape(s, psi);
282 
283  // Initialise value of u
284  double interpolated_w = 0.0;
285 
286  // Loop over the local nodes and sum
287  for (unsigned l = 0; l < n_node; l++)
288  {
289  interpolated_w += this->nodal_value(l, w_nodal_index) * psi[l];
290  }
291 
292  return (interpolated_w);
293  }
294 
295  /// Compute in-plane stresses.
297  double& sigma_r_r,
298  double& sigma_phi_phi) const;
299 
300 
301  /// Self-test: Return 0 for OK
302  unsigned self_test();
303 
304  /// Sets a flag to signify that we are solving the linear,
305  /// pure bending equations, and pin all the nodal values that will
306  /// not be used in this case
308  {
309  // Set the boolean flag
310  Linear_bending_model = true;
311 
312  // Get the index of the first FvK nodal value
313  unsigned first_fvk_nodal_index = nodal_index_fvk();
314 
315  // Get the total number of FvK nodal values (assuming they are stored
316  // contiguously) at node 0 (it's the same at all nodes anyway)
317  unsigned total_fvk_nodal_indices = 6;
318 
319  // Get the number of nodes in this element
320  unsigned n_node = nnode();
321 
322  // Loop over the appropriate nodal indices
323  for (unsigned index = first_fvk_nodal_index + 2;
324  index < first_fvk_nodal_index + total_fvk_nodal_indices;
325  index++)
326  {
327  // Loop over the nodes in the element
328  for (unsigned inod = 0; inod < n_node; inod++)
329  {
330  // Pin the nodal value at the current index
331  node_pt(inod)->pin(index);
332  }
333  }
334  }
335 
336 
337  protected:
338  /// Shape/test functions and derivs w.r.t. to global coords at
339  /// local coord. s; return Jacobian of mapping
341  const Vector<double>& s,
342  Shape& psi,
343  DShape& dpsidr,
344  Shape& test,
345  DShape& dtestdr) const = 0;
346 
347 
348  /// Shape/test functions and derivs w.r.t. to global coords at
349  /// integration point ipt; return Jacobian of mapping
351  const unsigned& ipt,
352  Shape& psi,
353  DShape& dpsidr,
354  Shape& test,
355  DShape& dtestdr) const = 0;
356 
357  /// Pointer to FvK parameter
358  double* Eta_pt;
359 
360  /// Pointer to pressure function:
362 
363  /// Pointer to Airy forcing function
365 
366  /// Default value for physical constants
367  static double Default_Physical_Constant_Value;
368 
369  /// Flag which stores whether we are using a linear,
370  /// pure bending model instead of the full non-linear Foeppl-von Karman
372  };
373 
374 
375  ///////////////////////////////////////////////////////////////////////////
376  ///////////////////////////////////////////////////////////////////////////
377  ///////////////////////////////////////////////////////////////////////////
378 
379 
380  //======================================================================
381  /// Axisym FoepplvonKarmanElement elements are 1D
382  /// Foeppl von Karman elements with isoparametric interpolation for the
383  /// function.
384  //======================================================================
385  template<unsigned NNODE_1D>
387  : public virtual QElement<1, NNODE_1D>,
388  public virtual AxisymFoepplvonKarmanEquations
389  {
390  private:
391  /// Static int that holds the number of variables at
392  /// nodes: always the same
393  static const unsigned Initial_Nvalue;
394 
395  public:
396  /// Constructor: Call constructors for QElement and
397  /// AxisymFoepplvonKarmanEquations
399  : QElement<1, NNODE_1D>(), AxisymFoepplvonKarmanEquations()
400  {
401  }
402 
403  /// Broken copy constructor
405  const AxisymFoepplvonKarmanElement<NNODE_1D>& dummy) = delete;
406 
407  /// Broken assignment operator
409 
410  /// Required # of `values' (pinned or dofs)
411  /// at node n
412  inline unsigned required_nvalue(const unsigned& n) const
413  {
414  return Initial_Nvalue;
415  }
416 
417 
418  /// Output function:
419  /// r,w,sigma_r_r,sigma_phi_phi
420  void output(std::ostream& outfile)
421  {
423  }
424 
425  /// Output function:
426  /// r,w,sigma_r_r,sigma_phi_phi at n_plot plot points
427  void output(std::ostream& outfile, const unsigned& n_plot)
428  {
430  }
431 
432  /// C-style output function:
433  /// r,w
434  void output(FILE* file_pt)
435  {
437  }
438 
439  /// C-style output function:
440  /// r,w at n_plot plot points
441  void output(FILE* file_pt, const unsigned& n_plot)
442  {
444  }
445 
446  /// Output function for an exact solution:
447  /// r,w_exact at n_plot plot points
448  void output_fct(std::ostream& outfile,
449  const unsigned& n_plot,
451  {
453  outfile, n_plot, exact_soln_pt);
454  }
455 
456  /// Output function for a time-dependent exact solution.
457  /// r,w_exact at n_plot plot points
458  /// (Calls the steady version)
459  void output_fct(std::ostream& outfile,
460  const unsigned& n_plot,
461  const double& time,
463  {
465  outfile, n_plot, time, exact_soln_pt);
466  }
467 
468 
469  protected:
470  /// Shape, test functions & derivs. w.r.t. to global coords.
471  /// Return Jacobian.
473  Shape& psi,
474  DShape& dpsidr,
475  Shape& test,
476  DShape& dtestdr) const;
477 
478  /// Shape, test functions & derivs. w.r.t. to global coords. at
479  /// integration point ipt. Return Jacobian.
481  const unsigned& ipt,
482  Shape& psi,
483  DShape& dpsidr,
484  Shape& test,
485  DShape& dtestdr) const;
486  };
487 
488 
489  // Inline functions:
490 
491  //======================================================================
492  /// Define the shape functions and test functions and derivatives
493  /// w.r.t. global coordinates and return Jacobian of mapping.
494  ///
495  /// Galerkin: Test functions = shape functions
496  //======================================================================
497  template<unsigned NNODE_1D>
499  NNODE_1D>::dshape_and_dtest_eulerian_axisym_fvk(const Vector<double>& s,
500  Shape& psi,
501  DShape& dpsidr,
502  Shape& test,
503  DShape& dtestdr) const
504 
505  {
506  // Call the geometrical shape functions and derivatives
507  const double J = this->dshape_eulerian(s, psi, dpsidr);
508 
509  // Set the test functions equal to the shape functions
510  test = psi;
511  dtestdr = dpsidr;
512 
513  // Return the jacobian
514  return J;
515  }
516 
517 
518  //======================================================================
519  /// Define the shape functions and test functions and derivatives
520  /// w.r.t. global coordinates and return Jacobian of mapping.
521  ///
522  /// Galerkin: Test functions = shape functions
523  //======================================================================
524  template<unsigned NNODE_1D>
527  Shape& psi,
528  DShape& dpsidr,
529  Shape& test,
530  DShape& dtestdr) const
531  {
532  // Call the geometrical shape functions and derivatives
533  const double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidr);
534 
535  // Set the pointers of the test functions
536  test = psi;
537  dtestdr = dpsidr;
538 
539  // Return the jacobian
540  return J;
541  }
542 
543 
544  ///////////////////////////////////////////////////////////////////////////
545  ///////////////////////////////////////////////////////////////////////////
546  ///////////////////////////////////////////////////////////////////////////
547 
548 
549  //======================================================================
550  /// FSI Axisym FoepplvonKarmanElement elements are 1D
551  /// Foeppl von Karman elements with isoparametric interpolation for the
552  /// function. Gets traction from adjacent fluid element(s) of type
553  /// FLUID_ELEMENT.
554  //======================================================================
555  template<unsigned NNODE_1D, class FLUID_ELEMENT>
557  : public virtual AxisymFoepplvonKarmanElement<NNODE_1D>,
558  public virtual ElementWithExternalElement
559  {
560  public:
561  /// Constructor
564  {
565  // Set source element storage: one interaction with an external
566  // element -- the fluid bulk element that provides the pressure
567  this->set_ninteraction(1);
568  }
569 
570  /// Empty virtual destructor
572 
573  /// Return the ratio of the stress scales used to non-dimensionalise
574  /// the fluid and elasticity equations.
575  const double& q() const
576  {
577  return *Q_pt;
578  }
579 
580  /// Return a pointer the ratio of stress scales used to
581  /// non-dimensionalise the fluid and solid equations.
582  double*& q_pt()
583  {
584  return Q_pt;
585  }
586 
587  /// How many items of Data does the shape of the object depend on?
588  /// All nodal data
589  virtual unsigned ngeom_data() const
590  {
591  return this->nnode();
592  }
593 
594  /// Return pointer to the j-th Data item that the object's
595  /// shape depends on.
596  virtual Data* geom_data_pt(const unsigned& j)
597  {
598  return this->node_pt(j);
599  }
600 
601  /// Overloaded position function: Return 2D position vector:
602  /// (r(zeta),z(zeta)) of material point whose "Lagrangian coordinate"
603  /// is given by zeta. Here r=zeta!
604  void position(const Vector<double>& zeta, Vector<double>& r) const
605  {
606  const unsigned t = 0;
607  this->position(t, zeta, r);
608  }
609 
610  /// Overloaded position function: Return 2D position vector:
611  /// (r(zeta),z(zeta)) of material point whose "Lagrangian coordinate"
612  /// is given by zeta.
613  void position(const unsigned& t,
614  const Vector<double>& zeta,
615  Vector<double>& r) const
616  {
617  // Find number of nodes
618  const unsigned n_node = this->nnode();
619 
620  // Get the index at which the poisson unknown is stored
621  const unsigned w_nodal_index = this->nodal_index_fvk(0);
622 
623  // Local shape function
624  Shape psi(n_node);
625 
626  // Find values of shape function
627  this->shape(zeta, psi);
628 
629  // Initialise
630  double interpolated_w = 0.0;
631  double interpolated_r = 0.0;
632 
633  // Loop over the local nodes and sum
634  for (unsigned l = 0; l < n_node; l++)
635  {
636  interpolated_w += this->nodal_value(t, l, w_nodal_index) * psi[l];
637  interpolated_r += this->node_pt(l)->x(t, 0) * psi[l];
638  }
639 
640  // Assign
641  r[0] = interpolated_r;
642  r[1] = interpolated_w;
643  }
644 
645  /// j-th time-derivative on object at current time:
646  /// \f$ \frac{d^{j} r(\zeta)}{dt^j} \f$.
647  void dposition_dt(const Vector<double>& zeta,
648  const unsigned& j,
649  Vector<double>& drdt)
650  {
651  // Find number of nodes
652  const unsigned n_node = this->nnode();
653 
654  // Get the index at which the poisson unknown is stored
655  const unsigned w_nodal_index = this->nodal_index_fvk(0);
656 
657  // Local shape function
658  Shape psi(n_node);
659 
660  // Find values of shape function
661  this->shape(zeta, psi);
662 
663  // Initialise
664  double interpolated_dwdt = 0.0;
665  double interpolated_drdt = 0.0;
666 
667  // Loop over the local nodes and sum
668  for (unsigned l = 0; l < n_node; l++)
669  {
670  // Get the timestepper
672 
673  // If we are doing an unsteady solve then calculate the derivative
674  if (!time_stepper_pt->is_steady())
675  {
676  // Get the number of values required to represent history
677  const unsigned n_time = time_stepper_pt->ntstorage();
678 
679  // Loop over history values
680  for (unsigned t = 0; t < n_time; t++)
681  {
682  // Add the contribution to the derivative
683  interpolated_dwdt += time_stepper_pt->weight(1, t) *
684  this->nodal_value(t, l, w_nodal_index) *
685  psi[l];
686  }
687  }
688  }
689 
690  // Assign
691  drdt[0] = interpolated_drdt;
692  drdt[1] = interpolated_dwdt;
693  }
694 
695 
696  /// Overload pressure term at (Eulerian) position r.
697  /// Adds fluid traction to pressure imposed by "pressure fct pointer"
698  /// (which can be regarded as applying an external (i.e.
699  /// "on the other side" of the fluid) pressure
700  inline virtual void get_pressure_fvk(const unsigned& ipt,
701  const double& r,
702  double& pressure) const
703  {
704  pressure = 0.0;
705 
706  // Get underlying version
708 
709  // Get FSI parameter
710  const double q_value = q();
711 
712  // Get fluid element
713  FLUID_ELEMENT* ext_el_pt =
714  dynamic_cast<FLUID_ELEMENT*>(external_element_pt(0, ipt));
716 
717  // Outer unit normal is vertically upward (in z direction)
718  // (within an essentiall flat) model for the wall)
719  Vector<double> normal(2);
720  normal[0] = 0.0;
721  normal[1] = 1.0;
722 
723  // Get traction
724  Vector<double> traction(3);
725  ext_el_pt->traction(s_ext, normal, traction);
726 
727  // Add z-component of traction
728  pressure -= q_value * traction[1];
729  }
730 
731 
732  /// Output integration points (for checking of fsi setup)
733  void output_integration_points(std::ostream& outfile)
734  {
735  // How many nodes do we have?
736  unsigned nnod = this->nnode();
737  Shape psi(nnod);
738 
739  // Get the index at which the unknown is stored
740  const unsigned w_nodal_index = this->nodal_index_fvk(0);
741 
742  // Loop over the integration points
743  const unsigned n_intpt = this->integral_pt()->nweight();
744  outfile << "ZONE I=" << n_intpt << std::endl;
745  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
746  {
747  // Get shape fct
748  Vector<double> s(1);
749  s[0] = this->integral_pt()->knot(ipt, 0);
750  shape(s, psi);
751 
752  // Initialise
753  double interpolated_w = 0.0;
754  double interpolated_r = 0.0;
755 
756  // Loop over the local nodes and sum
757  for (unsigned l = 0; l < nnod; l++)
758  {
759  interpolated_w += this->nodal_value(l, w_nodal_index) * psi[l];
760  interpolated_r += this->node_pt(l)->x(0) * psi[l];
761  }
762 
763  // Get fluid element
764  FLUID_ELEMENT* ext_el_pt =
765  dynamic_cast<FLUID_ELEMENT*>(external_element_pt(0, ipt));
767 
768  // Get veloc
769  Vector<double> veloc(3);
770  ext_el_pt->interpolated_u_axi_nst(s_ext, veloc);
771  Vector<double> x(2);
772  ext_el_pt->interpolated_x(s_ext, x);
773 
774  outfile << interpolated_r << " " << interpolated_w << " " << veloc[0]
775  << " " << veloc[1] << " " << x[0] << " " << x[1] << " "
776  << std::endl;
777  }
778  }
779 
780 
781  /// Output adjacent fluid elements (for checking of fsi setup)
782  void output_adjacent_fluid_elements(std::ostream& outfile,
783  const unsigned& nplot)
784  {
785  // Loop over the integration points
786  const unsigned n_intpt = this->integral_pt()->nweight();
787  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
788  {
789  // Get fluid element
790  FLUID_ELEMENT* ext_el_pt =
791  dynamic_cast<FLUID_ELEMENT*>(external_element_pt(0, ipt));
792 
793  // Dump it
794  ext_el_pt->output(outfile, nplot);
795  }
796  }
797 
798  /// Perform any auxiliary node update fcts of the adjacent
799  /// fluid elements
801  {
803  }
804 
805  /// Perform any auxiliary node update fcts of the adjacent
806  /// fluid elements
808  {
810  }
811 
812 
813  /// Perform any auxiliary node update fcts of the adjacent
814  /// fluid elements
816  {
818  }
819 
820 
821  /// Perform any auxiliary node update fcts of the adjacent
822  /// fluid elements
824  {
826  }
827 
828 
829  /// Update the nodal positions in all fluid elements that affect
830  /// the traction on this element
832  {
833  // Don't update elements repeatedly
834  std::map<FLUID_ELEMENT*, bool> done;
835 
836  // Number of integration points
837  unsigned n_intpt = integral_pt()->nweight();
838 
839  // Loop over all integration points in wall element
840  for (unsigned iint = 0; iint < n_intpt; iint++)
841  {
842  // Get fluid element that affects this integration point
843  FLUID_ELEMENT* el_f_pt =
844  dynamic_cast<FLUID_ELEMENT*>(external_element_pt(0, iint));
845 
846  // Is there an adjacent fluid element?
847  if (el_f_pt != 0)
848  {
849  // Have we updated its positions yet?
850  if (!done[el_f_pt])
851  {
852  // Update nodal positions
853  el_f_pt->node_update();
854  done[el_f_pt] = true;
855  }
856  }
857  }
858  }
859 
860 
861  /// Output FE representation of soln:
862  /// r,w,dwdt,sigma_r_r,sigma_phi_phi at n_plot plot points
863  void output(std::ostream& outfile, const unsigned& n_plot)
864  {
865  // Vector of local coordinates
866  Vector<double> s(1);
867 
868  // Tecplot header info
869  outfile << "ZONE\n";
870 
871  // Loop over plot points
872  unsigned num_plot_points = nplot_points(n_plot);
873  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
874  {
875  // Get local coordinates of plot point
876  get_s_plot(iplot, n_plot, s);
877 
878  // Get velocity
879  Vector<double> drdt(2);
880  dposition_dt(s, 1, drdt);
881 
882  // Get stress
883  double sigma_r_r = 0.0;
884  double sigma_phi_phi = 0.0;
885  this->interpolated_stress(s, sigma_r_r, sigma_phi_phi);
886  outfile << this->interpolated_x(s, 0) << " "
887  << this->interpolated_w_fvk(s) << " " << drdt[0] << " "
888  << drdt[1] << " " << sigma_r_r << " " << sigma_phi_phi
889  << std::endl;
890  }
891  }
892 
893 
894  protected:
895  /// Pointer to the ratio, \f$ Q \f$ , of the stress used to
896  /// non-dimensionalise the fluid stresses to the stress used to
897  /// non-dimensionalise the solid stresses.
898  double* Q_pt;
899  };
900 
901 
902 } // namespace oomph
903 
904 #endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,w at n_plot plot points.
double dshape_and_dtest_eulerian_axisym_fvk(const Vector< double > &s, Shape &psi, DShape &dpsidr, Shape &test, DShape &dtestdr) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
void output(FILE *file_pt)
C-style output function: r,w.
double dshape_and_dtest_eulerian_at_knot_axisym_fvk(const unsigned &ipt, Shape &psi, DShape &dpsidr, Shape &test, DShape &dtestdr) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt....
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. r,w_exact at n_plot plot points (Calls the stead...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: r,w,sigma_r_r,sigma_phi_phi at n_plot plot points.
AxisymFoepplvonKarmanElement(const AxisymFoepplvonKarmanElement< NNODE_1D > &dummy)=delete
Broken copy constructor.
void output(std::ostream &outfile)
Output function: r,w,sigma_r_r,sigma_phi_phi.
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
unsigned required_nvalue(const unsigned &n) const
Required # of ‘values’ (pinned or dofs) at node n.
void operator=(const AxisymFoepplvonKarmanElement< NNODE_1D > &)=delete
Broken assignment operator.
AxisymFoepplvonKarmanElement()
Constructor: Call constructors for QElement and AxisymFoepplvonKarmanEquations.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: r,w_exact at n_plot plot points.
A class for all isoparametric elements that solve the axisYm Foeppl von Karman equations in a displac...
AxisymFoepplvonKarmanPressureFctPt pressure_fct_pt() const
Access function: Pointer to pressure function. Const version.
AxisymFoepplvonKarmanPressureFctPt Pressure_fct_pt
Pointer to pressure function:
static double Default_Physical_Constant_Value
Default value for physical constants.
void operator=(const AxisymFoepplvonKarmanEquations &)=delete
Broken assignment operator.
virtual void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,w_exact at n_plot plot points (dummy time-dependent version to keep intel compil...
void output(std::ostream &outfile)
Output with default number of plot points.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals with this element's contribution.
void output(FILE *file_pt)
C_style output with default number of plot points.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
unsigned self_test()
Self-test: Return 0 for OK.
AxisymFoepplvonKarmanPressureFctPt Airy_forcing_fct_pt
Pointer to Airy forcing function.
virtual double dshape_and_dtest_eulerian_at_knot_axisym_fvk(const unsigned &ipt, Shape &psi, DShape &dpsidr, Shape &test, DShape &dtestdr) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
virtual unsigned nodal_index_fvk(const unsigned &i=0) const
Return the index at which the i-th unknown value is stored. The default value, i, is appropriate for ...
void(* AxisymFoepplvonKarmanPressureFctPt)(const double &r, double &f)
Function pointer to pressure function fct(r,f(r)) – r is a Vector!
void output(FILE *file_pt, const unsigned &n_plot)
C-style output FE representation of soln: r,w at n_plot plot points.
void get_gradient_of_deflection(const Vector< double > &s, Vector< double > &gradient) const
Get gradient of deflection: gradient[i] = dw/dr_i *‍/.
void interpolated_stress(const Vector< double > &s, double &sigma_r_r, double &sigma_phi_phi) const
Compute in-plane stresses.
AxisymFoepplvonKarmanPressureFctPt & airy_forcing_fct_pt()
Access function: Pointer to Airy forcing function.
AxisymFoepplvonKarmanPressureFctPt & pressure_fct_pt()
Access function: Pointer to pressure function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output FE representation of soln: r,w,sigma_r_r,sigma_phi_phi at n_plot plot points.
AxisymFoepplvonKarmanEquations()
Constructor (must initialise the Pressure_fct_pt and Airy_forcing_fct_pt to null)....
virtual double dshape_and_dtest_eulerian_axisym_fvk(const Vector< double > &s, Shape &psi, DShape &dpsidr, Shape &test, DShape &dtestdr) const =0
Shape/test functions and derivs w.r.t. to global coords at local coord. s; return Jacobian of mapping...
bool Linear_bending_model
Flag which stores whether we are using a linear, pure bending model instead of the full non-linear Fo...
const double & eta() const
FvK parameter.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,w_exact at n_plot plot points.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void use_linear_bending_model()
Sets a flag to signify that we are solving the linear, pure bending equations, and pin all the nodal ...
double interpolated_w_fvk(const Vector< double > &s) const
Return FE representation of vertical displacement, w_fvk(s) at local coordinate s.
virtual void get_pressure_fvk(const unsigned &ipt, const double &r, double &pressure) const
Get pressure term at (Eulerian) position r. This function is virtual to allow overloading in multi-ph...
double *& eta_pt()
Pointer to FvK parameter.
AxisymFoepplvonKarmanEquations(const AxisymFoepplvonKarmanEquations &dummy)=delete
Broken copy constructor.
virtual void get_airy_forcing_fvk(const unsigned &ipt, const double &r, double &airy_forcing) const
Get Airy forcing term at (Eulerian) position r. This function is virtual to allow overloading in mult...
AxisymFoepplvonKarmanPressureFctPt airy_forcing_fct_pt() const
Access function: Pointer to Airy forcing function. Const version.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
This is a base class for all elements that require external sources (e.g. FSI, multi-domain problems ...
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
void set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point.
FSI Axisym FoepplvonKarmanElement elements are 1D Foeppl von Karman elements with isoparametric inter...
void output(std::ostream &outfile, const unsigned &n_plot)
Output FE representation of soln: r,w,dwdt,sigma_r_r,sigma_phi_phi at n_plot plot points.
void reset_after_external_interaction_geometric_fd()
Perform any auxiliary node update fcts of the adjacent fluid elements.
void node_update_adjacent_fluid_elements()
Update the nodal positions in all fluid elements that affect the traction on this element.
virtual ~FSIAxisymFoepplvonKarmanElement()
Empty virtual destructor.
double *& q_pt()
Return a pointer the ratio of stress scales used to non-dimensionalise the fluid and solid equations.
virtual Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object's shape depends on.
void update_before_external_interaction_geometric_fd()
Perform any auxiliary node update fcts of the adjacent fluid elements.
double * Q_pt
Pointer to the ratio, , of the stress used to non-dimensionalise the fluid stresses to the stress us...
void output_adjacent_fluid_elements(std::ostream &outfile, const unsigned &nplot)
Output adjacent fluid elements (for checking of fsi setup)
void dposition_dt(const Vector< double > &zeta, const unsigned &j, Vector< double > &drdt)
j-th time-derivative on object at current time: .
void output_integration_points(std::ostream &outfile)
Output integration points (for checking of fsi setup)
virtual void get_pressure_fvk(const unsigned &ipt, const double &r, double &pressure) const
Overload pressure term at (Eulerian) position r. Adds fluid traction to pressure imposed by "pressure...
void position(const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
Overloaded position function: Return 2D position vector: (r(zeta),z(zeta)) of material point whose "L...
void position(const Vector< double > &zeta, Vector< double > &r) const
Overloaded position function: Return 2D position vector: (r(zeta),z(zeta)) of material point whose "L...
const double & q() const
Return the ratio of the stress scales used to non-dimensionalise the fluid and elasticity equations.
void update_in_external_interaction_geometric_fd(const unsigned &i)
Perform any auxiliary node update fcts of the adjacent fluid elements.
virtual unsigned ngeom_data() const
How many items of Data does the shape of the object depend on? All nodal data.
void reset_in_external_interaction_geometric_fd(const unsigned &i)
Perform any auxiliary node update fcts of the adjacent fluid elements.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2179
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2597
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...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3992
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2214
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1763
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1967
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3152
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3190
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:3328
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1769
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:192
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:389
//////////////////////////////////////////////////////////////////// ////////////////////////////////...