spine_channel2.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$
7 //LIC//
8 //LIC// $LastChangedDate$
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/channel_spine_mesh.h"
41 
42 using namespace std;
43 
44 using namespace oomph;
45 
46 //==start_of_namespace===================================================
47 /// Namespace for physical parameters
48 //=======================================================================
49 namespace Global_Physical_Variables2
50 {
51  /// Reynolds number
52  double Re=100;
53 } // end_of_namespace
54 
55 
56 ///////////////////////////////////////////////////////////////////////
57 ///////////////////////////////////////////////////////////////////////
58 // Deflected line as geometric object
59 ///////////////////////////////////////////////////////////////////////
60 ///////////////////////////////////////////////////////////////////////
61 
62 
63 
64 //=========================================================================
65 /// Steady, spiked 1D line in 2D space
66 /// \f[ x = \zeta \f]
67 /// \f[ y = \left\{
68 /// \begin{array}{cl}
69 /// H + 2A\left(\frac{\zeta - \zeta_{\mbox{min}}}
70 /// {\zeta_{\mbox{max}} - \zeta_{\mbox{min}}}\right) &
71 /// \mbox{for } \zeta \leq \frac{1}{2}
72 /// \left(\zeta_{\mbox{max}} + \zeta_{\mbox{min}}\right)\\H +
73 /// 2A\left(\frac{\zeta - \zeta_{\mbox{max}}}
74 /// {\zeta_{\mbox{min}} - \zeta_{\mbox{max}}}\right) &
75 /// \mbox{for } \zeta > \frac{1}{2}
76 /// \left(\zeta_{\mbox{max}}
77 /// + \zeta_{\mbox{min}}\right)\\ \end{array}
78 /// \right.\f]
79 //=========================================================================
80 class SpikedLine : public GeomObject
81 {
82 
83 public:
84 
85  /// \short Constructor: Four item of geometric data:
86  /// \code
87  /// Geom_data_pt[0]->value(0) = height
88  /// Geom_data_pt[0]->value(1) = amplitude
89  /// Geom_data_pt[0]->value(2) = zeta_min
90  /// Geom_data_pt[0]->value(3) = zeta_max
91  /// \endcode
92  SpikedLine(const Vector<Data*>& geom_data_pt) : GeomObject(1,2)
93  {
94 #ifdef PARANOID
95  if (geom_data_pt.size()!=1)
96  {
97  std::ostringstream error_stream;
98  error_stream
99  << "Wrong size for geom_data_pt " << geom_data_pt.size() << std::endl;
100  if (geom_data_pt[0]->nvalue()!=1)
101  {
102  error_stream << "Wrong geom_data_pt[0]->nvalue() "
103  << geom_data_pt[0]->nvalue() << std::endl;
104  }
105 
106  throw OomphLibError(error_stream.str(),
107  OOMPH_CURRENT_FUNCTION,
108  OOMPH_EXCEPTION_LOCATION);
109  }
110 #endif
111  Geom_data_pt.resize(1);
112  Geom_data_pt[0]=geom_data_pt[0];
113 
114  // Data has been created externally: Must not clean up
115  Must_clean_up=false;
116  }
117 
118 
119  /// \short Constructor: Pass height, amplitude, zeta min and zeta max
120  /// (all are pinned by default)
121  SpikedLine(const double& height, const double& amplitude,
122  const double& zeta_min, const double& zeta_max)
123  : GeomObject(1,2)
124  {
125  // Create Data for deflected-line object
126  Geom_data_pt.resize(1);
127 
128  // Create data: Four value, no timedependence, free by default
129  Geom_data_pt[0] = new Data(4);
130 
131  // I've created the data, I need to clean up
132  Must_clean_up=true;
133 
134  // Pin the data
135  for (unsigned i=0;i<4;i++) {Geom_data_pt[0]->pin(i);}
136 
137  // Give it a value: Initial height
138  Geom_data_pt[0]->set_value(0,height);
139  Geom_data_pt[0]->set_value(1,amplitude);
140  Geom_data_pt[0]->set_value(2,zeta_min);
141  Geom_data_pt[0]->set_value(3,zeta_max);
142  }
143 
144 
145  /// Destructor: Clean up if necessary
146  ~SpikedLine()
147  {
148  // Do I need to clean up?
149  if (Must_clean_up)
150  {
151  delete Geom_data_pt[0];
152  Geom_data_pt[0]=0;
153  }
154  }
155 
156 
157  /// \short Position Vector at Lagrangian coordinate zeta
158  void position(const Vector<double>& zeta, Vector<double>& r) const
159  {
160 #ifdef PARANOID
161  if (r.size()!=Ndim)
162  {
163  throw OomphLibError("The vector r has the wrong dimension\n",
164  OOMPH_CURRENT_FUNCTION,
165  OOMPH_EXCEPTION_LOCATION);
166  }
167 #endif
168  // Get parametres
169  double H = Geom_data_pt[0]->value(0);
170  double A = Geom_data_pt[0]->value(1);
171  double zeta_min = Geom_data_pt[0]->value(2);
172  double zeta_max = Geom_data_pt[0]->value(3);
173  double halfLz = (zeta_max-zeta_min)/2.0;
174 
175  // Position Vector
176  r[0] = zeta[0];
177  if (zeta[0]-zeta_min<=halfLz)
178  {
179  r[1] = H + A*(zeta[0]-zeta_min)/halfLz;
180  }
181  else
182  {
183  r[1] = H - A*(zeta[0]-zeta_max)/halfLz;
184  }
185  }
186 
187 
188  /// \short Parametrised position on object: r(zeta). Evaluated at
189  /// previous timestep. t=0: current time; t>0: previous
190  /// timestep.
191  void position(const unsigned& t, const Vector<double>& zeta,
192  Vector<double>& r) const
193  {
194 #ifdef PARANOID
195  if (t>Geom_data_pt[0]->time_stepper_pt()->nprev_values())
196  {
197  std::ostringstream error_stream;
198  error_stream
199  << "t > nprev_values() in SpikedLine: " << t << " "
200  << Geom_data_pt[0]->time_stepper_pt()->nprev_values() << std::endl;
201 
202  throw OomphLibError(error_stream.str(),
203  OOMPH_CURRENT_FUNCTION,
204  OOMPH_EXCEPTION_LOCATION);
205  }
206 #endif
207 
208  // Get parametres
209  double H = Geom_data_pt[0]->value(t,0);
210  double A = Geom_data_pt[0]->value(t,1);
211  double zeta_min = Geom_data_pt[0]->value(t,2);
212  double zeta_max = Geom_data_pt[0]->value(t,3);
213  double halfLz = (zeta_max-zeta_min)/2.0;
214 
215  // Position Vector
216  r[0] = zeta[0];
217  if (zeta[0]-zeta_min<=halfLz)
218  {
219  r[1] = H + A*(zeta[0]-zeta_min)/halfLz;
220  }
221  else
222  {
223  r[1] = H - A*(zeta[0]-zeta_max)/halfLz;
224  }
225  }
226 
227 
228  /// \short Derivative of position Vector w.r.t. to coordinates:
229  /// \f$ \frac{dR_i}{d \zeta_\alpha}\f$ = drdzeta(alpha,i).
230  /// Evaluated at current time.
231  virtual void dposition(const Vector<double>& zeta,
232  DenseMatrix<double> &drdzeta) const
233  {
234  // Get parametres
235  double A = Geom_data_pt[0]->value(1);
236  double zeta_min = Geom_data_pt[0]->value(2);
237  double zeta_max = Geom_data_pt[0]->value(3);
238  double halfLz = (zeta_max-zeta_min)/2.0;
239 
240  // Tangent vector
241  drdzeta(0,0)=1.0;
242  if (zeta[0]-zeta_min<=halfLz)
243  {
244  drdzeta(0,1)=A/halfLz;
245  }
246  else
247  {
248  drdzeta(0,1)=-A/halfLz;
249  }
250  }
251 
252 
253  /// \short 2nd derivative of position Vector w.r.t. to coordinates:
254  /// \f$ \frac{d^2R_i}{d \zeta_\alpha d \zeta_\beta}\f$ = ddrdzeta(alpha,beta,i).
255  /// Evaluated at current time.
256  virtual void d2position(const Vector<double>& zeta,
257  RankThreeTensor<double> &ddrdzeta) const
258  {
259  // Derivative of tangent vector
260  ddrdzeta(0,0,0)=0.0;
261  ddrdzeta(0,0,1)=0.0;
262  }
263 
264 
265  /// \short Posn Vector and its 1st & 2nd derivatives
266  /// w.r.t. to coordinates:
267  /// \f$ \frac{dR_i}{d \zeta_\alpha}\f$ = drdzeta(alpha,i).
268  /// \f$ \frac{d^2R_i}{d \zeta_\alpha d \zeta_\beta}\f$ = ddrdzeta(alpha,beta,i).
269  /// Evaluated at current time.
270  virtual void d2position(const Vector<double>& zeta, Vector<double>& r,
271  DenseMatrix<double> &drdzeta,
272  RankThreeTensor<double> &ddrdzeta) const
273  {
274  // Get parametres
275  double H = Geom_data_pt[0]->value(0);
276  double A = Geom_data_pt[0]->value(1);
277  double zeta_min = Geom_data_pt[0]->value(2);
278  double zeta_max = Geom_data_pt[0]->value(3);
279  double halfLz = (zeta_max-zeta_min)/2.0;
280 
281  // Position Vector
282  r[0] = zeta[0];
283  if (zeta[0]-zeta_min<=halfLz)
284  {
285  r[1] = H + A*(zeta[0]-zeta_min)/halfLz;
286  }
287  else
288  {
289  r[1] = H - A*(zeta[0]-zeta_max)/halfLz;
290  }
291 
292  // Tangent vector
293  drdzeta(0,0)=1.0;
294  if (zeta[0]-zeta_min<=halfLz)
295  {
296  drdzeta(0,1)=A/halfLz;
297  }
298  else
299  {
300  drdzeta(0,1)=-A/halfLz;
301  }
302 
303  // Derivative of tangent vector
304  ddrdzeta(0,0,0)=0.0;
305  ddrdzeta(0,0,1)=0.0;
306  }
307 
308  /// How many items of Data does the shape of the object depend on?
309  unsigned ngeom_data() const {return Geom_data_pt.size();}
310 
311  /// \short Return pointer to the j-th Data item that the object's
312  /// shape depends on
313  Data* geom_data_pt(const unsigned& j) {return Geom_data_pt[j];}
314 
315 
316 private:
317 
318  /// \short Vector of pointers to Data items that affects the object's shape
319  Vector<Data*> Geom_data_pt;
320 
321  /// Do I need to clean up?
322  bool Must_clean_up;
323 
324 };
325 
326 ///////////////////////////////////////////////////////////////////////
327 ///////////////////////////////////////////////////////////////////////
328 ///////////////////////////////////////////////////////////////////////
329 
330 
331 
332 //==start_of_problem_class============================================
333 /// Channel flow, through a non-uniform channel, using Spines.
334 //====================================================================
335 template<class ELEMENT>
336 class SpikedChannelSpineFlowProblem : public Problem
337 {
338 private:
339 
340  /// Width of channel
341  double Ly;
342 
343 public:
344 
345  /// Destructor: Empty
346  ~SpikedChannelSpineFlowProblem(){}
347 
348  /// Update the after solve (empty)
349  void actions_after_newton_solve() {}
350 
351  /// \short Update the problem specs before solve.
352  /// set velocity boundary conditions just to be on the safe side...
353  void actions_before_newton_solve()
354  {
355  // Update the mesh
356  mesh_pt()->node_update();
357 
358  // Setup inflow flow along boundary 3:
359  unsigned ibound=3;
360  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
361  for (unsigned inod=0;inod<num_nod;inod++)
362  {
363  double y=mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
364  // parabolic inflow
365  unsigned i=0;
366  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,y*(Ly-y));
367  // 0 Tangential flow
368  i=1;
369  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0);
370  }
371 
372  // Overwrite with no flow along top & bottom boundaries and
373  // parallel outflow
374  unsigned num_bound = mesh_pt()->nboundary();
375  for(unsigned ibound=0;ibound<num_bound-1;ibound++)
376  {
377  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
378  for (unsigned inod=0;inod<num_nod;inod++)
379  {
380  for (unsigned i=0;i<2;i++)
381  {
382  if (ibound!=1)
383  {
384  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0.0);
385  }
386  else
387  {
388  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
389  }
390  }
391  }
392  }
393  // Leave boundary 1 free!
394  } // end_of_actions_before_newton_solve
395 
396  /// Upcasted access function for the mesh
397  ChannelSpineMesh<ELEMENT>* mesh_pt()
398  {
399  return dynamic_cast<ChannelSpineMesh<ELEMENT>*>(Problem::mesh_pt());
400  }
401 
402  /// Constructor
403  SpikedChannelSpineFlowProblem();
404 
405  /// Doc the solution
406  void doc_solution(DocInfo& doc_info);
407 
408 }; // end_of_problem_class
409 
410 
411 
412 //==start_of_constructor==================================================
413 /// Constructor for SpikedChannelSpineFlow problem
414 ///
415 //========================================================================
416 template<class ELEMENT>
417 SpikedChannelSpineFlowProblem<ELEMENT>::SpikedChannelSpineFlowProblem()
418 {
419 
420  // Setup mesh
421 
422  // # of elements in x-direction in left region
423  unsigned Nx0=3;
424  // # of elements in x-direction in centre region
425  unsigned Nx1=12;
426  // # of elements in x-direction in right region
427  unsigned Nx2=8;
428 
429  // # of elements in y-direction
430  unsigned Ny=10;
431 
432  // Domain length in x-direction in left region
433  double Lx0=0.5;
434  // Domain length in x-direction in centre region
435  double Lx1=0.7;
436  // Domain length in x-direction in right region
437  double Lx2=1.5;
438 
439  // Domain length in y-direction
440  Ly=1.0;
441 
442  double amplitude_upper = -0.4*Ly;
443  double zeta_min=Lx0;
444  double zeta_max=Lx0+Lx1;
445  GeomObject* UpperWall =
446  new SpikedLine(Ly,amplitude_upper,zeta_min,zeta_max);
447 
448  // Build and assign mesh
449  Problem::mesh_pt() = new ChannelSpineMesh<ELEMENT>(Nx0,Nx1,Nx2,Ny,
450  Lx0,Lx1,Lx2,Ly,
451  UpperWall);
452 
453  // Set the boundary conditions for this problem: All nodes are
454  // free by default -- just pin the ones that have Dirichlet conditions
455  // here: All boundaries are Dirichlet boundaries, except on boundary 1
456  unsigned num_bound = mesh_pt()->nboundary();
457  for(unsigned ibound=0;ibound<num_bound;ibound++)
458  {
459  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
460  for (unsigned inod=0;inod<num_nod;inod++)
461  {
462  if (ibound!=1)
463  {
464  // Loop over values (u and v velocities)
465  for (unsigned i=0;i<2;i++)
466  {
467  mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
468  }
469  }
470  else
471  {
472  // Parallel outflow ==> no-slip
473  mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
474  }
475  }
476  } // end loop over boundaries
477 
478  // Find number of elements in mesh
479  unsigned n_element = mesh_pt()->nelement();
480 
481  // Loop over the elements to set up element-specific
482  // things that cannot be handled by constructor: Pass pointer to Reynolds
483  // number
484  for(unsigned e=0;e<n_element;e++)
485  {
486  // Upcast from GeneralisedElement to the present element
487  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
488  //Set the Reynolds number
489  el_pt->re_pt() = &Global_Physical_Variables2::Re;
490  } // end loop over elements
491 
492  // Setup equation numbering scheme
493  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
494 
495 } // end_of_constructor
496 
497 
498 
499 //==start_of_doc_solution=================================================
500 /// Doc the solution
501 //========================================================================
502 template<class ELEMENT>
503 void SpikedChannelSpineFlowProblem<ELEMENT>::doc_solution(DocInfo& doc_info)
504 {
505 
506  ofstream some_file;
507  char filename[100];
508 
509  // Number of plot points
510  unsigned npts=5;
511 
512  // Output solution
513  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
514  doc_info.number());
515  some_file.open(filename);
516  mesh_pt()->output(some_file,npts);
517  some_file.close();
518 
519 } // end_of_doc_solution
520 
521 
522 //==start_of_main======================================================
523 /// Driver for SpikedChannelSpineFlow test problem
524 //=====================================================================
525 int main()
526 {
527  // Set output directory
528  DocInfo doc_info;
529  doc_info.set_directory("RESLT");
530  doc_info.number()=0;
531 
532  // Solve problem with Taylor Hood elements
533  //---------------------------------------
534  {
535  //Build problem
536  SpikedChannelSpineFlowProblem<SpineElement<QTaylorHoodElement<2> > >
537  problem;
538 
539  // Solve the problem with automatic adaptation
540  problem.newton_solve();
541 
542  //Output solution
543  problem.doc_solution(doc_info);
544  // Step number
545  doc_info.number()++;
546 
547  } // end of Taylor Hood elements
548 
549 
550  // Solve problem with Crouzeix Raviart elements
551  //--------------------------------------------
552  {
553  // Build problem
554  SpikedChannelSpineFlowProblem<SpineElement<QCrouzeixRaviartElement<2> > >
555  problem;
556 
557  // Solve the problem with automatic adaptation
558  problem.newton_solve();
559 
560  //Output solution
561  problem.doc_solution(doc_info);
562  // Step number
563  doc_info.number()++;
564 
565  } // end of Crouzeix Raviart elements
566 
567 
568 } // end_of_main
569 
570 
571 
572 
573 
574 
575 
576 
577 
578 
579 
Global_Physical_Variables2
Namespace for physical parameters.
Definition: spine_channel2.cc:49
main
int main()
Driver for SpikedChannelSpineFlow test problem.
Definition: spine_channel2.cc:525
Global_Physical_Variables::height
double height(const double &x)
Height of domain.
Definition: simple_spine_channel.cc:433
SpikedLine::ngeom_data
unsigned ngeom_data() const
How many items of Data does the shape of the object depend on?
Definition: spine_channel2.cc:309
SpikedChannelSpineFlowProblem::actions_before_newton_solve
void actions_before_newton_solve()
Update the problem specs before solve. set velocity boundary conditions just to be on the safe side...
Definition: spine_channel2.cc:353
SpikedChannelSpineFlowProblem::~SpikedChannelSpineFlowProblem
~SpikedChannelSpineFlowProblem()
Destructor: Empty.
Definition: spine_channel2.cc:346
SpikedLine::~SpikedLine
~SpikedLine()
Destructor: Clean up if necessary.
Definition: spine_channel2.cc:146
SpikedLine::position
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...
Definition: spine_channel2.cc:191
SpikedLine::position
void position(const Vector< double > &zeta, Vector< double > &r) const
Position Vector at Lagrangian coordinate zeta.
Definition: spine_channel2.cc:158
SpikedLine::SpikedLine
SpikedLine(const Vector< Data *> &geom_data_pt)
Constructor: Four item of geometric data:
Definition: spine_channel2.cc:92
SpikedChannelSpineFlowProblem::doc_solution
void doc_solution(DocInfo &doc_info)
Doc the solution.
Definition: spine_channel2.cc:503
SpikedChannelSpineFlowProblem::SpikedChannelSpineFlowProblem
SpikedChannelSpineFlowProblem()
Constructor.
Definition: spine_channel2.cc:417
SpikedLine::d2position
virtual void d2position(const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
2nd derivative of position Vector w.r.t. to coordinates: = ddrdzeta(alpha,beta,i). Evaluated at current time.
Definition: spine_channel2.cc:256
SpikedLine::geom_data_pt
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object&#39;s shape depends on.
Definition: spine_channel2.cc:313
SpikedLine::SpikedLine
SpikedLine(const double &height, const double &amplitude, const double &zeta_min, const double &zeta_max)
Constructor: Pass height, amplitude, zeta min and zeta max (all are pinned by default) ...
Definition: spine_channel2.cc:121
Global_Physical_Variables2::Re
double Re
Reynolds number.
Definition: spine_channel2.cc:52
SpikedLine::Geom_data_pt
Vector< Data * > Geom_data_pt
Vector of pointers to Data items that affects the object&#39;s shape.
Definition: spine_channel2.cc:319
Global_Physical_Variables::H
double H
Undeformed height of domain.
Definition: simple_spine_channel.cc:427
SpikedLine::dposition
virtual void dposition(const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const
Derivative of position Vector w.r.t. to coordinates: = drdzeta(alpha,i). Evaluated at current time...
Definition: spine_channel2.cc:231
SpikedLine::Must_clean_up
bool Must_clean_up
Do I need to clean up?
Definition: spine_channel2.cc:322
SpikedChannelSpineFlowProblem::Ly
double Ly
Width of channel.
Definition: spine_channel2.cc:341
SpikedLine::d2position
virtual void d2position(const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
Posn Vector and its 1st & 2nd derivatives w.r.t. to coordinates: = drdzeta(alpha,i). = ddrdzeta(alpha,beta,i). Evaluated at current time.
Definition: spine_channel2.cc:270
SpikedChannelSpineFlowProblem
Channel flow, through a non-uniform channel, using Spines.
Definition: spine_channel2.cc:336
Global_Physical_Variables::A
double A
Amplitude of indentation.
Definition: simple_spine_channel.cc:430
SpikedChannelSpineFlowProblem::actions_after_newton_solve
void actions_after_newton_solve()
Update the after solve (empty)
Definition: spine_channel2.cc:349
SpikedChannelSpineFlowProblem::mesh_pt
ChannelSpineMesh< ELEMENT > * mesh_pt()
Upcasted access function for the mesh.
Definition: spine_channel2.cc:397