airy_cantilever.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: 1216 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-08-06 15:23:30 +0100 (Sat, 06 Aug 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 //Driver for Airy cantilever beam problem
31 
32 //Oomph-lib includes
33 #include "generic.h"
34 #include "solid.h"
35 #include "constitutive.h"
36 
37 //The mesh
38 #include "meshes/rectangular_quadmesh.h"
39 
40 using namespace std;
41 
42 using namespace oomph;
43 
44 namespace oomph
45 {
46 
47 //=================start_wrapper==================================
48 /// Wrapper class for solid elements to modify their output
49 /// functions.
50 //================================================================
51 template <class ELEMENT>
52 class MySolidElement : public virtual ELEMENT
53 {
54 
55 public:
56 
57  /// Constructor: Call constructor of underlying element
58  MySolidElement() : ELEMENT() {};
59 
60  /// Overload output function:
61  void output(std::ostream &outfile, const unsigned &n_plot)
62  {
63 
64  // Element dimension
65  unsigned el_dim = this->dim();
66 
67  Vector<double> s(el_dim);
68  Vector<double> x(el_dim);
69  DenseMatrix<double> sigma(el_dim,el_dim);
70 
71  switch(el_dim)
72  {
73 
74  case 2:
75 
76  //Tecplot header info
77  outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
78 
79  //Loop over element nodes
80  for(unsigned l2=0;l2<n_plot;l2++)
81  {
82  s[1] = -1.0 + l2*2.0/(n_plot-1);
83  for(unsigned l1=0;l1<n_plot;l1++)
84  {
85  s[0] = -1.0 + l1*2.0/(n_plot-1);
86 
87  // Get Eulerian coordinates and stress
88  this->interpolated_x(s,x);
89  this->get_stress(s,sigma);
90 
91  //Output the x,y,..
92  for(unsigned i=0;i<el_dim;i++)
93  {outfile << x[i] << " ";}
94 
95  // Output stress
96  outfile << sigma(0,0) << " "
97  << sigma(1,0) << " "
98  << sigma(1,1) << " "
99  << std::endl;
100  }
101  }
102 
103  break;
104 
105  default:
106 
107  std::ostringstream error_message;
108  error_message << "Output for dim !=2 not implemented" << std::endl;
109  throw OomphLibError(error_message.str(),
110  OOMPH_CURRENT_FUNCTION,
111  OOMPH_EXCEPTION_LOCATION);
112  }
113 
114  }
115 
116 };
117 
118 
119 
120 //===========start_face_geometry==============================================
121 /// FaceGeometry of wrapped element is the same as the underlying element
122 //============================================================================
123 template<class ELEMENT>
124 class FaceGeometry<MySolidElement<ELEMENT> > :
125  public virtual FaceGeometry<ELEMENT>
126 {
127 public:
128 
129  /// \short Constructor [this was only required explicitly
130  /// from gcc 4.5.2 onwards...]
131  FaceGeometry() : FaceGeometry<ELEMENT>() {}
132 
133 };
134 
135 
136 }
137 
138 
139 
140 
141 
142 
143 ///////////////////////////////////////////////////////////////////////
144 ///////////////////////////////////////////////////////////////////////
145 ///////////////////////////////////////////////////////////////////////
146 
147 
148 
149 //=======start_namespace==========================================
150 /// Global variables
151 //================================================================
152 namespace Global_Physical_Variables
153 {
154 
155  /// Half height of beam
156  double H=0.5;
157 
158  /// Length of beam
159  double L=10.0;
160 
161  /// Pointer to constitutive law
162  ConstitutiveLaw* Constitutive_law_pt;
163 
164  /// Elastic modulus
165  double E=1.0;
166 
167  /// Poisson's ratio
168  double Nu=0.3;
169 
170  /// Uniform pressure
171  double P = 0.0;
172 
173  /// \short Constant pressure load. The arguments to this function are imposed
174  /// on us by the SolidTractionElements which allow the traction to
175  /// depend on the Lagrangian and Eulerian coordinates x and xi, and on the
176  /// outer unit normal to the surface. Here we only need the outer unit
177  /// normal.
178  void constant_pressure(const Vector<double> &xi, const Vector<double> &x,
179  const Vector<double> &n, Vector<double> &traction)
180  {
181  unsigned dim = traction.size();
182  for(unsigned i=0;i<dim;i++)
183  {
184  traction[i] = -P*n[i];
185  }
186  } // end traction
187 
188 
189  /// Non-dim gravity
190  double Gravity=0.0;
191 
192  /// Non-dimensional gravity as body force
193  void gravity(const double& time,
194  const Vector<double> &xi,
195  Vector<double> &b)
196  {
197  b[0]=0.0;
198  b[1]=-Gravity;
199  }
200 
201 } //end namespace
202 
203 
204 
205 //=============begin_problem============================================
206 /// Problem class for the cantilever "beam" structure.
207 //======================================================================
208 template<class ELEMENT>
209 class CantileverProblem : public Problem
210 {
211 
212 public:
213 
214  /// Constructor:
216 
217  /// Update function (empty)
219 
220  /// Update function (empty)
222 
223  /// Access function for the solid mesh
224  ElasticRefineableRectangularQuadMesh<ELEMENT>*& solid_mesh_pt()
225  {return Solid_mesh_pt;}
226 
227  /// Access function to the mesh of surface traction elements
228  SolidMesh*& traction_mesh_pt(){return Traction_mesh_pt;}
229 
230  /// Actions before adapt: Wipe the mesh of traction elements
231  void actions_before_adapt();
232 
233  /// Actions after adapt: Rebuild the mesh of traction elements
234  void actions_after_adapt();
235 
236  /// Doc the solution
237  void doc_solution();
238 
239 private:
240 
241  /// \short Pass pointer to traction function to the
242  /// elements in the traction mesh
243  void set_traction_pt();
244 
245  /// \short Create traction elements
246  void create_traction_elements();
247 
248  /// Delete traction elements
249  void delete_traction_elements();
250 
251  /// Trace file
252  ofstream Trace_file;
253 
254  /// Pointers to node whose position we're tracing
256 
257  /// Pointer to solid mesh
258  ElasticRefineableRectangularQuadMesh<ELEMENT>* Solid_mesh_pt;
259 
260  /// Pointer to mesh of traction elements
261  SolidMesh* Traction_mesh_pt;
262 
263  /// DocInfo object for output
264  DocInfo Doc_info;
265 
266 };
267 
268 
269 //===========start_of_constructor=======================================
270 /// Constructor:
271 //======================================================================
272 template<class ELEMENT>
274 {
275 
276  // Create the mesh
277 
278  // # of elements in x-direction
279  unsigned n_x=20;
280 
281  // # of elements in y-direction
282  unsigned n_y=2;
283 
284  // Domain length in x-direction
285  double l_x= Global_Physical_Variables::L;
286 
287  // Domain length in y-direction
288  double l_y=2.0*Global_Physical_Variables::H;
289 
290  // Shift mesh downwards so that centreline is at y=0:
291  Vector<double> origin(2);
292  origin[0]=0.0;
293  origin[1]=-0.5*l_y;
294 
295  //Now create the mesh
296  solid_mesh_pt() = new ElasticRefineableRectangularQuadMesh<ELEMENT>(
297  n_x,n_y,l_x,l_y,origin);
298 
299  // Set error estimator
300  solid_mesh_pt()->spatial_error_estimator_pt()=new Z2ErrorEstimator;
301 
302  //Assign the physical properties to the elements before any refinement
303  //Loop over the elements in the main mesh
304  unsigned n_element =solid_mesh_pt()->nelement();
305  for(unsigned i=0;i<n_element;i++)
306  {
307  //Cast to a solid element
308  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(solid_mesh_pt()->element_pt(i));
309 
310  // Set the constitutive law
311  el_pt->constitutive_law_pt() =
313 
314  //Set the body force
315  el_pt->body_force_fct_pt() = Global_Physical_Variables::gravity;
316  }
317 
318 
319  // Choose a control node: The last node in the solid mesh
320  unsigned nnod=solid_mesh_pt()->nnode();
321  Trace_node_pt=solid_mesh_pt()->node_pt(nnod-1);
322 
323  // Refine the mesh uniformly
324  solid_mesh_pt()->refine_uniformly();
325 
326  // Construct the traction element mesh
327  Traction_mesh_pt=new SolidMesh;
328  create_traction_elements();
329 
330  // Pass pointer to traction function to the elements
331  // in the traction mesh
332  set_traction_pt();
333 
334  // Solid mesh is first sub-mesh
335  add_sub_mesh(solid_mesh_pt());
336 
337  // Add traction sub-mesh
338  add_sub_mesh(traction_mesh_pt());
339 
340  // Build combined "global" mesh
341  build_global_mesh();
342 
343  // Pin the left boundary (boundary 3) in both directions
344  unsigned n_side = mesh_pt()->nboundary_node(3);
345 
346  // Loop over the nodes
347  for(unsigned i=0;i<n_side;i++)
348  {
349  solid_mesh_pt()->boundary_node_pt(3,i)->pin_position(0);
350  solid_mesh_pt()->boundary_node_pt(3,i)->pin_position(1);
351  }
352 
353  // Pin the redundant solid pressures (if any)
354  PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
355  solid_mesh_pt()->element_pt());
356 
357  //Attach the boundary conditions to the mesh
358  cout << assign_eqn_numbers() << std::endl;
359 
360  // Set output directory
361  Doc_info.set_directory("RESLT");
362 
363  // Open trace file
364  char filename[100];
365  sprintf(filename,"%s/trace.dat",Doc_info.directory().c_str());
366  Trace_file.open(filename);
367 
368 
369 } //end of constructor
370 
371 
372 //=====================start_of_actions_before_adapt======================
373 /// Actions before adapt: Wipe the mesh of traction elements
374 //========================================================================
375 template<class ELEMENT>
377 {
378  // Kill the traction elements and wipe surface mesh
379  delete_traction_elements();
380 
381  // Rebuild the Problem's global mesh from its various sub-meshes
382  rebuild_global_mesh();
383 
384 }// end of actions_before_adapt
385 
386 
387 
388 //=====================start_of_actions_after_adapt=======================
389 /// Actions after adapt: Rebuild the mesh of traction elements
390 //========================================================================
391 template<class ELEMENT>
393 {
394  // Create traction elements from all elements that are
395  // adjacent to boundary 2 and add them to surface meshes
396  create_traction_elements();
397 
398  // Rebuild the Problem's global mesh from its various sub-meshes
399  rebuild_global_mesh();
400 
401  // Pin the redundant solid pressures (if any)
402  PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
403  solid_mesh_pt()->element_pt());
404 
405  // Set pointer to prescribed traction function for traction elements
406  set_traction_pt();
407 
408 }// end of actions_after_adapt
409 
410 
411 
412 //==================start_of_set_traction_pt==============================
413 /// Set pointer to traction function for the relevant
414 /// elements in the traction mesh
415 //========================================================================
416 template<class ELEMENT>
418 {
419  // Loop over the elements in the traction element mesh
420  // for elements on the top boundary (boundary 2)
421  unsigned n_element=traction_mesh_pt()->nelement();
422  for(unsigned i=0;i<n_element;i++)
423  {
424  //Cast to a solid traction element
425  SolidTractionElement<ELEMENT> *el_pt =
426  dynamic_cast<SolidTractionElement<ELEMENT>*>
427  (traction_mesh_pt()->element_pt(i));
428 
429  //Set the traction function
430  el_pt->traction_fct_pt() = Global_Physical_Variables::constant_pressure;
431  }
432 
433 }// end of set traction pt
434 
435 
436 
437 //============start_of_create_traction_elements==============================
438 /// Create traction elements
439 //=======================================================================
440 template<class ELEMENT>
442 {
443  // Traction elements are located on boundary 2:
444  unsigned b=2;
445 
446  // How many bulk elements are adjacent to boundary b?
447  unsigned n_element = solid_mesh_pt()->nboundary_element(b);
448 
449  // Loop over the bulk elements adjacent to boundary b
450  for(unsigned e=0;e<n_element;e++)
451  {
452  // Get pointer to the bulk element that is adjacent to boundary b
453  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
454  solid_mesh_pt()->boundary_element_pt(b,e));
455 
456  //Find the index of the face of element e along boundary b
457  int face_index = solid_mesh_pt()->face_index_at_boundary(b,e);
458 
459  // Create new element and add to mesh
460  Traction_mesh_pt->add_element_pt(new SolidTractionElement<ELEMENT>
461  (bulk_elem_pt,face_index));
462  }
463 
464  // Pass the pointer to the traction function to the traction elements
465  set_traction_pt();
466 
467 } // end of create_traction_elements
468 
469 
470 
471 
472 //============start_of_delete_traction_elements==============================
473 /// Delete traction elements and wipe the traction meshes
474 //=======================================================================
475 template<class ELEMENT>
477 {
478  // How many surface elements are in the surface mesh
479  unsigned n_element = Traction_mesh_pt->nelement();
480 
481  // Loop over the surface elements
482  for(unsigned e=0;e<n_element;e++)
483  {
484  // Kill surface element
485  delete Traction_mesh_pt->element_pt(e);
486  }
487 
488  // Wipe the mesh
489  Traction_mesh_pt->flush_element_and_node_storage();
490 
491 } // end of delete_traction_elements
492 
493 
494 
495 //==============start_doc===========================================
496 /// Doc the solution
497 //==================================================================
498 template<class ELEMENT>
500 {
501 
502  ofstream some_file;
503  char filename[100];
504 
505  // Number of plot points
506  unsigned n_plot = 5;
507 
508  // Output shape of and stress in deformed body
509  //--------------------------------------------
510  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
511  Doc_info.number());
512  some_file.open(filename);
513  solid_mesh_pt()->output(some_file,n_plot);
514  some_file.close();
515 
516 
517  // Output St. Venant solution
518  //---------------------------
519  sprintf(filename,"%s/exact_soln%i.dat",Doc_info.directory().c_str(),
520  Doc_info.number());
521  some_file.open(filename);
522 
523  // Element dimension
524  unsigned el_dim=2;
525 
526  Vector<double> s(el_dim);
527  Vector<double> x(el_dim);
528  Vector<double> xi(el_dim);
529  DenseMatrix<double> sigma(el_dim,el_dim);
530 
531  // Constants for exact (St. Venant) solution
532  double a=-1.0/4.0*Global_Physical_Variables::P;
534  double c=1.0/8.0*Global_Physical_Variables::P/
537 
538  // Loop over all elements to plot exact solution for stresses
539  unsigned nel=solid_mesh_pt()->nelement();
540  for (unsigned e=0;e<nel;e++)
541  {
542  // Get pointer to element
543  SolidFiniteElement* el_pt=dynamic_cast<SolidFiniteElement*>(
544  solid_mesh_pt()->element_pt(e));
545 
546  //Tecplot header info
547  some_file << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
548 
549  //Loop over plot points
550  for(unsigned l2=0;l2<n_plot;l2++)
551  {
552  s[1] = -1.0 + l2*2.0/(n_plot-1);
553  for(unsigned l1=0;l1<n_plot;l1++)
554  {
555  s[0] = -1.0 + l1*2.0/(n_plot-1);
556 
557  // Get Eulerian and Lagrangian coordinates
558  el_pt->interpolated_x(s,x);
559  el_pt->interpolated_xi(s,xi);
560 
561  //Output the x,y,..
562  for(unsigned i=0;i<el_dim;i++)
563  {some_file << x[i] << " ";}
564 
565  // Change orientation of coordinate system relative
566  // to solution in lecture notes
567  double xx=Global_Physical_Variables::L-xi[0];
568  double yy=xi[1];
569 
570  // Approximate analytical (St. Venant) solution
571  sigma(0,0)=c*(6.0*xx*xx*yy-4.0*yy*yy*yy)+
572  6.0*d*yy;
573  sigma(1,1)=2.0*(a+b*yy+c*yy*yy*yy);
574  sigma(1,0)=2.0*(b*xx+3.0*c*xx*yy*yy);
575  sigma(0,1)=sigma(1,0);
576 
577  // Output stress
578  some_file << sigma(0,0) << " "
579  << sigma(1,0) << " "
580  << sigma(1,1) << " "
581  << std::endl;
582  }
583  }
584  }
585  some_file.close();
586 
587  // Write trace file: Load/displacement characteristics
588  Trace_file << Global_Physical_Variables::P << " "
589  << Trace_node_pt->x(0) << " "
590  << Trace_node_pt->x(1) << " "
591  << std::endl;
592 
593  // Increment label for output files
594  Doc_info.number()++;
595 
596 } //end doc
597 
598 
599 
600 //=======start_of_main==================================================
601 /// Driver for cantilever beam loaded by surface traction and/or
602 /// gravity
603 //======================================================================
604 int main()
605 {
606 
607  // Create generalised Hookean constitutive equations
609  new GeneralisedHookean(&Global_Physical_Variables::Nu,
611 
612  //Set up the problem
614 
615 
616  // Uncomment these as an exercise
617 
618  // CantileverProblem<MySolidElement<
619  // RefineableQPVDElementWithContinuousPressure<2> > > problem;
620 
621  // CantileverProblem<MySolidElement<
622  // RefineableQPVDElementWithPressure<2> > > problem;
623 
624  // Initial values for parameter values
627 
628  // Max. number of adaptations per solve
629  unsigned max_adapt=3;
630 
631  //Parameter incrementation
632  unsigned nstep=5;
633  double p_increment=1.0e-5;
634  for(unsigned i=0;i<nstep;i++)
635  {
636  // Increment pressure load
637  Global_Physical_Variables::P+=p_increment;
638 
639  // Solve the problem with Newton's method, allowing
640  // up to max_adapt mesh adaptations after every solve.
641  problem.newton_solve(max_adapt);
642 
643  // Doc solution
644  problem.doc_solution();
645 
646  }
647 
648 } //end of main
649 
650 
651 
652 
653 
654 
655 
656 
ElasticRefineableRectangularQuadMesh< ELEMENT > *& solid_mesh_pt()
Access function for the solid mesh.
void delete_traction_elements()
Delete traction elements.
double Gravity
Non-dim gravity.
void set_traction_pt()
Pass pointer to traction function to the elements in the traction mesh.
Problem class for the cantilever "beam" structure.
ConstitutiveLaw * Constitutive_law_pt
Pointer to constitutive law.
SolidMesh * Traction_mesh_pt
Pointer to mesh of traction elements.
double P
Uniform pressure.
MySolidElement()
Constructor: Call constructor of underlying element.
void actions_before_newton_solve()
Update function (empty)
double L
Length of beam.
ElasticRefineableRectangularQuadMesh< ELEMENT > * Solid_mesh_pt
Pointer to solid mesh.
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of traction elements.
void gravity(const double &time, const Vector< double > &xi, Vector< double > &b)
Non-dimensional gravity as body force.
void constant_pressure(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Constant pressure load. The arguments to this function are imposed on us by the SolidTractionElements...
SolidMesh *& traction_mesh_pt()
Access function to the mesh of surface traction elements.
void output(std::ostream &outfile, const unsigned &n_plot)
Overload output function:
double H
Half height of beam.
void actions_after_newton_solve()
Update function (empty)
ofstream Trace_file
Trace file.
void doc_solution()
Doc the solution.
FaceGeometry()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
int main()
Node * Trace_node_pt
Pointers to node whose position we're tracing.
CantileverProblem()
Constructor:
void actions_before_adapt()
Actions before adapt: Wipe the mesh of traction elements.
void create_traction_elements()
Create traction elements.
double Nu
Poisson's ratio.
DocInfo Doc_info
DocInfo object for output.
double E
Elastic modulus.