foeppl_von_karman_displacement_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_DISPLACEMENT_ELEMENTS_HEADER
32 #define OOMPH_FOEPPLVONKARMAN_DISPLACEMENT_ELEMENTS_HEADER
33 
34 // Config header generated by autoconfig
35 #ifdef HAVE_CONFIG_H
36  #include <oomph-lib-config.h>
37 #endif
38 
39 #include<sstream>
40 
41 //OOMPH-LIB headers
42 #include "../generic/projection.h"
43 #include "../generic/nodes.h"
44 #include "../generic/oomph_utilities.h"
45 
46 
47 
48 namespace oomph
49 {
50 
51 //=============================================================
52 /// A class for all isoparametric elements that solve the
53 /// Foeppl von Karman equations.
54 /// \f[
55 /// \nabla^4 w - 12 \frac{\partial}{\partial x_{\beta}}
56 /// \left( \sigma^{\alpha \beta}
57 /// \frac{\partial w}{\partial x_{\alpha}}
58 /// \right) = p(x,y)
59 /// \f]
60 /// and
61 /// \f[
62 /// \frac{\sigma^{\alpha \beta}}{\partial x_{\beta}} = \tau_\alpha
63 /// \f]
64 /// This contains the generic maths. Shape functions, geometric
65 /// mapping etc. must get implemented in derived class.
66 //=============================================================
68 {
69 
70 public:
71 
72  /// \short Function pointer to pressure function fct(x,f(x)) --
73  /// x is a Vector!
75  const Vector<double>& x, double& f);
76 
77  /// \short Function pointer to traction function fct(x,f(x)) --
78  /// x is a Vector!
79  typedef void (*FoepplvonKarmanTractionFctPt)(const Vector<double>& x,
80  Vector<double>& f);
81 
82  /// \short Constructor (must initialise the Pressure_fct_pt and the
83  /// Traction_fct_pt. Also set physical parameters to their default
84  /// values.
87  {
88  //Set default
90 
91  // Linear bending model?
92  Linear_bending_model = false;
93  }
94 
95  /// Broken copy constructor
97  {
98  BrokenCopy::broken_copy("FoepplvonKarmanDisplacementEquations");
99  }
100 
101  /// Broken assignment operator
103  {
104  BrokenCopy::broken_assign("FoepplvonKarmanDisplacementEquations");
105  }
106 
107  /// Poisson's ratio
108  const double &nu() const {return *Nu_pt;}
109 
110  /// Pointer to Poisson's ratio
111  double* &nu_pt() {return Nu_pt;}
112 
113  /// \short Return the index at which the i-th unknown value
114  /// is stored. The default value, i, is appropriate for single-physics
115  /// problems. By default, these are:
116  /// 0: w
117  /// 1: laplacian w
118  /// 2: u_x
119  /// 3: u_y
120  /// In derived multi-physics elements, this function should be overloaded
121  /// to reflect the chosen storage scheme. Note that these equations require
122  /// that the unknown is always stored at the same index at each node.
123  virtual inline unsigned nodal_index_fvk(const unsigned& i=0) const {return i;}
124 
125  /// Output with default number of plot points
126  void output(std::ostream &outfile)
127  {
128  const unsigned n_plot=5;
129  output(outfile,n_plot);
130  }
131 
132  /// \short Output FE representation of soln: x,y,w at
133  /// n_plot^DIM plot points
134  void output(std::ostream &outfile, const unsigned &n_plot);
135 
136  /// C_style output with default number of plot points
137  void output(FILE* file_pt)
138  {
139  const unsigned n_plot=5;
140  output(file_pt,n_plot);
141  }
142 
143  /// \short C-style output FE representation of soln: x,y,w at
144  /// n_plot^DIM plot points
145  void output(FILE* file_pt, const unsigned &n_plot);
146 
147  /// Output exact soln: x,y,w_exact at n_plot^DIM plot points
148  void output_fct(std::ostream &outfile, const unsigned &n_plot,
150 
151  /// \short Output exact soln: x,y,w_exact at
152  /// n_plot^DIM plot points (dummy time-dependent version to
153  /// keep intel compiler happy)
154  virtual void output_fct(std::ostream &outfile, const unsigned &n_plot,
155  const double& time,
157  exact_soln_pt)
158  {
159  throw OomphLibError(
160  "There is no time-dependent output_fct() for Foeppl von Karman"
161  "elements ",
162  OOMPH_CURRENT_FUNCTION,
163  OOMPH_EXCEPTION_LOCATION);
164  }
165 
166 
167  /// Get error against and norm of exact solution
168  void compute_error(std::ostream &outfile,
170  double& error, double& norm);
171 
172 
173  /// Dummy, time dependent error checker
174  void compute_error(std::ostream &outfile,
176  const double& time, double& error, double& norm)
177  {
178  throw OomphLibError(
179  "There is no time-dependent compute_error() for Foeppl von Karman"
180  "elements",
181  OOMPH_CURRENT_FUNCTION,
182  OOMPH_EXCEPTION_LOCATION);
183  }
184 
185  /// Access function: Pointer to pressure function
187 
188  /// Access function: Pointer to pressure function. Const version
190  {return Pressure_fct_pt;}
191 
192  /// Access function: Pointer to in-plane traction function
194 
195  /// Access function: Pointer to in-plane traction function. Const version
197  {return Traction_fct_pt;}
198 
199  /// \short Get pressure term at (Eulerian) position x. This function
200  /// is virtual to allow overloading in multi-physics problems where
201  /// the strength of the pressure function might be determined by
202  /// another system of equations.
203  inline virtual void get_pressure_fvk(const unsigned& ipt,
204  const Vector<double>& x,
205  double& pressure) const
206  {
207  //If no pressure function has been set, return zero
208  if(Pressure_fct_pt==0) {pressure = 0.0;}
209  else
210  {
211  // Get pressure strength
212  (*Pressure_fct_pt)(x,pressure);
213  }
214  }
215 
216  /// \short Get in-plane traction term at (Eulerian) position x.
217  inline virtual void get_traction_fvk(Vector<double>& x,
218  Vector<double>& traction) const
219  {
220  //If no pressure function has been set, return zero
221  if(Traction_fct_pt==0)
222  {
223  traction[0] = 0.0;
224  traction[1] = 0.0;
225  }
226  else
227  {
228  // Get traction
229  (*Traction_fct_pt)(x,traction);
230  }
231  }
232 
233  /// Get gradient of deflection: gradient[i] = dw/dx_i
235  Vector<double>& gradient) const
236  {
237  //Find out how many nodes there are in the element
238  const unsigned n_node = nnode();
239 
240  //Get the index at which the unknown is stored
241  // Indexes for unknows are
242  // 0: W (vertical displacement)
243  // 1: L (laplacian of W)
244  // 2: Ux (In plane displacement x)
245  // 3: Uy (In plane displacement y)
246  unsigned w_nodal_index = nodal_index_fvk(0);
247 
248  //Set up memory for the shape and test functions
249  Shape psi(n_node);
250  DShape dpsidx(n_node,2);
251 
252  //Call the derivatives of the shape and test functions
253  dshape_eulerian(s,psi,dpsidx);
254 
255  //Initialise to zero
256  for(unsigned j=0;j<2;j++)
257  {
258  gradient[j] = 0.0;
259  }
260 
261  // Loop over nodes
262  for(unsigned l=0;l<n_node;l++)
263  {
264  //Loop over derivative directions
265  for(unsigned j=0;j<2;j++)
266  {
267  gradient[j] += this->nodal_value(l,w_nodal_index)*dpsidx(l,j);
268  }
269  }
270  }
271 
272  /// Get gradient of field: gradient[i] = d[.]/dx_i,
273  // Indices for fields are
274  // 0: W (vertical displacement)
275  // 1: L (laplacian of W)
276  // 2: Ux (In plane displacement x)
277  // 3: Uy (In plane displacement y)
279  Vector<double>& gradient,
280  const unsigned &index) const
281  {
282  //Find out how many nodes there are in the element
283  const unsigned n_node = nnode();
284 
285  //Get the index at which the unknown is stored
286  // Indexes for unknows are
287  // 0: W (vertical displacement)
288  // 1: L (laplacian of W)
289  // 2: Ux (In plane displacement x)
290  // 3: Uy (In plane displacement y)
291  const unsigned w_nodal_index = nodal_index_fvk(index);
292 
293  //Set up memory for the shape and test functions
294  Shape psi(n_node);
295  DShape dpsidx(n_node,2);
296 
297  //Call the derivatives of the shape and test functions
298  dshape_eulerian(s,psi,dpsidx);
299 
300  //Initialise to zero
301  for(unsigned j=0;j<2;j++)
302  {
303  gradient[j] = 0.0;
304  }
305 
306  // Loop over nodes
307  for(unsigned l=0;l<n_node;l++)
308  {
309  //Loop over derivative directions
310  for(unsigned j=0;j<2;j++)
311  {
312  gradient[j] += this->nodal_value(l,w_nodal_index)*dpsidx(l,j);
313  }
314  }
315  }
316 
317  /// Get the in-plane stress (sigma) as a fct of the pre=computed
318  /// displcement derivatives
320  const Vector<double> &interpolated_dwdx,
321  const Vector<double> &interpolated_duxdx,
322  const Vector<double> &interpolated_duydx)
323  {
324  // The strain tensor values
325  double e_xx = interpolated_duxdx[0];
326  double e_yy = interpolated_duydx[1];
327  double e_xy = 0.5 * (interpolated_duxdx[1] + interpolated_duydx[0]);
328  e_xx+= 0.5 * interpolated_dwdx[0] * interpolated_dwdx[0];
329  e_yy+= 0.5 * interpolated_dwdx[1] * interpolated_dwdx[1];
330  e_xy+= 0.5 * interpolated_dwdx[0] * interpolated_dwdx[1];
331 
332  // Get Poisson's ratio
333  const double _nu = nu();
334 
335  // The stress tensor values
336  // sigma_xx
337  sigma(0,0) = (e_xx + _nu*e_yy);
338 
339  // sigma_yy
340  sigma(1,1) = (e_yy + _nu*e_xx);
341 
342  // sigma_xy = sigma_yx
343  sigma(0,1) = sigma(1,0) = e_xy*(1.0 - _nu);
344 
345  }
346 
347  // Get the stress for output
349  DenseMatrix<double> &sigma,
350  DenseMatrix<double> &strain)
351  {
352  //Find out how many nodes there are
353  const unsigned n_node = nnode();
354 
355  //Set up memory for the shape and test functions
356  Shape psi(n_node);
357  DShape dpsidx(n_node,2);
358 
359  //Local shape function
360  dshape_eulerian(s,psi,dpsidx);
361 
362  // Get the derivatives of the shape and test functions
363  //double J = dshape_and_dtest_eulerian_fvk(s, psi, dpsidx, test, dtestdx);
364 
365  //Indices at which the unknowns are stored
366  const unsigned w_nodal_index = nodal_index_fvk(0);
367  const unsigned u_x_nodal_index = nodal_index_fvk(2);
368  const unsigned u_y_nodal_index = nodal_index_fvk(3);
369 
370  //Allocate and initialise to zero storage for the interpolated values
371 
372  // The variables
373  Vector<double> interpolated_dwdx(2,0.0);
374  Vector<double> interpolated_duxdx(2,0.0);
375  Vector<double> interpolated_duydx(2,0.0);
376 
377  //--------------------------------------------
378  //Calculate function values and derivatives:
379  //--------------------------------------------
381  // Loop over nodes
382  for(unsigned l=0;l<n_node;l++)
383  {
384  //Get the nodal values
385  nodal_value[0] = this->nodal_value(l,w_nodal_index);
386  nodal_value[2] = this->nodal_value(l,u_x_nodal_index);
387  nodal_value[3] = this->nodal_value(l,u_y_nodal_index);
388 
389  //Add contributions from current node/shape function
390 
391  // Loop over directions
392  for(unsigned j=0;j<2;j++)
393  {
394  interpolated_dwdx[j] += nodal_value[0]*dpsidx(l,j);
395  interpolated_duxdx[j] += nodal_value[2]*dpsidx(l,j);
396  interpolated_duydx[j] += nodal_value[3]*dpsidx(l,j);
397 
398  } // Loop over directions for (j<2)
399 
400  } // Loop over nodes for (l<n_node)
401 
402 
403  // Get in-plane stress
404  get_sigma(sigma, interpolated_dwdx,
405  interpolated_duxdx, interpolated_duydx);
406 
407 
408  // The strain tensor values
409  // E_xx
410  strain(0,0) = interpolated_duxdx[0];
411  strain(0,0)+= 0.5 * interpolated_dwdx[0] * interpolated_dwdx[0];
412 
413  // E_yy
414  strain(1,1) = interpolated_duydx[1];
415  strain(1,1)+= 0.5 * interpolated_dwdx[1] * interpolated_dwdx[1];
416 
417  // E_xy
418  strain(0,1) = 0.5 * (interpolated_duxdx[1] + interpolated_duydx[0]);
419  strain(0,1)+= 0.5 * interpolated_dwdx[0] * interpolated_dwdx[1];
420 
421  // E_yx
422  strain(1,0)=strain(0,1);
423 
424  }
425 
426 
427  /// hierher dummy
429  Vector<double> &residuals, DenseMatrix<double> &jacobian,
430  DenseMatrix<double> &mass_matrix)
431  {
432 
433  // Get Jacobian from base class (FD-ed)
435  jacobian);
436 
437  // Dummy diagonal (won't result in global unit matrix but
438  // doesn't matter for zero eigenvalue/eigenvector
439  unsigned ndof=mass_matrix.nrow();
440  for (unsigned i=0;i<ndof;i++)
441  {
442  mass_matrix(i,i)+=1.0;
443  }
444 
445  }
446 
447 
448 
449 
450  /// Fill in the residuals with this element's contribution
452  {
453  //Find out how many nodes there are
454  const unsigned n_node = nnode();
455 
456  //Set up memory for the shape and test functions
457  Shape psi(n_node), test(n_node);
458  DShape dpsidx(n_node,2), dtestdx(n_node,2);
459 
460  //Indices at which the unknowns are stored
461  const unsigned w_nodal_index = nodal_index_fvk(0);
462  const unsigned laplacian_w_nodal_index = nodal_index_fvk(1);
463  const unsigned u_x_nodal_index = nodal_index_fvk(2);
464  const unsigned u_y_nodal_index = nodal_index_fvk(3);
465 
466  //Set the value of n_intpt
467  const unsigned n_intpt = integral_pt()->nweight();
468 
469  //Integers to store the local equation numbers
470  int local_eqn=0;
471 
472  //Loop over the integration points
473  for(unsigned ipt=0;ipt<n_intpt;ipt++)
474  {
475  //Get the integral weight
476  double w = integral_pt()->weight(ipt);
477 
478  //Call the derivatives of the shape and test functions
479  double J = dshape_and_dtest_eulerian_at_knot_fvk(ipt,psi,dpsidx,
480  test,dtestdx);
481 
482  //Premultiply the weights and the Jacobian
483  double W = w*J;
484 
485  //Allocate and initialise to zero storage for the interpolated values
487 
488  // The variables
489  double interpolated_laplacian_w = 0;
490  Vector<double> interpolated_dwdx(2,0.0);
491  Vector<double> interpolated_dlaplacian_wdx(2,0.0);
492  Vector<double> interpolated_duxdx(2,0.0);
493  Vector<double> interpolated_duydx(2,0.0);
494 
495  //--------------------------------------------
496  //Calculate function values and derivatives:
497  //--------------------------------------------
499  // Loop over nodes
500  for(unsigned l=0;l<n_node;l++)
501  {
502  //Get the nodal values
503  nodal_value[0] = this->nodal_value(l,w_nodal_index);
504  nodal_value[1] = this->nodal_value(l,laplacian_w_nodal_index);
505  nodal_value[2] = this->nodal_value(l,u_x_nodal_index);
506  nodal_value[3] = this->nodal_value(l,u_y_nodal_index);
507 
508  //Add contributions from current node/shape function
509  interpolated_laplacian_w += nodal_value[1]*psi(l);
510 
511  // Loop over directions
512  for(unsigned j=0;j<2;j++)
513  {
514  interpolated_x[j] += nodal_position(l,j)*psi(l);
515  interpolated_dwdx[j] += nodal_value[0]*dpsidx(l,j);
516  interpolated_dlaplacian_wdx[j] += nodal_value[1]*dpsidx(l,j);
517  interpolated_duxdx[j] += nodal_value[2]*dpsidx(l,j);
518  interpolated_duydx[j] += nodal_value[3]*dpsidx(l,j);
519 
520  } // Loop over directions for (j<2)
521 
522  } // Loop over nodes for (l<n_node)
523 
524  // Get in-plane stress
525  DenseMatrix<double> sigma(2,2,0.0);
526 
527  // Stress not used if we have the linear bending model
529  {
530  // Get the value of in plane stress at the integration
531  // point
532  get_sigma(sigma, interpolated_dwdx,
533  interpolated_duxdx, interpolated_duydx);
534  }
535 
536  //Get pressure function
537  //-------------------
538  double pressure = 0.0;
539  get_pressure_fvk(ipt,interpolated_x,pressure);
540 
541  // Assemble residuals and Jacobian
542  //--------------------------------
543 
544  // Loop over the test functions
545  for(unsigned l=0;l<n_node;l++)
546  {
547  //Get the local equation (First equation)
548  local_eqn = nodal_local_eqn(l,w_nodal_index);
549 
550  // IF it's not a boundary condition
551  if(local_eqn >= 0)
552  {
553  residuals[local_eqn] += pressure*test(l)*W;
554 
555  // Reduced order biharmonic operator
556  for(unsigned k=0;k<2;k++)
557  {
558  residuals[local_eqn] +=
559  interpolated_dlaplacian_wdx[k]*dtestdx(l,k)*W;
560  }
561 
562  // Sigma_alpha_beta part
564  {
565  // Alpha loop
566  for (unsigned alpha=0;alpha<2;alpha++)
567  {
568  // Beta loop
569  for(unsigned beta=0;beta<2;beta++)
570  {
571  residuals[local_eqn]-=
572  12.0*sigma(alpha,beta)*
573  interpolated_dwdx[alpha]*dtestdx(l,beta)*W;
574  }
575  }
576  } // if(!Linear_bending_model)
577  }
578 
579  //Get the local equation (Second equation)
580  local_eqn = nodal_local_eqn(l,laplacian_w_nodal_index);
581 
582  // IF it's not a boundary condition
583  if(local_eqn >= 0)
584  {
585  // The coupled Poisson equations for the biharmonic operator
586  residuals[local_eqn] += interpolated_laplacian_w*test(l)*W;
587 
588  for(unsigned k=0;k<2;k++)
589  {
590  residuals[local_eqn] += interpolated_dwdx[k]*dtestdx(l,k)*W;
591  }
592  }
593 
594  // Get in plane traction
595  Vector<double> traction(2, 0.0);
596  get_traction_fvk(interpolated_x,traction);
597 
598  //Get the local equation (Third equation)
599  local_eqn = nodal_local_eqn(l,u_x_nodal_index);
600 
601  // IF it's not a boundary condition
602  if(local_eqn >= 0)
603  {
604  // tau_x
605  residuals[local_eqn] += traction[0]*test(l)*W;
606 
607  // r_{\alpha = x}
608  for(unsigned beta=0;beta<2;beta++)
609  {
610  residuals[local_eqn] += sigma(0,beta)*dtestdx(l,beta)*W;
611  }
612 
613  }
614 
615  //Get the local equation (Fourth equation)
616  local_eqn = nodal_local_eqn(l,u_y_nodal_index);
617 
618  // IF it's not a boundary condition
619  if(local_eqn >= 0)
620  {
621  // tau_y
622  residuals[local_eqn] += traction[1]*test(l)*W;
623 
624  // r_{\alpha = y}
625  for(unsigned beta=0;beta<2;beta++)
626  {
627  residuals[local_eqn] += sigma(1,beta)*dtestdx(l,beta)*W;
628  }
629 
630  }
631 
632  } // End loop over nodes or test functions (l<n_node)
633 
634  } // End of loop over integration points
635 
636  }
637 
638 
639 
640  /// \short Return FE representation of function value w_fvk(s)
641  /// at local coordinate s (by default - if index > 0, returns
642  /// FE representation of valued stored at index^th nodal index
643  inline double interpolated_w_fvk(const Vector<double> &s,
644  unsigned index=0) const
645  {
646  //Find number of nodes
647  const unsigned n_node = nnode();
648 
649  //Get the index at which the poisson unknown is stored
650  const unsigned w_nodal_index = nodal_index_fvk(index);
651 
652  //Local shape function
653  Shape psi(n_node);
654 
655  //Find values of shape function
656  shape(s,psi);
657 
658  //Initialise value of u
659  double interpolated_w = 0.0;
660 
661  //Loop over the local nodes and sum
662  for(unsigned l=0;l<n_node;l++)
663  {
664  interpolated_w += this->nodal_value(l,w_nodal_index)*psi[l];
665  }
666 
667  return(interpolated_w);
668  }
669 
670  /// \short Self-test: Return 0 for OK
671  unsigned self_test();
672 
673  /// \short Sets a flag to signify that we are solving the linear, pure bending
674  /// equations, and pin all the nodal values that will not be used in this case
676  {
677  // Set the boolean flag
678  Linear_bending_model = true;
679 
680  // Get the index of the first FvK nodal value
681  unsigned first_fvk_nodal_index = nodal_index_fvk();
682 
683  // Get the total number of FvK nodal values (assuming they are stored
684  // contiguously) at node 0 (it's the same at all nodes anyway)
685  unsigned total_fvk_nodal_indices = 4;
686 
687  // Get the number of nodes in this element
688  unsigned n_node = nnode();
689 
690  // Loop over the appropriate nodal indices
691  for(unsigned index=first_fvk_nodal_index+2; // because [2] is u_x
692  // and [3] is u_y
693  index<first_fvk_nodal_index+total_fvk_nodal_indices;
694  index++)
695  {
696  // Loop over the nodes in the element
697  for(unsigned inod=0;inod<n_node;inod++)
698  {
699  // Pin the nodal value at the current index
700  node_pt(inod)->pin(index);
701  }
702  }
703  }
704 
705 
706 protected:
707 
708  /// \short Shape/test functions and derivs w.r.t. to global coords at
709  /// local coord. s; return Jacobian of mapping
710  virtual double dshape_and_dtest_eulerian_fvk(const Vector<double> &s,
711  Shape &psi,
712  DShape &dpsidx, Shape &test,
713  DShape &dtestdx) const=0;
714 
715 
716  /// \short Shape/test functions and derivs w.r.t. to global coords at
717  /// integration point ipt; return Jacobian of mapping
718  virtual double dshape_and_dtest_eulerian_at_knot_fvk(const unsigned &ipt,
719  Shape &psi,
720  DShape &dpsidx,
721  Shape &test,
722  DShape &dtestdx)
723  const=0;
724 
725 
726  /// Pointer to global Poisson's ratio
727  double *Nu_pt;
728 
729  /// Pointer to pressure function:
731 
732  /// Pointer to traction function:
734 
735 private:
736 
737  /// Default value for Poisson's ratio
738  static double Default_Nu_Value;
739 
740  /// \short Flag which stores whether we are using a linear, pure
741  /// bending model instead of the full non-linear Foeppl-von Karman
743 
744 };
745 
746 
747 //==========================================================
748 /// Foeppl von Karman upgraded to become projectable
749 //==========================================================
750  template<class FVK_ELEMENT>
752  public virtual ProjectableElement<FVK_ELEMENT>
753  {
754 
755  public:
756 
757  /// \short Specify the values associated with field fld. The
758  /// information is returned in a vector of pairs which comprise the
759  /// Data object and the value within it, that correspond to field
760  /// fld.
762  {
763 
764 #ifdef PARANOID
765  if (fld > 3)
766  {
767  std::stringstream error_stream;
768  error_stream
769  << "Foeppl von Karman elements only store four fields so fld must be"
770  << "0 to 3 rather than " << fld << std::endl;
771  throw OomphLibError(
772  error_stream.str(),
773  OOMPH_CURRENT_FUNCTION,
774  OOMPH_EXCEPTION_LOCATION);
775  }
776 #endif
777 
778  // Create the vector
779  unsigned nnod=this->nnode();
780  Vector<std::pair<Data*,unsigned> > data_values(nnod);
781 
782  // Loop over all nodes
783  for (unsigned j=0;j<nnod;j++)
784  {
785  // Add the data value associated field: The node itself
786  data_values[j]=std::make_pair(this->node_pt(j),fld);
787  }
788 
789  // Return the vector
790  return data_values;
791  }
792 
793  /// \short Number of fields to be projected: Just two
795  {
796  return 4;
797  }
798 
799  /// \short Number of history values to be stored for fld-th field.
800  /// (Note: count includes current value!)
801  unsigned nhistory_values_for_projection(const unsigned &fld)
802  {
803 #ifdef PARANOID
804  if (fld > 3)
805  {
806  std::stringstream error_stream;
807  error_stream
808  << "Foeppl von Karman elements only store four fields so fld must be"
809  << "0 to 3 rather than " << fld << std::endl;
810  throw OomphLibError(
811  error_stream.str(),
812  OOMPH_CURRENT_FUNCTION,
813  OOMPH_EXCEPTION_LOCATION);
814  }
815 #endif
816  return this->node_pt(0)->ntstorage();
817  }
818 
819 
820  ///\short Number of positional history values
821  /// (Note: count includes current value!)
823  {
824  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
825  }
826 
827  /// \short Return Jacobian of mapping and shape functions of field fld
828  /// at local coordinate s
829  double jacobian_and_shape_of_field(const unsigned &fld,
830  const Vector<double> &s,
831  Shape &psi)
832  {
833 #ifdef PARANOID
834  if (fld > 3)
835  {
836  std::stringstream error_stream;
837  error_stream
838  << "Foeppl von Karman elements only store four fields so fld must be"
839  << "0 to 3 rather than " << fld << std::endl;
840  throw OomphLibError(
841  error_stream.str(),
842  OOMPH_CURRENT_FUNCTION,
843  OOMPH_EXCEPTION_LOCATION);
844  }
845 #endif
846  unsigned n_dim=this->dim();
847  unsigned n_node=this->nnode();
848  Shape test(n_node);
849  DShape dpsidx(n_node,n_dim), dtestdx(n_node,n_dim);
850  double J=this->dshape_and_dtest_eulerian_fvk(s,psi,dpsidx,
851  test,dtestdx);
852  return J;
853  }
854 
855 
856 
857  /// \short Return interpolated field fld at local coordinate s, at
858  /// time level t (t=0: present; t>0: history values)
859  double get_field(const unsigned &t,
860  const unsigned &fld,
861  const Vector<double>& s)
862  {
863 
864 #ifdef PARANOID
865  if (fld > 3)
866  {
867  std::stringstream error_stream;
868  error_stream
869  << "Foeppl von Karman elements only store four fields so fld must be"
870  << "0 to 3 rather than " << fld << std::endl;
871  throw OomphLibError(
872  error_stream.str(),
873  OOMPH_CURRENT_FUNCTION,
874  OOMPH_EXCEPTION_LOCATION);
875  }
876 #endif
877  //Find the index at which the variable is stored
878  unsigned w_nodal_index = this->nodal_index_fvk(fld);
879 
880  //Local shape function
881  unsigned n_node=this->nnode();
882  Shape psi(n_node);
883 
884  //Find values of shape function
885  this->shape(s,psi);
886 
887  //Initialise value of u
888  double interpolated_w = 0.0;
889 
890  //Sum over the local nodes
891  for(unsigned l=0;l<n_node;l++)
892  {
893  interpolated_w += this->nodal_value(t,l,w_nodal_index)*psi[l];
894  }
895  return interpolated_w;
896  }
897 
898 
899 
900 
901  ///Return number of values in field fld: One per node
902  unsigned nvalue_of_field(const unsigned &fld)
903  {
904 #ifdef PARANOID
905  if (fld > 3)
906  {
907  std::stringstream error_stream;
908  error_stream
909  << "Foeppl von Karman elements only store four fields so fld must be"
910  << "0 to 3 rather than " << fld << std::endl;
911  throw OomphLibError(
912  error_stream.str(),
913  OOMPH_CURRENT_FUNCTION,
914  OOMPH_EXCEPTION_LOCATION);
915  }
916 #endif
917  return this->nnode();
918  }
919 
920 
921  ///Return local equation number of field fld of node j.
922  int local_equation(const unsigned &fld,
923  const unsigned &j)
924  {
925 #ifdef PARANOID
926  if (fld > 3)
927  {
928  std::stringstream error_stream;
929  error_stream
930  << "Foeppl von Karman elements only store four fields so fld must be"
931  << "0 to 3 rather than " << fld << std::endl;
932  throw OomphLibError(
933  error_stream.str(),
934  OOMPH_CURRENT_FUNCTION,
935  OOMPH_EXCEPTION_LOCATION);
936  }
937 #endif
938  const unsigned w_nodal_index = this->nodal_index_fvk(fld);
939  return this->nodal_local_eqn(j,w_nodal_index);
940  }
941 
942  };
943 
944 //=======================================================================
945 /// Face geometry for element is the same as that for the underlying
946 /// wrapped element
947 //=======================================================================
948  template<class ELEMENT>
950  : public virtual FaceGeometry<ELEMENT>
951  {
952  public:
953  FaceGeometry() : FaceGeometry<ELEMENT>() {}
954  };
955 
956 
957 //=======================================================================
958 /// Face geometry of the Face Geometry for element is the same as
959 /// that for the underlying wrapped element
960 //=======================================================================
961  template<class ELEMENT>
963  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
964  {
965  public:
967  };
968 
969 }
970 
971 #endif
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
FoepplvonKarmanTractionFctPt & traction_fct_pt()
Access function: Pointer to in-plane traction function.
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 ...
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
virtual void get_traction_fvk(Vector< double > &x, Vector< double > &traction) const
Get in-plane traction term at (Eulerian) position x.
FoepplvonKarmanPressureFctPt & pressure_fct_pt()
Access function: Pointer to pressure function.
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 ...
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
FoepplvonKarmanDisplacementEquations(const FoepplvonKarmanDisplacementEquations &dummy)
Broken copy constructor.
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
void use_linear_bending_model()
Sets a flag to signify that we are solving the linear, pure bending equations, and pin all the nodal ...
cstr elem_len * i
Definition: cfortran.h:607
FoepplvonKarmanDisplacementEquations()
Constructor (must initialise the Pressure_fct_pt and the Traction_fct_pt. Also set physical parameter...
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 ...
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
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 get_sigma(DenseMatrix< double > &sigma, const Vector< double > &interpolated_dwdx, const Vector< double > &interpolated_duxdx, const Vector< double > &interpolated_duydx)
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
unsigned nfields_for_projection()
Number of fields to be projected: Just two.
FoepplvonKarmanTractionFctPt traction_fct_pt() const
Access function: Pointer to in-plane traction function. Const version.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
void(* FoepplvonKarmanPressureFctPt)(const Vector< double > &x, double &f)
Function pointer to pressure function fct(x,f(x)) – x is a Vector!
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:504
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...
void(* FoepplvonKarmanTractionFctPt)(const Vector< double > &x, Vector< double > &f)
Function pointer to traction function fct(x,f(x)) – x is a Vector!
void operator=(const FoepplvonKarmanDisplacementEquations &)
Broken assignment operator.
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 ...
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of field fld of node j.
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:3841
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
hierher dummy
bool Linear_bending_model
Flag which stores whether we are using a linear, pure bending model instead of the full non-linear Fo...
static char t char * s
Definition: cfortran.h:572
void output(FILE *file_pt)
C_style output with default number of plot points.
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...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
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(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
double 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.
FoepplvonKarmanPressureFctPt Pressure_fct_pt
Pointer to pressure function:
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2215
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
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
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals with this element's contribution.
void get_stress_and_strain_for_output(const Vector< double > &s, DenseMatrix< double > &sigma, DenseMatrix< double > &strain)
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.
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!) ...
FoepplvonKarmanTractionFctPt Traction_fct_pt
Pointer to traction function:
void output(std::ostream &outfile)
Output with default number of plot points.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the elemental contribution to the jacobian matrix. and the residuals vector. Note that this funct...
Definition: elements.h:1694
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 get_gradient_of_field(const Vector< double > &s, Vector< double > &gradient, const unsigned &index) const
Get gradient of field: gradient[i] = d[.]/dx_i,.
void get_gradient_of_deflection(const Vector< double > &s, Vector< double > &gradient) const
Get gradient of deflection: gradient[i] = dw/dx_i.
static double Default_Nu_Value
Default value for Poisson's ratio.