periodic_load.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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 // Driver for a periodically loaded elastic body
31 
32 // The oomphlib headers
33 #include "generic.h"
34 #include "linear_elasticity.h"
35 
36 // The mesh
37 #include "meshes/rectangular_quadmesh.h"
38 
39 using namespace std;
40 
41 using namespace oomph;
42 
43 //===start_of_namespace=================================================
44 /// Namespace for global parameters
45 //======================================================================
46 namespace Global_Parameters
47 {
48  /// Amplitude of traction applied
49  double Amplitude = 1.0;
50 
51  /// \short Specify problem to be solved (boundary conditons for finite or
52  /// infinite domain).
53  bool Finite=false;
54 
55  /// Define Poisson coefficient Nu
56  double Nu = 0.3;
57 
58  /// Length of domain in x direction
59  double Lx = 1.0;
60 
61  /// Length of domain in y direction
62  double Ly = 2.0;
63 
64  /// The elasticity tensor
65  IsotropicElasticityTensor E(Nu);
66 
67  /// The exact solution for infinite depth case
68  void exact_solution(const Vector<double> &x,
69  Vector<double> &u)
70  {
71  u[0] = -Amplitude*cos(2.0*MathematicalConstants::Pi*x[0]/Lx)*
72  exp(2.0*MathematicalConstants::Pi*(x[1]-Ly))/
73  (2.0/(1.0+Nu)*MathematicalConstants::Pi);
74  u[1] = -Amplitude*sin(2.0*MathematicalConstants::Pi*x[0]/Lx)*
75  exp(2.0*MathematicalConstants::Pi*(x[1]-Ly))/
76  (2.0/(1.0+Nu)*MathematicalConstants::Pi);
77  }
78 
79  /// The traction function
80 void periodic_traction(const double &time,
81  const Vector<double> &x,
82  const Vector<double> &n,
83  Vector<double> &result)
84  {
85  result[0] = -Amplitude*cos(2.0*MathematicalConstants::Pi*x[0]/Lx);
86  result[1] = -Amplitude*sin(2.0*MathematicalConstants::Pi*x[0]/Lx);
87  }
88 } // end_of_namespace
89 
90 
91 //===start_of_problem_class=============================================
92 /// Periodic loading problem
93 //======================================================================
94 template<class ELEMENT>
95 class PeriodicLoadProblem : public Problem
96 {
97 public:
98 
99  /// \short Constructor: Pass number of elements in x and y directions
100  /// and lengths
101  PeriodicLoadProblem(const unsigned &nx, const unsigned &ny,
102  const double &lx, const double &ly);
103 
104  /// Update before solve is empty
106 
107  /// Update after solve is empty
109 
110  /// Doc the solution
111  void doc_solution(DocInfo& doc_info);
112 
113 private:
114 
115  /// Allocate traction elements on the top surface
116  void assign_traction_elements();
117 
118  /// Pointer to the bulk mesh
120 
121  /// Pointer to the mesh of traction elements
123 
124 }; // end_of_problem_class
125 
126 
127 //===start_of_constructor=============================================
128 /// Problem constructor: Pass number of elements in coordinate
129 /// directions and size of domain.
130 //====================================================================
131 template<class ELEMENT>
133 (const unsigned &nx, const unsigned &ny,
134  const double &lx, const double& ly)
135 {
136  //Now create the mesh with periodic boundary conditions in x direction
137  bool periodic_in_x=true;
138  Bulk_mesh_pt =
139  new RectangularQuadMesh<ELEMENT>(nx,ny,lx,ly,periodic_in_x);
140 
141  //Create the surface mesh of traction elements
142  Surface_mesh_pt=new Mesh;
143  assign_traction_elements();
144 
145  // Set the boundary conditions for this problem: All nodes are
146  // free by default -- just pin & set the ones that have Dirichlet
147  // conditions here
148  unsigned ibound=0;
149  unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
150  for (unsigned inod=0;inod<num_nod;inod++)
151  {
152  // Get pointer to node
153  Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(ibound,inod);
154 
155  // Pinned in x & y at the bottom and set value
156  nod_pt->pin(0);
157  nod_pt->pin(1);
158 
159  // Check which boundary conditions to set and set them
161  {
162  // Set the displacements to zero
163  nod_pt->set_value(0,0);
164  nod_pt->set_value(1,0);
165  }
166  else
167  {
168  // Extract nodal coordinates from node:
169  Vector<double> x(2);
170  x[0]=nod_pt->x(0);
171  x[1]=nod_pt->x(1);
172 
173  // Compute the value of the exact solution at the nodal point
174  Vector<double> u(2);
176 
177  // Assign these values to the nodal values at this node
178  nod_pt->set_value(0,u[0]);
179  nod_pt->set_value(1,u[1]);
180  };
181  } // end_loop_over_boundary_nodes
182 
183  // Complete the problem setup to make the elements fully functional
184 
185  // Loop over the elements
186  unsigned n_el = Bulk_mesh_pt->nelement();
187  for(unsigned e=0;e<n_el;e++)
188  {
189  // Cast to a bulk element
190  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
191 
192  // Set the elasticity tensor
193  el_pt->elasticity_tensor_pt() = &Global_Parameters::E;
194  }// end loop over elements
195 
196  // Loop over the traction elements
197  unsigned n_traction = Surface_mesh_pt->nelement();
198  for(unsigned e=0;e<n_traction;e++)
199  {
200  // Cast to a surface element
201  LinearElasticityTractionElement<ELEMENT> *el_pt =
202  dynamic_cast<LinearElasticityTractionElement<ELEMENT>* >
203  (Surface_mesh_pt->element_pt(e));
204 
205  // Set the applied traction
206  el_pt->traction_fct_pt() = &Global_Parameters::periodic_traction;
207  }// end loop over traction elements
208 
209  // Add the submeshes to the problem
210  add_sub_mesh(Bulk_mesh_pt);
211  add_sub_mesh(Surface_mesh_pt);
212 
213  // Now build the global mesh
214  build_global_mesh();
215 
216  // Assign equation numbers
217  cout << assign_eqn_numbers() << " equations assigned" << std::endl;
218 } // end of constructor
219 
220 
221 //===start_of_traction===============================================
222 /// Make traction elements along the top boundary of the bulk mesh
223 //===================================================================
224 template<class ELEMENT>
226 {
227 
228  // How many bulk elements are next to boundary 2 (the top boundary)?
229  unsigned bound=2;
230  unsigned n_neigh = Bulk_mesh_pt->nboundary_element(bound);
231 
232  // Now loop over bulk elements and create the face elements
233  for(unsigned n=0;n<n_neigh;n++)
234  {
235  // Create the face element
236  FiniteElement *traction_element_pt
237  = new LinearElasticityTractionElement<ELEMENT>
238  (Bulk_mesh_pt->boundary_element_pt(bound,n),
239  Bulk_mesh_pt->face_index_at_boundary(bound,n));
240 
241  // Add to mesh
242  Surface_mesh_pt->add_element_pt(traction_element_pt);
243  }
244 
245 } // end of assign_traction_elements
246 
247 //==start_of_doc_solution=================================================
248 /// Doc the solution
249 //========================================================================
250 template<class ELEMENT>
252 {
253  ofstream some_file;
254  char filename[100];
255 
256  // Number of plot points
257  unsigned npts=5;
258 
259  // Output solution
260  sprintf(filename,"%s/soln.dat",doc_info.directory().c_str());
261  some_file.open(filename);
262  Bulk_mesh_pt->output(some_file,npts);
263  some_file.close();
264 
265  // Output exact solution
266  sprintf(filename,"%s/exact_soln.dat",doc_info.directory().c_str());
267  some_file.open(filename);
268  Bulk_mesh_pt->output_fct(some_file,npts,
270  some_file.close();
271 
272  // Doc error
273  double error=0.0;
274  double norm=0.0;
275  sprintf(filename,"%s/error.dat",doc_info.directory().c_str());
276  some_file.open(filename);
277  Bulk_mesh_pt->compute_error(some_file,
279  error,norm);
280  some_file.close();
281 
282 // Doc error norm:
283  cout << "\nNorm of error " << sqrt(error) << std::endl;
284  cout << "Norm of solution : " << sqrt(norm) << std::endl << std::endl;
285  cout << std::endl;
286 
287 
288 } // end_of_doc_solution
289 
290 
291 //===start_of_main======================================================
292 /// Driver code for PeriodicLoad linearly elastic problem
293 //======================================================================
294 int main(int argc, char* argv[])
295 {
296  // Number of elements in x-direction
297  unsigned nx=5;
298 
299  // Number of elements in y-direction (for (approximately) square elements)
300  unsigned ny=unsigned(double(nx)*Global_Parameters::Ly/Global_Parameters::Lx);
301 
302  // Set up doc info
303  DocInfo doc_info;
304 
305  // Set output directory
306  doc_info.set_directory("RESLT");
307 
308  // Set up problem
311 
312  // Solve
313  problem.newton_solve();
314 
315  // Output the solution
316  problem.doc_solution(doc_info);
317 
318 } // end_of_main
Mesh * Surface_mesh_pt
Pointer to the mesh of traction elements.
void doc_solution(DocInfo &doc_info)
Doc the solution.
PeriodicLoadProblem(const unsigned &nx, const unsigned &ny, const double &lx, const double &ly)
Constructor: Pass number of elements in x and y directions and lengths.
void actions_before_newton_solve()
Update before solve is empty.
bool Finite
Specify problem to be solved (boundary conditons for finite or infinite domain).
double Lx
Length of domain in x direction.
Periodic loading problem.
void periodic_traction(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
The traction function.
void assign_traction_elements()
Allocate traction elements on the top surface.
Mesh * Bulk_mesh_pt
Pointer to the bulk mesh.
double Amplitude
Amplitude of traction applied.
int main(int argc, char *argv[])
Driver code for PeriodicLoad linearly elastic problem.
double Nu
Define Poisson coefficient Nu.
void exact_solution(const Vector< double > &x, Vector< double > &u)
The exact solution for infinite depth case.
IsotropicElasticityTensor E(Nu)
The elasticity tensor.
void actions_after_newton_solve()
Update after solve is empty.
double Ly
Length of domain in y direction.