steady_ring.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-2022 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 ring problem with/without
27//diplacement control
28
29//Generic oomph-lib sources
30#include "generic.h"
31
32// Beam sources
33#include "beam.h"
34
35// The mesh
36#include "meshes/one_d_lagrangian_mesh.h"
37
38using namespace std;
39
40using namespace oomph;
41
42
43//=======start_of_namespace=========================
44/// Namespace for physical parameters
45//==================================================
47{
48
49 /// Nondim thickness
50 double H=0.05;
51
52 /// Prescribed position (only used for displacement control)
53 double Xprescr = 1.0;
54
55 /// Perturbation pressure
56 double Pcos=0.0;
57
58 /// Pointer to pressure load (stored in Data so it can
59 /// become an unknown in the problem when displacement control is used
61
62 /// Load function: Constant external pressure with cos variation to
63 /// induce buckling in n=2 mode
64 void press_load(const Vector<double>& xi,
65 const Vector<double> &x,
66 const Vector<double>& N,
67 Vector<double>& load)
68 {
69 for(unsigned i=0;i<2;i++)
70 {
71 load[i] = (Pext_data_pt->value(0)-Pcos*cos(2.0*xi[0]))*N[i];
72 }
73 }
74
75 /// Return a reference to the external pressure
76 /// load on the elastic ring.
77 /// A reference is obtained by de-referencing the pointer to the
78 /// data value that contains the external load
80 {return *Pext_data_pt->value_pt(0);}
81
82
83} // end of namespace
84
85
86
87/// //////////////////////////////////////////////////////////////////
88/// //////////////////////////////////////////////////////////////////
89/// //////////////////////////////////////////////////////////////////
90
91
92
93//========start_of_problem_class========================================
94/// Ring problem
95//======================================================================
96template<class ELEMENT>
97class ElasticRingProblem : public Problem
98{
99
100public:
101
102 /// Constructor: Number of elements and flags for displ control
103 /// and displacement control with existing data respectively.
104 ElasticRingProblem(const unsigned &n_element,
105 bool& displ_control,
106 bool& load_data_already_exists);
107
108 /// Access function for the specific mesh
109 OneDLagrangianMesh<ELEMENT>* mesh_pt()
110 {
111 return dynamic_cast<OneDLagrangianMesh<ELEMENT>*>(Problem::mesh_pt());
112 }
113
114 /// Update function is empty
116
117 /// Update function is empty
119
120 /// Doc solution
121 void doc_solution(DocInfo& doc_info, ofstream& trace_file);
122
123 /// Perform the parameter study
124 void parameter_study(DocInfo& doc_info);
125
126private:
127
128 /// Use displacement control?
130
131 /// Pointer to geometric object that represents the undeformed shape
132 GeomObject* Undef_geom_pt;
133
134 /// Number of elements in the beam mesh
136
137}; // end of problem class
138
139
140
141
142
143
144
145//======start_of_constructor============================================
146/// Constructor for elastic ring problem
147//======================================================================
148template<class ELEMENT>
150(const unsigned& n_element, bool& displ_control, bool& load_data_already_exists) :
151 Displ_control(displ_control),Nbeam_element(n_element)
152{
153
154 // Undeformed beam is an elliptical ring
155 Undef_geom_pt=new Ellipse(1.0,1.0);
156
157 // Length of the doamin (in terms of the Lagrangian coordinates)
158 double length=2.0*atan(1.0);
159
160 //Now create the (Lagrangian!) mesh
161 Problem::mesh_pt() =
162 new OneDLagrangianMesh<ELEMENT>(n_element,length,Undef_geom_pt);
163
164 // Boundary condition:
165
166 // Bottom:
167 unsigned ibound=0;
168 // No vertical displacement
169 mesh_pt()->boundary_node_pt(ibound,0)->pin_position(1);
170 // Infinite slope: Pin type 1 (slope) dof for displacement direction 0
171 mesh_pt()->boundary_node_pt(ibound,0)->pin_position(1,0);
172
173 // Top:
174 ibound=1;
175 // No horizontal displacement
176 mesh_pt()->boundary_node_pt(ibound,0)->pin_position(0);
177 // Zero slope: Pin type 1 (slope) dof for displacement direction 1
178 mesh_pt()->boundary_node_pt(ibound,0)->pin_position(1,1);
179
180
181 // Normal load incrementation
182 //===========================
183 if (!Displ_control)
184 {
185 // Create Data object whose one-and-only value contains the
186 // (in principle) adjustable load
188
189 //Pin the external pressure because it isn't actually adjustable.
191 }
192 // Displacement control
193 //=====================
194 else
195 {
196 // Choose element in which displacement control is applied: the last one
197 SolidFiniteElement* controlled_element_pt=
198 dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(Nbeam_element-1));
199
200 // Fix the displacement in the vertical (1) direction...
201 unsigned controlled_direction=1;
202
203 //... at right end of the control element
204 Vector<double> s_displ_control(1);
205 s_displ_control[0]=1.0;
206
207 // Pointer to displacement control element
208 DisplacementControlElement* displ_control_el_pt;
209
210 // Displacement control without previously existing load Data
211 //-----------------------------------------------------------
212 if (!load_data_already_exists)
213 {
214 // Build displacement control element
215 displ_control_el_pt=
216 new DisplacementControlElement(controlled_element_pt,
217 s_displ_control,
218 controlled_direction,
220
221 // The constructor of the DisplacementControlElement has created
222 // a new Data object whose one-and-only value contains the
223 // adjustable load: Use this Data object in the load function:
224 Global_Physical_Variables::Pext_data_pt=displ_control_el_pt->
225 displacement_control_load_pt();
226 }
227 // Demonstrate use of displacement control with some existing data
228 //----------------------------------------------------------------
229 else
230 {
231 // Create Data object whose one-and-only value contains the
232 // adjustable load
234
235 // Currently, nobody's "in charge of" this Data so it won't
236 // get included in any equation numbering schemes etc.
237 // --> declare it to be "global Data" for the Problem
238 // so the Problem is in charge and will perform such tasks.
240
241 // Build displacement control element and pass pointer to the
242 // already existing adjustable load Data.
243 displ_control_el_pt=
244 new DisplacementControlElement(controlled_element_pt,
245 s_displ_control,
246 controlled_direction,
249 }
250
251 // Add the displacement-control element to the mesh
252 mesh_pt()->add_element_pt(displ_control_el_pt);
253 }
254
255
256
257
258 //Loop over the elements to set physical parameters etc.
259 for(unsigned i=0;i<Nbeam_element;i++)
260 {
261 //Cast to proper element type
262 ELEMENT *elem_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
263
264 // Set wall thickness
265 elem_pt->h_pt() = &Global_Physical_Variables::H;
266
267 // Function that specifies load Vector
268 elem_pt->load_vector_fct_pt() = &Global_Physical_Variables::press_load;
269
270 //Assign the undeformed beam shape
271 elem_pt->undeformed_beam_pt() = Undef_geom_pt;
272
273 // Displacement control? If so, the load on *all* elements
274 // is affected by an unknown -- the external pressure, stored
275 // as the one-and-only value in a Data object: Add it to the
276 // elements' external Data.
277 if (Displ_control)
278 {
279 //The external pressure is external data for all elements
280 elem_pt->add_external_data(Global_Physical_Variables::Pext_data_pt);
281 }
282 }
283
284 // Do equation numbering
285 cout << "# of dofs " << assign_eqn_numbers() << std::endl;
286
287} // end of constructor
288
289
290//=======start_of_doc=====================================================
291/// Document solution
292//========================================================================
293template<class ELEMENT>
295 ofstream& trace_file)
296{
297 ofstream some_file;
298 char filename[100];
299
300 // Number of plot points
301 unsigned npts=5;
302
303 // Output solution
304 sprintf(filename,"%s/ring%i.dat",doc_info.directory().c_str(),
305 doc_info.number());
306 some_file.open(filename);
307 mesh_pt()->output(some_file,npts);
308 some_file.close();
309
310 // Local coordinates of plot points at left and right end of domain
311 Vector<double> s_left(1);
312 s_left[0]=-1.0;
313 Vector<double> s_right(1);
314 s_right[0]=1.0;
315
316 // Write trace file: Pressure, two radii
317 trace_file
319 (pow(Global_Physical_Variables::H,3)/12.0)
320 << " "
321 << dynamic_cast<ELEMENT*>(mesh_pt()->
322 element_pt(0))->interpolated_x(s_left,0)
323 << " "
324 << dynamic_cast<ELEMENT*>
325 (mesh_pt()->element_pt(Nbeam_element-1))->interpolated_x(s_right,1)
326 << std::endl;
327
328
329} // end of doc
330
331
332
333//========start_of_run=====================================================
334/// Solver loop to perform parameter study
335//=========================================================================
336template<class ELEMENT>
338{
339
340 //Open a trace file
341 char filename[100];
342 sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
343 ofstream trace_file(filename);
344 trace_file << "VARIABLES=\"p_e_x_t\",\"R_1\",\"R_2\"" << std::endl;
345 trace_file << "ZONE" << std::endl;
346
347 //Output initial data
348 doc_solution(doc_info,trace_file);
349
350 // Number of steps
351 unsigned nstep= 11; //51;
352
353 // Increments
354 double displ_increment=1.0/double(nstep-1);
355 double p_buckl=3.00*pow(Global_Physical_Variables::H,3)/12.0;
356 double p_owc =5.22*pow(Global_Physical_Variables::H,3)/12.0;
357 double pext_increment=(p_owc-p_buckl)/double(nstep-1);
358
359 // Set initial values for parameters that are to be incremented
360 Global_Physical_Variables::Xprescr=1.0+displ_increment;
361 Global_Physical_Variables::Pext_data_pt->set_value(0,p_buckl-pext_increment);
362
363
364
365 // Without displacement control the Newton method converges very slowly
366 // as we approach the axisymmetric state: Allow more iterations before
367 // calling it a day...
368 if (Displ_control)
369 {
370 Max_newton_iterations=100;
371 }
372
373
374 // Downward loop over parameter incrementation with pcos>0
375 //--------------------------------------------------------
376
377 /// Perturbation pressure
379
380
381 // Downward loop over parameter incrementation
382 //---------------------------------------------
383 for(unsigned i=1;i<=nstep;i++)
384 {
385
386 // Displacement control?
387 if (!Displ_control)
388 {
389 // Increment pressure
391 }
392 else
393 {
394 // Increment control displacement
395 Global_Physical_Variables::Xprescr-=displ_increment;
396 }
397
398 // Solve
399 newton_solve();
400
401 // Doc solution
402 doc_info.number()++;
403 doc_solution(doc_info,trace_file);
404
405 } // end of downward loop
406
407 // Reset perturbation pressure
408 //----------------------------
410
411 // Set initial values for parameters that are to be incremented
412 Global_Physical_Variables::Xprescr-=displ_increment;
414
415 // Start new zone for tecplot
416 trace_file << "ZONE" << std::endl;
417
418 // Upward loop over parameter incrementation
419 //------------------------------------------
420 for(unsigned i=nstep;i<2*nstep;i++)
421 {
422
423 // Displacement control?
424 if (!Displ_control)
425 {
426 // Increment pressure
428 }
429 else
430 {
431 // Increment control displacement
432 Global_Physical_Variables::Xprescr+=displ_increment;
433 }
434
435 // Solve
436 newton_solve();
437
438 // Doc solution
439 doc_info.number()++;
440 doc_solution(doc_info,trace_file);
441
442 }
443
444} // end of run
445
446
447
448/// /////////////////////////////////////////////////////////////////////
449/// /////////////////////////////////////////////////////////////////////
450/// /////////////////////////////////////////////////////////////////////
451
452
453//=======start_of_main=================================================
454/// Driver for ring test problem
455//=====================================================================
456int main()
457{
458 // Number of elements
459 unsigned n_element = 13;
460
461 // Displacement control?
462 bool displ_control=true;
463
464 // Label for output
465 DocInfo doc_info;
466
467 // Demonstrate how to use displacement control with already existing load Data
468 //----------------------------------------------------------------------------
469 {
470 bool load_data_already_exists=true;
471
472 //Set up the problem
474 problem(n_element,displ_control,load_data_already_exists);
475
476 // Output directory
477 doc_info.set_directory("RESLT_global");
478
479 // Do static run
480 problem.parameter_study(doc_info);
481 }
482
483 // Demonstrate how to use displacement control without existing load Data
484 //-----------------------------------------------------------------------
485 {
486 bool load_data_already_exists=false;
487
488 //Set up the problem
490 problem(n_element,displ_control,load_data_already_exists);
491
492 // Output directory
493 doc_info.set_directory("RESLT_no_global");
494
495 // Reset counter
496 doc_info.number()=0;
497
498 // Do static run
499 problem.parameter_study(doc_info);
500 }
501
502} // end of main
503
504
505
506
507
508
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: steady_ring.cc:98
GeomObject * Undef_geom_pt
Pointer to geometric object that represents the undeformed shape.
Definition: steady_ring.cc:132
OneDLagrangianMesh< ELEMENT > * mesh_pt()
Access function for the specific mesh.
Definition: steady_ring.cc:109
ElasticRingProblem(const unsigned &n_element, bool &displ_control, bool &load_data_already_exists)
Constructor: Number of elements and flags for displ control and displacement control with existing da...
Definition: steady_ring.cc:150
void actions_after_newton_solve()
Update function is empty.
Definition: steady_ring.cc:115
void actions_before_newton_solve()
Update function is empty.
Definition: steady_ring.cc:118
bool Displ_control
Use displacement control?
Definition: steady_ring.cc:129
unsigned Nbeam_element
Number of elements in the beam mesh.
Definition: steady_ring.cc:135
void doc_solution(DocInfo &doc_info, ofstream &trace_file)
Doc solution.
Definition: steady_ring.cc:294
void parameter_study(DocInfo &doc_info)
Perform the parameter study.
Definition: steady_ring.cc:337
Namespace for physical parameters.
Definition: steady_ring.cc:47
double Xprescr
Prescribed position (only used for displacement control)
Definition: steady_ring.cc:53
void press_load(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Load function: Constant external pressure with cos variation to induce buckling in n=2 mode.
Definition: steady_ring.cc:64
Data * Pext_data_pt
Pointer to pressure load (stored in Data so it can become an unknown in the problem when displacement...
Definition: steady_ring.cc:60
double Pcos
Perturbation pressure.
Definition: steady_ring.cc:56
double & external_pressure()
Return a reference to the external pressure load on the elastic ring. A reference is obtained by de-r...
Definition: steady_ring.cc:79
double H
Nondim thickness.
Definition: steady_ring.cc:50
int main()
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: steady_ring.cc:456