prescribed_displ_lagr_mult_precond.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 // The preconditioner
43 #include "multi_physics/pseudo_elastic_preconditioner.h"
44 
45 
46 using namespace std;
47 
48 using namespace oomph;
49 
50 namespace oomph {
51 
52 //=============start_wrapper==================================================
53 /// Pseudo-Elastic Solid element class to overload the block preconditioner
54 /// methods ndof_types() and get_dof_numbers_for_unknowns() to differentiate
55 /// between DOFs subject to Lagrange multiplier and those that are not.
56 //============================================================================
57 template <class ELEMENT>
59  public virtual ELEMENT
60 {
61 
62 public:
63 
64  /// Constructor
65  PseudoElasticBulkElement() : ELEMENT() {}
66 
67  /// \short Returns the number of DOF types associated with this element: Twice
68  /// the number of spatial dimensions (for the constrained and
69  /// unconstrained nodal positions).
70  unsigned ndof_types() const
71  {
72  return 2*ELEMENT::dim();
73  }
74 
75  /// \short Create a list of pairs for all unknowns in this element,
76  /// so that the first entry in each pair contains the global equation
77  /// number of the unknown, while the second one contains the number
78  /// of the "DOF" that this unknown is associated with.
79  /// (Function can obviously only be called if the equation numbering
80  /// scheme has been set up.)\n
81  /// E.g. in a 3D problem there are 6 types of DOF:\n
82  /// 0 - x displacement (without lagr mult traction)\n
83  /// 1 - y displacement (without lagr mult traction)\n
84  /// 2 - z displacement (without lagr mult traction)\n
85  /// 4 - x displacement (with lagr mult traction)\n
86  /// 5 - y displacement (with lagr mult traction)\n
87  /// 6 - z displacement (with lagr mult traction)\n
89  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
90  {
91  // temporary pair (used to store dof lookup prior to being added to list
92  std::pair<unsigned,unsigned> dof_lookup;
93 
94  // number of nodes
95  const unsigned n_node = this->nnode();
96 
97  //Get the number of position dofs and dimensions at the node
98  const unsigned n_position_type = ELEMENT::nnodal_position_type();
99  const unsigned nodal_dim = ELEMENT::nodal_dimension();
100 
101  //Integer storage for local unknown
102  int local_unknown=0;
103 
104  //Loop over the nodes
105  for(unsigned n=0;n<n_node;n++)
106  {
107  unsigned offset = 0;
108  if (this->node_pt(n)->nvalue() != this->required_nvalue(n))
109  {
110  offset = ELEMENT::dim();
111  }
112 
113  //Loop over position dofs
114  for(unsigned k=0;k<n_position_type;k++)
115  {
116  //Loop over dimension
117  for(unsigned i=0;i<nodal_dim;i++)
118  {
119  //If the variable is unpinned
120  local_unknown = ELEMENT::position_local_eqn(n,k,i);
121  if (local_unknown >= 0)
122  {
123  // store dof lookup in temporary pair: First entry in pair
124  // is global equation number; second entry is dof type
125  dof_lookup.first = this->eqn_number(local_unknown);
126  dof_lookup.second = offset+i;
127 
128  // add to list
129  dof_lookup_list.push_front(dof_lookup);
130  }
131  }
132  }
133  }
134  }
135 };
136 
137 
138 
139 
140 
141 //===========start_face_geometry==============================================
142 /// FaceGeometry of wrapped element is the same as the underlying element
143 //============================================================================
144 template<class ELEMENT>
145 class FaceGeometry<PseudoElasticBulkElement<ELEMENT> > :
146  public virtual FaceGeometry<ELEMENT>
147 {
148 
149 public:
150 
151  /// Constructor -- required for more recent versions of gcc
153 
154 };
155 
156 
157 
158 }
159 
160 
161 
162 
163 //////////////////////////////////////////////////////////////////////
164 //////////////////////////////////////////////////////////////////////
165 //////////////////////////////////////////////////////////////////////
166 
167 
168 
169 //================================================================
170 /// Namespace to "hide" global function that creates an instance
171 /// of oomph-lib's diagonal preconditioner.
172 //================================================================
173 namespace DiagonalPreconditionerHelper
174 {
175 
176 /// \short Create a matrix-based diagonal preconditioner for
177 /// subsidiary linear systems
178  Preconditioner* get_diagonal_preconditioner()
179  {
180  MatrixBasedDiagPreconditioner* prec_pt=
181  new MatrixBasedDiagPreconditioner;
182  return prec_pt;
183  }
184 
185 }
186 
187 
188 //////////////////////////////////////////////////////////////////////
189 //////////////////////////////////////////////////////////////////////
190 //////////////////////////////////////////////////////////////////////
191 
192 
193 
194 //================================================================
195 /// Function-type-object to compare finite elements based on
196 /// their x coordinate
197 //================================================================
198 class FiniteElementComp
199 {
200 
201 public:
202 
203  /// Comparison. Is x coordinate of el1_pt less than that of el2_pt?
204  bool operator()(FiniteElement* const& el1_pt, FiniteElement* const& el2_pt)
205  const
206  {
207  return el1_pt->node_pt(0)->x(0) < el2_pt->node_pt(0)->x(0);
208  }
209 
210 };
211 
212 
213 
214 //======Start_of_warped_line===============================================
215 /// Warped line in 2D space
216 //=========================================================================
217 class WarpedLine : public GeomObject
218 {
219 
220 public:
221 
222  /// Constructor: Specify amplitude of deflection from straight horizontal line
223  WarpedLine(const double& ampl) : GeomObject(1,2)
224  {
225  Ampl=ampl;
226  }
227 
228  /// Broken copy constructor
229  WarpedLine(const WarpedLine& dummy)
230  {
231  BrokenCopy::broken_copy("WarpedLine");
232  }
233 
234  /// Broken assignment operator
235  void operator=(const WarpedLine&)
236  {
237  BrokenCopy::broken_assign("WarpedLine");
238  }
239 
240 
241  /// Empty Destructor
243 
244  /// \short Position vector at Lagrangian coordinate zeta
245  void position(const Vector<double>& zeta, Vector<double>& r) const
246  {
247  // Position vector
248  r[0] = zeta[0]+5.0*Ampl*zeta[0]*(zeta[0]-1.0)*(zeta[0]-0.7);
249  r[1] = 1.0+Ampl*0.5*(1.0-cos(2.0*MathematicalConstants::Pi*zeta[0]));
250  }
251 
252  /// \short Parametrised position on object: r(zeta). Evaluated at
253  /// previous timestep. t=0: current time; t>0: previous
254  /// timestep. Forward to steady version
255  void position(const unsigned& t, const Vector<double>& zeta,
256  Vector<double>& r) const
257  {
258  position(zeta,r);
259  }
260 
261  /// Access to amplitude
262  double& ampl() {return Ampl;}
263 
264  /// \short How many items of Data does the shape of the object depend on?
265  /// None.
266  unsigned ngeom_data() const
267  {
268  return 0;
269  }
270 
271 private:
272 
273  /// Amplitude of perturbation
274  double Ampl;
275 
276 };
277 
278 
279 
280 ///////////////////////////////////////////////////////////////////////
281 ///////////////////////////////////////////////////////////////////////
282 ///////////////////////////////////////////////////////////////////////
283 
284 
285 //=======start_namespace==========================================
286 /// Global parameters
287 //================================================================
288 namespace Global_Physical_Variables
289 {
290 
291  /// GeomObject specifying the shape of the boundary: Initially it's flat.
293 
294  /// Poisson's ratio
295  double Nu=0.3;
296 
297  // Generalised Hookean constitutive equations
298  GeneralisedHookean Constitutive_law(&Global_Physical_Variables::Nu);
299 
300 } //end namespace
301 
302 
303 
304 //=============begin_problem============================================
305 /// Problem class for deformation of elastic DOF type by prescribed
306 /// boundary motion.
307 //======================================================================
308 template<class ELEMENT>
309 class PrescribedBoundaryDisplacementProblem : public Problem
310 {
311 
312 public:
313 
314  /// Constructor: Pass in number of elements along axes
315  PrescribedBoundaryDisplacementProblem(const unsigned& nel_1d);
316 
317  /// Update function (empty)
319 
320  /// \short Update function (empty)
322 
323  /// Access function for the solid mesh
324  ElasticRefineableRectangularQuadMesh<ELEMENT>*& solid_mesh_pt()
325  {return Solid_mesh_pt;}
326 
327  /// Actions before adapt: Wipe the mesh of Lagrange multiplier elements
328  void actions_before_adapt();
329 
330  /// Actions after adapt: Rebuild the mesh of Lagrange multiplier elements
331  void actions_after_adapt();
332 
333  /// Doc the solution
334  void doc_solution();
335 
336 private:
337 
338  /// \short Create elements that enforce prescribed boundary motion
339  /// by Lagrange multiplilers
340  void create_lagrange_multiplier_elements();
341 
342  /// Delete elements that enforce prescribed boundary motion
343  /// by Lagrange multiplilers
344  void delete_lagrange_multiplier_elements();
345 
346  /// Pointer to solid mesh
347  ElasticRefineableRectangularQuadMesh<ELEMENT>* Solid_mesh_pt;
348 
349  /// Pointers to meshes of Lagrange multiplier elements
350  SolidMesh* Lagrange_multiplier_mesh_pt;
351 
352  /// DocInfo object for output
353  DocInfo Doc_info;
354 
355 };
356 
357 
358 //===========start_of_constructor=======================================
359 /// Constructor: Pass in number of elements along axes
360 //======================================================================
361 template<class ELEMENT>
364 {
365 
366  // Create the mesh
367 
368  // # of elements in x-direction
369  unsigned n_x=nel_1d;
370 
371  // # of elements in y-direction
372  unsigned n_y=nel_1d;
373 
374  // Domain length in x-direction
375  double l_x=1.0;
376 
377  // Domain length in y-direction
378  double l_y=1.0;
379 
380  //Now create the mesh
381  solid_mesh_pt() = new ElasticRefineableRectangularQuadMesh<ELEMENT>(
382  n_x,n_y,l_x,l_y);
383 
384  // Set error estimator
385  solid_mesh_pt()->spatial_error_estimator_pt()=new Z2ErrorEstimator;
386 
387  //Assign the physical properties to the elements before any refinement
388  //Loop over the elements in the main mesh
389  unsigned n_element =solid_mesh_pt()->nelement();
390  for(unsigned i=0;i<n_element;i++)
391  {
392  //Cast to a solid element
393  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(solid_mesh_pt()->element_pt(i));
394 
395  // Set the constitutive law
396  el_pt->constitutive_law_pt()=&Global_Physical_Variables::Constitutive_law;
397  }
398 
399  // Refine the mesh uniformly
400  solid_mesh_pt()->refine_uniformly();
401 
402  // Construct the mesh of elements that enforce prescribed boundary motion
403  // by Lagrange multipliers
404  Lagrange_multiplier_mesh_pt=new SolidMesh;
405  create_lagrange_multiplier_elements();
406 
407  // Solid mesh is first sub-mesh
408  add_sub_mesh(solid_mesh_pt());
409 
410  // Add Lagrange multiplier sub-mesh
411  add_sub_mesh(Lagrange_multiplier_mesh_pt);
412 
413  // Build combined "global" mesh
414  build_global_mesh();
415 
416  // Pin nodal positions on all boundaries apart from the top one (2)
417  for (unsigned b=0;b<4;b++)
418  {
419  if (b!=2)
420  {
421  unsigned n_side = solid_mesh_pt()->nboundary_node(b);
422 
423  //Loop over the nodes
424  for(unsigned i=0;i<n_side;i++)
425  {
426  solid_mesh_pt()->boundary_node_pt(b,i)->pin_position(0);
427  solid_mesh_pt()->boundary_node_pt(b,i)->pin_position(1);
428  }
429  }
430  }
431 
432  // Pin the redundant solid pressures (if any)
433  PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
434  solid_mesh_pt()->element_pt());
435 
436  // Setup equation numbering scheme
437  cout << "Number of dofs: " << assign_eqn_numbers() << std::endl;
438 
439  // Set output directory
440  Doc_info.set_directory("RESLT");
441 
442  // Create the linear solver
443  IterativeLinearSolver* solver_pt=0;
444 
445  // If we have trilinos, use it
446 #ifdef OOMPH_HAS_TRILINOS
447 
448  // Create solver
449  solver_pt = new TrilinosAztecOOSolver;
450 
451  // Use GMRES
452  dynamic_cast<TrilinosAztecOOSolver*>(solver_pt)->solver_type()
453  = TrilinosAztecOOSolver::GMRES;
454 
455 #else
456 
457  // Use oomph-lib's own GMRES
458  solver_pt = new GMRES<CRDoubleMatrix>;
459 
460 #endif
461 
462  // Set solver
463  linear_solver_pt() = solver_pt;
464 
465  // Create the preconditioner
466  PseudoElasticPreconditioner * prec_pt = new PseudoElasticPreconditioner;
467 
468  // Set solid and lagrange multiplier meshes
469  prec_pt->set_elastic_mesh(Solid_mesh_pt);
470  prec_pt->set_lagrange_multiplier_mesh(Lagrange_multiplier_mesh_pt);
471 
472  // Set the preconditioner
473  solver_pt->preconditioner_pt() = prec_pt;
474 
475  // Use upper block triangular preconditioner for elastic block?
476  if (CommandLineArgs::command_line_flag_has_been_set
477  ("--block_upper_for_elastic_block"))
478  {
479  prec_pt->elastic_preconditioner_type()=
480  PseudoElasticPreconditioner::Block_upper_triangular_preconditioner;
481  }
482 
483 #ifdef OOMPH_HAS_HYPRE
484 
485  // Use Hypre as subsidiary preconditioner (inexact solver) for
486  // linear (sub-)systems to be solved in the elastic block?
487  if (CommandLineArgs::command_line_flag_has_been_set
488  ("--hypre_for_elastic_blocks"))
489  {
490  prec_pt->set_elastic_subsidiary_preconditioner
491  (Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper::
492  get_elastic_preconditioner_hypre);
493  }
494 
495 #endif
496 
497 
498 #ifdef OOMPH_HAS_TRILINOS
499 
500  // Use Trilinos CG as subsidiary preconditioner (inexact solver) for
501  // linear (sub-)systems to be solved in the Lagrange multiplier block?
502  if (CommandLineArgs::command_line_flag_has_been_set
503  ("--trilinos_cg_for_lagrange_multiplier_blocks"))
504  {
505  prec_pt->set_lagrange_multiplier_subsidiary_preconditioner
506  (Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper::
507  get_lagrange_multiplier_preconditioner);
508  }
509 
510 #endif
511 
512  // Use diagonal scaling as subsidiary preconditioner (inexact solver) for
513  // linear (sub-)systems to be solved in the elastic block?
514  if (CommandLineArgs::command_line_flag_has_been_set
515  ("--diagonal_scaling_for_elastic_blocks"))
516  {
517  prec_pt->set_elastic_subsidiary_preconditioner(
519  }
520 
521 } //end of constructor
522 
523 
524 //=====================start_of_actions_before_adapt======================
525 /// Actions before adapt: Wipe the mesh of elements that impose
526 /// the prescribed boundary displacements
527 //========================================================================
528 template<class ELEMENT>
530 {
531  // Kill the elements and wipe surface mesh
532  delete_lagrange_multiplier_elements();
533 
534  // Rebuild the Problem's global mesh from its various sub-meshes
535  rebuild_global_mesh();
536 
537 }// end of actions_before_adapt
538 
539 
540 
541 //=====================start_of_actions_after_adapt=======================
542 /// Actions after adapt: Rebuild the mesh of elements that impose
543 /// the prescribed boundary displacements
544 //========================================================================
545 template<class ELEMENT>
547 {
548  // Create the elements that impose the displacement constraint
549  // and attach them to the bulk elements that are
550  // adjacent to boundary 2
551  create_lagrange_multiplier_elements();
552 
553  // Rebuild the Problem's global mesh from its various sub-meshes
554  rebuild_global_mesh();
555 
556  // Pin the redundant solid pressures (if any)
557  PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
558  solid_mesh_pt()->element_pt());
559 
560 }// end of actions_after_adapt
561 
562 
563 
564 //============start_of_create_lagrange_multiplier_elements===============
565 /// Create elements that impose the prescribed boundary displacement
566 //=======================================================================
567 template<class ELEMENT>
570 {
571  // Lagrange multiplier elements are located on boundary 2:
572  unsigned b=2;
573 
574  // How many bulk elements are adjacent to boundary b?
575  unsigned n_element = solid_mesh_pt()->nboundary_element(b);
576 
577  // Loop over the bulk elements adjacent to boundary b?
578  for(unsigned e=0;e<n_element;e++)
579  {
580  // Get pointer to the bulk element that is adjacent to boundary b
581  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
582  solid_mesh_pt()->boundary_element_pt(b,e));
583 
584  //Find the index of the face of element e along boundary b
585  int face_index = solid_mesh_pt()->face_index_at_boundary(b,e);
586 
587  // Create new element and add to mesh
588  Lagrange_multiplier_mesh_pt->add_element_pt(
589  new ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>(
590  bulk_elem_pt,face_index));
591  }
592 
593 
594  // Loop over the elements in the Lagrange multiplier element mesh
595  // for elements on the top boundary (boundary 2)
596  n_element=Lagrange_multiplier_mesh_pt->nelement();
597  for(unsigned i=0;i<n_element;i++)
598  {
599  //Cast to a Lagrange multiplier element
600  ImposeDisplacementByLagrangeMultiplierElement<ELEMENT> *el_pt =
601  dynamic_cast<ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>*>
602  (Lagrange_multiplier_mesh_pt->element_pt(i));
603 
604  // Set the GeomObject that defines the boundary shape and
605  // specify which bulk boundary we are attached to (needed to extract
606  // the boundary coordinate from the bulk nodes)
607  el_pt->set_boundary_shape_geom_object_pt(
609 
610  // Loop over the nodes
611  unsigned nnod=el_pt->nnode();
612  for (unsigned j=0;j<nnod;j++)
613  {
614  Node* nod_pt = el_pt->node_pt(j);
615 
616  // Is the node also on boundary 1 or 3?
617  if ((nod_pt->is_on_boundary(1))||(nod_pt->is_on_boundary(3)))
618  {
619  // How many nodal values were used by the "bulk" element
620  // that originally created this node?
621  unsigned n_bulk_value=el_pt->nbulk_value(j);
622 
623  // The remaining ones are Lagrange multipliers and we pin them.
624  unsigned nval=nod_pt->nvalue();
625  for (unsigned j=n_bulk_value;j<nval;j++)
626  {
627  nod_pt->pin(j);
628  }
629  }
630  }
631  }
632 
633 } // end of create_lagrange_multiplier_elements
634 
635 
636 
637 
638 //====start_of_delete_lagrange_multiplier_elements=======================
639 /// Delete elements that impose the prescribed boundary displacement
640 /// and wipe the associated mesh
641 //=======================================================================
642 template<class ELEMENT>
644 {
645  // How many surface elements are in the surface mesh
646  unsigned n_element = Lagrange_multiplier_mesh_pt->nelement();
647 
648  // Loop over the surface elements
649  for(unsigned e=0;e<n_element;e++)
650  {
651  // Kill surface element
652  delete Lagrange_multiplier_mesh_pt->element_pt(e);
653  }
654 
655  // Wipe the mesh
656  Lagrange_multiplier_mesh_pt->flush_element_and_node_storage();
657 
658 } // end of delete_lagrange_multiplier_elements
659 
660 
661 
662 //==============start_doc===========================================
663 /// Doc the solution
664 //==================================================================
665 template<class ELEMENT>
667 {
668 
669  ofstream some_file;
670  char filename[100];
671 
672  // Number of plot points
673  unsigned n_plot = 5;
674 
675 
676  // Output shape of deformed body
677  //------------------------------
678  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
679  Doc_info.number());
680  some_file.open(filename);
681  solid_mesh_pt()->output(some_file,n_plot);
682  some_file.close();
683 
684  // Output Lagrange multipliers
685  //----------------------------
686  sprintf(filename,"%s/lagr%i.dat",Doc_info.directory().c_str(),
687  Doc_info.number());
688  some_file.open(filename);
689 
690  // This makes sure the elements are ordered in same way every time
691  // the code is run -- necessary for validation tests.
692  std::vector<FiniteElement*> el_pt;
693  unsigned nelem=Lagrange_multiplier_mesh_pt->nelement();
694  for (unsigned e=0;e<nelem;e++)
695  {
696  el_pt.push_back(Lagrange_multiplier_mesh_pt->finite_element_pt(e));
697  }
698  std::sort(el_pt.begin(),el_pt.end(),FiniteElementComp());
699  for (unsigned e=0;e<nelem;e++)
700  {
701  el_pt[e]->output(some_file);
702  }
703  some_file.close();
704 
705  // Increment label for output files
706  Doc_info.number()++;
707 
708 } //end doc
709 
710 
711 
712 //=======start_of_main==================================================
713 /// Driver code
714 //======================================================================
715 int main(int argc, char* argv[])
716 {
717 
718 
719 #ifdef OOMPH_HAS_MPI
720 
721  // Start up mpi if oomph-lib has been compiled with parallel support
722  // because parallel versions of trilinos and hypre which are then called
723  // by default) need it.
724  MPI_Helpers::init(argc,argv,false);
725 
726 #endif
727 
728  // Store command line arguments
729  CommandLineArgs::setup(argc,argv);
730 
731  // Define possible command line arguments and parse the ones that
732  // were actually specified
733 
734  // Number of elements along axes
735  unsigned nel_1d=5;
736  CommandLineArgs::specify_command_line_flag("--nel_1d",&nel_1d);
737 
738  // Suppress adaptation (for study of iteration count vs uniform mesh
739  // refinement, say)
740  CommandLineArgs::specify_command_line_flag("--no_adapt");
741 
742  // Use block upper triangular preconditioner for elasticity block
743  CommandLineArgs::specify_command_line_flag("--block_upper_for_elastic_block");
744 
745  // Use Hypre as subsidiary preconditioner (inexact solver) for
746  // linear (sub-)systems to be solved in the elastic block?
747  CommandLineArgs::specify_command_line_flag("--hypre_for_elastic_blocks");
748 
749  // Use Trilinos CG as subsidiary preconditioner (inexact solver) for
750  // linear (sub-)systems to be solved in the Lagrange multiplier block?
751  CommandLineArgs::specify_command_line_flag
752  ("--trilinos_cg_for_lagrange_multiplier_blocks");
753 
754  // Use diagonal scaling as subsidiary preconditioner (inexact solver) for
755  // linear (sub-)systems to be solved in the elastic block?
756  // [only used for exercise]
757  CommandLineArgs::specify_command_line_flag
758  ("--diagonal_scaling_for_elastic_blocks");
759 
760  // Parse command line
761  CommandLineArgs::parse_and_assign();
762 
763  // Doc what has actually been specified on the command line
764  CommandLineArgs::doc_specified_flags();
765 
766  //Set up the problem with specified number of elements along axes.
769 
770  // Doc initial domain shape
771  problem.doc_solution();
772 
773  // Max. number of adaptations per solve
774  unsigned max_adapt=1;
775  if (CommandLineArgs::command_line_flag_has_been_set("--no_adapt"))
776  {
777  max_adapt=0;
778  }
779 
780  //Parameter incrementation
781  unsigned nstep=2;
782  for(unsigned i=0;i<nstep;i++)
783  {
784  // Increment imposed boundary displacement
786 
787  // Solve the problem with Newton's method, allowing
788  // up to max_adapt mesh adaptations after every solve.
789  problem.newton_solve(max_adapt);
790 
791  // Doc solution
792  problem.doc_solution();
793 
794  // For maximum stability: Reset the current nodal positions to be
795  // the "stress-free" ones -- this assignment means that the
796  // parameter study no longer corresponds to a physical experiment
797  // but is what we'd do if we wanted to use the solid solve
798  // to update a fluid mesh in an FSI problem, say.
799  problem.solid_mesh_pt()->set_lagrangian_nodal_coordinates();
800 
801  }
802 
803 } //end of main
804 
805 
806 
807 
808 
809 
810 
811 
void operator=(const WarpedLine &)
Broken assignment operator.
WarpedLine Boundary_geom_object(0.0)
GeomObject specifying the shape of the boundary: Initially it's flat.
FaceGeometry()
Constructor – required for more recent versions of gcc.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of Lagrange multiplier elements.
Warped line in 2D space.
void position(const Vector< double > &zeta, Vector< double > &r) const
Position vector at Lagrangian coordinate zeta.
int main(int argc, char *argv[])
Driver code.
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.
Preconditioner * get_diagonal_preconditioner()
Create a matrix-based diagonal preconditioner for subsidiary linear systems.
unsigned ndof_types() const
Returns the number of DOF types associated with this element: Twice the number of spatial dimensions ...
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...
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
WarpedLine(const WarpedLine &dummy)
Broken copy constructor.