prescribed_displ_lagr_mult.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 via Lagrange multipliers
32 
33 
34 //Oomph-lib includes
35 #include "generic.h"
36 #include "solid.h"
37 #include "constitutive.h"
38 
39 //The mesh
40 #include "meshes/rectangular_quadmesh.h"
41 
42 using namespace std;
43 
44 using namespace oomph;
45 
46 
47 //================================================================
48 /// Function-type-object to compare finite elements based on
49 /// their x coordinate
50 //================================================================
52 {
53 
54 public:
55 
56  /// Comparison. Is x coordinate of el1_pt less than that of el2_pt?
57  bool operator()(FiniteElement* const& el1_pt, FiniteElement* const& el2_pt)
58  const
59  {
60  return el1_pt->node_pt(0)->x(0) < el2_pt->node_pt(0)->x(0);
61  }
62 
63 };
64 
65 
66 
67 //======Start_of_warped_line===============================================
68 /// Warped line in 2D space
69 //=========================================================================
70 class WarpedLine : public GeomObject
71 {
72 
73 public:
74 
75  /// Constructor: Specify amplitude of deflection from straight horizontal line
76  WarpedLine(const double& ampl) : GeomObject(1,2)
77  {
78  Ampl=ampl;
79  }
80 
81  /// Broken copy constructor
82  WarpedLine(const WarpedLine& dummy)
83  {
84  BrokenCopy::broken_copy("WarpedLine");
85  }
86 
87  /// Broken assignment operator
88  void operator=(const WarpedLine&)
89  {
90  BrokenCopy::broken_assign("WarpedLine");
91  }
92 
93 
94  /// Empty Destructor
96 
97  /// \short Position vector at Lagrangian coordinate zeta
98  void position(const Vector<double>& zeta, Vector<double>& r) const
99  {
100  // Position vector
101  r[0] = zeta[0]+5.0*Ampl*zeta[0]*(zeta[0]-1.0)*(zeta[0]-0.7);
102  r[1] = 1.0+Ampl*0.5*(1.0-cos(2.0*MathematicalConstants::Pi*zeta[0]));
103  }
104 
105  /// \short Parametrised position on object: r(zeta). Evaluated at
106  /// previous timestep. t=0: current time; t>0: previous
107  /// timestep. Forward to steady version
108  void position(const unsigned& t, const Vector<double>& zeta,
109  Vector<double>& r) const
110  {
111  position(zeta,r);
112  }
113 
114  /// Access to amplitude
115  double& ampl() {return Ampl;}
116 
117  /// \short How many items of Data does the shape of the object depend on?
118  /// None.
119  unsigned ngeom_data() const
120  {
121  return 0;
122  }
123 
124 private:
125 
126  /// Amplitude of perturbation
127  double Ampl;
128 
129 };
130 
131 
132 
133 ///////////////////////////////////////////////////////////////////////
134 ///////////////////////////////////////////////////////////////////////
135 ///////////////////////////////////////////////////////////////////////
136 
137 
138 //=======start_namespace==========================================
139 /// Global parameters
140 //================================================================
141 namespace Global_Physical_Variables
142 {
143 
144  /// GeomObject specifying the shape of the boundary: Initially it's flat.
146 
147  /// Poisson's ratio
148  double Nu=0.3;
149 
150  // Generalised Hookean constitutive equations
151  GeneralisedHookean Constitutive_law(&Global_Physical_Variables::Nu);
152 
153 } //end namespace
154 
155 
156 
157 //=============begin_problem============================================
158 /// Problem class for deformation of elastic block by prescribed
159 /// boundary motion.
160 //======================================================================
161 template<class ELEMENT>
163 {
164 
165 public:
166 
167  /// Constructor:
169 
170  /// Update function (empty)
172 
173  /// \short Update function (empty)
175 
176  /// Access function for the solid mesh
177  ElasticRefineableRectangularQuadMesh<ELEMENT>*& solid_mesh_pt()
178  {return Solid_mesh_pt;}
179 
180  /// Actions before adapt: Wipe the mesh of Lagrange multiplier elements
181  void actions_before_adapt();
182 
183  /// Actions after adapt: Rebuild the mesh of Lagrange multiplier elements
184  void actions_after_adapt();
185 
186  /// Doc the solution
187  void doc_solution();
188 
189 private:
190 
191  /// \short Create elements that enforce prescribed boundary motion
192  /// by Lagrange multiplilers
193  void create_lagrange_multiplier_elements();
194 
195  /// Delete elements that enforce prescribed boundary motion
196  /// by Lagrange multiplilers
197  void delete_lagrange_multiplier_elements();
198 
199  /// Pointer to solid mesh
200  ElasticRefineableRectangularQuadMesh<ELEMENT>* Solid_mesh_pt;
201 
202  /// Pointers to meshes of Lagrange multiplier elements
204 
205  /// DocInfo object for output
206  DocInfo Doc_info;
207 
208 };
209 
210 
211 //===========start_of_constructor=======================================
212 /// Constructor:
213 //======================================================================
214 template<class ELEMENT>
216 {
217 
218  // Create the mesh
219 
220  // # of elements in x-direction
221  unsigned n_x=5;
222 
223  // # of elements in y-direction
224  unsigned n_y=5;
225 
226  // Domain length in x-direction
227  double l_x=1.0;
228 
229  // Domain length in y-direction
230  double l_y=1.0;
231 
232  //Now create the mesh
233  solid_mesh_pt() = new ElasticRefineableRectangularQuadMesh<ELEMENT>(
234  n_x,n_y,l_x,l_y);
235 
236  // Set error estimator
237  solid_mesh_pt()->spatial_error_estimator_pt()=new Z2ErrorEstimator;
238 
239  //Assign the physical properties to the elements before any refinement
240  //Loop over the elements in the main mesh
241  unsigned n_element =solid_mesh_pt()->nelement();
242  for(unsigned i=0;i<n_element;i++)
243  {
244  //Cast to a solid element
245  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(solid_mesh_pt()->element_pt(i));
246 
247  // Set the constitutive law
248  el_pt->constitutive_law_pt()=&Global_Physical_Variables::Constitutive_law;
249  }
250 
251  // Refine the mesh uniformly
252  solid_mesh_pt()->refine_uniformly();
253 
254  // Construct the mesh of elements that enforce prescribed boundary motion
255  // by Lagrange multipliers
256  Lagrange_multiplier_mesh_pt=new SolidMesh;
257  create_lagrange_multiplier_elements();
258 
259  // Solid mesh is first sub-mesh
260  add_sub_mesh(solid_mesh_pt());
261 
262  // Add Lagrange multiplier sub-mesh
263  add_sub_mesh(Lagrange_multiplier_mesh_pt);
264 
265  // Build combined "global" mesh
266  build_global_mesh();
267 
268  // Pin nodal positions on all boundaries apart from the top one (2)
269  for (unsigned b=0;b<4;b++)
270  {
271  if (b!=2)
272  {
273  unsigned n_side = solid_mesh_pt()->nboundary_node(b);
274 
275  //Loop over the nodes
276  for(unsigned i=0;i<n_side;i++)
277  {
278  solid_mesh_pt()->boundary_node_pt(b,i)->pin_position(0);
279  solid_mesh_pt()->boundary_node_pt(b,i)->pin_position(1);
280  }
281  }
282  }
283 
284  // Pin the redundant solid pressures (if any)
285  PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
286  solid_mesh_pt()->element_pt());
287 
288  // Setup equation numbering scheme
289  cout << "Number of dofs: " << assign_eqn_numbers() << std::endl;
290 
291  // Set output directory
292  Doc_info.set_directory("RESLT");
293 
294 } //end of constructor
295 
296 
297 //=====================start_of_actions_before_adapt======================
298 /// Actions before adapt: Wipe the mesh of elements that impose
299 /// the prescribed boundary displacements
300 //========================================================================
301 template<class ELEMENT>
303 {
304  // Kill the elements and wipe surface mesh
305  delete_lagrange_multiplier_elements();
306 
307  // Rebuild the Problem's global mesh from its various sub-meshes
308  rebuild_global_mesh();
309 
310 }// end of actions_before_adapt
311 
312 
313 
314 //=====================start_of_actions_after_adapt=======================
315 /// Actions after adapt: Rebuild the mesh of elements that impose
316 /// the prescribed boundary displacements
317 //========================================================================
318 template<class ELEMENT>
320 {
321  // Create the elements that impose the displacement constraint
322  // and attach them to the bulk elements that are
323  // adjacent to boundary 2
324  create_lagrange_multiplier_elements();
325 
326  // Rebuild the Problem's global mesh from its various sub-meshes
327  rebuild_global_mesh();
328 
329  // Pin the redundant solid pressures (if any)
330  PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
331  solid_mesh_pt()->element_pt());
332 
333 }// end of actions_after_adapt
334 
335 
336 
337 //============start_of_create_lagrange_multiplier_elements===============
338 /// Create elements that impose the prescribed boundary displacement
339 //=======================================================================
340 template<class ELEMENT>
343 {
344  // Lagrange multiplier elements are located on boundary 2:
345  unsigned b=2;
346 
347  // How many bulk elements are adjacent to boundary b?
348  unsigned n_element = solid_mesh_pt()->nboundary_element(b);
349 
350  // Loop over the bulk elements adjacent to boundary b?
351  for(unsigned e=0;e<n_element;e++)
352  {
353  // Get pointer to the bulk element that is adjacent to boundary b
354  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
355  solid_mesh_pt()->boundary_element_pt(b,e));
356 
357  //Find the index of the face of element e along boundary b
358  int face_index = solid_mesh_pt()->face_index_at_boundary(b,e);
359 
360  // Create new element and add to mesh
361  Lagrange_multiplier_mesh_pt->add_element_pt(
362  new ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>(
363  bulk_elem_pt,face_index));
364  }
365 
366 
367  // Loop over the elements in the Lagrange multiplier element mesh
368  // for elements on the top boundary (boundary 2)
369  n_element=Lagrange_multiplier_mesh_pt->nelement();
370  for(unsigned i=0;i<n_element;i++)
371  {
372  //Cast to a Lagrange multiplier element
373  ImposeDisplacementByLagrangeMultiplierElement<ELEMENT> *el_pt =
374  dynamic_cast<ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>*>
375  (Lagrange_multiplier_mesh_pt->element_pt(i));
376 
377  // Set the GeomObject that defines the boundary shape and
378  // specify which bulk boundary we are attached to (needed to extract
379  // the boundary coordinate from the bulk nodes)
380  el_pt->set_boundary_shape_geom_object_pt(
382 
383  // Loop over the nodes
384  unsigned nnod=el_pt->nnode();
385  for (unsigned j=0;j<nnod;j++)
386  {
387  Node* nod_pt = el_pt->node_pt(j);
388 
389  // Is the node also on boundary 1 or 3?
390  if ((nod_pt->is_on_boundary(1))||(nod_pt->is_on_boundary(3)))
391  {
392  // How many nodal values were used by the "bulk" element
393  // that originally created this node?
394  unsigned n_bulk_value=el_pt->nbulk_value(j);
395 
396  // The remaining ones are Lagrange multipliers and we pin them.
397  unsigned nval=nod_pt->nvalue();
398  for (unsigned j=n_bulk_value;j<nval;j++)
399  {
400  nod_pt->pin(j);
401  }
402  }
403  }
404  }
405 
406 } // end of create_lagrange_multiplier_elements
407 
408 
409 
410 
411 //====start_of_delete_lagrange_multiplier_elements=======================
412 /// Delete elements that impose the prescribed boundary displacement
413 /// and wipe the associated mesh
414 //=======================================================================
415 template<class ELEMENT>
417 {
418  // How many surface elements are in the surface mesh
419  unsigned n_element = Lagrange_multiplier_mesh_pt->nelement();
420 
421  // Loop over the surface elements
422  for(unsigned e=0;e<n_element;e++)
423  {
424  // Kill surface element
425  delete Lagrange_multiplier_mesh_pt->element_pt(e);
426  }
427 
428  // Wipe the mesh
429  Lagrange_multiplier_mesh_pt->flush_element_and_node_storage();
430 
431 } // end of delete_lagrange_multiplier_elements
432 
433 
434 
435 //==============start_doc===========================================
436 /// Doc the solution
437 //==================================================================
438 template<class ELEMENT>
440 {
441 
442  ofstream some_file;
443  char filename[100];
444 
445  // Number of plot points
446  unsigned n_plot = 5;
447 
448 
449  // Output shape of deformed body
450  //------------------------------
451  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
452  Doc_info.number());
453  some_file.open(filename);
454  solid_mesh_pt()->output(some_file,n_plot);
455  some_file.close();
456 
457  // Output Lagrange multipliers
458  //----------------------------
459  sprintf(filename,"%s/lagr%i.dat",Doc_info.directory().c_str(),
460  Doc_info.number());
461  some_file.open(filename);
462 
463  // This makes sure the elements are ordered in same way every time
464  // the code is run -- necessary for validation tests.
465  std::vector<FiniteElement*> el_pt;
466  unsigned nelem=Lagrange_multiplier_mesh_pt->nelement();
467  for (unsigned e=0;e<nelem;e++)
468  {
469  el_pt.push_back(Lagrange_multiplier_mesh_pt->finite_element_pt(e));
470  }
471  std::sort(el_pt.begin(),el_pt.end(),FiniteElementComp());
472  for (unsigned e=0;e<nelem;e++)
473  {
474  el_pt[e]->output(some_file);
475  }
476  some_file.close();
477 
478  // Increment label for output files
479  Doc_info.number()++;
480 
481 } //end doc
482 
483 
484 
485 //=======start_of_main==================================================
486 /// Driver code
487 //======================================================================
488 int main()
489 {
490 
491  //Set up the problem
493 
494  // Doc initial domain shape
495  problem.doc_solution();
496 
497  // Max. number of adaptations per solve
498  unsigned max_adapt=1;
499 
500  //Parameter incrementation
501  unsigned nstep=2;
502  for(unsigned i=0;i<nstep;i++)
503  {
504  // Increment imposed boundary displacement
506 
507  // Solve the problem with Newton's method, allowing
508  // up to max_adapt mesh adaptations after every solve.
509  problem.newton_solve(max_adapt);
510 
511  // Doc solution
512  problem.doc_solution();
513 
514  // For maximum stability: Reset the current nodal positions to be
515  // the "stress-free" ones -- this assignment means that the
516  // parameter study no longer corresponds to a physical experiment
517  // but is what we'd do if we wanted to use the solid solve
518  // to update a fluid mesh in an FSI problem, say.
519  problem.solid_mesh_pt()->set_lagrangian_nodal_coordinates();
520 
521  }
522 
523 } //end of main
524 
525 
526 
527 
528 
529 
530 
531 
void operator=(const WarpedLine &)
Broken assignment operator.
double Ampl
Amplitude of perturbation.
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)
void actions_before_adapt()
Actions before adapt: Wipe the mesh of Lagrange multiplier elements.
Warped line in 2D space.
void actions_before_newton_solve()
Update function (empty)
int main()
Driver code.
void position(const Vector< double > &zeta, Vector< double > &r) const
Position vector at Lagrangian coordinate zeta.
unsigned ngeom_data() const
How many items of Data does the shape of the object depend on? None.
bool operator()(FiniteElement *const &el1_pt, FiniteElement *const &el2_pt) const
Comparison. Is x coordinate of el1_pt less than that of el2_pt?
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of Lagrange multiplier elements.
ElasticRefineableRectangularQuadMesh< ELEMENT > *& solid_mesh_pt()
Access function for the solid mesh.
WarpedLine(const double &ampl)
Constructor: Specify amplitude of deflection from straight horizontal line.
double & ampl()
Access to amplitude.
void create_lagrange_multiplier_elements()
Create elements that enforce prescribed boundary motion by Lagrange multiplilers. ...
~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...
ElasticRefineableRectangularQuadMesh< ELEMENT > * Solid_mesh_pt
Pointer to solid mesh.
WarpedLine(const WarpedLine &dummy)
Broken copy constructor.
SolidMesh * Lagrange_multiplier_mesh_pt
Pointers to meshes of Lagrange multiplier elements.
DocInfo Doc_info
DocInfo object for output.