unstructured_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-2023 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/triangle_mesh.h"
34 
35 using namespace std;
36 
37 using namespace oomph;
38 
39 
40 //===start_of_namespace=================================================
41 /// Namespace for global parameters
42 //======================================================================
43 namespace Global_Parameters
44 {
45  /// Define Poisson's ratio Nu
46  std::complex<double> Nu(0.3,0.05);
47 
48  /// Define the non-dimensional Young's modulus
49  std::complex<double> E(1.0,0.01);
50 
51  // Lame parameters
52  std::complex<double> lambda = E*Nu/(1.0+Nu)/(1.0-2.0*Nu);
53  std::complex<double> mu = E/2.0/(1.0+Nu);
54 
55  /// Define Fourier wavenumber
56  int Fourier_wavenumber = 3;
57 
58  /// Define the non-dimensional square angular frequency of
59  /// time-harmonic motion
60  std::complex<double> Omega_sq (10.0,5.0);
61 
62  /// Length of domain in r direction
63  double Lr = 1.0;
64 
65  /// Length of domain in z-direction
66  double Lz = 2.0;
67 
68  // Set up min & max (r,z) coordinates
69  double rmin = 0.1;
70  double zmin = 0.3;
71  double rmax = rmin+Lr;
72  double zmax = zmin+Lz;
73 
74  /// Define the imaginary unit
75  const std::complex<double> I(0.0,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  result[0] = -6.0*pow(x[0],2)*mu*cos(x[1])-
83  lambda*(I*double(Fourier_wavenumber)*pow(x[0],2)*pow(x[1],3)+
84  (4.0*pow(x[0],2)+pow(x[0],3))*cos(x[1]));
85  result[1] = -mu*(3.0*pow(x[0],2)-pow(x[0],3))*sin(x[1]);
86  result[2] = -mu*pow(x[0],2)*(2*pow(x[1],3)+I*double(Fourier_wavenumber)*
87  cos(x[1]));
88  }
89 
90 
91  /// The body force function; returns vector of complex doubles
92  /// in the order (b_r, b_z, b_theta)
93  void body_force(const Vector<double> &x,
94  Vector<std::complex<double> > &result)
95  {
96  result[0] =
97  x[0]*(-2.0*I*lambda*double(Fourier_wavenumber)*pow(x[1],3)-cos(x[1])*
98  (lambda*(8.0+3.0*x[0])-
99  mu*(pow(double(Fourier_wavenumber),2)
100  -16.0+x[0]*(x[0]-3.0))+pow(x[0],2)*Omega_sq));
101  result[1] =
102  x[0]*sin(x[1])*(mu*(pow(double(Fourier_wavenumber),2)-9.0)+
103  4.0*x[0]*(lambda+mu)+pow(x[0],2)*
104  (lambda+2.0*mu-Omega_sq))-
105  3.0*I*double(Fourier_wavenumber)*pow(x[0],2)*pow(x[1],2)*(lambda+mu);
106  result[2] =
107  -x[0]*(8.0*mu*pow(x[1],3)-pow(double(Fourier_wavenumber),2)*pow(x[1],3)*
108  (lambda+2.0*mu)+pow(x[0],2)*(pow(x[1],3)*Omega_sq+6.0*mu*x[1])+
109  I*cos(x[1])*double(Fourier_wavenumber)*
110  (lambda*(4.0+x[0])+mu*(6.0+x[0])));
111  }
112 
113  /// The exact solution in a flat-packed vector:
114  // 0: u_r[real], 1: u_z[real],..., 5: u_theta[imag]
115  void exact_solution(const Vector<double> &x,
116  Vector<double> &u)
117  {
118  u[0] = pow(x[0],3)*cos(x[1]);
119  u[1] = pow(x[0],3)*sin(x[1]);
120  u[2] = pow(x[0],3)*pow(x[1],3);
121  u[3] = 0.0;
122  u[4] = 0.0;
123  u[5] = 0.0;
124  }
125 
126 } // end_of_namespace
127 
128 
129 //===start_of_problem_class=============================================
130 /// Class to validate time harmonic linear elasticity (Fourier
131 /// decomposed)
132 //======================================================================
133 template<class ELEMENT>
135 {
136 public:
137 
138  /// Constructor: Pass boundary locations
140  const double& rmax,
141  const double &zmin,
142  const double& zmax);
143 
144  /// Update before solve is empty
146 
147  /// Update after solve is empty
149 
150 
151  /// Actions before adapt: Wipe the mesh of traction elements
153  {
154  // Kill the traction elements and wipe surface mesh
155  delete_traction_elements();
156 
157  // Rebuild the Problem's global mesh from its various sub-meshes
158  rebuild_global_mesh();
159  }
160 
161 
162  /// Actions after adapt: Rebuild the mesh of traction elements
164  {
165  // Create traction elements from all elements that are
166  // adjacent to FSI boundaries and add them to surface meshes
167  assign_traction_elements();
168 
169  // Rebuild the Problem's global mesh from its various sub-meshes
170  rebuild_global_mesh();
171 
172  // Complete problem setup
173  complete_problem_setup();
174  }
175 
176  /// Doc the solution
177  void doc_solution(DocInfo& doc_info);
178 
179 private:
180 
181  /// Allocate traction elements on the bottom surface
183 
184  /// Delete traction elements
186 
187  /// Helper function to complete problem setup
189 
190 #ifdef ADAPTIVE
191 
192  /// Pointer to the bulk mesh
193  RefineableTriangleMesh<ELEMENT>* Bulk_mesh_pt;
194 
195 #else
196 
197  /// Pointer to the bulk mesh
198  Mesh* Bulk_mesh_pt;
199 
200 #endif
201 
202  /// Pointer to the mesh of traction elements
203  Mesh* Surface_mesh_pt;
204 
205 }; // end_of_problem_class
206 
207 
208 //===start_of_constructor=============================================
209 /// Problem constructor: Pass size of domain.
210 //====================================================================
211 template<class ELEMENT>
214 (const double& rmin, const double& rmax, const double &zmin, const double& zmax)
215 {
216 
217  // The boundary is bounded by four distinct boundaries, each
218  // represented by its own polyline
219  Vector<TriangleMeshCurveSection*> boundary_polyline_pt(4);
220 
221  // Vertex coordinates on boundary
222  Vector<Vector<double> > bound_coords(2);
223  bound_coords[0].resize(2);
224  bound_coords[1].resize(2);
225 
226  // Horizontal bottom boundary
227  bound_coords[0][0]=rmin;
228  bound_coords[0][1]=zmin;
229  bound_coords[1][0]=rmax;
230  bound_coords[1][1]=zmin;
231 
232  // Build the boundary polyline
233  unsigned boundary_id=0;
234  boundary_polyline_pt[0]=new TriangleMeshPolyLine(bound_coords,boundary_id);
235 
236  // Vertical outer boundary
237  bound_coords[0][0]=rmax;
238  bound_coords[0][1]=zmin;
239  bound_coords[1][0]=rmax;
240  bound_coords[1][1]=zmax;
241 
242  // Build the boundary polyline
243  boundary_id=1;
244  boundary_polyline_pt[1]=new TriangleMeshPolyLine(bound_coords,boundary_id);
245 
246 
247  // Horizontal top boundary
248  bound_coords[0][0]=rmax;
249  bound_coords[0][1]=zmax;
250  bound_coords[1][0]=rmin;
251  bound_coords[1][1]=zmax;
252 
253  // Build the boundary polyline
254  boundary_id=2;
255  boundary_polyline_pt[2]=new TriangleMeshPolyLine(bound_coords,boundary_id);
256 
257  // Vertical inner boundary
258  bound_coords[0][0]=rmin;
259  bound_coords[0][1]=zmax;
260  bound_coords[1][0]=rmin;
261  bound_coords[1][1]=zmin;
262 
263  // Build the boundary polyline
264  boundary_id=3;
265  boundary_polyline_pt[3]=new TriangleMeshPolyLine(bound_coords,boundary_id);
266 
267  // Pointer to the closed curve that defines the outer boundary
268  TriangleMeshClosedCurve* closed_curve_pt=
269  new TriangleMeshPolygon(boundary_polyline_pt);
270 
271  // Use the TriangleMeshParameters object for helping on the manage of the
272  // TriangleMesh parameters
273  TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
274 
275  // Specify the maximum area element
276  double uniform_element_area=0.2;
277  triangle_mesh_parameters.element_area() = uniform_element_area;
278 
279 #ifdef ADAPTIVE
280 
281  // Create the mesh
282  Bulk_mesh_pt=new RefineableTriangleMesh<ELEMENT>(triangle_mesh_parameters);
283 
284  // Set error estimator
285  Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
286 
287 #else
288 
289  // Create the mesh
290  Bulk_mesh_pt=new TriangleMesh<ELEMENT>(triangle_mesh_parameters);
291 
292 #endif
293 
294  //Create the surface mesh of traction elements
295  Surface_mesh_pt=new Mesh;
296  assign_traction_elements();
297 
298  // Complete problem setup
299  complete_problem_setup();
300 
301  // Add the submeshes to the problem
302  add_sub_mesh(Bulk_mesh_pt);
303  add_sub_mesh(Surface_mesh_pt);
304 
305  // Now build the global mesh
306  build_global_mesh();
307 
308  // Assign equation numbers
309  cout << assign_eqn_numbers() << " equations assigned" << std::endl;
310 
311 } // end of constructor
312 
313 
314 
315 //===start_of_complete_problem_setup=================================
316 /// Complete problem setup
317 //===================================================================
318 template<class ELEMENT>
321 {
322 
323  // Set the boundary conditions for this problem: All nodes are
324  // free by default -- just pin & set the ones that have Dirichlet
325  // conditions here
326 
327  // storage for nodal position
328  Vector<double> x(2);
329 
330  // Storage for prescribed displacements
331  Vector<double> u(6);
332 
333  // Now set displacements on boundaries 0 (z=zmin),
334  //------------------------------------------------
335  // 1 (r=rmax) and 2 (z=zmax)
336  //--------------------------
337  for (unsigned ibound=0;ibound<=2;ibound++)
338  {
339  unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
340  for (unsigned inod=0;inod<num_nod;inod++)
341  {
342  // Get pointer to node
343  Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(ibound,inod);
344 
345  // get r and z coordinates
346  x[0]=nod_pt->x(0);
347  x[1]=nod_pt->x(1);
348 
349  // Pinned in r, z and theta
350  nod_pt->pin(0);nod_pt->pin(1);nod_pt->pin(2);
351  nod_pt->pin(3);nod_pt->pin(4);nod_pt->pin(5);
352 
353  // Compute the value of the exact solution at the nodal point
354  Vector<double> u(6);
356 
357  // Set the displacements
358  nod_pt->set_value(0,u[0]);
359  nod_pt->set_value(1,u[1]);
360  nod_pt->set_value(2,u[2]);
361  nod_pt->set_value(3,u[3]);
362  nod_pt->set_value(4,u[4]);
363  nod_pt->set_value(5,u[5]);
364  }
365  } // end_of_loop_over_boundary_nodes
366 
367 
368  // Complete the problem setup to make the elements fully functional
369 
370  // Loop over the elements
371  unsigned n_el = Bulk_mesh_pt->nelement();
372  for(unsigned e=0;e<n_el;e++)
373  {
374  // Cast to a bulk element
375  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
376 
377  // Set the body force
378  el_pt->body_force_fct_pt() = &Global_Parameters::body_force;
379 
380  // Set the pointer to Poisson's ratio
381  el_pt->nu_pt() = &Global_Parameters::Nu;
382 
383  // Set the pointer to Fourier wavenumber
384  el_pt->fourier_wavenumber_pt() = &Global_Parameters::Fourier_wavenumber;
385 
386  // Set the pointer to non-dim Young's modulus
387  el_pt->youngs_modulus_pt() = &Global_Parameters::E;
388 
389  // Set the pointer to square of the angular frequency
390  el_pt->omega_sq_pt() = &Global_Parameters::Omega_sq;
391 
392  }// end loop over elements
393 
394  // Loop over the traction elements
395  unsigned n_traction = Surface_mesh_pt->nelement();
396  for(unsigned e=0;e<n_traction;e++)
397  {
398  // Cast to a surface element
399  TimeHarmonicFourierDecomposedLinearElasticityTractionElement<ELEMENT>*
400  el_pt =
401  dynamic_cast<TimeHarmonicFourierDecomposedLinearElasticityTractionElement
402  <ELEMENT>* >(Surface_mesh_pt->element_pt(e));
403 
404  // Set the applied traction
405  el_pt->traction_fct_pt() = &Global_Parameters::boundary_traction;
406 
407  }// end loop over traction elements
408 }
409 
410 //===start_of_traction===============================================
411 /// Make traction elements along the boundary r=rmin
412 //===================================================================
413 template<class ELEMENT>
416 {
417  unsigned bound, n_neigh;
418 
419  // How many bulk elements are next to boundary 3
420  bound=3;
421  n_neigh = Bulk_mesh_pt->nboundary_element(bound);
422 
423  // Now loop over bulk elements and create the face elements
424  for(unsigned n=0;n<n_neigh;n++)
425  {
426  // Create the face element
427  FiniteElement *traction_element_pt
428  = new TimeHarmonicFourierDecomposedLinearElasticityTractionElement<ELEMENT>
429  (Bulk_mesh_pt->boundary_element_pt(bound,n),
430  Bulk_mesh_pt->face_index_at_boundary(bound,n));
431 
432  // Add to mesh
433  Surface_mesh_pt->add_element_pt(traction_element_pt);
434  }
435 
436 } // end of assign_traction_elements
437 
438 
439 
440 //===start_of_delete_traction========================================
441 /// Delete traction elements
442 //===================================================================
443 template<class ELEMENT>
446 {
447  // How many surface elements are in the surface mesh
448  unsigned n_element = Surface_mesh_pt->nelement();
449 
450  // Loop over the surface elements
451  for(unsigned e=0;e<n_element;e++)
452  {
453  // Kill surface element
454  delete Surface_mesh_pt->element_pt(e);
455  }
456 
457  // Wipe the mesh
458  Surface_mesh_pt->flush_element_and_node_storage();
459 
460 } // end of delete_traction_elements
461 
462 
463 //==start_of_doc_solution=================================================
464 /// Doc the solution
465 //========================================================================
466 template<class ELEMENT>
468 doc_solution(DocInfo& doc_info)
469 {
470  ofstream some_file;
471  char filename[100];
472 
473  // Number of plot points
474  unsigned npts=10;
475 
476  // Output solution
477  sprintf(filename,"%s/soln.dat",doc_info.directory().c_str());
478  some_file.open(filename);
479  Bulk_mesh_pt->output(some_file,npts);
480  some_file.close();
481 
482  // Output exact solution
483  sprintf(filename,"%s/exact_soln.dat",doc_info.directory().c_str());
484  some_file.open(filename);
485  Bulk_mesh_pt->output_fct(some_file,npts,
487  some_file.close();
488 
489  // Doc error
490  double error=0.0;
491  double norm=0.0;
492  sprintf(filename,"%s/error.dat",doc_info.directory().c_str());
493  some_file.open(filename);
494  Bulk_mesh_pt->compute_error(some_file,
496  error,norm);
497  some_file.close();
498 
499  // Doc error norm:
500  cout << "\nNorm of error: " << sqrt(error) << std::endl;
501  cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
502  cout << std::endl;
503 
504  // Output norm of solution (to allow validation of solution even
505  // if triangle generates a slightly different mesh)
506  sprintf(filename,"%s/norm.dat",doc_info.directory().c_str());
507  some_file.open(filename);
508  some_file << norm << std::endl;
509  some_file.close();
510 
511 } // end_of_doc_solution
512 
513 
514 //===start_of_main======================================================
515 /// Driver code
516 //======================================================================
517 int main(int argc, char* argv[])
518 {
519 
520  // Set up doc info
521  DocInfo doc_info;
522 
523  // Set output directory
524  doc_info.set_directory("RESLT");
525 
526 #ifdef ADAPTIVE
527 
528  // Set up problem
530  <ProjectableTimeHarmonicFourierDecomposedLinearElasticityElement
531  <TTimeHarmonicFourierDecomposedLinearElasticityElement<3> > >
534 
535  // Solve
536  unsigned max_adapt=3;
537  problem.newton_solve(max_adapt);
538 
539 #else
540 
541  // Set up problem
543  <TTimeHarmonicFourierDecomposedLinearElasticityElement<3> >
546 
547  // Solve
548  problem.newton_solve();
549 
550 #endif
551 
552 
553  // Output the solution
554  problem.doc_solution(doc_info);
555 
556 } // 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.
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
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of traction elements.
void assign_traction_elements()
Allocate traction elements on the bottom surface.
void actions_after_newton_solve()
Update after solve is empty.
void delete_traction_elements()
Delete traction elements.
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.
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
int main(int argc, char *argv[])
Driver code.