Toggle navigation
Documentation
Big picture
The finite element method
The data structure
Not-so-quick guide
Optimisation
Order of action functions
Example codes and tutorials
List of example codes and tutorials
Meshing
Solvers
MPI parallel processing
Post-processing/visualisation
Other
Change log
Creating documentation
Coding conventions
Index
FAQ
About
People
Contact/Get involved
Publications
Acknowledgements
Copyright
Picture show
Go
demo_drivers
solid
disk_compression
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
41
using namespace
std;
42
43
using namespace
oomph;
44
45
46
//============namespace_for_problem_parameters=====================
47
/// Global variables
48
//=================================================================
49
namespace
Global_Physical_Variables
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
//=====================================================================
111
template
<
class
ELEMENT>
112
class
ElasticRefineableQuarterCircleSectorMesh
:
113
public
virtual
RefineableQuarterCircleSectorMesh<ELEMENT>,
114
public
virtual
SolidMesh
115
{
116
117
118
public
:
119
120
/// Constructor: Build mesh and copy Eulerian coords to Lagrangian
121
/// ones so that the initial configuration is the stress-free one.
122
ElasticRefineableQuarterCircleSectorMesh<ELEMENT>
(
GeomObject
*
wall_pt
,
123
const
double
&
xi_lo
,
124
const
double
&
fract_mid
,
125
const
double
&
xi_hi
,
126
TimeStepper
*
time_stepper_pt
=
127
&
Mesh::Default_TimeStepper
) :
128
RefineableQuarterCircleSectorMesh<ELEMENT>
(
wall_pt
,
xi_lo
,
fract_mid
,
xi_hi
,
129
time_stepper_pt
)
130
{
131
/// Make the current configuration the undeformed one by
132
/// setting the nodal Lagrangian coordinates to their current
133
/// Eulerian ones
134
set_lagrangian_nodal_coordinates
();
135
}
136
137
138
/// Function to create mesh made of traction elements
139
void
make_traction_element_mesh
(SolidMesh*&
traction_mesh_pt
)
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:
151
FiniteElement
*
fe_pt
= this->
boundary_element_pt
(
b
,
e
);
152
153
// Find the index of the face of element e along boundary b
154
int
face_index
= this->
face_index_at_boundary
(
b
,
e
);
155
156
// Create new element
157
traction_mesh_pt
->add_element_pt(
new
SolidTractionElement<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
//======================================================================
169
template
<
class
ELEMENT>
170
class
StaticDiskCompressionProblem
:
public
Problem
171
{
172
173
public
:
174
175
/// Constructor:
176
StaticDiskCompressionProblem
();
177
178
/// Run simulation: Pass case number to label output files
179
void
parameter_study
(
const
unsigned
&
case_number
);
180
181
/// Doc the solution
182
void
doc_solution
(
DocInfo
&
doc_info
);
183
184
/// Update function (empty)
185
void
actions_after_newton_solve
() {}
186
187
/// Update function (empty)
188
void
actions_before_newton_solve
() {}
189
190
private
:
191
192
/// Trace file
193
ofstream
Trace_file
;
194
195
/// Vector of pointers to nodes whose position we're tracing
196
Vector<Node*>
Trace_node_pt
;
197
198
/// Pointer to solid mesh
199
ElasticRefineableQuarterCircleSectorMesh<ELEMENT>
*
Solid_mesh_pt
;
200
201
/// Pointer to mesh of traction elements
202
SolidMesh*
Traction_mesh_pt
;
203
204
};
205
206
//======================================================================
207
/// Constructor:
208
//======================================================================
209
template
<
class
ELEMENT>
210
StaticDiskCompressionProblem<ELEMENT>::StaticDiskCompressionProblem
()
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
226
Solid_mesh_pt =
new
ElasticRefineableQuarterCircleSectorMesh<ELEMENT>
(
227
curved_boundary_pt
,
xi_lo
,
fract_mid
,
xi_hi
);
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
249
build_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)
263
PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures
(
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() =
275
Global_Physical_Variables::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
286
SolidTractionElement<ELEMENT>
*
el_pt
=
287
dynamic_cast<
SolidTractionElement<ELEMENT>
*
>
288
(Traction_mesh_pt->element_pt(
i
));
289
290
//Set the traction function
291
el_pt
->traction_fct_pt() =
Global_Physical_Variables::constant_pressure
;
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
//==================================================================
302
template
<
class
ELEMENT>
303
void
StaticDiskCompressionProblem<ELEMENT>::doc_solution
(
DocInfo
&
doc_info
)
304
{
305
306
ofstream
some_file
;
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
328
double
nu
=
Global_Physical_Variables::Nu
;
329
double
exact_r
=
sqrt
(
Global_Physical_Variables::Uniform_gamma
)*
330
(1.0-
Global_Physical_Variables::P
/
Global_Physical_Variables::E
331
*((1.0+
nu
)*(1.0-2.0*
nu
)));
332
333
334
// Write trace file: Problem parameters
335
Trace_file <<
Global_Physical_Variables::P
<<
" "
336
<<
Global_Physical_Variables::Uniform_gamma
<<
" "
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
//==================================================================
355
template
<
class
ELEMENT>
356
void
StaticDiskCompressionProblem<ELEMENT>::parameter_study
(
357
const
unsigned
&
case_number
)
358
{
359
// Output
360
DocInfo
doc_info
;
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)
382
if
(
CommandLineArgs::Argc
!=1)
383
{
384
nstep
=3;
385
}
386
387
// Offset external pressure so that the computation sweeps
388
// over a range of positive and negative pressures
389
Global_Physical_Variables::P
=-
delta_p
*
double
(
nstep
-1)*0.5;
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
402
Global_Physical_Variables::P
+=
delta_p
;
403
}
404
405
}
// end of parameter study
406
407
408
//=====start_of_main====================================================
409
/// Driver code for disk-compression
410
//======================================================================
411
int
main
(
int
argc
,
char
*
argv
[])
412
{
413
414
// Store command line arguments
415
CommandLineArgs::setup
(
argc
,
argv
);
416
417
// Define a strain energy function: Generalised Mooney Rivlin
418
Global_Physical_Variables::Strain_energy_function_pt
=
419
new
GeneralisedMooneyRivlin
(&
Global_Physical_Variables::Nu
,
420
&
Global_Physical_Variables::C1
,
421
&
Global_Physical_Variables::E
);
422
423
// Define a constitutive law (based on strain energy function)
424
Global_Physical_Variables::Constitutive_law_pt
=
425
new
IsotropicStrainEnergyFunctionConstitutiveLaw
(
426
Global_Physical_Variables::Strain_energy_function_pt
);
427
428
// Case 0: No pressure, generalised Mooney Rivlin
429
//-----------------------------------------------
430
{
431
//Set up the problem
432
StaticDiskCompressionProblem<RefineableQPVDElement<2,3>
>
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
446
StaticDiskCompressionProblem
<
447
RefineableQPVDElementWithContinuousPressure<2>
>
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
463
StaticDiskCompressionProblem<RefineableQPVDElementWithPressure<2>
>
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
475
delete
Global_Physical_Variables::Constitutive_law_pt
;
476
477
// Create oomph-lib's generalised Hooke's law constitutive equation
478
Global_Physical_Variables::Constitutive_law_pt
=
479
new
GeneralisedHookean
(&
Global_Physical_Variables::Nu
,
480
&
Global_Physical_Variables::E
);
481
482
// Case 3: No pressure, generalised Hooke's law
483
//----------------------------------------------
484
{
485
//Set up the problem
486
StaticDiskCompressionProblem<RefineableQPVDElement<2,3>
>
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
500
StaticDiskCompressionProblem
<
501
RefineableQPVDElementWithContinuousPressure<2>
>
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
517
StaticDiskCompressionProblem<RefineableQPVDElementWithPressure<2>
>
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
527
delete
Global_Physical_Variables::Constitutive_law_pt
;
528
Global_Physical_Variables::Constitutive_law_pt
=0;
529
530
}
// end of main
531
532
533
534
535
536
537
538
ElasticRefineableQuarterCircleSectorMesh
Elastic quarter circle sector mesh with functionality to attach traction elements to the curved surfa...
Definition
disk_compression.cc:115
ElasticRefineableQuarterCircleSectorMesh::ElasticRefineableQuarterCircleSectorMesh
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...
Definition
disk_compression.cc:122
ElasticRefineableQuarterCircleSectorMesh::make_traction_element_mesh
void make_traction_element_mesh(SolidMesh *&traction_mesh_pt)
Function to create mesh made of traction elements.
Definition
disk_compression.cc:139
StaticDiskCompressionProblem
Uniform compression of a circular disk in a state of plane strain, subject to uniform growth.
Definition
disk_compression.cc:171
StaticDiskCompressionProblem::actions_before_newton_solve
void actions_before_newton_solve()
Update function (empty)
Definition
disk_compression.cc:188
StaticDiskCompressionProblem::StaticDiskCompressionProblem
StaticDiskCompressionProblem()
Constructor:
Definition
disk_compression.cc:210
StaticDiskCompressionProblem::actions_after_newton_solve
void actions_after_newton_solve()
Update function (empty)
Definition
disk_compression.cc:185
StaticDiskCompressionProblem::Traction_mesh_pt
SolidMesh * Traction_mesh_pt
Pointer to mesh of traction elements.
Definition
disk_compression.cc:202
StaticDiskCompressionProblem::Solid_mesh_pt
ElasticRefineableQuarterCircleSectorMesh< ELEMENT > * Solid_mesh_pt
Pointer to solid mesh.
Definition
disk_compression.cc:199
StaticDiskCompressionProblem::Trace_file
ofstream Trace_file
Trace file.
Definition
disk_compression.cc:193
StaticDiskCompressionProblem::Trace_node_pt
Vector< Node * > Trace_node_pt
Vector of pointers to nodes whose position we're tracing.
Definition
disk_compression.cc:196
StaticDiskCompressionProblem::doc_solution
void doc_solution(DocInfo &doc_info)
Doc the solution.
Definition
disk_compression.cc:303
StaticDiskCompressionProblem::parameter_study
void parameter_study(const unsigned &case_number)
Run simulation: Pass case number to label output files.
Definition
disk_compression.cc:356
main
int main(int argc, char *argv[])
Driver code for disk-compression.
Definition
disk_compression.cc:411
Global_Physical_Variables
Global variables.
Definition
disk_compression.cc:50
Global_Physical_Variables::E
double E
Elastic modulus.
Definition
disk_compression.cc:58
Global_Physical_Variables::constant_pressure
void constant_pressure(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Constant pressure load.
Definition
disk_compression.cc:70
Global_Physical_Variables::P
double P
Uniform pressure.
Definition
disk_compression.cc:67
Global_Physical_Variables::Constitutive_law_pt
ConstitutiveLaw * Constitutive_law_pt
Pointer to constitutive law.
Definition
disk_compression.cc:55
Global_Physical_Variables::Nu
double Nu
Poisson's ratio.
Definition
disk_compression.cc:61
Global_Physical_Variables::growth_function
void growth_function(const Vector< double > &xi, double &gamma)
Growth function.
Definition
disk_compression.cc:85
Global_Physical_Variables::Strain_energy_function_pt
StrainEnergyFunction * Strain_energy_function_pt
Pointer to strain energy function.
Definition
disk_compression.cc:52
Global_Physical_Variables::C1
double C1
"Mooney Rivlin" coefficient for generalised Mooney Rivlin law
Definition
disk_compression.cc:64
Global_Physical_Variables::Uniform_gamma
double Uniform_gamma
Uniform volumetric expansion.
Definition
disk_compression.cc:82