prescribed_displ_lagr_mult2.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 Solid deformation -- driven by boundary motion which
31 // is imposed directly (without Lagrange multipliers)
32 
33 //Oomph-lib includes
34 #include "generic.h"
35 #include "solid.h"
36 #include "constitutive.h"
37 
38 //The mesh
39 #include "meshes/rectangular_quadmesh.h"
40 
41 using namespace std;
42 
43 using namespace oomph;
44 
45 
46 //======Start_of_warped_line===============================================
47 /// Warped line in 2D space
48 //=========================================================================
49 class WarpedLine : public GeomObject
50 {
51 
52 public:
53 
54  /// Constructor: Specify amplitude of deflection from straight horizontal line
55  WarpedLine(const double& ampl) : GeomObject(1,2)
56  {
57  Ampl=ampl;
58  }
59 
60  /// Broken copy constructor
61  WarpedLine(const WarpedLine& dummy)
62  {
63  BrokenCopy::broken_copy("WarpedLine");
64  }
65 
66  /// Broken assignment operator
67  void operator=(const WarpedLine&)
68  {
69  BrokenCopy::broken_assign("WarpedLine");
70  }
71 
72 
73  /// Empty Destructor
75 
76  /// \short Position vector at Lagrangian coordinate zeta
77  void position(const Vector<double>& zeta, Vector<double>& r) const
78  {
79  // Position vector
80  r[0] = zeta[0]+5.0*Ampl*zeta[0]*(zeta[0]-1.0)*(zeta[0]-0.7);
81  r[1] = 1.0+Ampl*0.5*(1.0-cos(2.0*MathematicalConstants::Pi*zeta[0]));
82  }
83 
84  /// \short Parametrised position on object: r(zeta). Evaluated at
85  /// previous timestep. t=0: current time; t>0: previous
86  /// timestep. Forward to steady version
87  void position(const unsigned& t, const Vector<double>& zeta,
88  Vector<double>& r) const
89  {
90  position(zeta,r);
91  }
92 
93  /// Access to amplitude
94  double& ampl() {return Ampl;}
95 
96 private:
97 
98  /// Amplitude of perturbation
99  double Ampl;
100 
101 };
102 
103 
104 
105 ///////////////////////////////////////////////////////////////////////
106 ///////////////////////////////////////////////////////////////////////
107 ///////////////////////////////////////////////////////////////////////
108 
109 
110 
111 //=======start_namespace==========================================
112 /// Global parameters
113 //================================================================
114 namespace Global_Physical_Variables
115 {
116 
117  /// GeomObject specifying the shape of the boundary: Initially it's flat.
119 
120  /// Poisson's ratio
121  double Nu=0.3;
122 
123  // Generalised Hookean constitutive equations
124  GeneralisedHookean Constitutive_law(&Global_Physical_Variables::Nu);
125 
126 } //end namespace
127 
128 
129 
130 //=============begin_problem============================================
131 /// Problem class for deformation of elastic block by prescribed
132 /// boundary motion.
133 //======================================================================
134 template<class ELEMENT>
135 class PrescribedBoundaryDisplacementProblem : public Problem
136 {
137 
138 public:
139 
140  /// Constructor:
142 
143  /// Update function (empty)
145 
146  /// Update boundary position directly
148  {
149 
150  // Loop over all nodes on top boundary (boundary 2)
151  unsigned b=2;
152  unsigned n_nod = solid_mesh_pt()->nboundary_node(b);
153  for(unsigned i=0;i<n_nod;i++)
154  {
155  Node* nod_pt= solid_mesh_pt()->boundary_node_pt(b,i);
156 
157  // Get boundary coordinate associated with boundary 2
158  Vector<double> zeta(1);
159  nod_pt->get_coordinates_on_boundary(b,zeta);
160 
161  // Get prescribed position from GeomObject
162  Vector<double> r(2);
164 
165  // Update position
166  nod_pt->x(0)=r[0];
167  nod_pt->x(1)=r[1];
168  }
169 
170  } // end actions_before_newton_solve
171 
172  /// Access function for the solid mesh
173  ElasticRefineableRectangularQuadMesh<ELEMENT>*& solid_mesh_pt()
174  {return Solid_mesh_pt;}
175 
176  /// Actions after adapt: Pin the redundant solid pressures (if any)
178  {
179  PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
180  solid_mesh_pt()->element_pt());
181  }
182 
183  /// Doc the solution
184  void doc_solution();
185 
186 private:
187 
188  /// Pointer to solid mesh
189  ElasticRefineableRectangularQuadMesh<ELEMENT>* Solid_mesh_pt;
190 
191  /// DocInfo object for output
192  DocInfo Doc_info;
193 
194 };
195 
196 
197 //===========start_of_constructor=======================================
198 /// Constructor:
199 //======================================================================
200 template<class ELEMENT>
202 {
203 
204  // Create the mesh
205 
206  // # of elements in x-direction
207  unsigned n_x=5;
208 
209  // # of elements in y-direction
210  unsigned n_y=5;
211 
212  // Domain length in x-direction
213  double l_x=1.0;
214 
215  // Domain length in y-direction
216  double l_y=1.0;
217 
218  //Now create the mesh
219  solid_mesh_pt() = new ElasticRefineableRectangularQuadMesh<ELEMENT>(
220  n_x,n_y,l_x,l_y);
221 
222  // Copy across
223  Problem::mesh_pt()=solid_mesh_pt();
224 
225  // Set error estimator
226  solid_mesh_pt()->spatial_error_estimator_pt()=new Z2ErrorEstimator;
227 
228  //Assign the physical properties to the elements before any refinement
229  //Loop over the elements in the main mesh
230  unsigned n_element =solid_mesh_pt()->nelement();
231  for(unsigned i=0;i<n_element;i++)
232  {
233  //Cast to a solid element
234  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(solid_mesh_pt()->element_pt(i));
235 
236  // Set the constitutive law
237  el_pt->constitutive_law_pt()=&Global_Physical_Variables::Constitutive_law;
238  }
239 
240  // Refine the mesh uniformly
241  solid_mesh_pt()->refine_uniformly();
242 
243 
244 
245  // Pin nodal positions on all boundaries including the top one
246  for (unsigned b=0;b<4;b++)
247  {
248  unsigned n_side = solid_mesh_pt()->nboundary_node(b);
249 
250  //Loop over the nodes
251  for(unsigned i=0;i<n_side;i++)
252  {
253  solid_mesh_pt()->boundary_node_pt(b,i)->pin_position(0);
254  solid_mesh_pt()->boundary_node_pt(b,i)->pin_position(1);
255  }
256  }
257 
258 
259 // Pin the redundant solid pressures (if any)
260  PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
261  solid_mesh_pt()->element_pt());
262 
263  // Setup equation numbering scheme
264  cout << "Number of dofs: " << assign_eqn_numbers() << std::endl;
265 
266  // Set output directory
267  Doc_info.set_directory("RESLT");
268 
269 } //end of constructor
270 
271 
272 
273 //==============start_doc===========================================
274 /// Doc the solution
275 //==================================================================
276 template<class ELEMENT>
278 {
279 
280  ofstream some_file;
281  char filename[100];
282 
283  // Number of plot points
284  unsigned n_plot = 5;
285 
286  // Output shape of deformed body
287  //------------------------------
288  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
289  Doc_info.number());
290  some_file.open(filename);
291  solid_mesh_pt()->output(some_file,n_plot);
292  some_file.close();
293 
294  // Increment label for output files
295  Doc_info.number()++;
296 
297 } //end doc
298 
299 
300 
301 //=======start_of_main==================================================
302 /// Driver code
303 //======================================================================
304 int main()
305 {
306 
307  //Set up the problem
309 
310  // Doc initial domain shape
311  problem.doc_solution();
312 
313  // Max. number of adaptations per solve
314  unsigned max_adapt=1;
315 
316  //Parameter incrementation
317  unsigned nstep=2; // 16;
318  for(unsigned i=0;i<nstep;i++)
319  {
320  // Increment imposed boundary displacement
322 
323  // Solve the problem with Newton's method, allowing
324  // up to max_adapt mesh adaptations after every solve.
325  problem.newton_solve(max_adapt);
326 
327  // Doc solution
328  problem.doc_solution();
329 
330  // For maximum stability: Reset the current nodal positions to be
331  // the "stress-free" ones -- this assignment means that the
332  // parameter study no longer corresponds to a physical experiment
333  // but is what we'd do if we wanted to use the solid solve
334  // to update a fluid mesh in an FSI problem, say.
335  problem.solid_mesh_pt()->set_lagrangian_nodal_coordinates();
336 
337  }
338 
339 } //end of main
340 
341 
342 
343 
344 
345 
346 
347 
void operator=(const WarpedLine &)
Broken assignment operator.
WarpedLine Boundary_geom_object(0.0)
GeomObject specifying the shape of the boundary: Initially it's flat.
void actions_after_newton_solve()
Update function (empty)
Warped line in 2D space.
void actions_before_newton_solve()
Update boundary position directly.
void position(const Vector< double > &zeta, Vector< double > &r) const
Position vector at Lagrangian coordinate zeta.
void actions_after_adapt()
Actions after adapt: Pin the redundant solid pressures (if any)
ElasticRefineableRectangularQuadMesh< ELEMENT > *& solid_mesh_pt()
Access function for the solid mesh.
int main()
Driver code.
WarpedLine(const double &ampl)
Constructor: Specify amplitude of deflection from straight horizontal line.
double & ampl()
Access to amplitude.
~WarpedLine()
Empty Destructor.
void position(const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
Parametrised position on object: r(zeta). Evaluated at previous timestep. t=0: current time; t>0: pre...
WarpedLine(const WarpedLine &dummy)
Broken copy constructor.