airy_cantilever2.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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 //#define REFINE
45 
46 namespace oomph
47 {
48 
49 //=================start_wrapper==================================
50 /// Wrapper class for solid elements to modify their output
51 /// functions.
52 //================================================================
53 template <class ELEMENT>
54 class MySolidElement : public virtual ELEMENT
55 {
56 
57 public:
58 
59  /// Constructor: Call constructor of underlying element
60  MySolidElement() : ELEMENT() {};
61 
62  /// Overload output function:
63  void output(std::ostream &outfile, const unsigned &n_plot)
64  {
65 
66  // Element dimension
67  unsigned el_dim = this->dim();
68 
69  Vector<double> s(el_dim);
70  Vector<double> x(el_dim);
71  DenseMatrix<double> sigma(el_dim,el_dim);
72 
73  switch(el_dim)
74  {
75 
76  case 2:
77 
78  //Tecplot header info
79  outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
80 
81  //Loop over element nodes
82  for(unsigned l2=0;l2<n_plot;l2++)
83  {
84  s[1] = -1.0 + l2*2.0/(n_plot-1);
85  for(unsigned l1=0;l1<n_plot;l1++)
86  {
87  s[0] = -1.0 + l1*2.0/(n_plot-1);
88 
89  // Get Eulerian coordinates and stress
90  this->interpolated_x(s,x);
91  this->get_stress(s,sigma);
92 
93  //Output the x,y,..
94  for(unsigned i=0;i<el_dim;i++)
95  {outfile << x[i] << " ";}
96 
97  // Output stress
98  outfile << sigma(0,0) << " "
99  << sigma(1,0) << " "
100  << sigma(1,1) << " "
101  << std::endl;
102  }
103  }
104 
105  break;
106 
107  default:
108 
109  std::ostringstream error_message;
110  error_message << "Output for dim !=2 not implemented" << std::endl;
111  throw OomphLibError(error_message.str(),
112  OOMPH_CURRENT_FUNCTION,
113  OOMPH_EXCEPTION_LOCATION);
114  }
115 
116  }
117 
118 };
119 
120 
121 
122 //===========start_face_geometry==============================================
123 /// FaceGeometry of wrapped element is the same as the underlying element
124 //============================================================================
125 template<class ELEMENT>
126 class FaceGeometry<MySolidElement<ELEMENT> > :
127  public virtual FaceGeometry<ELEMENT>
128 {
129 
130 public:
131 
132  /// \short Constructor [this was only required explicitly
133  /// from gcc 4.5.2 onwards...]
134  FaceGeometry() : FaceGeometry<ELEMENT>() {}
135 
136 };
137 
138 
139 }
140 
141 
142 
143 
144 
145 
146 ///////////////////////////////////////////////////////////////////////
147 ///////////////////////////////////////////////////////////////////////
148 ///////////////////////////////////////////////////////////////////////
149 
150 
151 
152 //=======start_namespace==========================================
153 /// Global variables
154 //================================================================
155 namespace Global_Physical_Variables
156 {
157 
158  /// Pointer to strain energy function
159  StrainEnergyFunction*Strain_energy_function_pt;
160 
161  /// "Mooney Rivlin" coefficient for generalised Mooney Rivlin law
162  double C1=1.3;
163 
164  /// "Mooney Rivlin" coefficient for generalised Mooney Rivlin law
165  double C2=1.3;
166 
167  /// Half height of beam
168  double H=0.5;
169 
170  /// Length of beam
171  double L=10.0;
172 
173  /// Pointer to constitutive law
174  ConstitutiveLaw* Constitutive_law_pt;
175 
176  /// Elastic modulus
177  double E=1.0;
178 
179  /// Poisson's ratio
180  double Nu=0.3;
181 
182  /// Uniform pressure
183  double P = 0.0;
184 
185  /// \short Constant pressure load. The arguments to this function are imposed
186  /// on us by the SolidTractionElements which allow the traction to
187  /// depend on the Lagrangian and Eulerian coordinates x and xi, and on the
188  /// outer unit normal to the surface. Here we only need the outer unit
189  /// normal.
190  void constant_pressure(const Vector<double> &xi, const Vector<double> &x,
191  const Vector<double> &n, Vector<double> &traction)
192  {
193  unsigned dim = traction.size();
194  for(unsigned i=0;i<dim;i++)
195  {
196  traction[i] = -P*n[i];
197  }
198  } // end traction
199 
200 
201  /// Non-dim gravity
202  double Gravity=0.0;
203 
204  /// Non-dimensional gravity as body force
205  void gravity(const double& time,
206  const Vector<double> &xi,
207  Vector<double> &b)
208  {
209  b[0]=0.0;
210  b[1]=-Gravity;
211  }
212 
213 } //end namespace
214 
215 
216 
217 //=============begin_problem============================================
218 /// Problem class for the cantilever "beam" structure.
219 //======================================================================
220 template<class ELEMENT>
221 class CantileverProblem : public Problem
222 {
223 
224 public:
225 
226  /// Constructor:
227  CantileverProblem(const bool& incompress, const bool& use_fd);
228 
229  /// Update function (empty)
231 
232  /// Update function (empty)
234 
235 #ifdef REFINE
236 
237  /// Access function for the solid mesh
238  ElasticRefineableRectangularQuadMesh<ELEMENT>*& solid_mesh_pt()
239  {return Solid_mesh_pt;}
240 
241 #else
242 
243  /// Access function for the solid mesh
244  ElasticRectangularQuadMesh<ELEMENT>*& solid_mesh_pt()
245  {return Solid_mesh_pt;}
246 
247 #endif
248 
249 
250  /// Access function to the mesh of surface traction elements
251  SolidMesh*& traction_mesh_pt(){return Traction_mesh_pt;}
252 
253  /// Actions before adapt: Wipe the mesh of traction elements
254  void actions_before_adapt();
255 
256  /// Actions after adapt: Rebuild the mesh of traction elements
257  void actions_after_adapt();
258 
259  /// Doc the solution
260  void doc_solution();
261 
262  /// Run the job -- doc in RESLTi_case
263  void run_it(const unsigned& i_case);
264 
265 private:
266 
267  /// \short Pass pointer to traction function to the
268  /// elements in the traction mesh
269  void set_traction_pt();
270 
271  /// \short Create traction elements
272  void create_traction_elements();
273 
274  /// Delete traction elements
275  void delete_traction_elements();
276 
277  /// Trace file
278  ofstream Trace_file;
279 
280  /// Pointers to node whose position we're tracing
281  Node* Trace_node_pt;
282 
283 
284 #ifdef REFINE
285 
286  /// Pointer to solid mesh
287  ElasticRefineableRectangularQuadMesh<ELEMENT>* Solid_mesh_pt;
288 
289 #else
290 
291  /// Pointer to solid mesh
292  ElasticRectangularQuadMesh<ELEMENT>* Solid_mesh_pt;
293 
294 #endif
295 
296 
297  /// Pointers to meshes of traction elements
298  SolidMesh* Traction_mesh_pt;
299 
300  /// DocInfo object for output
301  DocInfo Doc_info;
302 
303 };
304 
305 
306 //===========start_of_constructor=======================================
307 /// Constructor:
308 //======================================================================
309 template<class ELEMENT>
311  const bool& use_fd)
312 {
313 
314  // Create the mesh
315 
316  // # of elements in x-direction
317  unsigned n_x=20;
318 
319  // # of elements in y-direction
320  unsigned n_y=2;
321 
322  // Domain length in x-direction
323  double l_x= Global_Physical_Variables::L;
324 
325  // Domain length in y-direction
326  double l_y=2.0*Global_Physical_Variables::H;
327 
328  // Shift mesh downwards so that centreline is at y=0:
329  Vector<double> origin(2);
330  origin[0]=0.0;
331  origin[1]=-0.5*l_y;
332 
333 #ifdef REFINE
334 
335  //Now create the mesh
336  solid_mesh_pt() = new ElasticRefineableRectangularQuadMesh<ELEMENT>(
337  n_x,n_y,l_x,l_y,origin);
338 
339  // Set error estimator
340  dynamic_cast<ElasticRefineableRectangularQuadMesh<ELEMENT>* >(
341  solid_mesh_pt())->spatial_error_estimator_pt()=new Z2ErrorEstimator;
342 
343 #else
344 
345  //Now create the mesh
346  solid_mesh_pt() = new ElasticRectangularQuadMesh<ELEMENT>(
347  n_x,n_y,l_x,l_y,origin);
348 
349 #endif
350 
351 
352  //Assign the physical properties to the elements before any refinement
353  //Loop over the elements in the main mesh
354  unsigned n_element =solid_mesh_pt()->nelement();
355  for(unsigned i=0;i<n_element;i++)
356  {
357  //Cast to a solid element
358  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(solid_mesh_pt()->element_pt(i));
359 
360  // Set the constitutive law
361  el_pt->constitutive_law_pt() =
363 
364  //Set the body force
365  el_pt->body_force_fct_pt() = Global_Physical_Variables::gravity;
366 
367  // Get Jacobian by FD?
368  if(use_fd)
369  {
370  el_pt->enable_evaluate_jacobian_by_fd();
371  }
372  else
373  {
374  el_pt->disable_evaluate_jacobian_by_fd();
375  }
376 
377  // Is it incompressible
378  if (incompress)
379  {
380  PVDEquationsWithPressure<2>* test_pt =
381  dynamic_cast<PVDEquationsWithPressure<2>*>(
382  solid_mesh_pt()->element_pt(i));
383  if (test_pt!=0)
384  {
385  test_pt->set_incompressible();
386  }
387  }
388  }
389 
390 
391  // Choose a control node: The last node in the solid mesh
392  unsigned nnod=solid_mesh_pt()->nnode();
393  Trace_node_pt=solid_mesh_pt()->node_pt(nnod-1);
394 
395 #ifdef REFINE
396 
397  // Refine the mesh uniformly
398  dynamic_cast<ElasticRefineableRectangularQuadMesh<ELEMENT>* >(
399  solid_mesh_pt())->refine_uniformly();
400 
401 #endif
402 
403  // Construct the traction element mesh
404  Traction_mesh_pt=new SolidMesh;
405  create_traction_elements();
406 
407  // Pass pointer to traction function to the elements
408  // in the traction mesh
409  set_traction_pt();
410 
411  // Solid mesh is first sub-mesh
412  add_sub_mesh(solid_mesh_pt());
413 
414  // Add traction sub-mesh
415  add_sub_mesh(traction_mesh_pt());
416 
417  // Build combined "global" mesh
418  build_global_mesh();
419 
420  // Pin the left boundary (boundary 3) in both directions
421  unsigned n_side = mesh_pt()->nboundary_node(3);
422 
423  //Loop over the nodes
424  for(unsigned i=0;i<n_side;i++)
425  {
426  solid_mesh_pt()->boundary_node_pt(3,i)->pin_position(0);
427  solid_mesh_pt()->boundary_node_pt(3,i)->pin_position(1);
428  }
429 
430 
431  // Pin the redundant solid pressures (if any)
432  PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
433  solid_mesh_pt()->element_pt());
434 
435  //Attach the boundary conditions to the mesh
436  cout << assign_eqn_numbers() << std::endl;
437 
438 
439 } //end of constructor
440 
441 
442 //=====================start_of_actions_before_adapt======================
443 /// Actions before adapt: Wipe the mesh of traction elements
444 //========================================================================
445 template<class ELEMENT>
447 {
448  // Kill the traction elements and wipe surface mesh
449  delete_traction_elements();
450 
451  // Rebuild the Problem's global mesh from its various sub-meshes
452  rebuild_global_mesh();
453 
454 }// end of actions_before_adapt
455 
456 
457 
458 //=====================start_of_actions_after_adapt=======================
459 /// Actions after adapt: Rebuild the mesh of traction elements
460 //========================================================================
461 template<class ELEMENT>
463 {
464  // Create traction elements from all elements that are
465  // adjacent to boundary 2 and add them to surface meshes
466  create_traction_elements();
467 
468  // Rebuild the Problem's global mesh from its various sub-meshes
469  rebuild_global_mesh();
470 
471  // Pin the redundant solid pressures (if any)
472  PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
473  solid_mesh_pt()->element_pt());
474 
475  // Set pointer to prescribed traction function for traction elements
476  set_traction_pt();
477 
478 }// end of actions_after_adapt
479 
480 
481 
482 //==================start_of_set_traction_pt==============================
483 /// Set pointer to traction function for the relevant
484 /// elements in the traction mesh
485 //========================================================================
486 template<class ELEMENT>
488 {
489  // Loop over the elements in the traction element mesh
490  // for elements on the top boundary (boundary 2)
491  unsigned n_element=traction_mesh_pt()->nelement();
492  for(unsigned i=0;i<n_element;i++)
493  {
494  //Cast to a solid traction element
495  SolidTractionElement<ELEMENT> *el_pt =
496  dynamic_cast<SolidTractionElement<ELEMENT>*>
497  (traction_mesh_pt()->element_pt(i));
498 
499  //Set the traction function
500  el_pt->traction_fct_pt() = Global_Physical_Variables::constant_pressure;
501  }
502 
503 }// end of set traction pt
504 
505 
506 
507 //============start_of_create_traction_elements==============================
508 /// Create traction elements
509 //=======================================================================
510 template<class ELEMENT>
512 {
513  // Traction elements are located on boundary 2:
514  unsigned b=2;
515 
516  // How many bulk elements are adjacent to boundary b?
517  unsigned n_element = solid_mesh_pt()->nboundary_element(b);
518 
519  // Loop over the bulk elements adjacent to boundary b?
520  for(unsigned e=0;e<n_element;e++)
521  {
522  // Get pointer to the bulk element that is adjacent to boundary b
523  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
524  solid_mesh_pt()->boundary_element_pt(b,e));
525 
526  //Find the index of the face of element e along boundary b
527  int face_index = solid_mesh_pt()->face_index_at_boundary(b,e);
528 
529  // Create new element and add to mesh
530  Traction_mesh_pt->add_element_pt(new SolidTractionElement<ELEMENT>
531  (bulk_elem_pt,face_index));
532  }
533 
534  // Pass the pointer to the traction function to the traction elements
535  set_traction_pt();
536 
537 } // end of create_traction_elements
538 
539 
540 
541 
542 //============start_of_delete_traction_elements==============================
543 /// Delete traction elements and wipe the traction meshes
544 //=======================================================================
545 template<class ELEMENT>
547 {
548  // How many surface elements are in the surface mesh
549  unsigned n_element = Traction_mesh_pt->nelement();
550 
551  // Loop over the surface elements
552  for(unsigned e=0;e<n_element;e++)
553  {
554  // Kill surface element
555  delete Traction_mesh_pt->element_pt(e);
556  }
557 
558  // Wipe the mesh
559  Traction_mesh_pt->flush_element_and_node_storage();
560 
561 } // end of delete_traction_elements
562 
563 
564 
565 //==============start_doc===========================================
566 /// Doc the solution
567 //==================================================================
568 template<class ELEMENT>
570 {
571 
572  ofstream some_file;
573  char filename[100];
574 
575  // Number of plot points
576  unsigned n_plot = 5;
577 
578  // Output shape of and stress in deformed body
579  //--------------------------------------------
580  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
581  Doc_info.number());
582  some_file.open(filename);
583  solid_mesh_pt()->output(some_file,n_plot);
584  some_file.close();
585 
586 
587  // Output St. Venant solution
588  //---------------------------
589  sprintf(filename,"%s/exact_soln%i.dat",Doc_info.directory().c_str(),
590  Doc_info.number());
591  some_file.open(filename);
592 
593  // Element dimension
594  unsigned el_dim=2;
595 
596  Vector<double> s(el_dim);
597  Vector<double> x(el_dim);
598  Vector<double> xi(el_dim);
599  DenseMatrix<double> sigma(el_dim,el_dim);
600 
601  // Constants for exact (St. Venant) solution
602  double a=-1.0/4.0*Global_Physical_Variables::P;
604  double c=1.0/8.0*Global_Physical_Variables::P/
607 
608  // Loop over all elements to plot exact solution for stresses
609  unsigned nel=solid_mesh_pt()->nelement();
610  for (unsigned e=0;e<nel;e++)
611  {
612  // Get pointer to element
613  ELEMENT* el_pt=dynamic_cast<ELEMENT*>(
614  solid_mesh_pt()->element_pt(e));
615 
616  //Tecplot header info
617  some_file << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
618 
619  //Loop over plot points
620  for(unsigned l2=0;l2<n_plot;l2++)
621  {
622  s[1] = -1.0 + l2*2.0/(n_plot-1);
623  for(unsigned l1=0;l1<n_plot;l1++)
624  {
625  s[0] = -1.0 + l1*2.0/(n_plot-1);
626 
627  // Get Eulerian and Lagrangian coordinates
628  el_pt->interpolated_x(s,x);
629  el_pt->interpolated_xi(s,xi);
630 
631  //Output the x,y,..
632  for(unsigned i=0;i<el_dim;i++)
633  {some_file << x[i] << " ";}
634 
635  // Change orientation of coordinate system relative
636  // to solution in lecture notes
637  double xx=Global_Physical_Variables::L-xi[0];
638  double yy=xi[1];
639 
640  // Approximate analytical (St. Venant) solution
641  sigma(0,0)=c*(6.0*xx*xx*yy-4.0*yy*yy*yy)+
642  6.0*d*yy;
643  sigma(1,1)=2.0*(a+b*yy+c*yy*yy*yy);
644  sigma(1,0)=2.0*(b*xx+3.0*c*xx*yy*yy);
645  sigma(0,1)=sigma(1,0);
646 
647  // Output stress
648  some_file << sigma(0,0) << " "
649  << sigma(1,0) << " "
650  << sigma(1,1) << " "
651  << std::endl;
652  }
653  }
654  }
655  some_file.close();
656 
657  // Write trace file: Load/displacement characteristics
658  Trace_file << Global_Physical_Variables::P << " "
659  << Trace_node_pt->x(0) << " "
660  << Trace_node_pt->x(1) << " "
661  << std::endl;
662 
663  // Increment label for output files
664  Doc_info.number()++;
665 
666 } //end doc
667 
668 
669 
670 
671 
672 //==============start_run_it========================================
673 /// Run it
674 //==================================================================
675 template<class ELEMENT>
676 void CantileverProblem<ELEMENT>::run_it(const unsigned& i_case)
677 {
678 
679 #ifdef TIME_SOLID_JAC
680  PVDEquationsBase<2>::Solid_timer.reset();
681 #endif
682 
683  // Set output directory
684  char dirname[100];
685 
686 #ifdef REFINE
687  sprintf(dirname,"RESLT_refine%i",i_case);
688 #else
689  sprintf(dirname,"RESLT_norefine%i",i_case);
690 #endif
691 
692  Doc_info.set_directory(dirname);
693 
694  // Open trace file
695  char filename[100];
696  sprintf(filename,"%s/trace.dat",Doc_info.directory().c_str());
697  Trace_file.open(filename);
698 
699 
700  // Doc solution
701  doc_solution();
702 
703  // Initial values for parameter values
706 
707  //Parameter incrementation
708  unsigned nstep=1;
709  double p_increment=1.0e-5;
710  for(unsigned i=0;i<nstep;i++)
711  {
712  // Increment pressure load
713  Global_Physical_Variables::P+=p_increment;
714 
715  // Solve the problem with Newton's method, allowing
716  // up to max_adapt mesh adaptations after every solve.
717 
718 #ifdef REFINE
719 
720  // Max. number of adaptations per solve
721  unsigned max_adapt=1;
722 
723  newton_solve(max_adapt);
724 
725 #else
726 
727  newton_solve();
728 
729 #endif
730 
731  // Doc solution
732  doc_solution();
733 
734  }
735 
736 #ifdef TIME_SOLID_JAC
737 
738  std::cout << "Total fill_in... : "
739  << PVDEquationsBase<2>::Solid_timer.cumulative_time(0)
740  << std::endl;
741 
742  std::cout << "Total d_stress_dG: "
743  << PVDEquationsBase<2>::Solid_timer.cumulative_time(1)
744  << std::endl;
745 
746  std::cout << "Total Jac : "
747  << PVDEquationsBase<2>::Solid_timer.cumulative_time(2)
748  << std::endl;
749 
750  PVDEquationsBase<2>::Solid_timer.reset();
751 
752 #endif
753 
754 }
755 
756 
757 //=======start_of_main==================================================
758 /// Driver for cantilever beam loaded by surface traction and/or
759 /// gravity
760 //======================================================================
761 int main()
762 {
763 
764  bool use_fd=false;
765 
766  // Is the material incomressible
767  bool incompress=false;
768 
769  // Number of cases per implementation
770  unsigned ncase=5;
771 
772 
773  // Loop over fd and analytical implementation
774  for (unsigned i=0;i<2;i++)
775  {
776 
777  // Generalised Hookean constitutive equations
778  //-------------------------------------------
779  {
781  new GeneralisedHookean(&Global_Physical_Variables::Nu,
783 
784  incompress=false;
785 
786 #ifdef REFINE
787  {
788  //Set up the problem with pure displacement based elements
790  problem(incompress,use_fd);
791 
792  problem.run_it(0+i*ncase);
793  }
794 #else
795  {
796  //Set up the problem with pure displacement based elements
798  problem(incompress,use_fd);
799 
800  problem.run_it(0+i*ncase);
801  }
802 #endif
803 
804 
805 #ifdef REFINE
806  {
807  //Set up the problem with continous pressure/displacement
809  RefineableQPVDElementWithContinuousPressure<2> > >
810  problem(incompress,use_fd);
811 
812  problem.run_it(1+i*ncase);
813  }
814 #else
815  {
816  //Set up the problem with continous pressure/displacement
818  QPVDElementWithContinuousPressure<2> > >
819  problem(incompress,use_fd);
820 
821  problem.run_it(1+i*ncase);
822  }
823 #endif
824 
825 
826 #ifdef REFINE
827  {
828  //Set up the problem with discontinous pressure/displacement
830  RefineableQPVDElementWithPressure<2> > > problem(incompress,use_fd);
831 
832  problem.run_it(2+i*ncase);
833  }
834 #else
835  {
836  //Set up the problem with discontinous pressure/displacement
838  QPVDElementWithPressure<2> > > problem(incompress,use_fd);
839 
840  problem.run_it(2+i*ncase);
841  }
842 #endif
843 
846  }
847 
848 
849 
850  // Incompressible Mooney Rivlin
851  //-----------------------------
852  {
854  new MooneyRivlin(&Global_Physical_Variables::C1,
856 
857  // Define a constitutive law (based on strain energy function)
859  new IsotropicStrainEnergyFunctionConstitutiveLaw(
861 
862  incompress=true;
863 
864 
865 #ifdef REFINE
866  {
867  //Set up the problem with continous pressure/displacement
869  RefineableQPVDElementWithContinuousPressure<2> > >
870  problem(incompress,use_fd);
871 
872  problem.run_it(3+i*ncase);
873  }
874 #else
875  {
876  //Set up the problem with continous pressure/displacement
878  QPVDElementWithContinuousPressure<2> > >
879  problem(incompress,use_fd);
880 
881  problem.run_it(3+i*ncase);
882  }
883 #endif
884 
885 
886 
887 #ifdef REFINE
888  {
889  //Set up the problem with discontinous pressure/displacement
891  RefineableQPVDElementWithPressure<2> > >
892  problem(incompress,use_fd);
893 
894  problem.run_it(4+i*ncase);
895  }
896 #else
897  {
898  //Set up the problem with discontinous pressure/displacement
900  QPVDElementWithPressure<2> > >
901  problem(incompress,use_fd);
902 
903  problem.run_it(4+i*ncase);
904  }
905 #endif
906 
909 
912  }
913 
914 
915  use_fd=true;
916  std::cout << "\n\n\n CR Total fill_in... : bla \n\n\n " << std::endl;
917 
918  }
919 
920 } //end of main
921 
922 
923 
924 
925 
926 
927 
928 
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.
ElasticRectangularQuadMesh< ELEMENT > * Solid_mesh_pt
Pointer to solid mesh.
ConstitutiveLaw * Constitutive_law_pt
Pointer to constitutive law.
double P
Uniform pressure.
int main()
MySolidElement()
Constructor: Call constructor of underlying element.
void actions_before_newton_solve()
Update function (empty)
double L
Length of beam.
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...
void run_it(const unsigned &i_case)
Run the job – doc in RESLTi_case.
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)
void doc_solution()
Doc the solution.
FaceGeometry()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
double C1
"Mooney Rivlin" coefficient for generalised Mooney Rivlin law
StrainEnergyFunction * Strain_energy_function_pt
Pointer to strain energy function.
CantileverProblem()
Constructor:
double C2
"Mooney Rivlin" coefficient for generalised Mooney Rivlin law
ElasticRectangularQuadMesh< ELEMENT > *& solid_mesh_pt()
Access function for the solid mesh.
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.
double E
Elastic modulus.