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-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// 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
41namespace 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)
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. Return boolean to indicate success
296 /// (false if attempt to evaluate stresses at zero radius)
298 double& sigma_r_r,
299 double& sigma_phi_phi);
300
301
302 /// Self-test: Return 0 for OK
303 unsigned self_test();
304
305 /// Sets a flag to signify that we are solving the linear,
306 /// pure bending equations, and pin all the nodal values that will
307 /// not be used in this case
309 {
310 // Set the boolean flag
312
313 // Get the index of the first FvK nodal value
314 unsigned first_fvk_nodal_index = nodal_index_fvk();
315
316 // Get the total number of FvK nodal values (assuming they are stored
317 // contiguously) at node 0 (it's the same at all nodes anyway)
318 unsigned total_fvk_nodal_indices = 6;
319
320 // Get the number of nodes in this element
321 unsigned n_node = nnode();
322
323 // Loop over the appropriate nodal indices
324 for (unsigned index = first_fvk_nodal_index + 2;
325 index < first_fvk_nodal_index + total_fvk_nodal_indices;
326 index++)
327 {
328 // Loop over the nodes in the element
329 for (unsigned inod = 0; inod < n_node; inod++)
330 {
331 // Pin the nodal value at the current index
332 node_pt(inod)->pin(index);
333 }
334 }
335 }
336
337
338 protected:
339 /// Shape/test functions and derivs w.r.t. to global coords at
340 /// local coord. s; return Jacobian of mapping
342 const Vector<double>& s,
343 Shape& psi,
344 DShape& dpsidr,
345 Shape& test,
346 DShape& dtestdr) const = 0;
347
348
349 /// Shape/test functions and derivs w.r.t. to global coords at
350 /// integration point ipt; return Jacobian of mapping
352 const unsigned& ipt,
353 Shape& psi,
354 DShape& dpsidr,
355 Shape& test,
356 DShape& dtestdr) const = 0;
357
358 /// Pointer to FvK parameter
359 double* Eta_pt;
360
361 /// Pointer to pressure function:
363
364 /// Pointer to Airy forcing function
366
367 /// Default value for physical constants
369
370 /// Flag which stores whether we are using a linear,
371 /// pure bending model instead of the full non-linear Foeppl-von Karman
373 };
374
375
376 /// ////////////////////////////////////////////////////////////////////////
377 /// ////////////////////////////////////////////////////////////////////////
378 /// ////////////////////////////////////////////////////////////////////////
379
380
381 //======================================================================
382 /// Axisym FoepplvonKarmanElement elements are 1D
383 /// Foeppl von Karman elements with isoparametric interpolation for the
384 /// function.
385 //======================================================================
386 template<unsigned NNODE_1D>
388 : public virtual QElement<1, NNODE_1D>,
389 public virtual AxisymFoepplvonKarmanEquations
390 {
391 private:
392 /// Static int that holds the number of variables at
393 /// nodes: always the same
394 static const unsigned Initial_Nvalue;
395
396 public:
397 /// Constructor: Call constructors for QElement and
398 /// AxisymFoepplvonKarmanEquations
400 : QElement<1, NNODE_1D>(), AxisymFoepplvonKarmanEquations()
401 {
402 }
403
404 /// Broken copy constructor
406 const AxisymFoepplvonKarmanElement<NNODE_1D>& dummy) = delete;
407
408 /// Broken assignment operator
410
411 /// Required # of `values' (pinned or dofs)
412 /// at node n
413 inline unsigned required_nvalue(const unsigned& n) const
414 {
415 return Initial_Nvalue;
416 }
417
418
419 /// Output function:
420 /// r,w,sigma_r_r,sigma_phi_phi
421 void output(std::ostream& outfile)
422 {
424 }
425
426 /// Output function:
427 /// r,w,sigma_r_r,sigma_phi_phi at n_plot plot points
428 void output(std::ostream& outfile, const unsigned& n_plot)
429 {
431 }
432
433 /// C-style output function:
434 /// r,w
435 void output(FILE* file_pt)
436 {
438 }
439
440 /// C-style output function:
441 /// r,w at n_plot plot points
442 void output(FILE* file_pt, const unsigned& n_plot)
443 {
445 }
446
447 /// Output function for an exact solution:
448 /// r,w_exact at n_plot plot points
449 void output_fct(std::ostream& outfile,
450 const unsigned& n_plot,
452 {
454 outfile, n_plot, exact_soln_pt);
455 }
456
457 /// Output function for a time-dependent exact solution.
458 /// r,w_exact at n_plot plot points
459 /// (Calls the steady version)
460 void output_fct(std::ostream& outfile,
461 const unsigned& n_plot,
462 const double& time,
464 {
466 outfile, n_plot, time, exact_soln_pt);
467 }
468
469
470 protected:
471 /// Shape, test functions & derivs. w.r.t. to global coords.
472 /// Return Jacobian.
474 Shape& psi,
475 DShape& dpsidr,
476 Shape& test,
477 DShape& dtestdr) const;
478
479 /// Shape, test functions & derivs. w.r.t. to global coords. at
480 /// integration point ipt. Return Jacobian.
482 const unsigned& ipt,
483 Shape& psi,
484 DShape& dpsidr,
485 Shape& test,
486 DShape& dtestdr) const;
487 };
488
489
490 // Inline functions:
491
492 //======================================================================
493 /// Define the shape functions and test functions and derivatives
494 /// w.r.t. global coordinates and return Jacobian of mapping.
495 ///
496 /// Galerkin: Test functions = shape functions
497 //======================================================================
498 template<unsigned NNODE_1D>
500 NNODE_1D>::dshape_and_dtest_eulerian_axisym_fvk(const Vector<double>& s,
501 Shape& psi,
502 DShape& dpsidr,
503 Shape& test,
504 DShape& dtestdr) const
505
506 {
507 // Call the geometrical shape functions and derivatives
508 const double J = this->dshape_eulerian(s, psi, dpsidr);
509
510 // Set the test functions equal to the shape functions
511 test = psi;
512 dtestdr = dpsidr;
513
514 // Return the jacobian
515 return J;
516 }
517
518
519 //======================================================================
520 /// Define the shape functions and test functions and derivatives
521 /// w.r.t. global coordinates and return Jacobian of mapping.
522 ///
523 /// Galerkin: Test functions = shape functions
524 //======================================================================
525 template<unsigned NNODE_1D>
528 Shape& psi,
529 DShape& dpsidr,
530 Shape& test,
531 DShape& dtestdr) const
532 {
533 // Call the geometrical shape functions and derivatives
534 const double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidr);
535
536 // Set the pointers of the test functions
537 test = psi;
538 dtestdr = dpsidr;
539
540 // Return the jacobian
541 return J;
542 }
543
544
545 /// ////////////////////////////////////////////////////////////////////////
546 /// ////////////////////////////////////////////////////////////////////////
547 /// ////////////////////////////////////////////////////////////////////////
548
549
550 //======================================================================
551 /// FSI Axisym FoepplvonKarmanElement elements are 1D
552 /// Foeppl von Karman elements with isoparametric interpolation for the
553 /// function. Gets traction from adjacent fluid element(s) of type
554 /// FLUID_ELEMENT.
555 //======================================================================
556 template<unsigned NNODE_1D, class FLUID_ELEMENT>
558 : public virtual AxisymFoepplvonKarmanElement<NNODE_1D>,
559 public virtual ElementWithExternalElement
560 {
561 public:
562 /// Constructor
565 {
566 // Set source element storage: one interaction with an external
567 // element -- the fluid bulk element that provides the pressure
568 this->set_ninteraction(1);
569 }
570
571 /// Empty virtual destructor
573
574 /// Return the ratio of the stress scales used to non-dimensionalise
575 /// the fluid and elasticity equations.
576 const double& q() const
577 {
578 return *Q_pt;
579 }
580
581 /// Return a pointer the ratio of stress scales used to
582 /// non-dimensionalise the fluid and solid equations.
583 double*& q_pt()
584 {
585 return Q_pt;
586 }
587
588 /// How many items of Data does the shape of the object depend on?
589 /// All nodal data
590 virtual unsigned ngeom_data() const
591 {
592 return this->nnode();
593 }
594
595 /// Return pointer to the j-th Data item that the object's
596 /// shape depends on.
597 virtual Data* geom_data_pt(const unsigned& j)
598 {
599 return this->node_pt(j);
600 }
601
602 /// Overloaded position function: Return 2D position vector:
603 /// (r(zeta),z(zeta)) of material point whose "Lagrangian coordinate"
604 /// is given by zeta. Here r=zeta!
605 void position(const Vector<double>& zeta, Vector<double>& r) const
606 {
607 const unsigned t = 0;
608 this->position(t, zeta, r);
609 }
610
611 /// Overloaded position function: Return 2D position vector:
612 /// (r(zeta),z(zeta)) of material point whose "Lagrangian coordinate"
613 /// is given by zeta.
614 void position(const unsigned& t,
615 const Vector<double>& zeta,
616 Vector<double>& r) const
617 {
618 // Find number of nodes
619 const unsigned n_node = this->nnode();
620
621 // Get the index at which the poisson unknown is stored
622 const unsigned w_nodal_index = this->nodal_index_fvk(0);
623
624 // Local shape function
625 Shape psi(n_node);
626
627 // Find values of shape function
628 this->shape(zeta, psi);
629
630 // Initialise
631 double interpolated_w = 0.0;
632 double interpolated_r = 0.0;
633
634 // Loop over the local nodes and sum
635 for (unsigned l = 0; l < n_node; l++)
636 {
637 interpolated_w += this->nodal_value(t, l, w_nodal_index) * psi[l];
638 interpolated_r += this->node_pt(l)->x(t, 0) * psi[l];
639 }
640
641 // Assign
642 r[0] = interpolated_r;
643 r[1] = interpolated_w;
644 }
645
646 /// j-th time-derivative on object at current time:
647 /// \f$ \frac{d^{j} r(\zeta)}{dt^j} \f$.
648 void dposition_dt(const Vector<double>& zeta,
649 const unsigned& j,
650 Vector<double>& drdt)
651 {
652 // Find number of nodes
653 const unsigned n_node = this->nnode();
654
655 // Get the index at which the poisson unknown is stored
656 const unsigned w_nodal_index = this->nodal_index_fvk(0);
657
658 // Local shape function
659 Shape psi(n_node);
660
661 // Find values of shape function
662 this->shape(zeta, psi);
663
664 // Initialise
665 double interpolated_dwdt = 0.0;
666 double interpolated_drdt = 0.0;
667
668 // Loop over the local nodes and sum
669 for (unsigned l = 0; l < n_node; l++)
670 {
671 // Get the timestepper
673
674 // If we are doing an unsteady solve then calculate the derivative
676 {
677 // Get the number of values required to represent history
678 const unsigned n_time = time_stepper_pt->ntstorage();
679
680 // Loop over history values
681 for (unsigned t = 0; t < n_time; t++)
682 {
683 // Add the contribution to the derivative
684 interpolated_dwdt += time_stepper_pt->weight(1, t) *
685 this->nodal_value(t, l, w_nodal_index) *
686 psi[l];
687 }
688 }
689 }
690
691 // Assign
692 drdt[0] = interpolated_drdt;
693 drdt[1] = interpolated_dwdt;
694 }
695
696
697 /// Overload pressure term at (Eulerian) position r.
698 /// Adds fluid traction to pressure imposed by "pressure fct pointer"
699 /// (which can be regarded as applying an external (i.e.
700 /// "on the other side" of the fluid) pressure
701 inline virtual void get_pressure_fvk(const unsigned& ipt,
702 const double& r,
703 double& pressure) const
704 {
705 pressure = 0.0;
706
707 // Get underlying version
709
710 // Get FSI parameter
711 const double q_value = q();
712
713 // Get fluid element
714 FLUID_ELEMENT* ext_el_pt =
715 dynamic_cast<FLUID_ELEMENT*>(external_element_pt(0, ipt));
717
718 // Outer unit normal is vertically upward (in z direction)
719 // (within an essentiall flat) model for the wall)
720 Vector<double> normal(2);
721 normal[0] = 0.0;
722 normal[1] = 1.0;
723
724 // Get traction
725 Vector<double> traction(3);
726 ext_el_pt->traction(s_ext, normal, traction);
727
728 // Add z-component of traction
729 pressure -= q_value * traction[1];
730 }
731
732
733 /// Output integration points (for checking of fsi setup)
734 void output_integration_points(std::ostream& outfile)
735 {
736 // How many nodes do we have?
737 unsigned nnod = this->nnode();
738 Shape psi(nnod);
739
740 // Get the index at which the unknown is stored
741 const unsigned w_nodal_index = this->nodal_index_fvk(0);
742
743 // Loop over the integration points
744 const unsigned n_intpt = this->integral_pt()->nweight();
745 outfile << "ZONE I=" << n_intpt << std::endl;
746 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
747 {
748 // Get shape fct
749 Vector<double> s(1);
750 s[0] = this->integral_pt()->knot(ipt, 0);
751 shape(s, psi);
752
753 // Initialise
754 double interpolated_w = 0.0;
755 double interpolated_r = 0.0;
756
757 // Loop over the local nodes and sum
758 for (unsigned l = 0; l < nnod; l++)
759 {
760 interpolated_w += this->nodal_value(l, w_nodal_index) * psi[l];
761 interpolated_r += this->node_pt(l)->x(0) * psi[l];
762 }
763
764 // Get fluid element
765 FLUID_ELEMENT* ext_el_pt =
766 dynamic_cast<FLUID_ELEMENT*>(external_element_pt(0, ipt));
768
769 // Get veloc
770 Vector<double> veloc(3);
771 ext_el_pt->interpolated_u_axi_nst(s_ext, veloc);
772 Vector<double> x(2);
773 ext_el_pt->interpolated_x(s_ext, x);
774
775 outfile << interpolated_r << " " << interpolated_w << " " << veloc[0]
776 << " " << veloc[1] << " " << x[0] << " " << x[1] << " "
777 << std::endl;
778 }
779 }
780
781
782 /// Output adjacent fluid elements (for checking of fsi setup)
783 void output_adjacent_fluid_elements(std::ostream& outfile,
784 const unsigned& nplot)
785 {
786 // Loop over the integration points
787 const unsigned n_intpt = this->integral_pt()->nweight();
788 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
789 {
790 // Get fluid element
791 FLUID_ELEMENT* ext_el_pt =
792 dynamic_cast<FLUID_ELEMENT*>(external_element_pt(0, ipt));
793
794 // Dump it
795 ext_el_pt->output(outfile, nplot);
796 }
797 }
798
799 /// Perform any auxiliary node update fcts of the adjacent
800 /// fluid elements
802 {
804 }
805
806 /// Perform any auxiliary node update fcts of the adjacent
807 /// fluid elements
809 {
811 }
812
813
814 /// Perform any auxiliary node update fcts of the adjacent
815 /// fluid elements
817 {
819 }
820
821
822 /// Perform any auxiliary node update fcts of the adjacent
823 /// fluid elements
825 {
827 }
828
829
830 /// Update the nodal positions in all fluid elements that affect
831 /// the traction on this element
833 {
834 // Don't update elements repeatedly
835 std::map<FLUID_ELEMENT*, bool> done;
836
837 // Number of integration points
838 unsigned n_intpt = integral_pt()->nweight();
839
840 // Loop over all integration points in wall element
841 for (unsigned iint = 0; iint < n_intpt; iint++)
842 {
843 // Get fluid element that affects this integration point
844 FLUID_ELEMENT* el_f_pt =
845 dynamic_cast<FLUID_ELEMENT*>(external_element_pt(0, iint));
846
847 // Is there an adjacent fluid element?
848 if (el_f_pt != 0)
849 {
850 // Have we updated its positions yet?
851 if (!done[el_f_pt])
852 {
853 // Update nodal positions
854 el_f_pt->node_update();
855 done[el_f_pt] = true;
856 }
857 }
858 }
859 }
860
861
862 /// Output FE representation of soln:
863 /// r,w,dwdt,sigma_r_r,sigma_phi_phi at n_plot plot points
864 void output(std::ostream& outfile, const unsigned& n_plot)
865 {
866 // Vector of local coordinates
867 Vector<double> s(1);
868
869 // Tecplot header info
870 outfile << "ZONE\n";
871
872 // Loop over plot points
873 unsigned num_plot_points = nplot_points(n_plot);
874 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
875 {
876 // Get local coordinates of plot point
877 get_s_plot(iplot, n_plot, s);
878
879 // Get velocity
880 Vector<double> drdt(2);
881 dposition_dt(s, 1, drdt);
882
883 // Get stress
884 double sigma_r_r = 0.0;
885 double sigma_phi_phi = 0.0;
886 bool success = this->interpolated_stress(s, sigma_r_r, sigma_phi_phi);
887 if (success)
888 {
889 outfile << this->interpolated_x(s, 0) << " "
890 << this->interpolated_w_fvk(s) << " " << drdt[0] << " "
891 << drdt[1] << " " << sigma_r_r << " " << sigma_phi_phi
892 << std::endl;
893 }
894 }
895 }
896
897
898 protected:
899 /// Pointer to the ratio, \f$ Q \f$ , of the stress used to
900 /// non-dimensionalise the fluid stresses to the stress used to
901 /// non-dimensionalise the solid stresses.
902 double* Q_pt;
903 };
904
905
906} // namespace oomph
907
908#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...
const double & eta() const
FvK parameter.
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.
AxisymFoepplvonKarmanPressureFctPt & airy_forcing_fct_pt()
Access function: Pointer to Airy forcing function.
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.
bool interpolated_stress(const Vector< double > &s, double &sigma_r_r, double &sigma_phi_phi) const
Compute in-plane stresses. Return boolean to indicate success (false if attempt to evaluate stresses ...
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.
double *& eta_pt()
Pointer to FvK parameter.
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 *‍/.
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...
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...
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
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
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.
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
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.
const double & q() const
Return the ratio of the stress scales used to non-dimensionalise the fluid and elasticity equations.
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...
virtual Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object's shape depends on.
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: .
double *& q_pt()
Return a pointer the ratio of stress scales used to non-dimensionalise the fluid and solid equations.
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...
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.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual void 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:3962
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
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:3148
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:3186
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3298
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1765
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...