simple_spine_channel.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 2-D Channel with changing width, using Taylor Hood
31 // and Crouzeix Raviart elements.
32 
33 // Generic oomph-lib header
34 #include "generic.h"
35 
36 // Navier Stokes headers
37 #include "navier_stokes.h"
38 
39 // The mesh
40 #include "meshes/simple_rectangular_quadmesh.h"
41 
42 using namespace std;
43 
44 using namespace oomph;
45 
46 ///////////////////////////////////////////////////////////////////////
47 ///////////////////////////////////////////////////////////////////////
48 
49 
50 
51 //=========================================================================
52 /// Geometric object representing a wavy wall, parametrised by
53 /// \f[ x = \zeta \f]
54 /// \f[ y = A (1 - \cos(2\pi\zeta/L) )\f]
55 //=========================================================================
56 class WavyWall : public GeomObject
57 {
58 
59 public:
60 
61  /// \short Constructor: Pass wavelength and amplitude
62  WavyWall(const double& l, const double& amplitude)
63  : GeomObject(1,2), L(l), A(amplitude) {}
64 
65  /// Position vector to wavy wall.
66  void position(const Vector<double>& zeta, Vector<double>& r) const
67  {
68  r[0]=zeta[0];
69  r[1]=A*(1.0-cos(2.0*MathematicalConstants::Pi*zeta[0]/L));
70  }
71 
72 protected:
73 
74  /// Wavelength
75  double L;
76 
77  /// Amplitude
78  double A;
79 
80 
81 };
82 
83 
84 ///////////////////////////////////////////////////////////////////////
85 ///////////////////////////////////////////////////////////////////////
86 
87 
88 //======================================================================
89 /// Simple spine mesh class derived from standard 2D mesh.
90 /// Vertical lines of nodes are located on spines.
91 //======================================================================
92 template <class ELEMENT>
93 class SimpleSpineMesh : public SimpleRectangularQuadMesh<ELEMENT >,
94  public SpineMesh
95 {
96 
97 public:
98 
99  /// \short Constructor: Pass number of elements in x-direction,
100  /// number of elements in y-direction, length in x direction, initial
101  /// height of mesh, and pointer to timestepper (defaults to
102  /// Steady timestepper)
103  SimpleSpineMesh(const unsigned &nx,
104  const unsigned &ny,
105  const double &lx,
106  const double &h,
107  GeomObject* substrate_pt,
108  TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper);
109 
110 
111  /// \short Function pointer to function, h(x), that may be used
112  /// to specify the "height" of the domain, by assigning the function
113  /// values to the spine heights.
114  typedef double (*HeightFctPt)(const double& x);
115 
116  /// Access function: Pointer to height function
117  HeightFctPt& height_fct_pt()
118  {
119  return Height_fct_pt;
120  }
121 
122 
123  /// \short Height function -- this is called by update_spine_heights()
124  /// when spine heights are assigned.
125  double height_fct(const double& x)
126  {
127  // Resolve function pointer if non-NULL
128  if (Height_fct_pt==0)
129  {
130  return Default_height;
131  }
132  else
133  {
134  return Height_fct_pt(x);
135  }
136  }
137 
138 
139  /// \short Update the spine heights according to the function specified
140  /// by height_fct_pt().
142  {
143  unsigned n_spine=nspine();
144  for (unsigned i=0;i<n_spine;i++)
145  {
146  // Zero-th geometric parameter of the spine is the x-coordinate
147  // of its basis
148  double x=spine_pt(i)->geom_parameter(0);
149 
150  // Set spine height
151  spine_pt(i)->height()=height_fct(x);
152  }
153  }
154 
155 
156  /// \short General node update function implements pure virtual function
157  /// defined in SpineMesh base class and performs specific node update
158  /// actions: Nodes are located along vertical "spines" that emanate
159  /// from the "substrate" (the lower wall) specified by a
160  /// GeomObject.
161  virtual void spine_node_update(SpineNode* spine_node_pt)
162  {
163 
164  //Get fraction along the spine
165  double w = spine_node_pt->fraction();
166 
167  // Get height of spine
168  double h = spine_node_pt->spine_pt()->height();
169 
170  // Get x-position of spine origin (as vector)
171  Vector<double> x_spine_origin(1);
172  x_spine_origin[0]=spine_node_pt->spine_pt()->geom_parameter(0);
173 
174  // Get position vector to substrate (lower wall) for this spine
175  Vector<double> spine_base(2);
176  spine_node_pt->spine_pt()->geom_object_pt(0)->
177  position(x_spine_origin,spine_base);
178 
179  //Update the nodal position
180  spine_node_pt->x(0) = spine_base[0];
181  spine_node_pt->x(1) = spine_base[1]+w*h;
182  }
183 
184 
185 private:
186 
187  /// Default height
189 
190  /// Pointer to height function
191  HeightFctPt Height_fct_pt;
192 
193  /// Pointer to GeomObject that specifies the "substrate" (the lower wall)
194  GeomObject* Substrate_pt;
195 
196 };
197 
198 
199 
200 //===========================================================================
201 /// Constructor for spine 2D mesh: Pass number of elements in x-direction,
202 /// number of elements in y-direction, axial length and height of layer,
203 /// pointer to geometric object that specifies the substrate (the lower wall)
204 /// and pointer to timestepper (defaults to Static timestepper).
205 ///
206 /// The mesh contains a layer of spine-ified fluid elements (of type ELEMENT;
207 /// e.g SpineElement<QCrouzeixRaviartElement<2>)
208 //===========================================================================
209 template<class ELEMENT>
211  const unsigned &ny,
212  const double &lx,
213  const double &h,
214  GeomObject* substrate_pt,
215  TimeStepper* time_stepper_pt) :
216  SimpleRectangularQuadMesh<ELEMENT>(nx,ny,lx,h,time_stepper_pt),
217  Default_height(h), Height_fct_pt(0), Substrate_pt(substrate_pt)
218 
219 {
220 
221  // We've already called the constructor for the SimpleRectangularQuadMesh
222  // which sets up nodal positons/topology etc. Now attach the spine
223  // information
224 
225 
226  // Allocate memory for the spines and fractions along spines
227 
228  // Read out number of linear points in the element
229  unsigned np = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
230 
231  // Number of spines
232  unsigned nspine = (np-1)*nx+1;
233  Spine_pt.reserve(nspine);
234 
235  // Set up x-increments between spine origins
236  double x_lo=0.0;
237  double dx=lx/double(nx);
238 
239  //FIRST SPINE
240  //-----------
241 
242  //Element 0
243  //Node 0
244  //Assign the new spine with specified height
245  Spine* new_spine_pt=new Spine(h);
246  Spine_pt.push_back(new_spine_pt);
247 
248  // Get pointer to node
249  SpineNode* nod_pt=element_node_pt(0,0);
250 
251  //Pass the pointer to the spine to the node
252  nod_pt->spine_pt() = new_spine_pt;
253 
254  //Set the node's fraction along the spine
255  nod_pt->fraction() = 0.0;
256 
257  // Set the pointer to the mesh that implements the update fct for this node
258  nod_pt->spine_mesh_pt() = this;
259 
260  // Set update fct id (not really needed here...)
261  nod_pt->node_update_fct_id()=0;
262 
263 
264  // Create vector that stores additional geometric parameters that
265  // are required during the execution of the remesh function:
266  // For our remesh function, we need the x-coordinate of the
267  // spine origin.
268 
269  // Vector with enough storage for one geometic parameter
270  Vector<double> parameters(1);
271 
272  // Store the x-coordinate in it
273  parameters[0]=x_lo;
274 
275  // Pass the parameter to the spine
276  new_spine_pt->set_geom_parameter(parameters);
277 
278 
279 
280  // The remesh function involving this spine only involves
281  // a single geometric object: The substrate
282  Vector<GeomObject*> geom_object_pt(1);
283  geom_object_pt[0] = Substrate_pt;
284 
285  // Pass geom object(s) to spine
286  new_spine_pt->set_geom_object_pt(geom_object_pt);
287 
288 
289  //Loop vertically along the first spine
290 
291  //Loop over the elements
292  for(unsigned i=0;i<ny;i++)
293  {
294  //Loop over the vertical nodes, apart from the first
295  for(unsigned l1=1;l1<np;l1++)
296  {
297  // Get pointer to node
298  SpineNode* nod_pt=element_node_pt(i*nx,l1*np);
299 
300  //Pass the pointer to the spine to the node
301  nod_pt->spine_pt() = new_spine_pt;
302 
303  //Set the fraction
304  nod_pt->fraction()=(double(i)+double(l1)/double(np-1))/double(ny);
305 
306  // Pointer to the mesh that implements the update fct
307  nod_pt->spine_mesh_pt() = this;
308 
309  // Set update fct id (not really needed here)
310  nod_pt->node_update_fct_id()=0;
311 
312  }
313  } // end loop over elements
314 
315 
316  //LOOP OVER OTHER SPINES
317  //----------------------
318 
319  //Now loop over the elements horizontally
320  for(unsigned j=0;j<nx;j++)
321  {
322  //Loop over the nodes in the elements horizontally, ignoring
323  //the first column
324  unsigned npmax=np;
325  for(unsigned l2=1;l2<npmax;l2++)
326  {
327  //Assign the new spine with specified height
328  new_spine_pt=new Spine(h);
329 
330  // Add to collection of spines
331  Spine_pt.push_back(new_spine_pt);
332 
333  // Get the node
334  SpineNode* nod_pt=element_node_pt(j,l2);
335 
336  //Pass the pointer to the spine to the node
337  nod_pt->spine_pt() = new_spine_pt;
338 
339  //Set the fraction
340  nod_pt->fraction() = 0.0;
341 
342  // Pointer to the mesh that implements the update fct
343  nod_pt->spine_mesh_pt() = this;
344 
345  // Set update fct id (not really needed here)
346  nod_pt->node_update_fct_id()=0;
347 
348  // Create vector that stores additional geometric parameters that
349  // are required during the execution of the remesh function:
350  // For our remesh function, we need the x-coordinate of the
351  // spine origin.
352 
353  // Vector with enough storage for one geometic parameter
354  Vector<double> parameters(1);
355 
356  // Store the x-coordinate in it
357  parameters[0]=x_lo+double(j)*dx + double(l2)*dx/double(np-1);
358 
359  // Pass the parameter to the spine
360  new_spine_pt->set_geom_parameter(parameters);
361 
362 
363  // The remesh function involving this spine only involves
364  // a single geometric object: The substrate
365  Vector<GeomObject*> geom_object_pt(1);
366  geom_object_pt[0] = Substrate_pt;
367 
368  // Pass geom object(s) to spine
369  new_spine_pt->set_geom_object_pt(geom_object_pt);
370 
371  //Loop vertically along the spine
372  //Loop over the elements
373  for(unsigned i=0;i<ny;i++)
374  {
375  //Loop over the vertical nodes, apart from the first
376  for(unsigned l1=1;l1<np;l1++)
377  {
378  // Get the node
379  SpineNode* nod_pt=element_node_pt(i*nx+j,l1*np+l2);
380 
381  //Set the pointer to the spine
382  nod_pt->spine_pt() = new_spine_pt;
383 
384  //Set the fraction
385  nod_pt->fraction()=(double(i)+double(l1)/double(np-1))/double(ny);
386 
387  // Pointer to the mesh that implements the update fct
388  nod_pt->spine_mesh_pt() = this;
389 
390  // Set update fct id (not really needed here)
391  nod_pt->node_update_fct_id()=0;
392  }
393  }
394  }
395  }
396 
397 } // end of constructor
398 
399 
400 
401 
402 
403 ///////////////////////////////////////////////////////////////////////
404 ///////////////////////////////////////////////////////////////////////
405 
406 
407 
408 //==start_of_namespace===================================================
409 /// Namespace for physical parameters
410 //=======================================================================
411 namespace Global_Physical_Variables
412 {
413 
414  /// Reynolds number
415  double Re=100;
416 
417  /// Start of indented region
418  double X_indent_start=0.5;
419 
420  /// Length of indented region
421  double L=1.0;
422 
423  /// Total length of domain
424  double L_total=4.0;
425 
426  /// Undeformed height of domain
427  double H=1.0;
428 
429  /// Amplitude of indentation
430  double A=-0.6;
431 
432  /// Height of domain
433  double height(const double& x)
434  {
435  if ((x>X_indent_start)&&(x<(X_indent_start+L)))
436  {
437  return H + A*sin(MathematicalConstants::Pi*(x-X_indent_start)/L);
438  }
439  else
440  {
441  return H;
442  }
443  }
444 
445 } // end_of_namespace
446 
447 
448 
449 
450 
451 ///////////////////////////////////////////////////////////////////////
452 ///////////////////////////////////////////////////////////////////////
453 ///////////////////////////////////////////////////////////////////////
454 
455 
456 
457 //==start_of_problem_class============================================
458 /// Channel flow through a non-uniform channel whose geometry is defined
459 /// by a spine mesh.
460 //====================================================================
461 template<class ELEMENT>
462 class ChannelSpineFlowProblem : public Problem
463 {
464 
465 public:
466 
467  /// Constructor
469 
470  /// Destructor: (empty)
472 
473  /// \short Update the problem specs before solve. Update the nodal
474  /// positions
476  {
477  // Update the mesh
478  mesh_pt()->node_update();
479 
480  } // end_of_actions_before_newton_solve
481 
482  /// Update the after solve (empty)
484 
485  /// Doc the solution
486  void doc_solution(DocInfo& doc_info);
487 
488  /// Overload access to mesh
490  {
491  return dynamic_cast<SimpleSpineMesh<ELEMENT>*>(Problem::mesh_pt());
492  }
493 
494 
495 private:
496 
497  /// Width of channel
498  double Ly;
499 
500 
501 }; // end_of_problem_class
502 
503 
504 
505 //==start_of_constructor==================================================
506 /// Constructor for ChannelSpineFlow problem
507 //========================================================================
508 template<class ELEMENT>
510 {
511 
512  // Setup mesh
513 
514  // Domain length in x-direction
516 
517  // Domain length in y-direction
519 
520  // # of elements in x-direction
521  unsigned Nx=40;
522 
523  // # of elements in y-direction
524  unsigned Ny=10;
525 
526  // Substrate (lower wall): A wavy wall
527  GeomObject* substrate_pt=
529 
530  // Build and assign mesh -- pass pointer to geometric object
531  // that represents the sinusoidal bump on the upper wall
532  Problem::mesh_pt() = new SimpleSpineMesh<ELEMENT>(Nx,Ny,Lx,Ly,substrate_pt);
533 
534  // Set function pointer for height function
536 
537  // Update spine heights according to the function just specified
538  mesh_pt()->update_spine_heights();
539 
540  // Update nodal positions
541  mesh_pt()->node_update();
542 
543  // Pin all spine heights
544  unsigned nspine=mesh_pt()->nspine();
545  for (unsigned i=0;i<nspine;i++)
546  {
547  mesh_pt()->spine_pt(i)->spine_height_pt()->pin(0);
548  }
549 
550 
551  // Set the boundary conditions for this problem: All nodes are
552  // free by default -- just pin the ones that have Dirichlet conditions
553  // here: All boundaries are Dirichlet boundaries, except on boundary 1
554  unsigned num_bound = mesh_pt()->nboundary();
555  for(unsigned ibound=0;ibound<num_bound;ibound++)
556  {
557  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
558  for (unsigned inod=0;inod<num_nod;inod++)
559  {
560  if (ibound!=1)
561  {
562  // Loop over values (u and v velocities)
563  for (unsigned i=0;i<2;i++)
564  {
565  mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
566  }
567  }
568  else
569  {
570  // Parallel outflow ==> no-slip
571  mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
572  }
573  }
574  } // end loop over boundaries
575 
576 
577  // No slip on stationary upper and lower walls (boundaries 0 and 2)
578  // and parallel outflow (boundary 1)
579  for (unsigned ibound=0;ibound<num_bound-1;ibound++)
580  {
581  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
582  for (unsigned inod=0;inod<num_nod;inod++)
583  {
584  if (ibound!=1)
585  {
586  for (unsigned i=0;i<2;i++)
587  {
588  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0.0);
589  }
590  }
591  else
592  {
593  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
594  }
595  }
596  }
597 
598 
599  // Setup parabolic inflow along boundary 3:
600  unsigned ibound=3;
601  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
602  for (unsigned inod=0;inod<num_nod;inod++)
603  {
604  double y=mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
605  // Parallel, parabolic inflow
606  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,y*(Ly-y));
607  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
608  }
609 
610 
611  // Loop over the elements to set up element-specific
612  // things that cannot be handled by constructor: Pass
613  // pointer to Reynolds number
614  unsigned n_element = mesh_pt()->nelement();
615  for(unsigned e=0;e<n_element;e++)
616  {
617  // Upcast from GeneralisedElement to the present element
618  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
619  //Set the Reynolds number
620  el_pt->re_pt() = &Global_Physical_Variables::Re;
621  }
622 
623  // Setup equation numbering scheme
624  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
625 
626 } // end_of_constructor
627 
628 
629 
630 
631 //==start_of_doc_solution=================================================
632 /// Doc the solution
633 //========================================================================
634 template<class ELEMENT>
636 {
637 
638  ofstream some_file;
639  char filename[100];
640 
641  // Number of plot points
642  unsigned npts=5;
643 
644  // Output solution
645  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
646  doc_info.number());
647  some_file.open(filename);
648  mesh_pt()->output(some_file,npts);
649  some_file.close();
650 
651 } // end_of_doc_solution
652 
653 
654 
655 //==start_of_main======================================================
656 /// Driver for channel flow problem with spine mesh.
657 //=====================================================================
658 int main()
659 {
660 
661  // Set output directory
662  DocInfo doc_info;
663  doc_info.set_directory("RESLT");
664  doc_info.number()=0;
665 
666  // Solve problem with Taylor Hood elements
667  //---------------------------------------
668  {
669  //Build problem
671 
672  // Solve the problem
673  problem.newton_solve();
674 
675  //Output solution
676  problem.doc_solution(doc_info);
677 
678  // Step number
679  doc_info.number()++;
680 
681  } // end of Taylor Hood elements
682 
683 
684  // Solve problem with Crouzeix Raviart elements
685  //--------------------------------------------
686  {
687  // Build problem
689  problem;
690 
691  // Solve the problem with automatic adaptation
692  problem.newton_solve();
693 
694  //Output solution
695  problem.doc_solution(doc_info);
696  // Step number
697  doc_info.number()++;
698 
699  } // end of Crouzeix Raviart elements
700 
701 
702 
703 } // end_of_main
704 
705 
706 
707 
708 // //====================================================================
709 // /// Create the files to illustrate the sparse node update operation
710 // //====================================================================
711 // void doc_sparse_node_update()
712 // {
713 
714 // // Set output directory
715 // DocInfo doc_info;
716 // doc_info.set_directory("RESLT");
717 
718 // // Setup mesh
719 
720 // // # of elements in x-direction
721 // unsigned Nx=5;
722 
723 // // # of elements in y-direction
724 // unsigned Ny=5;
725 
726 // // Domain length in x-direction
727 // double Lx=5.0;
728 
729 // // Domain length in y-direction
730 // double Ly=1.0;
731 
732 // // Build and assign mesh
733 // SimpleSpineMesh<SpineElement<QTaylorHoodElement<2> > >* mesh_pt =
734 // new SimpleSpineMesh<SpineElement<QTaylorHoodElement<2> > >(Nx,Ny,Lx,Ly);
735 
736 // // Update *all* nodal positions
737 // mesh_pt->node_update();
738 
739 // unsigned count=0;
740 // ofstream some_file;
741 // char filename[100];
742 
743 // // Number of plot points
744 // unsigned npts=5;
745 
746 // // Output solution
747 // sprintf(filename,"%s/mesh_update%i.dat",doc_info.directory().c_str(),
748 // count);
749 // count++;
750 
751 // // Loop over spines
752 // unsigned n_node = mesh_pt->nnode();
753 // for (unsigned inode=0;inode<n_node;inode++)
754 // {
755 // SpineNode* node_pt=dynamic_cast<SpineNode*>(mesh_pt->node_pt(inode));
756 // node_pt->node_update();
757 // // Output solution
758 // some_file.open(filename);
759 // sprintf(filename,"%s/mesh_update%i.dat",doc_info.directory().c_str(),
760 // count);
761 // count++;
762 // mesh_pt->output(some_file,npts);
763 // some_file.close();
764 // }
765 // }
766 
767 
768 
769 
770 
771 
772 
773 
774 
775 
776 
777 
778 
779 
SimpleSpineMesh(const unsigned &nx, const unsigned &ny, const double &lx, const double &h, GeomObject *substrate_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in x-direction, number of elements in y-direction, length in x direction, initial height of mesh, and pointer to timestepper (defaults to Steady timestepper)
void doc_solution(DocInfo &doc_info)
Doc the solution.
double height(const double &x)
Height of domain.
virtual void spine_node_update(SpineNode *spine_node_pt)
General node update function implements pure virtual function defined in SpineMesh base class and per...
void actions_after_newton_solve()
Update the after solve (empty)
GeomObject * Substrate_pt
Pointer to GeomObject that specifies the "substrate" (the lower wall)
HeightFctPt & height_fct_pt()
Access function: Pointer to height function.
ChannelSpineFlowProblem()
Constructor.
double Default_height
Default height.
HeightFctPt Height_fct_pt
Pointer to height function.
double Ly
Width of channel.
double L
Length of indented region.
~ChannelSpineFlowProblem()
Destructor: (empty)
double A
Amplitude.
double Re
Reynolds number.
double height_fct(const double &x)
Height function – this is called by update_spine_heights() when spine heights are assigned...
double L_total
Total length of domain.
double H
Undeformed height of domain.
void update_spine_heights()
Update the spine heights according to the function specified by height_fct_pt().
double X_indent_start
Start of indented region.
void position(const Vector< double > &zeta, Vector< double > &r) const
Position vector to wavy wall.
int main()
Driver for channel flow problem with spine mesh.
WavyWall(const double &l, const double &amplitude)
Constructor: Pass wavelength and amplitude.
void actions_before_newton_solve()
Update the problem specs before solve. Update the nodal positions.
double A
Amplitude of indentation.
double L
Wavelength.
SimpleSpineMesh< ELEMENT > * mesh_pt()
Overload access to mesh.