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
35using namespace std;
36
37using namespace oomph;
38
39//===start_of_namespace=================================================
40/// Namespace for global parameters
41//======================================================================
43{
44 /// Define Poisson's ratio Nu
45 std::complex<double> Nu(0.3,0.05);
46
47 /// Define the non-dimensional Young's modulus
48 std::complex<double> E(1.0,0.01);
49
50 // Lame parameters
51 std::complex<double> lambda = E*Nu/(1.0+Nu)/(1.0-2.0*Nu);
52 std::complex<double> mu = E/2.0/(1.0+Nu);
53
54 /// Define Fourier wavenumber
56
57 /// Define the non-dimensional square angular frequency of
58 /// time-harmonic motion
59 std::complex<double> Omega_sq (10.0,5.0);
60
61 /// Length of domain in r direction
62 double Lr = 1.0;
63
64 /// Length of domain in z-direction
65 double Lz = 2.0;
66
67 // Set up min & max (r,z) coordinates
68 double rmin = 0.1;
69 double zmin = 0.3;
70 double rmax = rmin+Lr;
71 double zmax = zmin+Lz;
72
73 /// Define the imaginary unit
74 const std::complex<double> I(0.0,1.0);
75
76 /// The traction function at r=rmin: (t_r, t_z, t_theta)
77 void boundary_traction(const Vector<double> &x,
78 const Vector<double> &n,
79 Vector<std::complex<double> > &result)
80 {
81 result[0] = -6.0*pow(x[0],2)*mu*cos(x[1])-
82 lambda*(I*double(Fourier_wavenumber)*pow(x[0],2)*pow(x[1],3)+
83 (4.0*pow(x[0],2)+pow(x[0],3))*cos(x[1]));
84 result[1] = -mu*(3.0*pow(x[0],2)-pow(x[0],3))*sin(x[1]);
85 result[2] = -mu*pow(x[0],2)*(2*pow(x[1],3)+I*double(Fourier_wavenumber)*
86 cos(x[1]));
87 }
88
89
90 /// The body force function; returns vector of complex doubles
91 /// in the order (b_r, b_z, b_theta)
92 void body_force(const Vector<double> &x,
93 Vector<std::complex<double> > &result)
94 {
95 result[0] =
96 x[0]*(-2.0*I*lambda*double(Fourier_wavenumber)*pow(x[1],3)-cos(x[1])*
97 (lambda*(8.0+3.0*x[0])-
98 mu*(pow(double(Fourier_wavenumber),2)
99 -16.0+x[0]*(x[0]-3.0))+pow(x[0],2)*Omega_sq));
100 result[1] =
101 x[0]*sin(x[1])*(mu*(pow(double(Fourier_wavenumber),2)-9.0)+
102 4.0*x[0]*(lambda+mu)+pow(x[0],2)*
103 (lambda+2.0*mu-Omega_sq))-
104 3.0*I*double(Fourier_wavenumber)*pow(x[0],2)*pow(x[1],2)*(lambda+mu);
105 result[2] =
106 -x[0]*(8.0*mu*pow(x[1],3)-pow(double(Fourier_wavenumber),2)*pow(x[1],3)*
107 (lambda+2.0*mu)+pow(x[0],2)*(pow(x[1],3)*Omega_sq+6.0*mu*x[1])+
108 I*cos(x[1])*double(Fourier_wavenumber)*
109 (lambda*(4.0+x[0])+mu*(6.0+x[0])));
110 }
111
112 /// The exact solution in a flat-packed vector:
113 // 0: u_r[real], 1: u_z[real],..., 5: u_theta[imag]
114 void exact_solution(const Vector<double> &x,
115 Vector<double> &u)
116 {
117 u[0] = pow(x[0],3)*cos(x[1]);
118 u[1] = pow(x[0],3)*sin(x[1]);
119 u[2] = pow(x[0],3)*pow(x[1],3);
120 u[3] = 0.0;
121 u[4] = 0.0;
122 u[5] = 0.0;
123 }
124
125} // end_of_namespace
126
127
128//===start_of_problem_class=============================================
129/// Class to validate time harmonic linear elasticity (Fourier
130/// decomposed)
131//======================================================================
132template<class ELEMENT>
134{
135public:
136
137 /// Constructor: Pass number of elements in r and z directions
138 /// and boundary locations
140 const unsigned &nr, const unsigned &nz,
141 const double &rmin, const double& rmax,
142 const double &zmin, const double& zmax);
143
144
145 /// Update before solve is empty
147
148 /// Update after solve is empty
150
151 /// Doc the solution
152 void doc_solution(DocInfo& doc_info);
153
154private:
155
156 /// Allocate traction elements on the bottom surface
158
159 /// Pointer to the bulk mesh
161
162 /// Pointer to the mesh of traction elements
164}; // end_of_problem_class
165
166
167//===start_of_constructor=============================================
168/// Problem constructor: Pass number of elements in coordinate
169/// directions and size of domain.
170//====================================================================
171template<class ELEMENT>
174(const unsigned &nr, const unsigned &nz,
175 const double &rmin, const double& rmax,
176 const double &zmin, const double& zmax)
177{
178 //Now create the mesh
179 Bulk_mesh_pt = new RectangularQuadMesh<ELEMENT>(nr,nz,rmin,rmax,zmin,zmax);
180
181 //Create the surface mesh of traction elements
182 Surface_mesh_pt=new Mesh;
183 assign_traction_elements();
184
185 // Set the boundary conditions for this problem: All nodes are
186 // free by default -- just pin & set the ones that have Dirichlet
187 // conditions here
188
189 // storage for nodal position
190 Vector<double> x(2);
191
192 // Storage for prescribed displacements
193 Vector<double> u(6);
194
195 // Now set displacements on boundaries 0 (z=zmin),
196 //------------------------------------------------
197 // 1 (r=rmax) and 2 (z=zmax)
198 //--------------------------
199 for (unsigned ibound=0;ibound<=2;ibound++)
200 {
201 unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
202 for (unsigned inod=0;inod<num_nod;inod++)
203 {
204 // Get pointer to node
205 Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(ibound,inod);
206
207 // get r and z coordinates
208 x[0]=nod_pt->x(0);
209 x[1]=nod_pt->x(1);
210
211 // Pinned in r, z and theta
212 nod_pt->pin(0);nod_pt->pin(1);nod_pt->pin(2);
213 nod_pt->pin(3);nod_pt->pin(4);nod_pt->pin(5);
214
215 // Compute the value of the exact solution at the nodal point
216 Vector<double> u(6);
218
219 // Set the displacements
220 nod_pt->set_value(0,u[0]);
221 nod_pt->set_value(1,u[1]);
222 nod_pt->set_value(2,u[2]);
223 nod_pt->set_value(3,u[3]);
224 nod_pt->set_value(4,u[4]);
225 nod_pt->set_value(5,u[5]);
226 }
227 } // end_of_loop_over_boundary_nodes
228
229
230 // Complete the problem setup to make the elements fully functional
231
232 // Loop over the elements
233 unsigned n_el = Bulk_mesh_pt->nelement();
234 for(unsigned e=0;e<n_el;e++)
235 {
236 // Cast to a bulk element
237 ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
238
239 // Set the body force
240 el_pt->body_force_fct_pt() = &Global_Parameters::body_force;
241
242 // Set the pointer to Poisson's ratio
243 el_pt->nu_pt() = &Global_Parameters::Nu;
244
245 // Set the pointer to Fourier wavenumber
246 el_pt->fourier_wavenumber_pt() = &Global_Parameters::Fourier_wavenumber;
247
248 // Set the pointer to non-dim Young's modulus
249 el_pt->youngs_modulus_pt() = &Global_Parameters::E;
250
251 // Set the pointer to square of the angular frequency
252 el_pt->omega_sq_pt() = &Global_Parameters::Omega_sq;
253
254 }// end loop over elements
255
256 // Loop over the traction elements
257 unsigned n_traction = Surface_mesh_pt->nelement();
258 for(unsigned e=0;e<n_traction;e++)
259 {
260 // Cast to a surface element
261 TimeHarmonicFourierDecomposedLinearElasticityTractionElement<ELEMENT>*
262 el_pt =
263 dynamic_cast<TimeHarmonicFourierDecomposedLinearElasticityTractionElement
264 <ELEMENT>* >(Surface_mesh_pt->element_pt(e));
265
266 // Set the applied traction
267 el_pt->traction_fct_pt() = &Global_Parameters::boundary_traction;
268
269 }// end loop over traction elements
270
271 // Add the submeshes to the problem
272 add_sub_mesh(Bulk_mesh_pt);
273 add_sub_mesh(Surface_mesh_pt);
274
275 // Now build the global mesh
276 build_global_mesh();
277
278 // Assign equation numbers
279 cout << assign_eqn_numbers() << " equations assigned" << std::endl;
280
281} // end of constructor
282
283
284//===start_of_traction===============================================
285/// Make traction elements along the boundary r=rmin
286//===================================================================
287template<class ELEMENT>
290{
291 unsigned bound, n_neigh;
292
293 // How many bulk elements are next to boundary 3
294 bound=3;
295 n_neigh = Bulk_mesh_pt->nboundary_element(bound);
296
297 // Now loop over bulk elements and create the face elements
298 for(unsigned n=0;n<n_neigh;n++)
299 {
300 // Create the face element
301 FiniteElement *traction_element_pt
302 = new TimeHarmonicFourierDecomposedLinearElasticityTractionElement<ELEMENT>
303 (Bulk_mesh_pt->boundary_element_pt(bound,n),
304 Bulk_mesh_pt->face_index_at_boundary(bound,n));
305
306 // Add to mesh
307 Surface_mesh_pt->add_element_pt(traction_element_pt);
308 }
309
310} // end of assign_traction_elements
311
312
313//==start_of_doc_solution=================================================
314/// Doc the solution
315//========================================================================
316template<class ELEMENT>
318doc_solution(DocInfo& doc_info)
319{
320 ofstream some_file;
321 char filename[100];
322
323 // Number of plot points
324 unsigned npts=5;
325
326 // Output solution
327 sprintf(filename,"%s/soln.dat",doc_info.directory().c_str());
328 some_file.open(filename);
329 Bulk_mesh_pt->output(some_file,npts);
330 some_file.close();
331
332 // Output exact solution
333 sprintf(filename,"%s/exact_soln.dat",doc_info.directory().c_str());
334 some_file.open(filename);
335 Bulk_mesh_pt->output_fct(some_file,npts,
337 some_file.close();
338
339 // Doc error
340 double error=0.0;
341 double norm=0.0;
342 sprintf(filename,"%s/error.dat",doc_info.directory().c_str());
343 some_file.open(filename);
344 Bulk_mesh_pt->compute_error(some_file,
346 error,norm);
347 some_file.close();
348
349 // Doc error norm:
350 cout << "\nNorm of error: " << sqrt(error) << std::endl;
351 cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
352 cout << std::endl;
353
354} // end_of_doc_solution
355
356
357//===start_of_main======================================================
358/// Driver code
359//======================================================================
360int main(int argc, char* argv[])
361{
362 // Number of elements in r-direction
363 unsigned nr=5;
364
365 // Number of elements in z-direction (for (approximately) square elements)
366 unsigned nz=unsigned(double(nr)*Global_Parameters::Lz/Global_Parameters::Lr);
367
368 // Set up doc info
369 DocInfo doc_info;
370
371 // Set output directory
372 doc_info.set_directory("RESLT");
373
374 // Set up problem
376 <QTimeHarmonicFourierDecomposedLinearElasticityElement<3> >
379
380 // Solve
381 problem.newton_solve();
382
383 // Output the solution
384 problem.doc_solution(doc_info);
385
386} // end_of_main
Class to validate time harmonic linear elasticity (Fourier decomposed)
Definition: cylinder.cc:134
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.
Definition: cylinder.cc:174
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_after_newton_solve()
Update after solve is empty.
Definition: cylinder.cc:149
void actions_before_newton_solve()
Update before solve is empty.
Definition: cylinder.cc:146
void doc_solution(DocInfo &doc_info)
Doc the solution.
Definition: cylinder.cc:318
int main(int argc, char *argv[])
Driver code.
Definition: cylinder.cc:360
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.
std::complex< double > lambda
Definition: cylinder.cc:51
double Lr
Length of domain in r direction.
Definition: cylinder.cc:62
std::complex< double > mu
Definition: cylinder.cc:52
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.
void exact_solution(const Vector< double > &x, Vector< double > &u)
The exact solution in a flat-packed vector:
Definition: cylinder.cc:114
void body_force(const Vector< double > &x, Vector< std::complex< double > > &result)
The body force function; returns vector of complex doubles in the order (b_r, b_z,...
Definition: cylinder.cc:92
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