oomph-lib

00001 //LIC// ====================================================================
00002 //LIC// This file forms part of oomph-lib, the object-oriented, 
00003 //LIC// multi-physics finite-element library, available 
00004 //LIC// at http://www.oomph-lib.org.
00005 //LIC// 
00006 //LIC//           Version 0.90. August 3, 2009.
00007 //LIC// 
00008 //LIC// Copyright (C) 2006-2009 Matthias Heil and Andrew Hazel
00009 //LIC// 
00010 //LIC// This library is free software; you can redistribute it and/or
00011 //LIC// modify it under the terms of the GNU Lesser General Public
00012 //LIC// License as published by the Free Software Foundation; either
00013 //LIC// version 2.1 of the License, or (at your option) any later version.
00014 //LIC// 
00015 //LIC// This library is distributed in the hope that it will be useful,
00016 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 //LIC// Lesser General Public License for more details.
00019 //LIC// 
00020 //LIC// You should have received a copy of the GNU Lesser General Public
00021 //LIC// License along with this library; if not, write to the Free Software
00022 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
00023 //LIC// 02110-1301  USA.
00024 //LIC// 
00025 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
00026 //LIC// 
00027 //LIC//====================================================================
00028 //Driver for 2D Navier-Stokes problem in moving domain
00029 
00030 //Generic routines
00031 #include "generic.h"
00032 
00033 // The Navier Stokes equations
00034 #include "navier_stokes.h"
00035 
00036 // Mesh
00037 #include "meshes/quarter_circle_sector_mesh.h"
00038 
00039 using namespace std;
00040 
00041 using namespace oomph;
00042 
00043 //============start_of_MyEllipse===========================================
00044 /// \short Oscillating ellipse
00045 /// \f[ x = (A + \widehat{A} \cos(2\pi t/T)) \cos(\xi)  \f]
00046 /// \f[ y = \frac{\sin(\xi)}{A + \widehat{A} \cos(2\pi t/T)}   \f]
00047 /// Note that cross-sectional area is conserved.
00048 //=========================================================================
00049 class MyEllipse : public GeomObject
00050 {
00051 
00052 public:
00053 
00054  /// \short Constructor:  Pass initial x-half axis, amplitude of x-variation, 
00055  /// period of oscillation and pointer to time object.
00056  MyEllipse(const double& a, const double& a_hat,
00057            const double& period, Time* time_pt) : 
00058   GeomObject(1,2), A(a), A_hat(a_hat), T(period), Time_pt(time_pt) {}
00059 
00060  /// Destructor: Empty
00061  virtual ~MyEllipse() {}
00062 
00063  /// \short Current position vector to material point at 
00064  /// Lagrangian coordinate xi 
00065  void position(const Vector<double>& xi, Vector<double>& r) const
00066   {
00067    // Get current time:
00068    double time=Time_pt->time();
00069 
00070    // Position vector
00071    double axis=A+A_hat*cos(2.0*MathematicalConstants::Pi*time/T);
00072    r[0] = axis*cos(xi[0]);
00073    r[1] = (1.0/axis)*sin(xi[0]);
00074   } 
00075 
00076  /// \short Parametrised position on object: r(xi). Evaluated at
00077  /// previous time level. t=0: current time; t>0: previous
00078  /// time level.
00079  void position(const unsigned& t, const Vector<double>& xi,
00080                Vector<double>& r) const
00081   {
00082    // Get current time:
00083    double time=Time_pt->time(t);
00084    
00085    // Position vector
00086    double axis=A+A_hat*cos(2.0*MathematicalConstants::Pi*time/T);
00087    r[0] = axis*cos(xi[0]);
00088    r[1] = (1.0/axis)*sin(xi[0]);
00089   } 
00090 
00091 private:
00092 
00093  /// x-half axis
00094  double A;
00095 
00096  /// Amplitude of variation in x-half axis
00097  double A_hat;
00098 
00099  /// Period of oscillation
00100  double T;
00101 
00102  /// Pointer to time object
00103  Time* Time_pt;
00104 
00105 }; // end of MyEllipse
00106 
00107 
00108 
00109 /////////////////////////////////////////////////////////////////////// 
00110 ///////////////////////////////////////////////////////////////////////
00111 ////////////////////////////////////////////////////////////////////////
00112 
00113 
00114 
00115 //===start_of_namespace=================================================
00116 /// Namepspace for global parameters
00117 //======================================================================
00118 namespace Global_Physical_Variables
00119 {
00120 
00121  /// Reynolds number
00122  double Re=100.0;
00123 
00124  /// Womersley = Reynolds times Strouhal
00125  double ReSt=100.0;
00126 
00127  /// x-Half axis length
00128  double A=1.0;
00129 
00130  /// x-Half axis amplitude
00131  double A_hat=0.1;
00132 
00133  /// Period of oscillations
00134  double T=1.0;
00135 
00136  /// Exact solution of the problem as a vector containing u,v,p
00137  void get_exact_u(const double& t, const Vector<double>& x, Vector<double>& u)
00138  {
00139   using namespace MathematicalConstants;
00140 
00141   // Strouhal number
00142   double St = ReSt/Re;
00143 
00144   // Half axis
00145   double a=A+A_hat*cos(2.0*Pi*t/T);
00146   double adot=-2.0*A_hat*Pi*sin(2.0*Pi*t/T)/T; 
00147   u.resize(3);
00148 
00149   // Velocity solution
00150   u[0]=adot*x[0]/a;
00151   u[1]=-adot*x[1]/a;
00152 
00153   // Pressure solution
00154   u[2]=(2.0*A_hat*Pi*Pi*Re*(x[0]*x[0]*St*cos(2.0*Pi*t/T)*A + 
00155                             x[0]*x[0]*St*A_hat - x[0]*x[0]*A_hat +
00156                             x[0]*x[0]*A_hat*cos(2.0*Pi*t/T)*cos(2.0*Pi*t/T) -
00157                             x[1]*x[1]*St*cos(2.0*Pi*t/T)*A - 
00158                             x[1]*x[1]*St*A_hat - x[1]*x[1]*A_hat +
00159                             x[1]*x[1]*A_hat*cos(2.0*Pi*t/T)*cos(2.0*Pi*t/T) ))
00160    /(T*T*(A*A + 2.0*A*A_hat*cos(2.0*Pi*t/T) + 
00161           A_hat*A_hat*cos(2.0*Pi*t/T)*cos(2.0*Pi*t/T) ));
00162  }
00163 
00164 
00165 } // end of namespace
00166 
00167 
00168 
00169 
00170 /////////////////////////////////////////////////////////////////////// 
00171 ///////////////////////////////////////////////////////////////////////
00172 ////////////////////////////////////////////////////////////////////////
00173 
00174 
00175 
00176 //=====start_of_problem_class=========================================
00177 /// Navier-Stokes problem in an oscillating ellipse domain.
00178 //====================================================================
00179 template<class ELEMENT, class TIMESTEPPER>
00180 class OscEllipseProblem : public Problem
00181 {
00182 
00183 public:
00184 
00185  /// Constructor
00186  OscEllipseProblem();
00187 
00188  /// Destructor (empty)
00189  ~OscEllipseProblem() {}
00190 
00191  /// Update the problem specs after solve (empty)
00192  void actions_after_newton_solve(){}
00193 
00194  /// \short Update problem specs before solve (empty)
00195  void actions_before_newton_solve() {} 
00196  
00197  /// Actions before adapt (empty)
00198  void actions_before_adapt(){}
00199 
00200  /// Actions after adaptation, pin relevant pressures
00201  void actions_after_adapt()
00202   {
00203    // Unpin all pressure dofs
00204    RefineableNavierStokesEquations<2>::
00205     unpin_all_pressure_dofs(mesh_pt()->element_pt());
00206 
00207    // Pin redundant pressure dofs
00208    RefineableNavierStokesEquations<2>::
00209     pin_redundant_nodal_pressures(mesh_pt()->element_pt());
00210    
00211    // Now set the first pressure dof in the first element to 0.0
00212    fix_pressure(0,0,0.0);
00213 
00214   } // end of actions_after_adapt
00215 
00216 
00217  /// \short Update the problem specs before next timestep
00218  void actions_before_implicit_timestep()
00219   {
00220    // Update the domain shape
00221    mesh_pt()->node_update();
00222 
00223    // Ring boundary: No slip; this implies that the velocity needs
00224    // to be updated in response to wall motion
00225    unsigned ibound=1;
00226    unsigned num_nod=mesh_pt()->nboundary_node(ibound);
00227    for (unsigned inod=0;inod<num_nod;inod++)
00228     {
00229      // Which node are we dealing with?
00230      Node* node_pt=mesh_pt()->boundary_node_pt(ibound,inod);
00231      
00232      // Apply no slip
00233      FSI_functions::apply_no_slip_on_moving_wall(node_pt);
00234     }
00235   } 
00236 
00237  /// Update the problem specs after timestep (empty)
00238  void actions_after_implicit_timestep(){}
00239  
00240  /// Doc the solution
00241  void doc_solution(DocInfo& doc_info);
00242 
00243  /// Timestepping loop
00244  void unsteady_run(DocInfo& doc_info);
00245 
00246  /// \short Set initial condition
00247  void set_initial_condition();
00248 
00249 private:
00250 
00251  ///Fix pressure in element e at pressure dof pdof and set to pvalue
00252  void fix_pressure(const unsigned &e, const unsigned &pdof, 
00253                    const double &pvalue)
00254   {
00255    //Cast to proper element and fix pressure
00256    dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
00257     fix_pressure(pdof,pvalue);
00258   } // end_of_fix_pressure
00259 
00260  /// Pointer to GeomObject that specifies the domain bondary
00261  GeomObject* Wall_pt;
00262 
00263 }; // end of problem_class
00264 
00265 
00266 
00267 //========start_of_constructor============================================
00268 /// Constructor for Navier-Stokes problem on an oscillating ellipse domain.
00269 //========================================================================
00270 template<class ELEMENT, class TIMESTEPPER>
00271 OscEllipseProblem<ELEMENT,TIMESTEPPER>::OscEllipseProblem()
00272 { 
00273 
00274  //Create the timestepper and add it to the problem
00275  add_time_stepper_pt(new TIMESTEPPER);
00276 
00277  // Setup mesh
00278  //-----------
00279 
00280  // Build geometric object that forms the curvilinear domain boundary:
00281  // an oscillating ellipse
00282 
00283  // Half axes
00284  double a=Global_Physical_Variables::A;
00285 
00286  // Variations of half axes
00287  double a_hat=Global_Physical_Variables::A_hat;
00288 
00289  // Period of the oscillation
00290  double period=Global_Physical_Variables::T;
00291 
00292  // Create GeomObject that specifies the domain bondary
00293  Wall_pt=new MyEllipse(a,a_hat,period,Problem::time_pt()); 
00294 
00295 
00296  // Start and end coordinates of curvilinear domain boundary on ellipse
00297  double xi_lo=0.0;
00298  double xi_hi=MathematicalConstants::Pi/2.0;
00299 
00300  // Now create the mesh. Separating line between the two 
00301  // elements next to the curvilinear boundary is located half-way
00302  // along the boundary.
00303  double fract_mid=0.5;
00304  Problem::mesh_pt() = new RefineableQuarterCircleSectorMesh<ELEMENT>(
00305   Wall_pt,xi_lo,fract_mid,xi_hi,time_stepper_pt());
00306 
00307  // Set error estimator NOT NEEDED IN CURRENT PROBLEM SINCE 
00308  // WE'RE ONLY REFINING THE MESH UNIFORMLY
00309  //Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
00310  //mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
00311 
00312 
00313  // Fluid boundary conditions
00314  //--------------------------
00315  // Ring boundary: No slip; this also implies that the velocity needs
00316  // to be updated in response to wall motion
00317  unsigned ibound=1;
00318  {
00319   unsigned num_nod= mesh_pt()->nboundary_node(ibound);
00320   for (unsigned inod=0;inod<num_nod;inod++)
00321    {
00322      
00323     // Pin both velocities
00324     for (unsigned i=0;i<2;i++)
00325      {
00326       mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
00327      }
00328    }
00329  } // end boundary 1
00330 
00331 //   // hierher keep for algebraic element version
00332 //   unsigned num_nod= mesh_pt()->nboundary_node(ibound);
00333 //   for (unsigned inod=0;inod<num_nod;inod++)
00334 //    {
00335 //     // Pin both velocities
00336 //     for (unsigned i=0;i<2;i++)
00337 //      {
00338 //      // Which node are we dealing with?
00339 //      Node* node_pt=mesh_pt()->boundary_node_pt(ibound,inod);
00340  
00341 //      // Set auxiliary update function pointer
00342 //      node_pt->set_auxiliary_node_update_fct_pt(
00343 //       FSI_functions::apply_no_slip_on_moving_wall);
00344 
00345 //      mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
00346 //      }
00347 //    }
00348  
00349  // Bottom boundary: 
00350  ibound=0;
00351  {
00352   unsigned num_nod= mesh_pt()->nboundary_node(ibound);
00353   for (unsigned inod=0;inod<num_nod;inod++)
00354    {
00355     // Pin vertical velocity
00356     {
00357      mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
00358     }
00359    }
00360  } // end boundary 0
00361 
00362  // Left boundary:
00363  ibound=2;
00364  {
00365   unsigned num_nod= mesh_pt()->nboundary_node(ibound);
00366   for (unsigned inod=0;inod<num_nod;inod++)
00367    {
00368     // Pin horizontal velocity
00369     {
00370      mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
00371     }
00372    }
00373  } // end boundary 2
00374 
00375  
00376  // Complete the build of all elements so they are fully functional
00377  //----------------------------------------------------------------
00378 
00379  // Find number of elements in mesh
00380  unsigned n_element = mesh_pt()->nelement();
00381 
00382  // Loop over the elements to set up element-specific 
00383  // things that cannot be handled by constructor
00384  for(unsigned i=0;i<n_element;i++)
00385   {
00386    // Upcast from FiniteElement to the present element
00387    ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
00388 
00389    //Set the Reynolds number, etc
00390    el_pt->re_pt() = &Global_Physical_Variables::Re;
00391    el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
00392 
00393    // Set pointer to continous time
00394    el_pt->time_pt()=time_pt();
00395   }
00396 
00397  // Pin redundant pressure dofs
00398  RefineableNavierStokesEquations<2>::
00399   pin_redundant_nodal_pressures(mesh_pt()->element_pt());
00400  
00401  // Now set the first pressure dof in the first element to 0.0
00402  fix_pressure(0,0,0.0);
00403 
00404  // Do equation numbering
00405  cout << "Number of equations: " << assign_eqn_numbers() << std::endl; 
00406 
00407 } // end of constructor
00408 
00409 
00410 //======================start_of_set_initial_condition====================
00411 /// \short Set initial condition: Assign previous and current values
00412 /// from exact solution.
00413 //========================================================================
00414 template<class ELEMENT,class TIMESTEPPER>
00415 void OscEllipseProblem<ELEMENT,TIMESTEPPER>::set_initial_condition()
00416 { 
00417  // Backup time in global timestepper
00418  double backed_up_time=time_pt()->time();
00419    
00420  // Past history for velocities must be established for t=time0-deltat, ...
00421  // Then provide current values (at t=time0) which will also form
00422  // the initial guess for first solve at t=time0+deltat
00423  
00424  // Vector of exact solution value
00425  Vector<double> soln(3);
00426  Vector<double> x(2);
00427  
00428  //Find number of nodes in mesh
00429  unsigned num_nod = mesh_pt()->nnode();
00430  
00431  // Get continuous times at previous timesteps
00432  int nprev_steps=time_stepper_pt()->nprev_values();
00433  Vector<double> prev_time(nprev_steps+1);
00434  for (int itime=nprev_steps;itime>=0;itime--)
00435   {
00436    prev_time[itime]=time_pt()->time(unsigned(itime));
00437   }
00438  
00439  // Loop over current & previous timesteps (in outer loop because
00440  // the mesh also moves!)
00441  for (int itime=nprev_steps;itime>=0;itime--)
00442   {
00443    double time=prev_time[itime];
00444    
00445    // Set global time (because this is how the geometric object refers 
00446    // to continous time )
00447    time_pt()->time()=time;
00448    
00449    cout << "setting IC at time =" << time << std::endl;
00450    
00451    // Update the mesh for this value of the continuous time
00452    // (The wall object reads the continous time from the 
00453    // global time object)
00454    mesh_pt()->node_update(); 
00455    
00456    // Loop over the nodes to set initial guess everywhere
00457    for (unsigned jnod=0;jnod<num_nod;jnod++)
00458     {
00459      // Get nodal coordinates
00460      x[0]=mesh_pt()->node_pt(jnod)->x(0);
00461      x[1]=mesh_pt()->node_pt(jnod)->x(1);
00462      
00463      // Get exact solution (unsteady stagnation point flow)
00464      Global_Physical_Variables::get_exact_u(time,x,soln);
00465      
00466      // Assign solution
00467      mesh_pt()->node_pt(jnod)->set_value(itime,0,soln[0]);
00468      mesh_pt()->node_pt(jnod)->set_value(itime,1,soln[1]);
00469        
00470      // Loop over coordinate directions
00471      for (unsigned i=0;i<2;i++)
00472       {
00473        mesh_pt()->node_pt(jnod)->x(itime,i)=x[i];
00474       }
00475     } 
00476   } // end of loop over previous timesteps
00477  
00478  // Reset backed up time for global timestepper
00479  time_pt()->time()=backed_up_time;
00480  
00481 } // end of set_initial_condition
00482 
00483 
00484 
00485 //=======start_of_doc_solution============================================
00486 /// Doc the solution
00487 //========================================================================
00488 template<class ELEMENT, class TIMESTEPPER>
00489 void OscEllipseProblem<ELEMENT,TIMESTEPPER>::doc_solution(DocInfo& doc_info)
00490 { 
00491  ofstream some_file;
00492  char filename[100];
00493 
00494  // Number of plot points
00495  unsigned npts;
00496  npts=5;
00497 
00498  // Output solution 
00499  //-----------------
00500  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
00501          doc_info.number());
00502  some_file.open(filename);
00503  mesh_pt()->output(some_file,npts);
00504  some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = " 
00505            << time_pt()->time() << "\"";
00506  some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
00507  some_file << "1" << std::endl;
00508  some_file << "2" << std::endl;
00509  some_file << " 0 0" << std::endl;
00510  some_file << time_pt()->time()*20.0 << " 0" << std::endl;
00511 
00512  // Write dummy zones that force tecplot to keep the axis limits constant
00513  // while the domain is moving.
00514  some_file << "ZONE I=2,J=2" << std::endl;
00515  some_file << "0.0 0.0 -0.65 -0.65 -200.0" << std::endl;
00516  some_file << "1.15 0.0 -0.65 -0.65 -200.0" << std::endl;
00517  some_file << "0.0 1.15 -0.65 -0.65 -200.0" << std::endl;
00518  some_file << "1.15 1.15 -0.65 -0.65 -200.0" << std::endl;
00519  some_file << "ZONE I=2,J=2" << std::endl;
00520  some_file << "0.0 0.0 0.65 0.65 300.0" << std::endl;
00521  some_file << "1.15 0.0 0.65 0.65 300.0" << std::endl;
00522  some_file << "0.0 1.15 0.65 0.65 300.0" << std::endl;
00523  some_file << "1.15 1.15 0.65 0.65 300.0" << std::endl;
00524 
00525  some_file.close();
00526 
00527  // Output exact solution 
00528  //----------------------
00529  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
00530          doc_info.number());
00531  some_file.open(filename);
00532  mesh_pt()->output_fct(some_file,npts,time_pt()->time(),
00533                        Global_Physical_Variables::get_exact_u); 
00534  some_file.close();
00535 
00536  // Doc error
00537  //----------
00538  double error,norm;
00539  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
00540          doc_info.number());
00541  some_file.open(filename);
00542  mesh_pt()->compute_error(some_file,
00543                           Global_Physical_Variables::get_exact_u,
00544                           time_pt()->time(),
00545                           error,norm); 
00546  some_file.close();
00547 
00548 
00549  // Doc solution and error
00550  //-----------------------
00551  cout << "error: " << error << std::endl; 
00552  cout << "norm : " << norm << std::endl << std::endl;
00553 
00554 
00555  // Plot wall posn
00556  //---------------
00557  sprintf(filename,"%s/Wall%i.dat",doc_info.directory().c_str(),
00558          doc_info.number());
00559  some_file.open(filename);
00560  
00561  unsigned nplot=100;
00562  for (unsigned iplot=0;iplot<nplot;iplot++)
00563   {
00564    Vector<double> xi_wall(1), r_wall(2);
00565    xi_wall[0]=0.5*MathematicalConstants::Pi*double(iplot)/double(nplot-1);
00566    Wall_pt->position(xi_wall,r_wall);
00567    some_file << r_wall[0] << " " << r_wall[1] << std::endl;
00568   }
00569  some_file.close();
00570  
00571  // Increment number of doc
00572  doc_info.number()++;
00573 
00574 } // end of doc_solution
00575 
00576 
00577 //=======start_of_unsteady_run============================================
00578 /// Unsteady run
00579 //========================================================================
00580 template<class ELEMENT, class TIMESTEPPER>
00581 void OscEllipseProblem<ELEMENT,TIMESTEPPER>::unsteady_run(DocInfo& doc_info)
00582 {
00583 
00584  // Specify duration of the simulation
00585  double t_max=3.0;
00586 
00587  // Initial timestep
00588  double dt=0.025;
00589 
00590  // Initialise timestep
00591  initialise_dt(dt);
00592 
00593  // Set initial conditions.
00594  set_initial_condition();
00595 
00596  // Alternative initial conditions: impulsive start; see exercise.
00597  //assign_initial_values_impulsive(); 
00598 
00599  // find number of steps
00600  unsigned nstep = unsigned(t_max/dt);
00601 
00602  // If validation: Reduce number of timesteps performed and 
00603  // use coarse-ish mesh
00604  if (CommandLineArgs::Argc>1)
00605   {
00606    nstep=2;
00607    refine_uniformly();
00608    cout << "validation run" << std::endl;
00609   }
00610  else
00611   {
00612    // Refine the mesh three times, to resolve the pressure distribution
00613    // (the velocities could be represented accurately on a much coarser mesh).
00614    refine_uniformly();
00615    refine_uniformly();
00616    refine_uniformly();
00617   }
00618 
00619  // Output solution initial 
00620  doc_solution(doc_info);
00621 
00622  // Timestepping loop
00623  for (unsigned istep=0;istep<nstep;istep++)
00624   {
00625    cout << "TIMESTEP " << istep << std::endl;
00626    cout << "Time is now " << time_pt()->time() << std::endl;
00627 
00628    // Take timestep 
00629    unsteady_newton_solve(dt);
00630      
00631    //Output solution
00632    doc_solution(doc_info);
00633   }
00634 
00635 } // end of unsteady_run
00636 
00637 
00638 ////////////////////////////////////////////////////////////////////////
00639 ////////////////////////////////////////////////////////////////////////
00640 
00641 
00642 //======start_of_main=================================================
00643 /// Driver code for unsteady Navier-Stokes flow, driven by
00644 /// oscillating ellipse. If the code is executed with command line
00645 /// arguments, a validation run is performed. 
00646 //====================================================================
00647 int main(int argc, char* argv[])
00648 {
00649 
00650  // Store command line arguments
00651  CommandLineArgs::setup(argc,argv);
00652 
00653 
00654  // Solve with Crouzeix-Raviart elements
00655  {
00656   // Create DocInfo object with suitable directory name for output
00657   DocInfo doc_info;
00658   doc_info.set_directory("RESLT_CR");
00659   
00660   //Set up problem
00661   OscEllipseProblem<RefineableQCrouzeixRaviartElement<2>,BDF<2> > problem;
00662   
00663   // Run the unsteady simulation
00664   problem.unsteady_run(doc_info);
00665  }
00666 
00667  // Solve with Taylor-Hood elements
00668  {
00669   // Create DocInfo object with suitable directory name for output
00670   DocInfo doc_info;
00671   doc_info.set_directory("RESLT_TH");
00672 
00673   //Set up problem
00674   OscEllipseProblem<RefineableQTaylorHoodElement<2>,BDF<2> > problem;
00675   
00676   // Run the unsteady simulation
00677   problem.unsteady_run(doc_info);
00678  }
00679  
00680 
00681 }; // end of main
00682 
00683 
00684
osc_quarter_ellipse.cc
