axisym_linear_elasticity_elements.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1104 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-01-09 08:06:50 +0000 (Sat, 09 Jan 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 // Non-inline functions for elements that solve the equations of linear
31 // elasticity in cartesian coordinates
32 
34 
35 
36 namespace oomph
37 {
38 
39  /// \short Static default value for Young's modulus (1.0 -- for natural
40  /// scaling, i.e. all stresses have been non-dimensionalised by
41  /// the same reference Young's modulus. Setting the "non-dimensional"
42  /// Young's modulus (obtained by de-referencing Youngs_modulus_pt)
43  /// to a number larger than one means that the material is stiffer
44  /// than assumed in the non-dimensionalisation.
47 
48  /// Static default value for timescale ratio (1.0 -- for natural scaling)
50  = 1.0;
51 
52 //////////////////////////////////////////////////////////////////
53 //////////////////////////////////////////////////////////////////
54 //////////////////////////////////////////////////////////////////
55 
56 //=======================================================================
57 /// Get strain (3x3 entries; r, z, phi)
58 //=======================================================================
61  {
62  //Find out how many nodes there are
63  unsigned n_node = this->nnode();
64 
65 
66 #ifdef PARANOID
67  //Find out how many positional dofs there are
68  unsigned n_position_type = this->nnodal_position_type();
69 
70  if(n_position_type != 1)
71  {
72  throw OomphLibError(
73  "AxisymmetricLinearElasticity is not yet implemented for more than one position type",
74  OOMPH_CURRENT_FUNCTION,
75  OOMPH_EXCEPTION_LOCATION);
76  }
77 #endif
78 
79  //Find the indices at which the local displacements are stored
80  unsigned u_nodal_index[3];
81  for(unsigned i=0;i<3;i++)
82  {
83  u_nodal_index[i] = this->
85  }
86 
87  //Set up memory for the shape functions
88  Shape psi(n_node);
89 
90  // Derivs only w.r.t. r [0] and z [1]
91  DShape dpsidx(n_node,2);
92 
93  //Storage for Eulerian coordinates (r,z; initialised to zero)
95 
96  // Displacements u_r,u_z,u_theta
97  Vector<double> interpolated_u(3,0.0);
98 
99  //Calculate interpolated values of the derivatives w.r.t.
100  //Eulerian coordinates
101  DenseMatrix<double> interpolated_dudx(3,2,0.0);
102 
103 
104  //Calculate displacements and derivatives
105  for(unsigned l=0;l<n_node;l++)
106  {
107  //Calculate the coordinates
108  for(unsigned i=0;i<2;i++)
109  {
110  interpolated_x[i] += this->raw_nodal_position(l,i)*psi(l);
111  }
112 
113  //Get the nodal displacements
114  for(unsigned i=0;i<3;i++)
115  {
116  const double u_value
117  = this->raw_nodal_value(l,u_nodal_index[i]);
118  interpolated_u[i]+=u_value*psi(l);
119 
120  //Loop over derivative directions
121  for(unsigned j=0;j<2;j++)
122  {
123  interpolated_dudx(i,j) += u_value*dpsidx(l,j);
124  }
125  }
126  }
127 
128 
129  // define shorthand notation for regularly-occurring terms
130  double r = interpolated_x[0];
131 
132  // r component of displacement
133  double ur = interpolated_u[0];
134 
135  // theta component of displacement
136  double uth = interpolated_u[2];
137 
138  // derivatives w.r.t. r and z:
139  double durdr = interpolated_dudx(0,0);
140  double durdz = interpolated_dudx(0,1);
141  double duzdr = interpolated_dudx(1,0);
142  double duzdz = interpolated_dudx(1,1);
143  double duthdr = interpolated_dudx(2,0);
144  double duthdz = interpolated_dudx(2,1);
145 
146 
147  // e_rr
148  strain(0,0)=durdr;
149  // e_rz
150  strain(0,1)=durdz+duzdr;
151  strain(1,0)=durdz+duzdr;
152  // e_rphi
153  strain(0,2)=duthdr-uth/r;
154  strain(2,0)=duthdr-uth/r;
155  // e_zz
156  strain(1,1)=duzdz;
157  // e_zphi
158  strain(1,2)=duthdz;
159  strain(2,1)=duthdz;
160  // e_phiphi
161  strain(2,2)=ur/r;
162 
163  }
164 
165 //=======================================================================
166 /// Compute the residuals for the axisymmetric (in cyl. polars)
167 /// linear elasticity equations in. Flag indicates if we want Jacobian too.
168 //=======================================================================
171  Vector<double> &residuals, DenseMatrix<double> &jacobian, unsigned flag)
172  {
173  //Find out how many nodes there are
174  unsigned n_node = this->nnode();
175 
176  // Get continuous time from timestepper of first node
177  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
178 
179 #ifdef PARANOID
180  //Find out how many positional dofs there are
181  unsigned n_position_type = this->nnodal_position_type();
182 
183  if(n_position_type != 1)
184  {
185  throw OomphLibError(
186  "AxisymmetricLinearElasticity is not yet implemented for more than one position type",
187  OOMPH_CURRENT_FUNCTION,
188  OOMPH_EXCEPTION_LOCATION);
189  }
190 #endif
191 
192  //Find the indices at which the local displacements are stored
193  unsigned u_nodal_index[3];
194  for(unsigned i=0;i<3;i++)
195  {
196  u_nodal_index[i] = this->
198  }
199 
200  // Get elastic parameters
201  double nu_local = this->nu();
202  double youngs_modulus_local = this->youngs_modulus();
203 
204  // Obtain Lame parameters from Young's modulus and Poisson's ratio
205  double lambda =
206  youngs_modulus_local*nu_local/(1.0+nu_local)/(1.0-2.0*nu_local);
207  double mu = youngs_modulus_local/2.0/(1.0+nu_local);
208 
209 
210  // Lambda squared --- time scaling, NOT sqaure of Lame parameter lambda
211  const double lambda_sq = this->lambda_sq();
212 
213  //Set up memory for the shape functions
214  Shape psi(n_node);
215 
216  // Derivs only w.r.t. r [0] and z [1]
217  DShape dpsidx(n_node,2);
218 
219  //Set the value of Nintpt -- the number of integration points
220  unsigned n_intpt = this->integral_pt()->nweight();
221 
222  //Set the vector to hold the local coordinates in the element
223  Vector<double> s(2);
224 
225  //Integers to store the local equation numbers
226  int local_eqn=0, local_unknown=0;
227 
228  //Loop over the integration points
229  for(unsigned ipt=0;ipt<n_intpt;ipt++)
230  {
231  //Assign the values of s
232  for(unsigned i=0;i<2;++i)
233  {
234  s[i] = this->integral_pt()->knot(ipt,i);
235  }
236 
237  //Get the integral weight
238  double w = this->integral_pt()->weight(ipt);
239 
240  //Call the derivatives of the shape functions (and get Jacobian)
241  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
242 
243  //Storage for Eulerian coordinates (r,z; initialised to zero)
245 
246  // Displacements u_r,u_z,u_theta
248  interpolated_u(3,0.0);
249 
250  //Calculate interpolated values of the derivatives w.r.t.
251  //Eulerian coordinates
253  interpolated_dudx(3,2,0.0);
254 
255  Vector<double> d2u_dt2(3,0.0);
256 
257  //Calculate displacements and derivatives
258  for(unsigned l=0;l<n_node;l++)
259  {
260  //Calculate the coordinates
261  for(unsigned i=0;i<2;i++)
262  {
263  interpolated_x[i] += this->raw_nodal_position(l,i)*psi(l);
264  }
265  //Get the nodal displacements
266  for(unsigned i=0;i<3;i++)
267  {
268  const double u_value
269  = this->raw_nodal_value(l,u_nodal_index[i]);
270 
271  interpolated_u[i]+=u_value*psi(l);
272 
273  d2u_dt2[i]+=d2u_dt2_axisymmetric_linear_elasticity(l,i)*psi(l);
274 
275  //Loop over derivative directions
276  for(unsigned j=0;j<2;j++)
277  {
278  interpolated_dudx(i,j) += u_value*dpsidx(l,j);
279  }
280  }
281  }
282 
283  //Get body force
284  Vector<double> b(3);
285  this->body_force(time,interpolated_x,b);
286 
287  //Premultiply the weights and the Jacobian
288  double W = w*J;
289 
290 //=====EQUATIONS OF AXISYMMETRIC LINEAR ELASTICITY ========
291 
292  // define shorthand notation for regularly-occurring terms
293  double r = interpolated_x[0];
294  double rsq = pow(r,2);
295 
296  // r component of displacement
297  double ur = interpolated_u[0];
298 
299  // theta component of displacement
300  double uth = interpolated_u[2];
301 
302  // derivatives w.r.t. r and z:
303  double durdr = interpolated_dudx(0,0);
304  double durdz = interpolated_dudx(0,1);
305  double duzdr = interpolated_dudx(1,0);
306  double duzdz = interpolated_dudx(1,1);
307  double duthdr = interpolated_dudx(2,0);
308  double duthdz = interpolated_dudx(2,1);
309 
310  // storage for terms required for analytic Jacobian
311  double G_r, G_z, G_theta;
312 
313  //Loop over the test functions, nodes of the element
314  for(unsigned l=0;l<n_node;l++)
315  {
316  //Loop over the displacement components
317  for(unsigned a=0;a<3;a++)
318  {
319  //Get the equation number
320  local_eqn = this->nodal_local_eqn(l,u_nodal_index[a]);
321 
322  /*IF it's not a boundary condition*/
323  if(local_eqn >= 0)
324  {
325  // Acceleration and body force
326  residuals[local_eqn] +=
327  (lambda_sq*d2u_dt2[a] - b[a])*psi(l)*r*W;
328 
329  // Three components of the stress divergence term:
330  // a=0: r; a=1: z; a=2: theta
331 
332  // r-equation
333  if(a==0)
334  {
335  residuals[local_eqn] +=
336  (
337  mu*(
338  2.0*durdr*dpsidx(l,0)+
339  +dpsidx(l,1)*(durdz+duzdr)+
340  2.0*psi(l)/pow(r,2)*(ur)
341  )
342  +lambda*(
343  durdr+ur/r+duzdz
344  )*(dpsidx(l,0)+psi(l)/r)
345  )*r*W;
346  }
347  // z-equation
348  else if(a==1)
349  {
350  residuals[local_eqn] +=
351  (
352  mu*(
353  dpsidx(l,0)*(durdz+duzdr)
354  +
355  2.0*duzdz*dpsidx(l,1)
356  )
357  +lambda*(
358  durdr+ur/r+duzdz
359  )*dpsidx(l,1)
360  )*r*W;
361  }
362  // theta-equation
363  else if(a==2)
364  {
365  residuals[local_eqn] +=
366  (
367  mu*(
368  (duthdr-uth/r)*
369  (dpsidx(l,0)-psi(l)/r)
370  +
371  dpsidx(l,1)*(duthdz)
372  )
373  )*r*W;
374  }
375  // error: a should be 0, 1 or 2
376  else
377  {
378  throw OomphLibError(
379  "a should equal 0, 1 or 2",
380  OOMPH_CURRENT_FUNCTION,
381  OOMPH_EXCEPTION_LOCATION);
382  }
383 
384  //Jacobian entries
385  if(flag)
386  {
387  //Loop over the displacement basis functions again
388  for(unsigned l2=0;l2<n_node;l2++)
389  {
390  // define terms used to obtain entries of current row in the Jacobian:
391 
392  // terms for rows of Jacobian matrix corresponding to r-equation
393  if(a==0)
394  {
395  G_r = (mu*(2.0*dpsidx(l2,0)*dpsidx(l,0)+
396  2.0/rsq*psi(l2)*psi(l)+
397  dpsidx(l2,1)*dpsidx(l,1))+
398  lambda*(dpsidx(l2,0)+psi(l2)/r)*
399  (dpsidx(l,0)+psi(l)/r)+
400  lambda_sq*node_pt(l2)->
401  time_stepper_pt()->weight(2,0)*psi(l2)*psi(l)
402  )*r*W;
403 
404  G_z=(mu*dpsidx(l2,0)*dpsidx(l,1)+
405  lambda*dpsidx(l2,1)*(dpsidx(l,0)+psi(l)/r))*r*W;
406 
407  G_theta=0;
408  }
409 
410  // terms for rows of Jacobian matrix corresponding to z-equation
411  else if(a==1)
412  {
413  G_r=(mu*dpsidx(l2,1)*dpsidx(l,0)+
414  lambda*(dpsidx(l2,0)+psi(l2)/r)*dpsidx(l,1))*r*W;
415 
416  G_z=(mu*(dpsidx(l2,0)*dpsidx(l,0)
417  +2.0*dpsidx(l2,1)*dpsidx(l,1))
418  +lambda*dpsidx(l2,1)*dpsidx(l,1)
419  +lambda_sq*node_pt(l2)->
420  time_stepper_pt()->weight(2,0)*psi(l2)*psi(l)
421  )*r*W;
422 
423  G_theta=0;
424  }
425 
426  // terms for rows of Jacobian matrix corresponding to theta-equation
427  else if(a==2)
428  {
429  G_r=0;
430 
431  G_z=0;
432 
433  G_theta=(mu*((dpsidx(l2,0)-psi(l2)/r)*(dpsidx(l,0)-psi(l)/r)
434  +dpsidx(l2,1)*dpsidx(l,1))
435  +lambda_sq*node_pt(l2)->
436  time_stepper_pt()->weight(2,0)*psi(l2)*psi(l)
437  )*r*W;
438  }
439 
440  //Loop over the displacement components
441  for(unsigned c=0;c<3;c++)
442  {
443  // Get local unknown
444  local_unknown=this->nodal_local_eqn(l2,u_nodal_index[c]);
445 
446  //If the local unknown is not pinned
447  if(local_unknown >= 0)
448  {
449  if(c==0)
450  {
451  jacobian(local_eqn,local_unknown) += G_r;
452  }
453  else if(c==1)
454  {
455  jacobian(local_eqn,local_unknown) += G_z;
456  }
457  else if(c==2)
458  {
459  jacobian(local_eqn,local_unknown) += G_theta;
460  }
461  }
462  }
463  }
464  } //End of jacobian calculation
465 
466  } //End of if not boundary condition
467  } //End of loop over coordinate directions
468  } //End of loop over shape functions
469  } //End of loop over integration points
470 
471  }
472 
473 //=======================================================================
474 /// Output exact solution r,z, u_r, u_z, u_theta
475 //=======================================================================
477  std::ostream &outfile,
478  const unsigned &nplot,
480  {
481  //Vector of local coordinates
482  Vector<double> s(2);
483 
484  // Vector for coordintes
485  Vector<double> x(2);
486 
487  // Tecplot header info
488  outfile << this->tecplot_zone_string(nplot);
489 
490  // Exact solution Vector
491  Vector<double> exact_soln(9);
492 
493  // Loop over plot points
494  unsigned num_plot_points=this->nplot_points(nplot);
495  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
496  {
497 
498  // Get local coordinates of plot point
499  this->get_s_plot(iplot,nplot,s);
500 
501  // Get x position as Vector
502  this->interpolated_x(s,x);
503 
504  // Get exact solution at this point
505  (*exact_soln_pt)(x,exact_soln);
506 
507  //Output r,z,...,u_exact,...
508  for(unsigned i=0;i<2;i++)
509  {
510  outfile << x[i] << " ";
511  }
512 
513  for(unsigned i=0;i<3;i++)
514  {
515  outfile << exact_soln[i] << " ";
516  }
517  outfile << std::endl;
518  }
519 
520  // Write tecplot footer (e.g. FE connectivity lists)
521  this->write_tecplot_zone_footer(outfile,nplot);
522 
523  }
524 
525 //=======================================================================
526 /// Output exact solution r,z, u_r, u_z, u_theta
527 /// Time dependent version
528 //=======================================================================
530  std::ostream &outfile,
531  const unsigned &nplot,
532  const double &time,
534  {
535  //Vector of local coordinates
536  Vector<double> s(2);
537 
538  // Vector for coordintes
539  Vector<double> x(2);
540 
541  // Tecplot header info
542  outfile << this->tecplot_zone_string(nplot);
543 
544  // Exact solution Vector
545  Vector<double> exact_soln(9);
546 
547  // Loop over plot points
548  unsigned num_plot_points=this->nplot_points(nplot);
549  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
550  {
551 
552  // Get local coordinates of plot point
553  this->get_s_plot(iplot,nplot,s);
554 
555  // Get x position as Vector
556  this->interpolated_x(s,x);
557 
558  // Get exact solution at this point
559  (*exact_soln_pt)(time,x,exact_soln);
560 
561  //Output r,z,...,u_exact,...
562  for(unsigned i=0;i<2;i++)
563  {
564  outfile << x[i] << " ";
565  }
566 
567  for(unsigned i=0;i<9;i++)
568  {
569  outfile << exact_soln[i] << " ";
570  }
571  outfile << std::endl;
572  }
573 
574  // Write tecplot footer (e.g. FE connectivity lists)
575  this->write_tecplot_zone_footer(outfile,nplot);
576 
577  }
578 
579 //=======================================================================
580 /// Output: r,z, u_r, u_z, u_theta
581 //=======================================================================
583  output(std::ostream &outfile,
584  const unsigned &nplot)
585  {
586  //Set output Vector
587  Vector<double> s(2);
588  Vector<double> x(2);
589  Vector<double> u(3);
590  Vector<double> du_dt(3);
591  Vector<double> d2u_dt2(3);
592 
593 
594  // Tecplot header info
595  outfile << this->tecplot_zone_string(nplot);
596 
597  // Loop over plot points
598  unsigned num_plot_points=this->nplot_points(nplot);
599  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
600  {
601 
602  // Get local coordinates of plot point
603  this->get_s_plot(iplot,nplot,s);
604 
605  // Get Eulerian coordinates and displacements
606  this->interpolated_x(s,x);
610 
611  //Output the r,z,..
612  for(unsigned i=0;i<2;i++)
613  {outfile << x[i] << " ";}
614 
615  // Output displacement
616  for(unsigned i=0;i<3;i++)
617  {outfile << u[i] << " ";}
618 
619  // Output veloc
620  for(unsigned i=0;i<3;i++)
621  {outfile << du_dt[i] << " ";}
622 
623  // Output accel
624  for(unsigned i=0;i<3;i++)
625  {outfile << d2u_dt2[i] << " ";}
626 
627  outfile << std::endl;
628  }
629 
630  // Write tecplot footer (e.g. FE connectivity lists)
631  this->write_tecplot_zone_footer(outfile,nplot);
632  }
633 
634 
635 //=======================================================================
636 /// C-style output:r,z, u_r, u_z, u_theta
637 //=======================================================================
639  FILE* file_pt,
640  const unsigned &nplot)
641 {
642  //Vector of local coordinates
643  Vector<double> s(2);
644 
645  // Tecplot header info
646  fprintf(file_pt,"%s",this->tecplot_zone_string(nplot).c_str());
647 
648  // Loop over plot points
649  unsigned num_plot_points=this->nplot_points(nplot);
650  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
651  {
652  // Get local coordinates of plot point
653  this->get_s_plot(iplot,nplot,s);
654 
655  // Coordinates
656  for(unsigned i=0;i<2;i++)
657  {
658  fprintf(file_pt,"%g ",this->interpolated_x(s,i));
659  }
660 
661  // Displacement
662  for(unsigned i=0;i<3;i++)
663  {
664  fprintf(file_pt,"%g ",
665  this->
667  }
668  }
669  fprintf(file_pt,"\n");
670 
671  // Write tecplot footer (e.g. FE connectivity lists)
672  this->write_tecplot_zone_footer(file_pt,nplot);
673 
674 }
675 
676 //======================================================================
677 /// Validate against exact solution
678 /// Solution is provided via function pointer.
679 /// Plot at a given number of plot points and compute L2 error
680 /// and L2 norm of velocity solution over element.
681 //=======================================================================
683  std::ostream &outfile,
685  double& error, double& norm)
686 {
687 
688  error=0.0;
689  norm=0.0;
690 
691  //Vector of local coordinates
692  Vector<double> s(2);
693 
694  // Vector for coordinates
695  Vector<double> x(2);
696 
697  //Set the value of n_intpt
698  unsigned n_intpt = this->integral_pt()->nweight();
699 
700  outfile << "ZONE" << std::endl;
701 
702  // Exact solution Vector (u_r, u_z, u_theta)
703  Vector<double> exact_soln(9);
704 
705  //Loop over the integration points
706  for(unsigned ipt=0;ipt<n_intpt;ipt++)
707  {
708 
709  //Assign values of s
710  for(unsigned i=0;i<2;i++)
711  {
712  s[i] = this->integral_pt()->knot(ipt,i);
713  }
714 
715  //Get the integral weight
716  double w = this->integral_pt()->weight(ipt);
717 
718  // Get jacobian of mapping
719  double J= this->J_eulerian(s);
720 
721  //Premultiply the weights and the Jacobian
722  double W = w*J;
723 
724  // Get x position as Vector
725  this->interpolated_x(s,x);
726 
727  // Get exact solution at this point
728  (*exact_soln_pt)(x,exact_soln);
729 
730  // Displacement error
731  for(unsigned i=0;i<3;i++)
732  {
733  norm+=(exact_soln[i]*exact_soln[i])*W;
734  error+=(
735  (exact_soln[i]-
736  this->
738  (exact_soln[i]-
739  this->
741  )*W;
742  }
743 
744 
745  //Output r,z coordinates
746  for(unsigned i=0;i<2;i++)
747  {
748  outfile << x[i] << " ";
749  }
750 
751  //Output ur_error, uz_error, uth_error
752  for(unsigned i=0;i<3;i++)
753  {
754  outfile << exact_soln[i]-
755  this->
757  << " ";
758  }
759  outfile << std::endl;
760  }
761 }
762 
763 //======================================================================
764 /// Validate against exact solution
765 /// Solution is provided via function pointer.
766 /// Plot at a given number of plot points and compute L2 error
767 /// and L2 norm of velocity solution over element.
768 //=======================================================================
770  std::ostream &outfile,
772  const double& time, double& error, double& norm)
773 {
774 
775  error=0.0;
776  norm=0.0;
777 
778  //Vector of local coordinates
779  Vector<double> s(2);
780 
781  // Vector for coordinates
782  Vector<double> x(2);
783 
784  //Set the value of n_intpt
785  unsigned n_intpt = this->integral_pt()->nweight();
786 
787  outfile << "ZONE" << std::endl;
788 
789  // Exact solution Vector (u_r, u_z, u_theta)
790  Vector<double> exact_soln(9);
791 
792  //Loop over the integration points
793  for(unsigned ipt=0;ipt<n_intpt;ipt++)
794  {
795 
796  //Assign values of s
797  for(unsigned i=0;i<2;i++)
798  {
799  s[i] = this->integral_pt()->knot(ipt,i);
800  }
801 
802  //Get the integral weight
803  double w = this->integral_pt()->weight(ipt);
804 
805  // Get jacobian of mapping
806  double J= this->J_eulerian(s);
807 
808  //Premultiply the weights and the Jacobian
809  double W = w*J;
810 
811  // Get x position as Vector
812  this->interpolated_x(s,x);
813 
814  // Get exact solution at this point
815  (*exact_soln_pt)(time,x,exact_soln);
816 
817  // Displacement error
818  for(unsigned i=0;i<3;i++)
819  {
820  norm+=(exact_soln[i]*exact_soln[i])*W;
821  error+=(
822  (exact_soln[i]-
823  this->
825  (exact_soln[i]-
826  this->
828  )*W;
829  }
830 
831 
832  //Output r,z coordinates
833  for(unsigned i=0;i<2;i++)
834  {
835  outfile << x[i] << " ";
836  }
837 
838  //Output ur_error, uz_error, uth_error
839  for(unsigned i=0;i<3;i++)
840  {
841  outfile << exact_soln[i]-
842  this->
844  << " ";
845  }
846  outfile << std::endl;
847  }
848 }
849 
850 // Instantiate required elements
854 
855 
856 }
const double & lambda_sq() const
Access function for timescale ratio (nondim density)
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
void body_force(const double &time, const Vector< double > &x, Vector< double > &b) const
Evaluate body force at Eulerian coordinate x at present time (returns zero vector if no body force fu...
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
static double Default_lambda_sq_value
Static default value for timescale ratio (1.0 for natural scaling)
static double Default_youngs_modulus_value
Static default value for Young's modulus (1.0 – for natural scaling, i.e. all stresses have been non-...
double d2u_dt2_axisymmetric_linear_elasticity(const unsigned &n, const unsigned &i) const
d^2u/dt^2 at local node n
double youngs_modulus() const
Access function to Young's modulus.
cstr elem_len * i
Definition: cfortran.h:607
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
virtual unsigned u_index_axisymmetric_linear_elasticity(const unsigned &i) const
Return the index at which the i-th (i=0: r, i=1: z; i=2: theta) unknown displacement component is sto...
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain)
Get strain (3x3 entries; r, z, phi)
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2340
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:3981
virtual void fill_in_generic_contribution_to_residuals_axisymmetric_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...
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:557
void output(std::ostream &outfile)
Output: r,z, u_r, u_z, u_theta.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution: r,z, u_r, u_z, u_theta.
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation. NOTE: Moved to cc file because of a possible compiler bug in gcc (yes, really!). The move to the cc file avoids inlining which appears to cause problems (only) when compiled with gcc and -O3; offensive "illegal read" is in optimised-out section of code and data that is allegedly illegal is readily readable (by other means) just before this function is called so I can't really see how we could possibly be responsible for this...
Definition: elements.cc:1657
void interpolated_du_dt_axisymmetric_linear_elasticity(const Vector< double > &s, Vector< double > &du_dt) const
Compute vector of FE interpolated velocity du/dt at local coordinate s.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
void interpolated_d2u_dt2_axisymmetric_linear_elasticity(const Vector< double > &s, Vector< double > &d2u_dt2) const
Compute vector of FE interpolated accel d2u/dt2 at local coordinate s.
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
static char t char * s
Definition: cfortran.h:572
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
void interpolated_u_axisymmetric_linear_elasticity(const Vector< double > &s, Vector< double > &disp) const
Compute vector of FE interpolated displacement u at local coordinate s.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
double & nu() const
Access function for Poisson's ratio.
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2446
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3252