time_harmonic_fourier_decomposed_linear_elasticity_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 
27 // Include guards to prevent multiple inclusion of the header
28 #ifndef OOMPH_FOURIER_DECOMPOSED_TIME_HARMONIC_LINEAR_ELASTICITY_ELEMENTS_HEADER
29 #define OOMPH_FOURIER_DECOMPOSED_TIME_HARMONIC_LINEAR_ELASTICITY_ELEMENTS_HEADER
30 
31 // Config header generated by autoconfig
32 #ifdef HAVE_CONFIG_H
33 #include <oomph-lib-config.h>
34 #endif
35 
36 
37 #ifdef OOMPH_HAS_MPI
38 #include "mpi.h"
39 #endif
40 
41 #include <complex>
42 
43 
44 // OOMPH-LIB headers
45 #include "../generic/Qelements.h"
46 #include "../generic/Telements.h"
47 #include "../generic/projection.h"
48 #include "../generic/error_estimator.h"
49 
50 
51 namespace oomph
52 {
53  //=======================================================================
54  /// A base class for elements that solve the Fourier decomposed (in
55  /// cylindrical polars) equations of time-harmonic linear elasticity.
56  //=======================================================================
58  : public virtual FiniteElement
59  {
60  public:
61  /// Return the index at which the i-th (i=0: r, i=1: z; i=2: theta)
62  /// real or imag unknown displacement component is stored at the nodes.
63  /// The default assignment here (u_r_real, u_z_real, ..., u_theta_imag)
64  /// is appropriate for single-physics problems.
65  virtual inline std::complex<unsigned> u_index_time_harmonic_fourier_decomposed_linear_elasticity(
66  const unsigned i) const
67  {
68  return std::complex<unsigned>(i, i + 3);
69  }
70 
71  /// Compute vector of FE interpolated displacement u at local coordinate s
73  const Vector<double>& s, Vector<std::complex<double>>& disp) const
74  {
75  // Find number of nodes
76  unsigned n_node = nnode();
77 
78  // Local shape function
79  Shape psi(n_node);
80 
81  // Find values of shape function
82  shape(s, psi);
83 
84  for (unsigned i = 0; i < 3; i++)
85  {
86  // Index at which the nodal value is stored
87  std::complex<unsigned> u_nodal_index =
89 
90  // Initialise value of u
91  disp[i] = std::complex<double>(0.0, 0.0);
92 
93  // Loop over the local nodes and sum
94  for (unsigned l = 0; l < n_node; l++)
95  {
96  const std::complex<double> u_value(
97  nodal_value(l, u_nodal_index.real()),
98  nodal_value(l, u_nodal_index.imag()));
99 
100  disp[i] += u_value * psi[l];
101  }
102  }
103  }
104 
105  /// Return FE interpolated displacement u[i] (i=0: r, i=1: z; i=2:
106  /// theta) at local coordinate s
108  const Vector<double>& s, const unsigned& i) const
109  {
110  // Find number of nodes
111  unsigned n_node = nnode();
112 
113  // Local shape function
114  Shape psi(n_node);
115 
116  // Find values of shape function
117  shape(s, psi);
118 
119  // Get nodal index at which i-th velocity is stored
120  std::complex<unsigned> u_nodal_index =
122 
123  // Initialise value of u
124  std::complex<double> interpolated_u(0.0, 0.0);
125 
126  // Loop over the local nodes and sum
127  for (unsigned l = 0; l < n_node; l++)
128  {
129  const std::complex<double> u_value(
130  nodal_value(l, u_nodal_index.real()),
131  nodal_value(l, u_nodal_index.imag()));
132 
133  interpolated_u += u_value * psi[l];
134  }
135 
136  return (interpolated_u);
137  }
138 
139 
140  /// Function pointer to function that specifies the body force
141  /// as a function of the Cartesian coordinates and time FCT(x,b) --
142  /// x and b are Vectors!
143  typedef void (*BodyForceFctPt)(const Vector<double>& x,
145 
146  /// Constructor: Set null pointers for constitutive law.
147  /// Set physical parameter values to
148  /// default values, and set body force to zero.
152  Nu_pt(0),
155  {
156  }
157 
158  /// Access function for square of non-dim frequency
159  const std::complex<double>& omega_sq() const
160  {
161  return *Omega_sq_pt;
162  }
163 
164  /// Access function for square of non-dim frequency
165  std::complex<double>*& omega_sq_pt()
166  {
167  return Omega_sq_pt;
168  }
169 
170  /// Return the pointer to Young's modulus
171  std::complex<double>*& youngs_modulus_pt()
172  {
173  return Youngs_modulus_pt;
174  }
175 
176  /// Access function to Young's modulus
177  inline std::complex<double> youngs_modulus() const
178  {
179  return (*Youngs_modulus_pt);
180  }
181 
182  /// Access function for Poisson's ratio
183  std::complex<double>& nu() const
184  {
185 #ifdef PARANOID
186  if (Nu_pt == 0)
187  {
188  std::ostringstream error_message;
189  error_message << "No pointer to Poisson's ratio set. Please set one!\n";
190  throw OomphLibError(error_message.str(),
191  OOMPH_CURRENT_FUNCTION,
192  OOMPH_EXCEPTION_LOCATION);
193  }
194 #endif
195  return *Nu_pt;
196  }
197 
198  /// Access function for pointer to Poisson's ratio
199  std::complex<double>*& nu_pt()
200  {
201  return Nu_pt;
202  }
203 
204  /// Access function for Fourier wavenumber
205  int& fourier_wavenumber() const
206  {
207 #ifdef PARANOID
208  if (Fourier_wavenumber_pt == 0)
209  {
210  std::ostringstream error_message;
211  error_message
212  << "No pointer to Fourier wavenumber set. Please set one!\n";
213  throw OomphLibError(error_message.str(),
214  OOMPH_CURRENT_FUNCTION,
215  OOMPH_EXCEPTION_LOCATION);
216  }
217 #endif
218  return *Fourier_wavenumber_pt;
219  }
220 
221  /// Access function for pointer to Fourier wavenumber
223  {
224  return Fourier_wavenumber_pt;
225  }
226 
227  /// Access function: Pointer to body force function
229  {
230  return Body_force_fct_pt;
231  }
232 
233  /// Access function: Pointer to body force function (const version)
235  {
236  return Body_force_fct_pt;
237  }
238 
239  /// Evaluate body force at Eulerian coordinate x at present time
240  /// (returns zero vector if no body force function pointer has been set)
241  inline void body_force(const Vector<double>& x,
242  Vector<std::complex<double>>& b) const
243  {
244  // If no function has been set, return zero vector
245  if (Body_force_fct_pt == 0)
246  {
247  // Get spatial dimension of element
248  unsigned n = dim();
249  for (unsigned i = 0; i < n; i++)
250  {
251  b[i] = std::complex<double>(0.0, 0.0);
252  }
253  }
254  else
255  {
256  (*Body_force_fct_pt)(x, b);
257  }
258  }
259 
260  /// The number of "DOF types" that degrees of freedom in this element
261  /// are sub-divided into: for now lump them all into one DOF type.
262  /// Can be adjusted later
263  unsigned ndof_types() const
264  {
265  return 1;
266  }
267 
268  /// Create a list of pairs for all unknowns in this element,
269  /// so that the first entry in each pair contains the global equation
270  /// number of the unknown, while the second one contains the number
271  /// of the "DOF type" that this unknown is associated with.
272  /// (Function can obviously only be called if the equation numbering
273  /// scheme has been set up.)
275  std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
276  {
277  // temporary pair (used to store dof lookup prior to being added
278  // to list)
279  std::pair<unsigned long, unsigned> dof_lookup;
280 
281  // number of nodes
282  const unsigned n_node = this->nnode();
283 
284  // Integer storage for local unknown
285  int local_unknown = 0;
286 
287  // Loop over the nodes
288  for (unsigned n = 0; n < n_node; n++)
289  {
290  // Loop over dimension (real and imag displacement components)
291  for (unsigned i = 0; i < 6; i++)
292  {
293  // If the variable is free
294  local_unknown = nodal_local_eqn(n, i);
295 
296  // ignore pinned values
297  if (local_unknown >= 0)
298  {
299  // store dof lookup in temporary pair: First entry in pair
300  // is global equation number; second entry is DOF type
301  dof_lookup.first = this->eqn_number(local_unknown);
302  dof_lookup.second = 0;
303 
304  // add to list
305  dof_lookup_list.push_front(dof_lookup);
306  }
307  }
308  }
309  }
310 
311 
312  protected:
313  /// Square of nondim frequency
314  std::complex<double>* Omega_sq_pt;
315 
316  /// Pointer to the Young's modulus
317  std::complex<double>* Youngs_modulus_pt;
318 
319  /// Pointer to Poisson's ratio
320  std::complex<double>* Nu_pt;
321 
322  /// Pointer to Fourier wavenumber
324 
325  /// Pointer to body force function
327 
328  /// Static default value for squared frequency
329  static std::complex<double> Default_omega_sq_value;
330 
331  /// Static default value for Young's modulus (1.0 -- for natural
332  /// scaling, i.e. all stresses have been non-dimensionalised by
333  /// the same reference Young's modulus. Setting the "non-dimensional"
334  /// Young's modulus (obtained by de-referencing Youngs_modulus_pt)
335  /// to a number larger than one means that the material is stiffer
336  /// than assumed in the non-dimensionalisation.
337  static std::complex<double> Default_youngs_modulus_value;
338  };
339 
340 
341  /// ////////////////////////////////////////////////////////////////////
342  /// ////////////////////////////////////////////////////////////////////
343  /// ////////////////////////////////////////////////////////////////////
344 
345 
346  //=======================================================================
347  /// A class for elements that solve the Fourier decomposed (in cylindrical
348  /// polars) equations of time-harmonic linear elasticity
349  //=======================================================================
352  {
353  public:
354  /// Constructor
356 
357  /// Number of values required at node n.
358  unsigned required_nvalue(const unsigned& n) const
359  {
360  return 6;
361  }
362 
363  /// Return the residuals for the equations (the discretised
364  /// principle of virtual displacements)
366  {
368  residuals, GeneralisedElement::Dummy_matrix, 0);
369  }
370 
371 
372  /// The jacobian is calculated by finite differences by default,
373  /// We need only to take finite differences w.r.t. positional variables
374  /// For this element
376  DenseMatrix<double>& jacobian)
377  {
378  // Add the contribution to the residuals
379  this
380  ->fill_in_generic_contribution_to_residuals_fourier_decomp_time_harmonic_linear_elasticity(
381  residuals, jacobian, 1);
382  }
383 
384 
385  /// Get strain tensor
386  void get_strain(const Vector<double>& s,
387  DenseMatrix<std::complex<double>>& strain);
388 
389  /// Compute norm of solution: square of the L2 norm
390  void compute_norm(double& norm);
391 
392  /// Output exact solution: r,z, u_r_real, u_z_real, ..., u_theta_imag
393  void output_fct(std::ostream& outfile,
394  const unsigned& nplot,
396 
397  /// Output: r,z, u_r_real, u_z_real, ..., u_theta_imag
398  void output(std::ostream& outfile)
399  {
400  unsigned n_plot = 5;
401  output(outfile, n_plot);
402  }
403 
404  /// Output: r,z, u_r_real, u_z_real, ..., u_theta_imag
405  void output(std::ostream& outfile, const unsigned& n_plot);
406 
407  /// C-style output: r,z, u_r_real, u_z_real, ..., u_theta_imag
408  void output(FILE* file_pt)
409  {
410  unsigned n_plot = 5;
411  output(file_pt, n_plot);
412  }
413 
414  /// Output: r,z, u_r_real, u_z_real, ..., u_theta_imag
415  void output(FILE* file_pt, const unsigned& n_plot);
416 
417  /// Validate against exact solution.
418  /// Solution is provided via function pointer.
419  /// Plot at a given number of plot points and compute L2 error
420  /// and L2 norm of displacement solution over element
421  void compute_error(std::ostream& outfile,
423  double& error,
424  double& norm);
425 
426 
427  private:
428  /// Private helper function to compute residuals and (if requested
429  /// via flag) also the Jacobian matrix.
431  Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
432  };
433 
434 
435  /// /////////////////////////////////////////////////////////////////////
436  /// /////////////////////////////////////////////////////////////////////
437  /// /////////////////////////////////////////////////////////////////////
438 
439 
440  //===========================================================================
441  /// An Element that solves the equations of Fourier decomposed (in cylindrical
442  /// polars) time-harmonic linear elasticity, using QElements for the geometry.
443  //============================================================================
444  template<unsigned NNODE_1D>
446  : public virtual QElement<2, NNODE_1D>,
448  {
449  public:
450  /// Constructor
452  : QElement<2, NNODE_1D>(),
454  {
455  }
456 
457  /// Output function
458  void output(std::ostream& outfile)
459  {
461  }
462 
463  /// Output function
464  void output(std::ostream& outfile, const unsigned& n_plot)
465  {
467  n_plot);
468  }
469 
470 
471  /// C-style output function
472  void output(FILE* file_pt)
473  {
475  }
476 
477  /// C-style output function
478  void output(FILE* file_pt, const unsigned& n_plot)
479  {
481  n_plot);
482  }
483  };
484 
485 
486  //============================================================================
487  /// FaceGeometry of a linear
488  /// QTimeHarmonicFourierDecomposedLinearElasticityElement element
489  //============================================================================
490  template<unsigned NNODE_1D>
493  : public virtual QElement<1, NNODE_1D>
494  {
495  public:
496  /// Constructor must call the constructor of the underlying element
497  FaceGeometry() : QElement<1, NNODE_1D>() {}
498  };
499 
500 
501  /// /////////////////////////////////////////////////////////////////////
502  /// /////////////////////////////////////////////////////////////////////
503  /// /////////////////////////////////////////////////////////////////////
504 
505 
506  //===========================================================================
507  /// An Element that solves the equations of Fourier decomposed (in cylindrical
508  /// polars) time-harmonic linear elasticity, using TElements for the geometry.
509  //============================================================================
510  template<unsigned NNODE_1D>
512  : public virtual TElement<2, NNODE_1D>,
514  public virtual ElementWithZ2ErrorEstimator
515  {
516  public:
517  /// Constructor
519  : TElement<2, NNODE_1D>(),
521  {
522  }
523 
524  /// Output function
525  void output(std::ostream& outfile)
526  {
528  }
529 
530  /// Output function
531  void output(std::ostream& outfile, const unsigned& n_plot)
532  {
534  n_plot);
535  }
536 
537  /// C-style output function
538  void output(FILE* file_pt)
539  {
541  }
542 
543  /// C-style output function
544  void output(FILE* file_pt, const unsigned& n_plot)
545  {
547  n_plot);
548  }
549 
550 
551  /// Number of vertex nodes in the element
552  unsigned nvertex_node() const
553  {
555  }
556 
557  /// Pointer to the j-th vertex node in the element
558  Node* vertex_node_pt(const unsigned& j) const
559  {
561  }
562 
563  /// Order of recovery shape functions for Z2 error estimation:
564  /// Same order as shape functions.
565  unsigned nrecovery_order()
566  {
567  return NNODE_1D - 1;
568  }
569 
570  /// Number of 'flux' terms for Z2 error estimation
571  unsigned num_Z2_flux_terms()
572  {
573  // 3 Diagonal strain rates and 3 off diagonal terms for real and imag part
574  return 12;
575  }
576 
577  /// Get 'flux' for Z2 error recovery: Upper triangular entries
578  /// in strain tensor.
580  {
581 #ifdef PARANOID
582  unsigned num_entries = 12;
583  if (flux.size() != num_entries)
584  {
585  std::ostringstream error_message;
586  error_message << "The flux vector has the wrong number of entries, "
587  << flux.size() << ", whereas it should be " << num_entries
588  << std::endl;
589  throw OomphLibError(error_message.str(),
590  OOMPH_CURRENT_FUNCTION,
591  OOMPH_EXCEPTION_LOCATION);
592  }
593 #endif
594 
595  // Get strain matrix
597  this->get_strain(s, strain);
598 
599  // Pack into flux Vector
600  unsigned icount = 0;
601 
602  // Start with diagonal terms
603  for (unsigned i = 0; i < 3; i++)
604  {
605  flux[icount] = strain(i, i).real();
606  icount++;
607  flux[icount] = strain(i, i).imag();
608  icount++;
609  }
610 
611  // Off diagonals row by row
612  for (unsigned i = 0; i < 3; i++)
613  {
614  for (unsigned j = i + 1; j < 3; j++)
615  {
616  flux[icount] = strain(i, j).real();
617  icount++;
618  flux[icount] = strain(i, j).imag();
619  icount++;
620  }
621  }
622  }
623  };
624 
625 
626  //============================================================================
627  /// FaceGeometry of a linear
628  /// TTimeHarmonicFourierDecomposedLinearElasticityElement element
629  //============================================================================
630  template<unsigned NNODE_1D>
633  : public virtual TElement<1, NNODE_1D>
634  {
635  public:
636  /// Constructor must call the constructor of the underlying element
637  FaceGeometry() : TElement<1, NNODE_1D>() {}
638  };
639 
640 
641  /// /////////////////////////////////////////////////////////////////
642  /// /////////////////////////////////////////////////////////////////
643  /// /////////////////////////////////////////////////////////////////
644 
645 
646  //==========================================================
647  /// Fourier-decomposed time-harmonic linear elasticity
648  /// upgraded to become projectable
649  //==========================================================
650  template<class TIME_HARMONIC_LINEAR_ELAST_ELEMENT>
652  : public virtual ProjectableElement<TIME_HARMONIC_LINEAR_ELAST_ELEMENT>
653  {
654  public:
655  /// Constructor [this was only required explicitly
656  /// from gcc 4.5.2 onwards...]
658 
659  /// Specify the values associated with field fld.
660  /// The information is returned in a vector of pairs which comprise
661  /// the Data object and the value within it, that correspond to field fld.
662  /// In the underlying time-harmonic linear elasticity elemements the
663  /// real and complex parts of the displacements are stored
664  /// at the nodal values
666  {
667  // Create the vector
669 
670  // Loop over all nodes and extract the fld-th nodal value
671  unsigned nnod = this->nnode();
672  for (unsigned j = 0; j < nnod; j++)
673  {
674  // Add the data value associated with the velocity components
675  data_values.push_back(std::make_pair(this->node_pt(j), fld));
676  }
677 
678  // Return the vector
679  return data_values;
680  }
681 
682  /// Number of fields to be projected: 3*dim, corresponding to
683  /// real and imag parts of the displacement components
685  {
686  return 3 * this->dim();
687  }
688 
689  /// Number of history values to be stored for fld-th field.
690  /// (includes present value!)
691  unsigned nhistory_values_for_projection(const unsigned& fld)
692  {
693 #ifdef PARANOID
694  if (fld > 5)
695  {
696  std::stringstream error_stream;
697  error_stream << "Elements only store six fields so fld can't be"
698  << " " << fld << std::endl;
699  throw OomphLibError(
700  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
701  }
702 #endif
703  return this->node_pt(0)->ntstorage();
704  }
705 
706  /// Number of positional history values: Read out from
707  /// positional timestepper
708  /// (Note: count includes current value!)
710  {
711  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
712  }
713 
714  /// Return Jacobian of mapping and shape functions of field fld
715  /// at local coordinate s
716  double jacobian_and_shape_of_field(const unsigned& fld,
717  const Vector<double>& s,
718  Shape& psi)
719  {
720  unsigned n_dim = this->dim();
721  unsigned n_node = this->nnode();
722  DShape dpsidx(n_node, n_dim);
723 
724  // Call the derivatives of the shape functions and return
725  // the Jacobian
726  return this->dshape_eulerian(s, psi, dpsidx);
727  }
728 
729 
730  /// Return interpolated field fld at local coordinate s, at time
731  /// level t (t=0: present; t>0: history values)
732  double get_field(const unsigned& t,
733  const unsigned& fld,
734  const Vector<double>& s)
735  {
736  unsigned n_node = this->nnode();
737 
738 #ifdef PARANOID
739  unsigned n_dim = this->node_pt(0)->ndim();
740 #endif
741 
742  // Local shape function
743  Shape psi(n_node);
744 
745  // Find values of shape function
746  this->shape(s, psi);
747 
748  // Initialise value of u
749  double interpolated_u = 0.0;
750 
751  // Sum over the local nodes at that time
752  for (unsigned l = 0; l < n_node; l++)
753  {
754 #ifdef PARANOID
755  unsigned nvalue = this->node_pt(l)->nvalue();
756  if (nvalue != 3 * n_dim)
757  {
758  std::stringstream error_stream;
759  error_stream
760  << "Current implementation only works for non-resized nodes\n"
761  << "but nvalue= " << nvalue << "!= 3 dim = " << 3 * n_dim
762  << std::endl;
763  throw OomphLibError(error_stream.str(),
764  OOMPH_CURRENT_FUNCTION,
765  OOMPH_EXCEPTION_LOCATION);
766  }
767 #endif
768  interpolated_u += this->nodal_value(t, l, fld) * psi[l];
769  }
770  return interpolated_u;
771  }
772 
773 
774  /// Return number of values in field fld
775  unsigned nvalue_of_field(const unsigned& fld)
776  {
777  return this->nnode();
778  }
779 
780 
781  /// Return local equation number of value j in field fld.
782  int local_equation(const unsigned& fld, const unsigned& j)
783  {
784 #ifdef PARANOID
785  unsigned n_dim = this->node_pt(0)->ndim();
786  unsigned nvalue = this->node_pt(j)->nvalue();
787  if (nvalue != 3 * n_dim)
788  {
789  std::stringstream error_stream;
790  error_stream
791  << "Current implementation only works for non-resized nodes\n"
792  << "but nvalue= " << nvalue << "!= 3 dim = " << 3 * n_dim
793  << std::endl;
794  throw OomphLibError(
795  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
796  }
797 #endif
798  return this->nodal_local_eqn(j, fld);
799  }
800  };
801 
802 
803  //=======================================================================
804  /// Face geometry for element is the same as that for the underlying
805  /// wrapped element
806  //=======================================================================
807  template<class ELEMENT>
810  : public virtual FaceGeometry<ELEMENT>
811  {
812  public:
813  FaceGeometry() : FaceGeometry<ELEMENT>() {}
814  };
815 
816 
817  //=======================================================================
818  /// Face geometry of the Face Geometry for element is the same as
819  /// that for the underlying wrapped element
820  //=======================================================================
821  template<class ELEMENT>
824  : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
825  {
826  public:
828  };
829 
830 
831 } // namespace oomph
832 
833 
834 #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
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
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
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
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...
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 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
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Definition: elements.h:704
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:227
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
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....
Wrapper class for projectable elements. Adds "projectability" to the underlying ELEMENT.
Definition: projection.h:183
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (includes present value!)
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.
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 nhistory_values_for_coordinate_projection()
Number of positional history values: Read out from positional timestepper (Note: count includes curre...
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
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 nfields_for_projection()
Number of fields to be projected: 3*dim, corresponding to real and imag parts of the displacement com...
ProjectableTimeHarmonicFourierDecomposedLinearElasticityElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
General TElement class.
Definition: Telements.h:1208
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain tensor.
A base class for elements that solve the Fourier decomposed (in cylindrical polars) equations of time...
std::complex< double > interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(const Vector< double > &s, const unsigned &i) const
Return FE interpolated displacement u[i] (i=0: r, i=1: z; i=2: theta) at local coordinate s.
static std::complex< double > Default_youngs_modulus_value
Static default value for Young's modulus (1.0 – for natural scaling, i.e. all stresses have been non-...
static std::complex< double > Default_omega_sq_value
Static default value for squared frequency.
virtual std::complex< unsigned > u_index_time_harmonic_fourier_decomposed_linear_elasticity(const unsigned i) const
Return the index at which the i-th (i=0: r, i=1: z; i=2: theta) real or imag unknown displacement com...
void interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(const Vector< double > &s, Vector< std::complex< double >> &disp) const
Compute vector of FE interpolated displacement u at local coordinate s.
void body_force(const Vector< double > &x, Vector< std::complex< double >> &b) const
Evaluate body force at Eulerian coordinate x at present time (returns zero vector if no body force fu...
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: for now lump ...
std::complex< double > *& nu_pt()
Access function for pointer to Poisson's ratio.
TimeHarmonicFourierDecomposedLinearElasticityEquationsBase()
Constructor: Set null pointers for constitutive law. Set physical parameter values to default values,...
std::complex< double > *& omega_sq_pt()
Access function for square of non-dim frequency.
void(* BodyForceFctPt)(const Vector< double > &x, Vector< std::complex< double >> &b)
Function pointer to function that specifies the body force as a function of the Cartesian coordinates...
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
BodyForceFctPt body_force_fct_pt() const
Access function: Pointer to body force function (const version)
const std::complex< double > & omega_sq() const
Access function for square of non-dim frequency.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void output(std::ostream &outfile)
Output: r,z, u_r_real, u_z_real, ..., u_theta_imag.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals for the equations (the discretised principle of virtual displacements)
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Validate against exact solution. Solution is provided via function pointer. Plot at a given number of...
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution: r,z, u_r_real, u_z_real, ..., u_theta_imag.
unsigned required_nvalue(const unsigned &n) const
Number of values required at node n.
void get_strain(const Vector< double > &s, DenseMatrix< std::complex< double >> &strain)
Get strain tensor.
virtual void fill_in_generic_contribution_to_residuals_fourier_decomp_time_harmonic_linear_elasticity(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Private helper function to compute residuals and (if requested via flag) also the Jacobian matrix.
void output(FILE *file_pt)
C-style output: r,z, u_r_real, u_z_real, ..., u_theta_imag.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
The jacobian is calculated by finite differences by default, We need only to take finite differences ...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
//////////////////////////////////////////////////////////////////// ////////////////////////////////...