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-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for 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
43namespace 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!
60 double& f);
61
62 /// Function pointer to traction function fct(x,f(x)) --
63 /// x is a Vector!
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?
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 //--------------------------------------------
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 //--------------------------------------------
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);
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
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)
FoepplvonKarmanTractionFctPt & traction_fct_pt()
Access function: Pointer to in-plane traction 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!
FoepplvonKarmanPressureFctPt & pressure_fct_pt()
Access function: Pointer to pressure function.
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.
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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
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
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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1765
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...
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 ...
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...