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