foeppl_von_karman_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 //Header file for FoepplvonKarman elements
31 #ifndef OOMPH_FOEPPLVONKARMAN_ELEMENTS_HEADER
32 #define OOMPH_FOEPPLVONKARMAN_ELEMENTS_HEADER
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37  #include <oomph-lib-config.h>
38 #endif
39 
40 #include<sstream>
41 
42 //OOMPH-LIB headers
43 #include "../generic/projection.h"
44 #include "../generic/nodes.h"
45 #include "../generic/Qelements.h"
46 #include "../generic/oomph_utilities.h"
47 
48 
49 
50 namespace oomph
51 {
52 
53 //=============================================================
54 /// A class for all isoparametric elements that solve the
55 /// Foeppl von Karman equations.
56 /// \f[
57 /// \nabla^4 w - \eta Diamond^4(w,\phi) = p(x,y)
58 /// \f]
59 /// and
60 /// \f[
61 /// \nabla^4 \phi + \frac{1}{2} Diamond^4(w,w) = 0
62 /// \f]
63 /// This contains the generic maths. Shape functions, geometric
64 /// mapping etc. must get implemented in derived class.
65 //=============================================================
67 {
68 
69 public:
70 
71  /// \short Function pointer to pressure function fct(x,f(x)) --
72  /// x is a Vector!
74  const Vector<double>& x,
75  double& f);
76 
77  /// \short Constructor (must initialise the Pressure_fct_pt and
78  /// Airy_forcing_fct_pt to null). Also set physical parameters to their
79  /// default values. No volume constraint applied by default.
82  {
83  //Set all the physical constants to the default value (zero)
85  Linear_bending_model = false;
86 
87  // No volume constraint
89  }
90 
91  /// Broken copy constructor
93  {
94  BrokenCopy::broken_copy("FoepplvonKarmanEquations");
95  }
96 
97  /// Broken assignment operator
99  {
100  BrokenCopy::broken_assign("FoepplvonKarmanEquations");
101  }
102 
103  // Access functions for the physical constants
104 
105  /// Eta
106  const double &eta() const {return *Eta_pt;}
107 
108  /// Pointer to eta
109  double* &eta_pt() {return Eta_pt;}
110 
111 
112  /// \short Set Data value containing a single value which represents the
113  /// volume control pressure as external data for this element.
114  /// Only used for volume controlled problems in conjunction with
115  /// FoepplvonKarmanVolumeConstraintElement.
117  {
118 #ifdef PARANOID
119  if (data_pt->nvalue()!=1)
120  {
121  throw OomphLibError(
122  "Data object that contains volume control pressure should only contain a single value. ",
123  OOMPH_CURRENT_FUNCTION,
124  OOMPH_EXCEPTION_LOCATION);
125  }
126 #endif
127 
128  // Add as external data and remember the index in the element's storage
129  // scheme
131  }
132 
133  /// \short Return the index at which the i-th unknown value
134  /// is stored. The default value, i, is appropriate for single-physics
135  /// problems. By default, these are:
136  /// 0: w
137  /// 1: laplacian w
138  /// 2: phi
139  /// 3: laplacian phi
140  /// 4-8: smooth first derivatives
141  /// In derived multi-physics elements, this function should be overloaded
142  /// to reflect the chosen storage scheme. Note that these equations require
143  /// that the unknown is always stored at the same index at each node.
144  virtual inline unsigned nodal_index_fvk(const unsigned& i=0) const {return i;}
145 
146  /// Output with default number of plot points
147  void output(std::ostream &outfile)
148  {
149  const unsigned n_plot=5;
150  output(outfile,n_plot);
151  }
152 
153  /// \short Output FE representation of soln: x,y,w at
154  /// n_plot^DIM plot points
155  void output(std::ostream &outfile, const unsigned &n_plot);
156 
157  /// C_style output with default number of plot points
158  void output(FILE* file_pt)
159  {
160  const unsigned n_plot=5;
161  output(file_pt,n_plot);
162  }
163 
164  /// \short C-style output FE representation of soln: x,y,w at
165  /// n_plot^DIM plot points
166  void output(FILE* file_pt, const unsigned &n_plot);
167 
168  /// Output exact soln: x,y,w_exact at n_plot^DIM plot points
169  void output_fct(std::ostream &outfile, const unsigned &n_plot,
171 
172  /// \short Output exact soln: x,y,w_exact at
173  /// n_plot^DIM plot points (dummy time-dependent version to
174  /// keep intel compiler happy)
175  virtual void output_fct(std::ostream &outfile, const unsigned &n_plot,
176  const double& time,
178  exact_soln_pt)
179  {
180  throw OomphLibError(
181  "There is no time-dependent output_fct() for Foeppl von Karman"
182  "elements ",
183  OOMPH_CURRENT_FUNCTION,
184  OOMPH_EXCEPTION_LOCATION);
185  }
186 
187 
188  /// Get error against and norm of exact solution
189  void compute_error(std::ostream &outfile,
191  double& error, double& norm);
192 
193 
194  /// Dummy, time dependent error checker
195  void compute_error(std::ostream &outfile,
197  const double& time, double& error, double& norm)
198  {
199  throw OomphLibError(
200  "There is no time-dependent compute_error() for Foeppl von Karman"
201  "elements",
202  OOMPH_CURRENT_FUNCTION,
203  OOMPH_EXCEPTION_LOCATION);
204  }
205 
206  /// Access function: Pointer to pressure function
208 
209  /// Access function: Pointer to pressure function. Const version
211  {return Pressure_fct_pt;}
212 
213  /// Access function: Pointer to Airy forcing function
215  {return Airy_forcing_fct_pt;}
216 
217  /// Access function: Pointer to Airy forcing function. Const version
219  const {return Airy_forcing_fct_pt;}
220 
221  /// \short Get pressure term at (Eulerian) position x. This function is
222  /// virtual to allow overloading in multi-physics problems where
223  /// the strength of the pressure function might be determined by
224  /// another system of equations.
225  inline virtual void get_pressure_fvk(const unsigned& ipt,
226  const Vector<double>& x,
227  double& pressure) const
228  {
229  //If no pressure function has been set, return zero
230  if(Pressure_fct_pt==0) {pressure = 0.0;}
231  else
232  {
233  // Get pressure strength
234  (*Pressure_fct_pt)(x,pressure);
235  }
236  }
237 
238  /// \short Get Airy forcing term at (Eulerian) position x. This function is
239  /// virtual to allow overloading in multi-physics problems where
240  /// the strength of the pressure function might be determined by
241  /// another system of equations.
242  inline virtual void get_airy_forcing_fvk(const unsigned& ipt,
243  const Vector<double>& x,
244  double& airy_forcing) const
245  {
246  //If no pressure function has been set, return zero
247  if(Airy_forcing_fct_pt==0) {airy_forcing = 0.0;}
248  else
249  {
250  // Get pressure strength
251  (*Airy_forcing_fct_pt)(x,airy_forcing);
252  }
253  }
254 
255  /// Get gradient of deflection: gradient[i] = dw/dx_i
257  Vector<double>& gradient) const
258  {
259  //Find out how many nodes there are in the element
260  const unsigned n_node = nnode();
261 
262  //Get the index at which the unknown is stored
263  const unsigned w_nodal_index = nodal_index_fvk(0);
264 
265  //Set up memory for the shape and test functions
266  Shape psi(n_node);
267  DShape dpsidx(n_node,2);
268 
269  //Call the derivatives of the shape and test functions
270  dshape_eulerian(s,psi,dpsidx);
271 
272  //Initialise to zero
273  for(unsigned j=0;j<2;j++)
274  {
275  gradient[j] = 0.0;
276  }
277 
278  // Loop over nodes
279  for(unsigned l=0;l<n_node;l++)
280  {
281  //Loop over derivative directions
282  for(unsigned j=0;j<2;j++)
283  {
284  gradient[j] += this->nodal_value(l,w_nodal_index)*dpsidx(l,j);
285  }
286  }
287  }
288 
289  /// Fill in the residuals with this element's contribution
291 
292  //void fill_in_contribution_to_jacobian(Vector<double> &residuals,
293  // DenseMatrix<double> &jacobian);
294 
295  /// \short Return FE representation of function value w_fvk(s)
296  /// at local coordinate s (by default - if index > 0, returns
297  /// FE representation of valued stored at index^th nodal index
298  inline double interpolated_w_fvk(const Vector<double> &s,
299  unsigned index=0) const
300  {
301  //Find number of nodes
302  const unsigned n_node = nnode();
303 
304  //Get the index at which the poisson unknown is stored
305  const unsigned w_nodal_index = nodal_index_fvk(index);
306 
307  //Local shape function
308  Shape psi(n_node);
309 
310  //Find values of shape function
311  shape(s,psi);
312 
313  //Initialise value of u
314  double interpolated_w = 0.0;
315 
316  //Loop over the local nodes and sum
317  for(unsigned l=0;l<n_node;l++)
318  {
319  interpolated_w += this->nodal_value(l,w_nodal_index)*psi[l];
320  }
321 
322  return(interpolated_w);
323  }
324 
325  /// Compute in-plane stresses
327  double& sigma_xx,
328  double& sigma_yy,
329  double& sigma_xy);
330 
331  /// \short Return the integral of the displacement over the current
332  /// element, effectively calculating its contribution to the volume under
333  /// the membrane. Virtual so it can be overloaded in multi-physics
334  /// where the volume may incorporate an offset, say.
335  virtual double get_bounded_volume() const
336  {
337 
338  //Number of nodes and integration points for the current element
339  const unsigned n_node = nnode();
340  const unsigned n_intpt = integral_pt()->nweight();
341 
342  //Shape functions and their derivatives
343  Shape psi(n_node);
344  DShape dpsidx(n_node,2);
345 
346  //The nodal index at which the displacement is stored
347  const unsigned w_nodal_index = nodal_index_fvk();
348 
349  //Initalise the integral variable
350  double integral_w = 0;
351 
352  //Loop over the integration points
353  for(unsigned ipt=0;ipt<n_intpt;ipt++)
354  {
355  //Get the integral weight
356  double w = integral_pt()->weight(ipt);
357 
358  //Get determinant of the Jacobian of the mapping
359  double J = dshape_eulerian_at_knot(ipt,psi,dpsidx);
360 
361  //Premultiply the weight and Jacobian
362  double W = w*J;
363 
364  //Initialise storage for the w value and nodal value
365  double interpolated_w = 0;
366  double w_nodal_value;
367 
368  //Loop over the shape functions/nodes
369  for(unsigned l=0;l<n_node;l++)
370  {
371  //Get the current nodal value
372  w_nodal_value = raw_nodal_value(l,w_nodal_index);
373  //Add the contribution to interpolated w
374  interpolated_w += w_nodal_value*psi(l);
375  }
376 
377  //Add the contribution from the current integration point
378  integral_w += interpolated_w*W;
379  }
380 
381  //Return the calculated integral
382  return integral_w;
383  }
384 
385  /// \short Self-test: Return 0 for OK
386  unsigned self_test();
387 
388  /// \short Sets a flag to signify that we are solving the linear, pure bending
389  /// equations, and pin all the nodal values that will not be used in this case
391  {
392  // Set the boolean flag
393  Linear_bending_model = true;
394 
395  // Get the index of the first FvK nodal value
396  unsigned first_fvk_nodal_index = nodal_index_fvk();
397 
398  // Get the total number of FvK nodal values (assuming they are stored
399  // contiguously) at node 0 (it's the same at all nodes anyway)
400  unsigned total_fvk_nodal_indicies = 8;
401 
402  // Get the number of nodes in this element
403  unsigned n_node = nnode();
404 
405  // Loop over the appropriate nodal indices
406  for(unsigned index=first_fvk_nodal_index+2;
407  index<first_fvk_nodal_index+total_fvk_nodal_indicies;
408  index++)
409  {
410  // Loop over the nodes in the element
411  for(unsigned inod=0;inod<n_node;inod++)
412  {
413  // Pin the nodal value at the current index
414  node_pt(inod)->pin(index);
415  }
416  }
417  }
418 
419 
420 protected:
421 
422  /// \short Shape/test functions and derivs w.r.t. to global coords at
423  /// local coord. s; return Jacobian of mapping
424  virtual double dshape_and_dtest_eulerian_fvk(const Vector<double> &s,
425  Shape &psi,
426  DShape &dpsidx, Shape &test,
427  DShape &dtestdx) const=0;
428 
429 
430  /// \short Shape/test functions and derivs w.r.t. to global coords at
431  /// integration point ipt; return Jacobian of mapping
432  virtual double dshape_and_dtest_eulerian_at_knot_fvk(const unsigned &ipt,
433  Shape &psi,
434  DShape &dpsidx,
435  Shape &test,
436  DShape &dtestdx)
437  const=0;
438 
439  // Pointers to global physical constants
440 
441  /// Pointer to global eta
442  double *Eta_pt;
443 
444  /// Pointer to pressure function:
446 
447  //mjr Ridicu-hack for made-up second pressure, q
449 
450 private:
451 
452  /// Default value for physical constants
454 
455  /// \short Flag which stores whether we are using a linear, pure bending model
456  /// instead of the full non-linear Foeppl-von Karman
458 
459  /// \short Index of the external Data object that represents the volume
460  /// constraint pressure (initialised to -1 indicating no such constraint)
461  /// Gets overwritten when calling
462  /// set_volume_constraint_pressure_data_as_external_data(...)
464 
465 };
466 
467 
468 
469 
470 
471 
472 ///////////////////////////////////////////////////////////////////////////
473 ///////////////////////////////////////////////////////////////////////////
474 ///////////////////////////////////////////////////////////////////////////
475 
476 
477 // mjr QFoepplvonKarmanElements are untested!
478 
479 //======================================================================
480 /// QFoepplvonKarmanElement elements are linear/quadrilateral/brick-shaped
481 /// Foeppl von Karman elements with isoparametric interpolation for the
482 /// function.
483 //======================================================================
484 
485 template <unsigned NNODE_1D>
486  class QFoepplvonKarmanElement : public virtual QElement<2,NNODE_1D>,
487  public virtual FoepplvonKarmanEquations
488 {
489 
490 private:
491 
492  /// \short Static int that holds the number of variables at
493  /// nodes: always the same
494  static const unsigned Initial_Nvalue;
495 
496 public:
497 
498  ///\short Constructor: Call constructors for QElement and
499  /// FoepplvonKarmanEquations
502  {}
503 
504  /// Broken copy constructor
506  {
507  BrokenCopy::broken_copy("QFoepplvonKarmanElement");
508  }
509 
510  /// Broken assignment operator
512  {
513  BrokenCopy::broken_assign("QFoepplvonKarmanElement");
514  }
515 
516 
517  /// \short Required # of `values' (pinned or dofs)
518  /// at node n
519  inline unsigned required_nvalue(const unsigned &n) const
520  {return Initial_Nvalue;}
521 
522  /// \short Output function:
523  /// x,y,w
524  void output(std::ostream &outfile)
526 
527 
528  /// \short Output function:
529  /// x,y,w at n_plot^DIM plot points
530  void output(std::ostream &outfile, const unsigned &n_plot)
531  {FoepplvonKarmanEquations::output(outfile,n_plot);}
532 
533 
534  /// \short C-style output function:
535  /// x,y,w
536  void output(FILE* file_pt)
538 
539 
540  /// \short C-style output function:
541  /// x,y,w at n_plot^DIM plot points
542  void output(FILE* file_pt, const unsigned &n_plot)
543  {FoepplvonKarmanEquations::output(file_pt,n_plot);}
544 
545 
546  /// \short Output function for an exact solution:
547  /// x,y,w_exact at n_plot^DIM plot points
548  void output_fct(std::ostream &outfile, const unsigned &n_plot,
550  {FoepplvonKarmanEquations::output_fct(outfile,n_plot,exact_soln_pt);}
551 
552 
553 
554  /// \short Output function for a time-dependent exact solution.
555  /// x,y,w_exact at n_plot^DIM plot points
556  /// (Calls the steady version)
557  void output_fct(std::ostream &outfile, const unsigned &n_plot,
558  const double& time,
560  {FoepplvonKarmanEquations::output_fct(outfile,n_plot,time,exact_soln_pt);}
561 
562 
563 protected:
564 
565 /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
566  inline double dshape_and_dtest_eulerian_fvk(
567  const Vector<double> &s, Shape &psi, DShape &dpsidx,
568  Shape &test, DShape &dtestdx) const;
569 
570 
571  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
572  /// integration point ipt. Return Jacobian.
573  inline double dshape_and_dtest_eulerian_at_knot_fvk(const unsigned& ipt,
574  Shape &psi,
575  DShape &dpsidx,
576  Shape &test,
577  DShape &dtestdx)
578  const;
579 
580 };
581 
582 
583 
584 
585 //Inline functions:
586 
587 
588 //======================================================================
589 /// Define the shape functions and test functions and derivatives
590 /// w.r.t. global coordinates and return Jacobian of mapping.
591 ///
592 /// Galerkin: Test functions = shape functions
593 //======================================================================
594 template<unsigned NNODE_1D>
596  const Vector<double> &s,
597  Shape &psi,
598  DShape &dpsidx,
599  Shape &test,
600  DShape &dtestdx) const
601 {
602  //Call the geometrical shape functions and derivatives
603  const double J = this->dshape_eulerian(s,psi,dpsidx);
604 
605  //Set the test functions equal to the shape functions
606  test = psi;
607  dtestdx= dpsidx;
608 
609  //Return the jacobian
610  return J;
611 }
612 
613 
614 
615 
616 //======================================================================
617 /// Define the shape functions and test functions and derivatives
618 /// w.r.t. global coordinates and return Jacobian of mapping.
619 ///
620 /// Galerkin: Test functions = shape functions
621 //======================================================================
622 template<unsigned NNODE_1D>
625  const unsigned &ipt,
626  Shape &psi,
627  DShape &dpsidx,
628  Shape &test,
629  DShape &dtestdx) const
630 {
631  //Call the geometrical shape functions and derivatives
632  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
633 
634  //Set the pointers of the test functions
635  test = psi;
636  dtestdx = dpsidx;
637 
638  //Return the jacobian
639  return J;
640 }
641 
642 
643 ////////////////////////////////////////////////////////////////////////
644 ////////////////////////////////////////////////////////////////////////
645 ////////////////////////////////////////////////////////////////////////
646 
647 
648 //=======================================================================
649 /// Face geometry for the QFoepplvonKarmanElement elements: The spatial
650 /// dimension of the face elements is one lower than that of the
651 /// bulk element but they have the same number of points
652 /// along their 1D edges.
653 //=======================================================================
654 template<unsigned NNODE_1D>
656  public virtual QElement<1,NNODE_1D>
657 {
658 
659  public:
660 
661  /// \short Constructor: Call the constructor for the
662  /// appropriate lower-dimensional QElement
663  FaceGeometry() : QElement<1,NNODE_1D>() {}
664 
665 };
666 
667 ////////////////////////////////////////////////////////////////////////
668 ////////////////////////////////////////////////////////////////////////
669 ////////////////////////////////////////////////////////////////////////
670 
671 
672 //==========================================================
673 /// Foeppl von Karman upgraded to become projectable
674 //==========================================================
675  template<class FVK_ELEMENT>
677  public virtual ProjectableElement<FVK_ELEMENT>
678  {
679 
680  public:
681 
682  /// \short Specify the values associated with field fld.
683  /// The information is returned in a vector of pairs which comprise
684  /// the Data object and the value within it, that correspond to field fld.
686  {
687 
688 #ifdef PARANOID
689  if (fld > 7)
690  {
691  std::stringstream error_stream;
692  error_stream
693  << "Foeppl von Karman elements only store eight fields so fld must be"
694  << "0 to 7 rather than " << fld << std::endl;
695  throw OomphLibError(
696  error_stream.str(),
697  OOMPH_CURRENT_FUNCTION,
698  OOMPH_EXCEPTION_LOCATION);
699  }
700 #endif
701 
702  // Create the vector
703  unsigned nnod=this->nnode();
704  Vector<std::pair<Data*,unsigned> > data_values(nnod);
705 
706  // Loop over all nodes
707  for (unsigned j=0;j<nnod;j++)
708  {
709  // Add the data value associated field: The node itself
710  data_values[j]=std::make_pair(this->node_pt(j),fld);
711  }
712 
713  // Return the vector
714  return data_values;
715  }
716 
717  /// \short Number of fields to be projected: Just two
719  {
720  return 8;
721  }
722 
723  /// \short Number of history values to be stored for fld-th field.
724  /// (Note: count includes current value!)
725  unsigned nhistory_values_for_projection(const unsigned &fld)
726  {
727 #ifdef PARANOID
728  if (fld > 7)
729  {
730  std::stringstream error_stream;
731  error_stream
732  << "Foeppl von Karman elements only store eight fields so fld must be"
733  << "0 to 7 rather than " << fld << std::endl;
734  throw OomphLibError(
735  error_stream.str(),
736  OOMPH_CURRENT_FUNCTION,
737  OOMPH_EXCEPTION_LOCATION);
738  }
739 #endif
740  return this->node_pt(0)->ntstorage();
741  }
742 
743 
744  ///\short Number of positional history values
745  /// (Note: count includes current value!)
747  {
748  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
749  }
750 
751  /// \short Return Jacobian of mapping and shape functions of field fld
752  /// at local coordinate s
753  double jacobian_and_shape_of_field(const unsigned &fld,
754  const Vector<double> &s,
755  Shape &psi)
756  {
757 #ifdef PARANOID
758  if (fld > 7)
759  {
760  std::stringstream error_stream;
761  error_stream
762  << "Foeppl von Karman elements only store eight fields so fld must be"
763  << "0 to 7 rather than " << fld << std::endl;
764  throw OomphLibError(
765  error_stream.str(),
766  OOMPH_CURRENT_FUNCTION,
767  OOMPH_EXCEPTION_LOCATION);
768  }
769 #endif
770  unsigned n_dim=this->dim();
771  unsigned n_node=this->nnode();
772  Shape test(n_node);
773  DShape dpsidx(n_node,n_dim), dtestdx(n_node,n_dim);
774  double J=this->dshape_and_dtest_eulerian_fvk(s,psi,dpsidx,
775  test,dtestdx);
776  return J;
777  }
778 
779 
780 
781  /// \short Return interpolated field fld at local coordinate s, at time level
782  /// t (t=0: present; t>0: history values)
783  double get_field(const unsigned &t,
784  const unsigned &fld,
785  const Vector<double>& s)
786  {
787 
788 #ifdef PARANOID
789  if (fld > 7)
790  {
791  std::stringstream error_stream;
792  error_stream
793  << "Foeppl von Karman elements only store eight fields so fld must be"
794  << "0 to 7 rather than " << fld << std::endl;
795  throw OomphLibError(
796  error_stream.str(),
797  OOMPH_CURRENT_FUNCTION,
798  OOMPH_EXCEPTION_LOCATION);
799  }
800 #endif
801  //Find the index at which the variable is stored
802  unsigned w_nodal_index = this->nodal_index_fvk(fld);
803 
804  //Local shape function
805  unsigned n_node=this->nnode();
806  Shape psi(n_node);
807 
808  //Find values of shape function
809  this->shape(s,psi);
810 
811  //Initialise value of u
812  double interpolated_w = 0.0;
813 
814  //Sum over the local nodes
815  for(unsigned l=0;l<n_node;l++)
816  {
817  interpolated_w += this->nodal_value(t,l,w_nodal_index)*psi[l];
818  }
819  return interpolated_w;
820  }
821 
822 
823 
824 
825  ///Return number of values in field fld: One per node
826  unsigned nvalue_of_field(const unsigned &fld)
827  {
828 #ifdef PARANOID
829  if (fld > 7)
830  {
831  std::stringstream error_stream;
832  error_stream
833  << "Foeppl von Karman elements only store eight fields so fld must be"
834  << "0 to 7 rather than " << fld << std::endl;
835  throw OomphLibError(
836  error_stream.str(),
837  OOMPH_CURRENT_FUNCTION,
838  OOMPH_EXCEPTION_LOCATION);
839  }
840 #endif
841  return this->nnode();
842  }
843 
844 
845  ///Return local equation number of field fld of node j.
846  int local_equation(const unsigned &fld,
847  const unsigned &j)
848  {
849 #ifdef PARANOID
850  if (fld > 7)
851  {
852  std::stringstream error_stream;
853  error_stream
854  << "Foeppl von Karman elements only store eight fields so fld must be"
855  << "0 to 7 rather than " << fld << std::endl;
856  throw OomphLibError(
857  error_stream.str(),
858  OOMPH_CURRENT_FUNCTION,
859  OOMPH_EXCEPTION_LOCATION);
860  }
861 #endif
862  const unsigned w_nodal_index = this->nodal_index_fvk(fld);
863  return this->nodal_local_eqn(j,w_nodal_index);
864  }
865 
866  };
867 
868 //=======================================================================
869 /// Face geometry for element is the same as that for the underlying
870 /// wrapped element
871 //=======================================================================
872  template<class ELEMENT>
874  : public virtual FaceGeometry<ELEMENT>
875  {
876  public:
877  FaceGeometry() : FaceGeometry<ELEMENT>() {}
878  };
879 
880 
881 //=======================================================================
882 /// Face geometry of the Face Geometry for element is the same as
883 /// that for the underlying wrapped element
884 //=======================================================================
885  template<class ELEMENT>
887  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
888  {
889  public:
891  };
892 
893 }
894 
895 #endif
896 
FoepplvonKarmanPressureFctPt & pressure_fct_pt()
Access function: Pointer to pressure function.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,w_exact at n_plot^DIM plot points.
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
static double Default_Physical_Constant_Value
Default value for physical constants.
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!) ...
FoepplvonKarmanPressureFctPt pressure_fct_pt() const
Access function: Pointer to pressure function. Const version.
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:847
bool Linear_bending_model
Flag which stores whether we are using a linear, pure bending model instead of the full non-linear Fo...
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 operator=(const FoepplvonKarmanEquations &)
Broken assignment operator.
void output(std::ostream &outfile)
Output function: x,y,w.
FoepplvonKarmanEquations()
Constructor (must initialise the Pressure_fct_pt and Airy_forcing_fct_pt to null). Also set physical parameters to their default values. No volume constraint applied by default.
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
A general Finite Element class.
Definition: elements.h:1271
char t
Definition: cfortran.h:572
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
unsigned nfields_for_projection()
Number of fields to be projected: Just two.
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
FoepplvonKarmanPressureFctPt & airy_forcing_fct_pt()
Access function: Pointer to Airy forcing function.
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
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
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Definition: elements.cc:293
int Volume_constraint_pressure_external_data_index
Index of the external Data object that represents the volume constraint pressure (initialised to -1 i...
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
FoepplvonKarmanPressureFctPt Pressure_fct_pt
Pointer to pressure function:
QFoepplvonKarmanElement(const QFoepplvonKarmanElement< NNODE_1D > &dummy)
Broken copy constructor.
void interpolated_stress(const Vector< double > &s, double &sigma_xx, double &sigma_yy, double &sigma_xy)
Compute in-plane stresses.
double * Eta_pt
Pointer to global eta.
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
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.
virtual void get_airy_forcing_fvk(const unsigned &ipt, const Vector< double > &x, double &airy_forcing) const
Get Airy forcing term at (Eulerian) position x. This function is virtual to allow overloading in mult...
void use_linear_bending_model()
Sets a flag to signify that we are solving the linear, pure bending equations, and pin all the nodal ...
FoepplvonKarmanEquations(const FoepplvonKarmanEquations &dummy)
Broken copy constructor.
static char t char * s
Definition: cfortran.h:572
unsigned self_test()
Self-test: Return 0 for OK.
void(* FoepplvonKarmanPressureFctPt)(const Vector< double > &x, double &f)
Function pointer to pressure function fct(x,f(x)) – x is a Vector!
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
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 fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals with this element's contribution.
virtual double dshape_and_dtest_eulerian_at_knot_fvk(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 ...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void operator=(const QFoepplvonKarmanElement< NNODE_1D > &)
Broken assignment operator.
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
FoepplvonKarmanPressureFctPt Airy_forcing_fct_pt
void output(std::ostream &outfile)
Output with default number of plot points.
QFoepplvonKarmanElement()
Constructor: Call constructors for QElement and FoepplvonKarmanEquations.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,w at n_plot^DIM plot points.
double dshape_and_dtest_eulerian_at_knot_fvk(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt. Return Jacobian.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of field fld of node j.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
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 ...
FoepplvonKarmanPressureFctPt airy_forcing_fct_pt() const
Access function: Pointer to Airy forcing function. Const version.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2446
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
Foeppl von Karman upgraded to become projectable.
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
virtual double get_bounded_volume() const
Return the integral of the displacement over the current element, effectively calculating its contrib...
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,w_exact at n_plot^DIM plot points (Calls the...
void get_gradient_of_deflection(const Vector< double > &s, Vector< double > &gradient) const
Get gradient of deflection: gradient[i] = dw/dx_i.
double dshape_and_dtest_eulerian_fvk(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.
virtual void get_pressure_fvk(const unsigned &ipt, const Vector< double > &x, double &pressure) const
Get pressure term at (Eulerian) position x. This function is virtual to allow overloading in multi-ph...
void set_volume_constraint_pressure_data_as_external_data(Data *data_pt)
Set Data value containing a single value which represents the volume control pressure as external dat...
virtual void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,w_exact at n_plot^DIM plot points (dummy time-dependent version to keep intel ...
virtual double dshape_and_dtest_eulerian_fvk(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...
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...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: x,y,w at n_plot^DIM plot points.
void output(FILE *file_pt)
C-style output function: x,y,w.
double interpolated_w_fvk(const Vector< double > &s, unsigned index=0) const
Return FE representation of function value w_fvk(s) at local coordinate s (by default - if index > 0...
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: x,y,w_exact at n_plot^DIM plot points.
unsigned required_nvalue(const unsigned &n) const
Required # of `values' (pinned or dofs) at node n.
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3252