pressure_loaded_cylinder.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
27
28// The oomphlib headers
29#include "generic.h"
30#include "time_harmonic_fourier_decomposed_linear_elasticity.h"
31
32// The mesh
33#include "meshes/rectangular_quadmesh.h"
34#include "meshes/triangle_mesh.h"
35
36using namespace std;
37
38using namespace oomph;
39
40
41//===start_of_namespace=================================================
42/// Namespace for global parameters
43//======================================================================
44namespace Global_Parameters
45{
46 /// Define Poisson's ratio Nu
47 std::complex<double> Nu(0.3,0.0);
48
49 /// Define the non-dimensional Young's modulus
50 std::complex<double> E(1.0,0.0);
51
52 /// Define Fourier wavenumber
53 int Fourier_wavenumber = 0;
54
55 /// Define the non-dimensional square angular frequency of
56 /// time-harmonic motion
57 std::complex<double> Omega_sq (10.0,0.0);
58
59 /// Length of domain in r direction
60 double Lr = 1.0;
61
62 /// Length of domain in z-direction
63 double Lz = 2.0;
64
65 // Set up min & max (r,z) coordinates
66 double rmin = 0.1;
67 double zmin = 0.3;
68 double rmax = rmin+Lr;
69 double zmax = zmin+Lz;
70
71 /// Define the imaginary unit
72 const std::complex<double> I(0.0,1.0);
73
74 // Pressure load
75 double P=1.0;
76
77 /// The traction function at r=rmin: (t_r, t_z, t_theta)
78 void boundary_traction(const Vector<double> &x,
79 const Vector<double> &n,
80 Vector<std::complex<double> > &result)
81 {
82 // Radial traction
83 result[0] = P;
84 // Axial traction
85 result[1] = 0.0;
86 // Azimuthal traction
87 result[2] = 0.0;
88 }
89
90
91} // end_of_namespace
92
93
94//===start_of_problem_class=============================================
95/// Class to validate time harmonic linear elasticity (Fourier
96/// decomposed)
97//======================================================================
98template<class ELEMENT>
100{
101public:
102
103 /// Constructor: Pass number of elements in r and z directions
104 /// and boundary locations
106 const unsigned &nr, const unsigned &nz,
107 const double &rmin, const double& rmax,
108 const double &zmin, const double& zmax);
109
110
111 /// Update before solve is empty
113
114 /// Update after solve is empty
116
117 /// Delete traction elements
119
120 /// Helper function to complete problem setup
122
123 /// Actions before adapt: Wipe the mesh of traction elements
125 {
126 // Kill the traction elements and wipe surface mesh
128
129 // Rebuild the Problem's global mesh from its various sub-meshes
130 rebuild_global_mesh();
131 }
132
133
134 /// Actions after adapt: Rebuild the mesh of traction elements
136 {
137 // Create traction elements from all elements that are
138 // adjacent to FSI boundaries and add them to surface meshes
140
141 // Rebuild the Problem's global mesh from its various sub-meshes
142 rebuild_global_mesh();
143
144 // Complete problem setup
146 }
147
148 /// Doc the solution
149 void doc_solution(DocInfo& doc_info);
150
151private:
152
153 /// Allocate traction elements on the bottom surface
155
156#ifdef ADAPTIVE
157
158 /// Pointer to the bulk mesh
159 RefineableTriangleMesh<ELEMENT>* Bulk_mesh_pt;
160
161#else
162
163 /// Pointer to the bulk mesh
164 Mesh* Bulk_mesh_pt;
165
166#endif
167
168 /// Pointer to the mesh of traction elements
169 Mesh* Surface_mesh_pt;
170}; // end_of_problem_class
171
172
173//===start_of_constructor=============================================
174/// Problem constructor: Pass number of elements in coordinate
175/// directions and size of domain.
176//====================================================================
177template<class ELEMENT>
180(const unsigned &nr, const unsigned &nz,
181 const double &rmin, const double& rmax,
182 const double &zmin, const double& zmax)
183{
184
185#ifdef ADAPTIVE
186
187 // The boundary is bounded by four distinct boundaries, each
188 // represented by its own polyline
189 Vector<TriangleMeshCurveSection*> boundary_polyline_pt(4);
190
191 // Vertex coordinates on boundary
192 Vector<Vector<double> > bound_coords(2);
193 bound_coords[0].resize(2);
194 bound_coords[1].resize(2);
195
196 // Horizontal bottom boundary
197 bound_coords[0][0]=rmin;
198 bound_coords[0][1]=zmin;
199 bound_coords[1][0]=rmax;
200 bound_coords[1][1]=zmin;
201
202 // Build the boundary polyline
203 unsigned boundary_id=0;
204 boundary_polyline_pt[0]=new TriangleMeshPolyLine(bound_coords,boundary_id);
205
206 // Vertical outer boundary
207 bound_coords[0][0]=rmax;
208 bound_coords[0][1]=zmin;
209 bound_coords[1][0]=rmax;
210 bound_coords[1][1]=zmax;
211
212 // Build the boundary polyline
213 boundary_id=1;
214 boundary_polyline_pt[1]=new TriangleMeshPolyLine(bound_coords,boundary_id);
215
216
217 // Horizontal top boundary
218 bound_coords[0][0]=rmax;
219 bound_coords[0][1]=zmax;
220 bound_coords[1][0]=rmin;
221 bound_coords[1][1]=zmax;
222
223 // Build the boundary polyline
224 boundary_id=2;
225 boundary_polyline_pt[2]=new TriangleMeshPolyLine(bound_coords,boundary_id);
226
227 // Vertical inner boundary
228 bound_coords[0][0]=rmin;
229 bound_coords[0][1]=zmax;
230 bound_coords[1][0]=rmin;
231 bound_coords[1][1]=zmin;
232
233 // Build the boundary polyline
234 boundary_id=3;
235 boundary_polyline_pt[3]=new TriangleMeshPolyLine(bound_coords,boundary_id);
236
237 // Pointer to the closed curve that defines the outer boundary
238 TriangleMeshClosedCurve* closed_curve_pt=
239 new TriangleMeshPolygon(boundary_polyline_pt);
240
241 // Use the TriangleMeshParameters object for helping on the manage of the
242 // TriangleMesh parameters
243 TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
244
245 // Specify the maximum area element
246 double uniform_element_area=0.2;
247 triangle_mesh_parameters.element_area() = uniform_element_area;
248
249 // Create the mesh
250 Bulk_mesh_pt=new RefineableTriangleMesh<ELEMENT>(triangle_mesh_parameters);
251
252 // Set error estimator
253 Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
254
255#else
256
257 //Now create the mesh
258 Bulk_mesh_pt = new RectangularQuadMesh<ELEMENT>(nr,nz,rmin,rmax,zmin,zmax);
259
260#endif
261
262 //Create the surface mesh of traction elements
263 Surface_mesh_pt=new Mesh;
264 assign_traction_elements();
265
266 // Complete problem setup
267 complete_problem_setup();
268
269 // Add the submeshes to the problem
270 add_sub_mesh(Bulk_mesh_pt);
271 add_sub_mesh(Surface_mesh_pt);
272
273 // Now build the global mesh
274 build_global_mesh();
275
276 // Assign equation numbers
277 cout << assign_eqn_numbers() << " equations assigned" << std::endl;
278
279} // end of constructor
280
281
282
283//===start_of_complete_problem_setup=================================
284/// Complete problem setup
285//===================================================================
286template<class ELEMENT>
289{
290 // Set the boundary conditions for this problem: All nodes are
291 // free by default -- just pin & set the ones that have Dirichlet
292 // conditions here
293
294 // Pin displacements everywhere apart from boundaries 1 and 3
295 //-----------------------------------------------------------
296 for (unsigned ibound=0;ibound<3;ibound=ibound+2)
297 {
298 unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
299 for (unsigned inod=0;inod<num_nod;inod++)
300 {
301 // Get pointer to node
302 Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(ibound,inod);
303
304 // Pinned in r, z and theta
305 nod_pt->pin(0);nod_pt->pin(1);nod_pt->pin(2);
306 nod_pt->pin(3);nod_pt->pin(4);nod_pt->pin(5);
307
308 // Set the displacements
309 nod_pt->set_value(0,0.0);
310 nod_pt->set_value(1,0.0);
311 nod_pt->set_value(2,0.0);
312 nod_pt->set_value(3,0.0);
313 nod_pt->set_value(4,0.0);
314 nod_pt->set_value(5,0.0);
315 }
316 }
317
318
319 // Complete the problem setup to make the elements fully functional
320
321 // Loop over the elements
322 unsigned n_el = Bulk_mesh_pt->nelement();
323 for(unsigned e=0;e<n_el;e++)
324 {
325 // Cast to a bulk element
326 ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
327
328 // Set the pointer to Poisson's ratio
329 el_pt->nu_pt() = &Global_Parameters::Nu;
330
331 // Set the pointer to Fourier wavenumber
332 el_pt->fourier_wavenumber_pt() = &Global_Parameters::Fourier_wavenumber;
333
334 // Set the pointer to non-dim Young's modulus
335 el_pt->youngs_modulus_pt() = &Global_Parameters::E;
336
337 // Set the pointer to square of the angular frequency
338 el_pt->omega_sq_pt() = &Global_Parameters::Omega_sq;
339
340 }// end loop over elements
341
342 // Loop over the traction elements
343 unsigned n_traction = Surface_mesh_pt->nelement();
344 for(unsigned e=0;e<n_traction;e++)
345 {
346 // Cast to a surface element
347 TimeHarmonicFourierDecomposedLinearElasticityTractionElement<ELEMENT>*
348 el_pt =
349 dynamic_cast<TimeHarmonicFourierDecomposedLinearElasticityTractionElement
350 <ELEMENT>* >(Surface_mesh_pt->element_pt(e));
351
352 // Set the applied traction
353 el_pt->traction_fct_pt() = &Global_Parameters::boundary_traction;
354
355 }// end loop over traction elements
356
357}
358
359//===start_of_traction===============================================
360/// Make traction elements along the boundary r=rmin
361//===================================================================
362template<class ELEMENT>
365{
366 unsigned bound, n_neigh;
367
368 // How many bulk elements are next to boundary 3
369 bound=3;
370 n_neigh = Bulk_mesh_pt->nboundary_element(bound);
371
372 // Now loop over bulk elements and create the face elements
373 for(unsigned n=0;n<n_neigh;n++)
374 {
375 // Create the face element
376 FiniteElement *traction_element_pt
377 = new TimeHarmonicFourierDecomposedLinearElasticityTractionElement<ELEMENT>
378 (Bulk_mesh_pt->boundary_element_pt(bound,n),
379 Bulk_mesh_pt->face_index_at_boundary(bound,n));
380
381 // Add to mesh
382 Surface_mesh_pt->add_element_pt(traction_element_pt);
383 }
384
385} // end of assign_traction_elements
386
387
388//===start_of_delete_traction========================================
389/// Delete traction elements
390//===================================================================
391template<class ELEMENT>
394{
395 // How many surface elements are in the surface mesh
396 unsigned n_element = Surface_mesh_pt->nelement();
397
398 // Loop over the surface elements
399 for(unsigned e=0;e<n_element;e++)
400 {
401 // Kill surface element
402 delete Surface_mesh_pt->element_pt(e);
403 }
404
405 // Wipe the mesh
406 Surface_mesh_pt->flush_element_and_node_storage();
407
408} // end of delete_traction_elements
409
410
411//==start_of_doc_solution=================================================
412/// Doc the solution
413//========================================================================
414template<class ELEMENT>
416doc_solution(DocInfo& doc_info)
417{
418 ofstream some_file;
419 char filename[100];
420
421 // Number of plot points
422 unsigned npts=5;
423
424 // Output solution
425 sprintf(filename,"%s/soln.dat",doc_info.directory().c_str());
426 some_file.open(filename);
427 Bulk_mesh_pt->output(some_file,npts);
428 some_file.close();
429
430 // Output norm of solution (to allow validation of solution even
431 // if triangle generates a slightly different mesh)
432 sprintf(filename,"%s/norm.dat",doc_info.directory().c_str());
433 some_file.open(filename);
434 double norm=0.0;
435 unsigned nel=Bulk_mesh_pt->nelement();
436 for (unsigned e=0;e<nel;e++)
437 {
438 double el_norm=0.0;
439 Bulk_mesh_pt->compute_norm(el_norm);
440 norm+=el_norm;
441 }
442 some_file << norm << std::endl;
443
444
445} // end_of_doc_solution
446
447
448//===start_of_main======================================================
449/// Driver code
450//======================================================================
451int main(int argc, char* argv[])
452{
453 // Number of elements in r-direction
454 unsigned nr=10;
455
456 // Number of elements in z-direction (for (approximately) square elements)
457 unsigned nz=unsigned(double(nr)*Global_Parameters::Lz/Global_Parameters::Lr);
458
459 // Set up doc info
460 DocInfo doc_info;
461
462 // Set output directory
463 doc_info.set_directory("RESLT");
464
465#ifdef ADAPTIVE
466
467 // Set up problem
469 <ProjectableTimeHarmonicFourierDecomposedLinearElasticityElement
470 <TTimeHarmonicFourierDecomposedLinearElasticityElement<3> > >
473
474 // Solve
475 unsigned max_adapt=3;
476 problem.newton_solve(max_adapt);
477
478#else
479
480 // Set up problem
482 <QTimeHarmonicFourierDecomposedLinearElasticityElement<3> >
485
486 // Solve
487 problem.newton_solve();
488
489#endif
490
491 // Output the solution
492 problem.doc_solution(doc_info);
493
494} // end_of_main
Class to validate time harmonic linear elasticity (Fourier decomposed)
Definition: cylinder.cc:134
void actions_before_adapt()
Actions before adapt: Wipe the mesh of traction elements.
RefineableTriangleMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the bulk mesh.
FourierDecomposedTimeHarmonicLinearElasticityProblem(const unsigned &nr, const unsigned &nz, const double &rmin, const double &rmax, const double &zmin, const double &zmax)
Constructor: Pass number of elements in r and z directions and boundary locations.
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of traction elements.
Mesh * Surface_mesh_pt
Pointer to the mesh of traction elements.
Definition: cylinder.cc:163
void assign_traction_elements()
Allocate traction elements on the bottom surface.
Definition: cylinder.cc:289
Mesh * Bulk_mesh_pt
Pointer to the bulk mesh.
Definition: cylinder.cc:160
void actions_before_newton_solve()
Update before solve is empty.
void complete_problem_setup()
Helper function to complete problem setup.
void doc_solution(DocInfo &doc_info)
Doc the solution.
Namespace for global parameters.
Definition: cylinder.cc:43
double Lz
Length of domain in z-direction.
Definition: cylinder.cc:65
const std::complex< double > I(0.0, 1.0)
Define the imaginary unit.
double Lr
Length of domain in r direction.
Definition: cylinder.cc:62
void boundary_traction(const Vector< double > &x, const Vector< double > &n, Vector< std::complex< double > > &result)
The traction function at r=rmin: (t_r, t_z, t_theta)
Definition: cylinder.cc:77
std::complex< double > Nu(0.3, 0.05)
Define Poisson's ratio Nu.
std::complex< double > Omega_sq(10.0, 5.0)
Define the non-dimensional square angular frequency of time-harmonic motion.
std::complex< double > E(1.0, 0.01)
Define the non-dimensional Young's modulus.
int Fourier_wavenumber
Define Fourier wavenumber.
Definition: cylinder.cc:55
int main(int argc, char *argv[])
Driver code.