disk_compression.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// Copyright (C) 2006-2025 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//Driver function for a simple test elasticity problem: the
27//compression of an axisymmetric disk. We also demonstrate how
28//to incorporate isotropic growth into the model and how to
29//switch between different constitutive equations.
30#include <iostream>
31#include <fstream>
32#include <cmath>
33
34//My own includes
35#include "generic.h"
36#include "solid.h"
37
38//Need to instantiate templated mesh
39#include "meshes/quarter_circle_sector_mesh.h"
40
41using namespace std;
42
43using namespace oomph;
44
45
46//============namespace_for_problem_parameters=====================
47/// Global variables
48//=================================================================
50{
51 /// Pointer to strain energy function
52 StrainEnergyFunction* Strain_energy_function_pt;
53
54 /// Pointer to constitutive law
55 ConstitutiveLaw* Constitutive_law_pt;
56
57 /// Elastic modulus
58 double E=1.0;
59
60 /// Poisson's ratio
61 double Nu=0.3;
62
63 /// "Mooney Rivlin" coefficient for generalised Mooney Rivlin law
64 double C1=1.3;
65
66 /// Uniform pressure
67 double P = 0.0;
68
69 /// Constant pressure load
70 void constant_pressure(const Vector<double> &xi,const Vector<double> &x,
71 const Vector<double> &n, Vector<double> &traction)
72 {
73 unsigned dim = traction.size();
74 for(unsigned i=0;i<dim;i++)
75 {
76 traction[i] = -P*n[i];
77 }
78 } // end of pressure load
79
80
81 /// Uniform volumetric expansion
82 double Uniform_gamma=1.1;
83
84 /// Growth function
85 void growth_function(const Vector<double>& xi, double& gamma)
86 {
87 gamma = Uniform_gamma;
88 }
89
90} // end namespace
91
92
93///////////////////////////////////////////////////////////////////////
94///////////////////////////////////////////////////////////////////////
95///////////////////////////////////////////////////////////////////////
96
97
98
99//==========start_mesh=================================================
100/// Elastic quarter circle sector mesh with functionality to
101/// attach traction elements to the curved surface. We "upgrade"
102/// the RefineableQuarterCircleSectorMesh to become an
103/// SolidMesh and equate the Eulerian and Lagrangian coordinates,
104/// thus making the domain represented by the mesh the stress-free
105/// configuration.
106/// \n\n
107/// The member function \c make_traction_element_mesh() creates
108/// a separate mesh of SolidTractionElements that are attached to the
109/// mesh's curved boundary (boundary 1).
110//=====================================================================
111template <class ELEMENT>
113 public virtual RefineableQuarterCircleSectorMesh<ELEMENT>,
114 public virtual SolidMesh
115{
116
117
118public:
119
120 /// Constructor: Build mesh and copy Eulerian coords to Lagrangian
121 /// ones so that the initial configuration is the stress-free one.
123 const double& xi_lo,
124 const double& fract_mid,
125 const double& xi_hi,
130 {
131 /// Make the current configuration the undeformed one by
132 /// setting the nodal Lagrangian coordinates to their current
133 /// Eulerian ones
135 }
136
137
138 /// Function to create mesh made of traction elements
140 {
141
142 // Make new mesh
143 traction_mesh_pt = new SolidMesh;
144
145 // Loop over all elements on boundary 1:
146 unsigned b=1;
147 unsigned n_element = this->nboundary_element(b);
148 for (unsigned e=0;e<n_element;e++)
149 {
150 // The element itself:
152
153 // Find the index of the face of element e along boundary b
155
156 // Create new element
158 (fe_pt,face_index));
159 }
160 }
161
162};
163
164
165//======================================================================
166/// Uniform compression of a circular disk in a state of plane strain,
167/// subject to uniform growth.
168//======================================================================
169template<class ELEMENT>
170class StaticDiskCompressionProblem : public Problem
171{
172
173public:
174
175 /// Constructor:
177
178 /// Run simulation: Pass case number to label output files
179 void parameter_study(const unsigned& case_number);
180
181 /// Doc the solution
183
184 /// Update function (empty)
186
187 /// Update function (empty)
189
190private:
191
192 /// Trace file
194
195 /// Vector of pointers to nodes whose position we're tracing
197
198 /// Pointer to solid mesh
200
201 /// Pointer to mesh of traction elements
203
204};
205
206//======================================================================
207/// Constructor:
208//======================================================================
209template<class ELEMENT>
211{
212 // Build the geometric object that describes the curvilinear
213 // boundary of the quarter circle domain
214 Ellipse* curved_boundary_pt = new Ellipse(1.0,1.0);
215
216 // The curved boundary of the mesh is defined by the geometric object
217 // What follows are the start and end coordinates on the geometric object:
218 double xi_lo=0.0;
219 double xi_hi=2.0*atan(1.0);
220
221 // Fraction along geometric object at which the radial dividing line
222 // is placed
223 double fract_mid=0.5;
224
225 //Now create the mesh using the geometric object
228
229 // Setup trace nodes as the nodes on boundary 1 (=curved boundary)
230 // in the original mesh.
231 unsigned n_boundary_node = Solid_mesh_pt->nboundary_node(1);
232 Trace_node_pt.resize(n_boundary_node);
233 for(unsigned j=0;j<n_boundary_node;j++)
234 {Trace_node_pt[j]=Solid_mesh_pt->boundary_node_pt(1,j);}
235
236 // Refine the mesh uniformly
237 Solid_mesh_pt->refine_uniformly();
238
239 // Now construct the traction element mesh
240 Solid_mesh_pt->make_traction_element_mesh(Traction_mesh_pt);
241
242 // Solid mesh is first sub-mesh
243 add_sub_mesh(Solid_mesh_pt);
244
245 // Traction mesh is second sub-mesh
246 add_sub_mesh(Traction_mesh_pt);
247
248 // Build combined "global" mesh
250
251
252 // Pin the left edge in the horizontal direction
253 unsigned n_side = mesh_pt()->nboundary_node(2);
254 for(unsigned i=0;i<n_side;i++)
255 {Solid_mesh_pt->boundary_node_pt(2,i)->pin_position(0);}
256
257 // Pin the bottom in the vertical direction
258 unsigned n_bottom = mesh_pt()->nboundary_node(0);
259 for(unsigned i=0;i<n_bottom;i++)
260 {Solid_mesh_pt->boundary_node_pt(0,i)->pin_position(1);}
261
262 // Pin the redundant solid pressures (if any)
264 Solid_mesh_pt->element_pt());
265
266 //Complete the build process for elements in "bulk" solid mesh
267 unsigned n_element =Solid_mesh_pt->nelement();
268 for(unsigned i=0;i<n_element;i++)
269 {
270 //Cast to a solid element
271 ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Solid_mesh_pt->element_pt(i));
272
273 // Set the constitutive law
274 el_pt->constitutive_law_pt() =
276
277 // Set the isotropic growth function pointer
278 el_pt->isotropic_growth_fct_pt()=Global_Physical_Variables::growth_function;
279 }
280
281 // Complete build process for SolidTractionElements
282 n_element=Traction_mesh_pt->nelement();
283 for(unsigned i=0;i<n_element;i++)
284 {
285 //Cast to a solid traction element
287 dynamic_cast<SolidTractionElement<ELEMENT>*>
288 (Traction_mesh_pt->element_pt(i));
289
290 //Set the traction function
292 }
293
294 //Set up equation numbering scheme
295 cout << "Number of equations: " << assign_eqn_numbers() << std::endl;
296}
297
298
299//==================================================================
300/// Doc the solution
301//==================================================================
302template<class ELEMENT>
304{
305
307 char filename[100];
308
309 // Number of plot points
310 unsigned npts = 5;
311
312 // Output shape of deformed body
313 snprintf(filename, sizeof(filename), "%s/soln%i.dat",doc_info.directory().c_str(),
314 doc_info.number());
315 some_file.open(filename);
316 Solid_mesh_pt->output(some_file,npts);
317 some_file.close();
318
319 //Find number of solid elements
320 unsigned nelement = Solid_mesh_pt->nelement();
321
322 // Work out volume
323 double volume = 0.0;
324 for(unsigned e=0;e<nelement;e++)
325 {volume+= Solid_mesh_pt->finite_element_pt(e)->size();}
326
327 // Exact outer radius for linear elasticity
331 *((1.0+nu)*(1.0-2.0*nu)));
332
333
334 // Write trace file: Problem parameters
335 Trace_file << Global_Physical_Variables::P << " "
337 << volume << " "
338 << exact_r << " ";
339
340 // Write radii of trace nodes
341 unsigned ntrace_node=Trace_node_pt.size();
342 for (unsigned j=0;j<ntrace_node;j++)
343 {
344 Trace_file << sqrt(pow(Trace_node_pt[j]->x(0),2)+
345 pow(Trace_node_pt[j]->x(1),2)) << " ";
346 }
347 Trace_file << std::endl;
348
349} // end of doc_solution
350
351
352//==================================================================
353/// Run the paramter study
354//==================================================================
355template<class ELEMENT>
357 const unsigned& case_number)
358{
359 // Output
361
362 char dirname[100];
363 snprintf(dirname, sizeof(dirname), "RESLT%i",case_number);
364
365 // Set output directory
366 doc_info.set_directory(dirname);
367
368 // Step number
369 doc_info.number()=0;
370
371 // Open trace file
372 char filename[100];
373 snprintf(filename, sizeof(filename), "%s/trace.dat",doc_info.directory().c_str());
374 Trace_file.open(filename);
375
376 //Parameter incrementation
377 double delta_p=0.0125;
378 unsigned nstep=21;
379
380 // Perform fewer steps if run as self-test (indicated by nonzero number
381 // of command line arguments)
383 {
384 nstep=3;
385 }
386
387 // Offset external pressure so that the computation sweeps
388 // over a range of positive and negative pressures
390
391 // Do the parameter study
392 for(unsigned i=0;i<nstep;i++)
393 {
394 //Solve the problem for current load
395 newton_solve();
396
397 // Doc solution
398 doc_solution(doc_info);
399 doc_info.number()++;
400
401 // Increment pressure load
403 }
404
405} // end of parameter study
406
407
408//=====start_of_main====================================================
409/// Driver code for disk-compression
410//======================================================================
411int main(int argc, char* argv[])
412{
413
414 // Store command line arguments
416
417 // Define a strain energy function: Generalised Mooney Rivlin
422
423 // Define a constitutive law (based on strain energy function)
427
428 // Case 0: No pressure, generalised Mooney Rivlin
429 //-----------------------------------------------
430 {
431 //Set up the problem
433
434 cout << "gen. M.R.: RefineableQPVDElement<2,3>" << std::endl;
435
436 //Run the simulation
437 problem.parameter_study(0);
438
439 } // done case 0
440
441
442 // Case 1: Continuous pressure formulation with generalised Mooney Rivlin
443 //------------------------------------------------------------------------
444 {
445 //Set up the problem
448
449 cout << "gen. M.R.: RefineableQPVDElementWithContinuousPressure<2> "
450 << std::endl;
451
452 //Run the simulation
453 problem.parameter_study(1);
454
455 } // done case 1
456
457
458
459 // Case 2: Discontinuous pressure formulation with generalised Mooney Rivlin
460 //--------------------------------------------------------------------------
461 {
462 //Set up the problem
464 problem;
465
466 cout << "gen. M.R.: RefineableQPVDElementWithPressure<2>" << std::endl;
467
468 //Run the simulation
469 problem.parameter_study(2);
470
471 } // done case 2
472
473
474 // Change the consitutive law: Delete the old one
476
477 // Create oomph-lib's generalised Hooke's law constitutive equation
481
482 // Case 3: No pressure, generalised Hooke's law
483 //----------------------------------------------
484 {
485 //Set up the problem
487
488 cout << "gen. Hooke: RefineableQPVDElement<2,3> " << std::endl;
489
490 //Run the simulation
491 problem.parameter_study(3);
492
493 } // done case 3
494
495 // Case 4: Continuous pressure formulation with generalised Hooke's law
496 //---------------------------------------------------------------------
497 {
498
499 //Set up the problem
502
503 cout << "gen. Hooke: RefineableQPVDElementWithContinuousPressure<2> "
504 << std::endl;
505
506 //Run the simulation
507 problem.parameter_study(4);
508
509 } // done case 4
510
511
512 // Case 5: Discontinous pressure formulation with generalised Hooke's law
513 //------------------------------------------------------------------------
514 {
515
516 //Set up the problem
518
519 cout << "gen. Hooke: RefineableQPVDElementWithPressure<2> " << std::endl;
520
521 //Run the simulation
522 problem.parameter_study(5);
523
524 } // done case 5
525
526 // Clean up
529
530} // end of main
531
532
533
534
535
536
537
538
Elastic quarter circle sector mesh with functionality to attach traction elements to the curved surfa...
ElasticRefineableQuarterCircleSectorMesh(GeomObject *wall_pt, const double &xi_lo, const double &fract_mid, const double &xi_hi, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Build mesh and copy Eulerian coords to Lagrangian ones so that the initial configuration...
void make_traction_element_mesh(SolidMesh *&traction_mesh_pt)
Function to create mesh made of traction elements.
Uniform compression of a circular disk in a state of plane strain, subject to uniform growth.
void actions_before_newton_solve()
Update function (empty)
void actions_after_newton_solve()
Update function (empty)
SolidMesh * Traction_mesh_pt
Pointer to mesh of traction elements.
ElasticRefineableQuarterCircleSectorMesh< ELEMENT > * Solid_mesh_pt
Pointer to solid mesh.
Vector< Node * > Trace_node_pt
Vector of pointers to nodes whose position we're tracing.
void doc_solution(DocInfo &doc_info)
Doc the solution.
void parameter_study(const unsigned &case_number)
Run simulation: Pass case number to label output files.
int main(int argc, char *argv[])
Driver code for disk-compression.
double E
Elastic modulus.
void constant_pressure(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Constant pressure load.
double P
Uniform pressure.
ConstitutiveLaw * Constitutive_law_pt
Pointer to constitutive law.
double Nu
Poisson's ratio.
void growth_function(const Vector< double > &xi, double &gamma)
Growth function.
StrainEnergyFunction * Strain_energy_function_pt
Pointer to strain energy function.
double C1
"Mooney Rivlin" coefficient for generalised Mooney Rivlin law
double Uniform_gamma
Uniform volumetric expansion.