refineable_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 -- spatially adaptive version
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  /// 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 RefineablePeriodicLoadProblem : public Problem
96 {
97 public:
98 
99  /// \short Constructor: Pass number of elements in x and y directions
100  /// and lengths.
101  RefineablePeriodicLoadProblem(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  /// Actions before adapt: Wipe the mesh of traction elements
112  {
113  // Kill the traction elements and wipe surface mesh
114  delete_traction_elements();
115 
116  // Rebuild the Problem's global mesh from its various sub-meshes
117  rebuild_global_mesh();
118  }
119 
120  /// Actions after adapt: Rebuild the mesh of traction elements
122  {
123  // Create traction elements
124  assign_traction_elements();
125 
126  // Rebuild the Problem's global mesh from its various sub-meshes
127  rebuild_global_mesh();
128  }
129 
130  /// Doc the solution
131  void doc_solution(DocInfo& doc_info);
132 
133 private:
134 
135  /// Allocate traction elements on the top surface
136  void assign_traction_elements();
137 
138  /// Kill traction elements on the top surface
140  {
141  // How many surface elements are in the surface mesh
142  unsigned n_element = Surface_mesh_pt->nelement();
143 
144  // Loop over the traction elements
145  for(unsigned e=0;e<n_element;e++)
146  {
147  // Kill surface element
148  delete Surface_mesh_pt->element_pt(e);
149  }
150 
151  // Wipe the mesh
152  Surface_mesh_pt->flush_element_and_node_storage();
153  }
154 
155  /// Pointer to the (refineable!) bulk mesh
156  TreeBasedRefineableMeshBase* Bulk_mesh_pt;
157 
158  /// Pointer to the mesh of traction elements
160 
161 }; // end_of_problem_class
162 
163 
164 //===start_of_constructor=============================================
165 /// Problem constructor: Pass number of elements in the coordinate
166 /// directions and the domain sizes.
167 //====================================================================
168 template<class ELEMENT>
170 (const unsigned &nx, const unsigned &ny,
171  const double &lx, const double& ly)
172 {
173  // Create the mesh
174  Bulk_mesh_pt = new RefineableRectangularQuadMesh<ELEMENT>(nx,ny,lx,ly);
175 
176  // Create/set error estimator
177  Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
178 
179  // Make the mesh periodic in the x-direction by setting the nodes on
180  // right boundary (boundary 1) to be the periodic counterparts of
181  // those on the left one (boundary 3).
182  unsigned n_node = Bulk_mesh_pt->nboundary_node(1);
183  for(unsigned n=0;n<n_node;n++)
184  {
185  Bulk_mesh_pt->boundary_node_pt(1,n)
186  ->make_periodic(Bulk_mesh_pt->boundary_node_pt(3,n));
187  }
188 
189 
190  // Now establish the new neighbours (created by "wrapping around"
191  // the domain) in the TreeForst representation of the mesh
192 
193  // Get pointers to tree roots associated with elements on the
194  // left and right boundaries
195  Vector<TreeRoot*> left_root_pt(ny);
196  Vector<TreeRoot*> right_root_pt(ny);
197  for(unsigned i=0;i<ny;i++)
198  {
199  left_root_pt[i] =
200  dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(i*nx))->
201  tree_pt()->root_pt();
202  right_root_pt[i] =
203  dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(nx-1+i*nx))->
204  tree_pt()->root_pt();
205  }
206 
207  // Switch on QuadTreeNames for enumeration of directions
208  using namespace QuadTreeNames;
209 
210  //Set the neighbour and periodicity
211  for(unsigned i=0;i<ny;i++)
212  {
213  // The western neighbours of the elements on the left
214  // boundary are those on the right
215  left_root_pt[i]->neighbour_pt(W) = right_root_pt[i];
216  left_root_pt[i]->set_neighbour_periodic(W);
217 
218  // The eastern neighbours of the elements on the right
219  // boundary are those on the left
220  right_root_pt[i]->neighbour_pt(E) = left_root_pt[i];
221  right_root_pt[i]->set_neighbour_periodic(E);
222  } // done
223 
224 
225  //Create the surface mesh of traction elements
226  Surface_mesh_pt=new Mesh;
227  assign_traction_elements();
228 
229  // Set the boundary conditions for this problem: All nodes are
230  // free by default -- just pin & set the ones that have Dirichlet
231  // conditions here
232  unsigned ibound=0;
233  unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
234  for (unsigned inod=0;inod<num_nod;inod++)
235  {
236  // Get pointer to node
237  Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(ibound,inod);
238 
239  // Pinned in x & y at the bottom and set value
240  nod_pt->pin(0);
241  nod_pt->pin(1);
242 
243  // Check which boundary conditions to set and set them
245  {
246  // Set both displacements to zero
247  nod_pt->set_value(0,0);
248  nod_pt->set_value(1,0);
249  }
250  else
251  {
252  // Extract nodal coordinates from node:
253  Vector<double> x(2);
254  x[0]=nod_pt->x(0);
255  x[1]=nod_pt->x(1);
256 
257  // Compute the value of the exact solution at the nodal point
258  Vector<double> u(2);
260 
261  // Assign these values to the nodal values at this node
262  nod_pt->set_value(0,u[0]);
263  nod_pt->set_value(1,u[1]);
264  };
265  } // end_loop_over_boundary_nodes
266 
267  // Complete the problem setup to make the elements fully functional
268 
269  // Loop over the elements
270  unsigned n_el = Bulk_mesh_pt->nelement();
271  for(unsigned e=0;e<n_el;e++)
272  {
273  // Cast to a bulk element
274  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
275 
276  // Set the elasticity tensor
277  el_pt->elasticity_tensor_pt() = &Global_Parameters::E;
278  }// end loop over elements
279 
280 
281  // Do selective refinement of one element so that we can test
282  // whether periodic hanging nodes work: Choose a single element
283  // (the zero-th one) as the to-be-refined element.
284  // This creates a hanging node on the periodic boundary
285  Vector<unsigned> refine_pattern(1,0);
286  Bulk_mesh_pt->refine_selected_elements(refine_pattern);
287 
288  // Add the submeshes to the problem
289  add_sub_mesh(Bulk_mesh_pt);
290  add_sub_mesh(Surface_mesh_pt);
291 
292  // Now build the global mesh
293  build_global_mesh();
294 
295  // Assign equation numbers
296  cout << assign_eqn_numbers() << " equations assigned" << std::endl;
297 } // end of constructor
298 
299 
300 //===start_of_traction===============================================
301 /// Make traction elements along the top boundary of the bulk mesh
302 //===================================================================
303 template<class ELEMENT>
305 {
306 
307  // How many bulk elements are next to boundary 2 (the top boundary)?
308  unsigned bound=2;
309  unsigned n_neigh = Bulk_mesh_pt->nboundary_element(bound);
310 
311  // Now loop over bulk elements and create the face elements
312  for(unsigned n=0;n<n_neigh;n++)
313  {
314  // Create the face element
315  FiniteElement *traction_element_pt
316  = new LinearElasticityTractionElement<ELEMENT>
317  (Bulk_mesh_pt->boundary_element_pt(bound,n),
318  Bulk_mesh_pt->face_index_at_boundary(bound,n));
319 
320  // Add to mesh
321  Surface_mesh_pt->add_element_pt(traction_element_pt);
322  }
323 
324  // Now set function pointer to applied traction
325  unsigned n_traction = Surface_mesh_pt->nelement();
326  for(unsigned e=0;e<n_traction;e++)
327  {
328  // Cast to a surface element
329  LinearElasticityTractionElement<ELEMENT> *el_pt =
330  dynamic_cast<LinearElasticityTractionElement<ELEMENT>* >
331  (Surface_mesh_pt->element_pt(e));
332 
333  // Set the applied traction
334  el_pt->traction_fct_pt() = &Global_Parameters::periodic_traction;
335  }// end loop over traction elements
336 
337 
338 } // end of assign_traction_elements
339 
340 //==start_of_doc_solution=================================================
341 /// Doc the solution
342 //========================================================================
343 template<class ELEMENT>
345 {
346  ofstream some_file;
347  char filename[100];
348 
349  // Number of plot points
350  unsigned npts=5;
351 
352  // Output solution
353  sprintf(filename,"%s/soln.dat",doc_info.directory().c_str());
354  some_file.open(filename);
355  Bulk_mesh_pt->output(some_file,npts);
356  some_file.close();
357 
358  // Output exact solution
359  sprintf(filename,"%s/exact_soln.dat",doc_info.directory().c_str());
360  some_file.open(filename);
361  Bulk_mesh_pt->output_fct(some_file,npts,
363  some_file.close();
364 
365  // Doc error
366  double error=0.0;
367  double norm=0.0;
368  sprintf(filename,"%s/error.dat",doc_info.directory().c_str());
369  some_file.open(filename);
370  Bulk_mesh_pt->compute_error(some_file,
372  error,norm);
373  some_file.close();
374 
375 // Doc error norm:
376  cout << "\nNorm of error " << sqrt(error) << std::endl;
377  cout << "Norm of solution : " << sqrt(norm) << std::endl << std::endl;
378  cout << std::endl;
379 
380 
381 } // end_of_doc_solution
382 
383 
384 //===start_of_main======================================================
385 /// Driver code for PeriodicLoad linearly elastic problem
386 //======================================================================
387 int main(int argc, char* argv[])
388 {
389  // Number of elements in x-direction
390  unsigned nx=2;
391 
392  // Number of elements in y-direction (for (approximately) square elements)
393  unsigned ny=unsigned(double(nx)*Global_Parameters::Ly/Global_Parameters::Lx);
394 
395  // Set up doc info
396  DocInfo doc_info;
397 
398  // Set output directory
399  doc_info.set_directory("RESLT");
400 
401  // Set up problem
404 
405  // Run the simulation, allowing for up to 4 spatial adaptations
406  unsigned max_adapt=4;
407  problem.newton_solve(max_adapt);
408 
409  // Output the solution
410  problem.doc_solution(doc_info);
411 
412 } // end_of_main
void actions_after_newton_solve()
Update after solve is empty.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of traction elements.
RefineablePeriodicLoadProblem(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).
void assign_traction_elements()
Allocate traction elements on the top surface.
double Lx
Length of domain in x direction.
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of traction elements.
TreeBasedRefineableMeshBase * Bulk_mesh_pt
Pointer to the (refineable!) bulk mesh.
void periodic_traction(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
The traction function.
Mesh * Surface_mesh_pt
Pointer to the mesh of traction elements.
int main(int argc, char *argv[])
Driver code for PeriodicLoad linearly elastic problem.
void doc_solution(DocInfo &doc_info)
Doc the solution.
double Amplitude
Amplitude of traction applied.
double Nu
Define Poisson coefficient Nu.
void exact_solution(const Vector< double > &x, Vector< double > &u)
The exact solution for infinite depth case.
void delete_traction_elements()
Kill traction elements on the top surface.
IsotropicElasticityTensor E(Nu)
The elasticity tensor.
double Ly
Length of domain in y direction.