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// 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 for a periodically loaded elastic body
27 
28 // The oomphlib headers
29 #include "generic.h"
30 #include "linear_elasticity.h"
31 
32 // The mesh
33 #include "meshes/rectangular_quadmesh.h"
34 
35 using namespace std;
36 
37 using namespace oomph;
38 
39 //===start_of_namespace=================================================
40 /// Namespace for global parameters
41 //======================================================================
43 {
44  /// Amplitude of traction applied
45  double Amplitude = 1.0;
46 
47  /// Specify problem to be solved (boundary conditons for finite or
48  /// infinite domain).
49  bool Finite=false;
50 
51  /// Define Poisson coefficient Nu
52  double Nu = 0.3;
53 
54  /// Length of domain in x direction
55  double Lx = 1.0;
56 
57  /// Length of domain in y direction
58  double Ly = 2.0;
59 
60  /// The elasticity tensor
61  IsotropicElasticityTensor E(Nu);
62 
63  /// The exact solution for infinite depth case
64  void exact_solution(const Vector<double> &x,
65  Vector<double> &u)
66  {
67  u[0] = -Amplitude*cos(2.0*MathematicalConstants::Pi*x[0]/Lx)*
68  exp(2.0*MathematicalConstants::Pi*(x[1]-Ly))/
69  (2.0/(1.0+Nu)*MathematicalConstants::Pi);
70  u[1] = -Amplitude*sin(2.0*MathematicalConstants::Pi*x[0]/Lx)*
71  exp(2.0*MathematicalConstants::Pi*(x[1]-Ly))/
72  (2.0/(1.0+Nu)*MathematicalConstants::Pi);
73  }
74 
75  /// The traction function
76 void periodic_traction(const double &time,
77  const Vector<double> &x,
78  const Vector<double> &n,
79  Vector<double> &result)
80  {
81  result[0] = -Amplitude*cos(2.0*MathematicalConstants::Pi*x[0]/Lx);
82  result[1] = -Amplitude*sin(2.0*MathematicalConstants::Pi*x[0]/Lx);
83  }
84 } // end_of_namespace
85 
86 
87 //===start_of_problem_class=============================================
88 /// Periodic loading problem
89 //======================================================================
90 template<class ELEMENT>
91 class PeriodicLoadProblem : public Problem
92 {
93 public:
94 
95  /// Constructor: Pass number of elements in x and y directions
96  /// and lengths
97  PeriodicLoadProblem(const unsigned &nx, const unsigned &ny,
98  const double &lx, const double &ly);
99 
100  /// Update before solve is empty
102 
103  /// Update after solve is empty
105 
106  /// Doc the solution
107  void doc_solution(DocInfo& doc_info);
108 
109 private:
110 
111  /// Allocate traction elements on the top surface
112  void assign_traction_elements();
113 
114  /// Pointer to the bulk mesh
116 
117  /// Pointer to the mesh of traction elements
119 
120 }; // end_of_problem_class
121 
122 
123 //===start_of_constructor=============================================
124 /// Problem constructor: Pass number of elements in coordinate
125 /// directions and size of domain.
126 //====================================================================
127 template<class ELEMENT>
129 (const unsigned &nx, const unsigned &ny,
130  const double &lx, const double& ly)
131 {
132  //Now create the mesh with periodic boundary conditions in x direction
133  bool periodic_in_x=true;
134  Bulk_mesh_pt =
135  new RectangularQuadMesh<ELEMENT>(nx,ny,lx,ly,periodic_in_x);
136 
137  //Create the surface mesh of traction elements
138  Surface_mesh_pt=new Mesh;
139  assign_traction_elements();
140 
141  // Set the boundary conditions for this problem: All nodes are
142  // free by default -- just pin & set the ones that have Dirichlet
143  // conditions here
144  unsigned ibound=0;
145  unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
146  for (unsigned inod=0;inod<num_nod;inod++)
147  {
148  // Get pointer to node
149  Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(ibound,inod);
150 
151  // Pinned in x & y at the bottom and set value
152  nod_pt->pin(0);
153  nod_pt->pin(1);
154 
155  // Check which boundary conditions to set and set them
157  {
158  // Set the displacements to zero
159  nod_pt->set_value(0,0);
160  nod_pt->set_value(1,0);
161  }
162  else
163  {
164  // Extract nodal coordinates from node:
165  Vector<double> x(2);
166  x[0]=nod_pt->x(0);
167  x[1]=nod_pt->x(1);
168 
169  // Compute the value of the exact solution at the nodal point
170  Vector<double> u(2);
172 
173  // Assign these values to the nodal values at this node
174  nod_pt->set_value(0,u[0]);
175  nod_pt->set_value(1,u[1]);
176  };
177  } // end_loop_over_boundary_nodes
178 
179  // Complete the problem setup to make the elements fully functional
180 
181  // Loop over the elements
182  unsigned n_el = Bulk_mesh_pt->nelement();
183  for(unsigned e=0;e<n_el;e++)
184  {
185  // Cast to a bulk element
186  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
187 
188  // Set the elasticity tensor
189  el_pt->elasticity_tensor_pt() = &Global_Parameters::E;
190  }// end loop over elements
191 
192  // Loop over the traction elements
193  unsigned n_traction = Surface_mesh_pt->nelement();
194  for(unsigned e=0;e<n_traction;e++)
195  {
196  // Cast to a surface element
197  LinearElasticityTractionElement<ELEMENT> *el_pt =
198  dynamic_cast<LinearElasticityTractionElement<ELEMENT>* >
199  (Surface_mesh_pt->element_pt(e));
200 
201  // Set the applied traction
202  el_pt->traction_fct_pt() = &Global_Parameters::periodic_traction;
203  }// end loop over traction elements
204 
205  // Add the submeshes to the problem
206  add_sub_mesh(Bulk_mesh_pt);
207  add_sub_mesh(Surface_mesh_pt);
208 
209  // Now build the global mesh
210  build_global_mesh();
211 
212  // Assign equation numbers
213  cout << assign_eqn_numbers() << " equations assigned" << std::endl;
214 } // end of constructor
215 
216 
217 //===start_of_traction===============================================
218 /// Make traction elements along the top boundary of the bulk mesh
219 //===================================================================
220 template<class ELEMENT>
222 {
223 
224  // How many bulk elements are next to boundary 2 (the top boundary)?
225  unsigned bound=2;
226  unsigned n_neigh = Bulk_mesh_pt->nboundary_element(bound);
227 
228  // Now loop over bulk elements and create the face elements
229  for(unsigned n=0;n<n_neigh;n++)
230  {
231  // Create the face element
232  FiniteElement *traction_element_pt
233  = new LinearElasticityTractionElement<ELEMENT>
234  (Bulk_mesh_pt->boundary_element_pt(bound,n),
235  Bulk_mesh_pt->face_index_at_boundary(bound,n));
236 
237  // Add to mesh
238  Surface_mesh_pt->add_element_pt(traction_element_pt);
239  }
240 
241 } // end of assign_traction_elements
242 
243 //==start_of_doc_solution=================================================
244 /// Doc the solution
245 //========================================================================
246 template<class ELEMENT>
248 {
249  ofstream some_file;
250  char filename[100];
251 
252  // Number of plot points
253  unsigned npts=5;
254 
255  // Output solution
256  sprintf(filename,"%s/soln.dat",doc_info.directory().c_str());
257  some_file.open(filename);
258  Bulk_mesh_pt->output(some_file,npts);
259  some_file.close();
260 
261  // Output exact solution
262  sprintf(filename,"%s/exact_soln.dat",doc_info.directory().c_str());
263  some_file.open(filename);
264  Bulk_mesh_pt->output_fct(some_file,npts,
266  some_file.close();
267 
268  // Doc error
269  double error=0.0;
270  double norm=0.0;
271  sprintf(filename,"%s/error.dat",doc_info.directory().c_str());
272  some_file.open(filename);
273  Bulk_mesh_pt->compute_error(some_file,
275  error,norm);
276  some_file.close();
277 
278 // Doc error norm:
279  cout << "\nNorm of error " << sqrt(error) << std::endl;
280  cout << "Norm of solution : " << sqrt(norm) << std::endl << std::endl;
281  cout << std::endl;
282 
283 
284 } // end_of_doc_solution
285 
286 
287 //===start_of_main======================================================
288 /// Driver code for PeriodicLoad linearly elastic problem
289 //======================================================================
290 int main(int argc, char* argv[])
291 {
292  // Number of elements in x-direction
293  unsigned nx=5;
294 
295  // Number of elements in y-direction (for (approximately) square elements)
296  unsigned ny=unsigned(double(nx)*Global_Parameters::Ly/Global_Parameters::Lx);
297 
298  // Set up doc info
299  DocInfo doc_info;
300 
301  // Set output directory
302  doc_info.set_directory("RESLT");
303 
304  // Set up problem
307 
308  // Solve
309  problem.newton_solve();
310 
311  // Output the solution
312  problem.doc_solution(doc_info);
313 
314 } // end_of_main
Periodic loading problem.
Mesh * Bulk_mesh_pt
Pointer to the bulk mesh.
Mesh * Surface_mesh_pt
Pointer to the mesh of traction elements.
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.
void actions_after_newton_solve()
Update after solve is empty.
void doc_solution(DocInfo &doc_info)
Doc the solution.
void assign_traction_elements()
Allocate traction elements on the top surface.
Namespace for global parameters.
void periodic_traction(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
The traction function.
double Amplitude
Amplitude of traction applied.
double Nu
Define Poisson coefficient Nu.
double Ly
Length of domain in y direction.
IsotropicElasticityTensor E(Nu)
The elasticity tensor.
bool Finite
Specify problem to be solved (boundary conditons for finite or infinite domain).
void exact_solution(const Vector< double > &x, Vector< double > &u)
The exact solution for infinite depth case.
double Lx
Length of domain in x direction.
int main(int argc, char *argv[])
Driver code for PeriodicLoad linearly elastic problem.