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 #ifndef OOMPH_CHANNEL_WITH_LEAFLET_MESH_HEADER
00029 #define OOMPH_CHANNEL_WITH_LEAFLET_MESH_HEADER
00030 
00031 // Generic includes
00032 #include "../generic/refineable_quad_mesh.h"
00033 #include "../generic/macro_element.h"
00034 #include "../generic/domain.h"
00035 #include "../generic/quad_mesh.h"
00036 
00037 // Mesh is based on simple_rectangular_quadmesh
00038 #include "simple_rectangular_quadmesh.template.h"
00039 #include "simple_rectangular_quadmesh.template.cc"
00040 
00041 //Include macro elements 
00042 #include "../generic/macro_element_node_update_element.h"
00043 
00044 //and algebraic elements
00045 #include "../generic/algebraic_elements.h"
00046 
00047 //Include the headers file for domain
00048 #include "channel_with_leaflet_domain.h"
00049 
00050 namespace oomph
00051 {
00052 
00053 //===================================================================
00054 /// Channel with leaflet mesh
00055 //===================================================================
00056 template <class ELEMENT> 
00057 class ChannelWithLeafletMesh : public SimpleRectangularQuadMesh<ELEMENT>
00058 {
00059   public:
00060 
00061 
00062  ///\short Constructor: Pass pointer to GeomObject that represents the leaflet,
00063  /// the length of the domain to left and right of the leaflet, the
00064  /// height of the leaflet and the overall height of the channel,
00065  /// the number of element columns to the left and right of the leaflet,
00066  /// the number of rows of elements from the bottom of the channel to 
00067  /// the end of the leaflet, the number of rows of elements above the
00068  /// end of the leaflet. Timestepper defaults to Steady default 
00069  /// Timestepper defined in the Mesh base class
00070  ChannelWithLeafletMesh(GeomObject* leaflet_pt, const double& lleft,
00071                         const double& lright, const double& hleaflet,
00072                         const double& htot,
00073                         const unsigned& nleft, const unsigned& nright,
00074                         const unsigned& ny1, const unsigned&  ny2,
00075                         TimeStepper* time_stepper_pt=
00076                         &Mesh::Default_TimeStepper);
00077 
00078  ///Destructor : empty
00079  virtual ~ChannelWithLeafletMesh(){}
00080  
00081  /// Access function to domain
00082  ChannelWithLeafletDomain * domain_pt(){return Domain_pt;}
00083 
00084 protected:
00085 
00086  /// Pointer to domain
00087  ChannelWithLeafletDomain * Domain_pt;
00088 
00089  /// Pointer to GeomObject that represents the leaflet
00090  GeomObject * Leaflet_pt;
00091 
00092 }; 
00093 
00094 
00095 /////////////////////////////////////////////////////////////////////
00096 /////////////////////////////////////////////////////////////////////
00097 /////////////////////////////////////////////////////////////////////
00098 
00099 
00100 
00101 //===================================================================
00102 /// Refineable version of ChannelWithLeafletMesh
00103 //===================================================================
00104 template <class ELEMENT> 
00105 class RefineableChannelWithLeafletMesh :
00106 public virtual ChannelWithLeafletMesh<ELEMENT>,
00107  public RefineableQuadMesh<ELEMENT>
00108 {
00109   public:
00110  
00111  ///\short Constructor: Pass pointer to GeomObject that represents the leaflet,
00112  /// the length of the domain to left and right of the leaflet, the
00113  /// height of the leaflet and the overall height of the channel,
00114  /// the number of element columns to the left and right of the leaflet,
00115  /// the number of rows of elements from the bottom of the channel to 
00116  /// the end of the leaflet, the number of rows of elements above the
00117  /// end of the leaflet. Timestepper defaults to Steady default 
00118  /// Timestepper defined in the Mesh base class
00119  RefineableChannelWithLeafletMesh(GeomObject* leaflet_pt,
00120                                   const double& lleft,const double& lright,
00121                                   const double& hleaflet,const double& htot,
00122                                   const unsigned& nleft,const unsigned& nright,
00123                                   const unsigned& ny1,const unsigned&  ny2,
00124                                   TimeStepper* time_stepper_pt=
00125                                   &Mesh::Default_TimeStepper):
00126   ChannelWithLeafletMesh<ELEMENT>(leaflet_pt,lleft,lright,hleaflet,htot,
00127                                   nleft,nright,ny1,ny2,
00128                                   time_stepper_pt)
00129   {  
00130    // Build quadtree forest
00131    this->setup_quadtree_forest();                 
00132   }
00133  
00134  /// Destructor (empty)
00135  virtual ~RefineableChannelWithLeafletMesh(){}
00136  
00137 };
00138 
00139 
00140 
00141 /////////////////////////////////////////////////////////////////////
00142 /////////////////////////////////////////////////////////////////////
00143 /////////////////////////////////////////////////////////////////////
00144 
00145 
00146 //=====start_of_mesh=======================================================
00147 /// Channel with leaflet mesh with MacroElement-based node update.
00148 /// The leaflet is represented by the specified geometric object.
00149 /// Some or all of the geometric Data in that geometric object
00150 /// may contain unknowns in the global Problem. The dependency
00151 /// on these unknowns is taken into account when setting up
00152 /// the Jacobian matrix of the elements. For this purpose,
00153 /// the element (whose type is specified by the template parameter)
00154 /// must inherit from MacroElementNodeUpdateElementBase. 
00155 //========================================================================
00156 template<class ELEMENT>
00157 class MacroElementNodeUpdateChannelWithLeafletMesh : 
00158 public virtual MacroElementNodeUpdateMesh, 
00159  public virtual ChannelWithLeafletMesh<ELEMENT>
00160 {
00161  
00162   public: 
00163  
00164  ///\short Constructor: Pass pointer to GeomObject that represents the leaflet,
00165  /// the length of the domain to left and right of the leaflet, the
00166  /// height of the leaflet and the overall height of the channel,
00167  /// the number of element columns to the left and right of the leaflet,
00168  /// the number of rows of elements from the bottom of the channel to 
00169  /// the end of the leaflet, the number of rows of elements above the
00170  /// end of the leaflet. Timestepper defaults to Steady default 
00171  /// Timestepper defined in the Mesh base class
00172  MacroElementNodeUpdateChannelWithLeafletMesh(
00173   GeomObject* leaflet_pt, const double& lleft,
00174   const double& lright, const double& hleaflet,
00175   const double& htot,
00176   const unsigned& nleft, const unsigned& nright,
00177   const unsigned& ny1, const unsigned&  ny2,
00178   TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper) :
00179   ChannelWithLeafletMesh<ELEMENT>(leaflet_pt,lleft,lright,hleaflet,htot,
00180                                   nleft,nright,ny1,ny2,
00181                                   time_stepper_pt)
00182   {
00183 #ifdef PARANOID
00184    ELEMENT* el_pt=new ELEMENT;
00185    if (dynamic_cast<MacroElementNodeUpdateElementBase*>(el_pt)==0)
00186     {
00187      std::ostringstream error_message;
00188      error_message 
00189       << "Base class for ELEMENT in "
00190       << "MacroElementNodeUpdateChannelWithLeafletMesh needs" 
00191       << "to be of type MacroElementNodeUpdateElement!\n";
00192      error_message << "Whereas it is: typeid(el_pt).name()" 
00193                    << typeid(el_pt).name() 
00194                    << std::endl;
00195      
00196      std::string function_name =
00197       "MacroElementNodeUpdateChannelWithLeafletMesh::\n";
00198      function_name += 
00199       "MacroElementNodeUpdateChannelWithLeafletMesh()";
00200      
00201      throw OomphLibError(error_message.str(),function_name,
00202                          OOMPH_EXCEPTION_LOCATION);
00203     }
00204    delete el_pt;
00205 #endif
00206    
00207    // Setup all the information that's required for MacroElement-based
00208    // node update: Tell the elements that their geometry depends on the
00209    // wall geometric object
00210    unsigned n_element = this->nelement();
00211    for(unsigned i=0;i<n_element;i++)
00212     {
00213      // Upcast from FiniteElement to the present element
00214      ELEMENT *el_pt = dynamic_cast<ELEMENT*>(this->element_pt(i));
00215      
00216 #ifdef PARANOID
00217      // Check if cast is successful
00218      MacroElementNodeUpdateElementBase* m_el_pt=dynamic_cast<
00219       MacroElementNodeUpdateElementBase*>(el_pt);
00220      if (m_el_pt==0)
00221       {
00222        std::ostringstream error_message;
00223        error_message 
00224         << "Failed to upcast to MacroElementNodeUpdateElementBase\n";
00225        error_message
00226         << "Element must be derived from MacroElementNodeUpdateElementBase\n";
00227        error_message << "but it is of type " << typeid(el_pt).name();
00228        
00229        std::string function_name =
00230         "MacroElementNodeUpdateChannelWithLeafletMesh::\n";
00231        function_name += 
00232         "MacroElementNodeUpdateChannelWithLeafletMesh()";
00233        
00234        throw OomphLibError(error_message.str(),function_name,
00235                            OOMPH_EXCEPTION_LOCATION);
00236       }
00237 #endif       
00238      
00239      // There's just one GeomObject
00240      Vector<GeomObject*> geom_object_pt(1);
00241      geom_object_pt[0] = this->Leaflet_pt;
00242      
00243      // Tell the element which geom objects its macro-element-based
00244      // node update depends on     
00245      el_pt->set_node_update_info(geom_object_pt);
00246     }
00247 
00248    // Add the geometric object(s) for the wall to the mesh's storage
00249    Vector<GeomObject*> geom_object_pt(1);
00250    geom_object_pt[0] = this->Leaflet_pt;
00251    MacroElementNodeUpdateMesh::set_geom_object_vector_pt(geom_object_pt);
00252 
00253    // Fill in the domain pointer to the mesh's storage in the base class
00254    MacroElementNodeUpdateMesh::macro_domain_pt()=this->domain_pt();
00255 
00256   } // end of constructor
00257  
00258  
00259  /// \short Destructor: empty
00260  virtual ~MacroElementNodeUpdateChannelWithLeafletMesh(){}
00261  
00262 
00263 }; //end of mesh
00264 
00265 
00266 
00267 
00268 ////////////////////////////////////////////////////////////////////////////
00269 ////////////////////////////////////////////////////////////////////////////
00270 ////////////////////////////////////////////////////////////////////////////
00271 
00272 
00273 //=====start_of_mesh=======================================================
00274 /// Refineable  mesh with MacroElement-based node update.
00275 //========================================================================
00276 template<class ELEMENT>
00277 class MacroElementNodeUpdateRefineableChannelWithLeafletMesh : 
00278  public virtual MacroElementNodeUpdateChannelWithLeafletMesh<ELEMENT>, 
00279  public virtual RefineableQuadMesh<ELEMENT>
00280 {
00281   
00282 public: 
00283 
00284  ///\short Constructor: Pass pointer to GeomObject that represents the leaflet,
00285  /// the length of the domain to left and right of the leaflet, the
00286  /// height of the leaflet and the overall height of the channel,
00287  /// the number of element columns to the left and right of the leaflet,
00288  /// the number of rows of elements from the bottom of the channel to 
00289  /// the end of the leaflet, the number of rows of elements above the
00290  /// end of the leaflet. Timestepper defaults to Steady default 
00291  /// Timestepper defined in the Mesh base class
00292  MacroElementNodeUpdateRefineableChannelWithLeafletMesh(
00293   GeomObject* leaflet_pt, const double& lleft,
00294   const double& lright, const double& hleaflet,
00295   const double& htot,
00296   const unsigned& nleft, const unsigned& nright,
00297   const unsigned& ny1, const unsigned&  ny2,
00298   TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper) :
00299   ChannelWithLeafletMesh<ELEMENT>(leaflet_pt,lleft,lright,hleaflet,htot,
00300                                   nleft,nright,ny1,ny2,
00301                                   time_stepper_pt),
00302   MacroElementNodeUpdateChannelWithLeafletMesh<ELEMENT>(
00303    leaflet_pt,lleft,lright,hleaflet,htot,
00304    nleft,nright,ny1,ny2,
00305    time_stepper_pt)
00306   {
00307    // Build quadtree forest
00308    this->setup_quadtree_forest();   
00309   } 
00310  
00311 
00312  /// \short Destructor: empty
00313  virtual ~MacroElementNodeUpdateRefineableChannelWithLeafletMesh(){}
00314  
00315  }; //end of mesh
00316 
00317 
00318 ///////////////////////////////////////////////////////////////////////
00319 //////////////////////////////////////////////////////////////////////
00320 ///////////////////////////////////////////////////////////////////////
00321 
00322 
00323 //=================================================================
00324 /// Algebraic version of ChannelWithLeafletMesh. Leaflet is
00325 /// assumed to be in its undeformed (straight vertical) position
00326 /// when the algebraic node update is set up.
00327 //=================================================================
00328 template<class ELEMENT>
00329 class AlgebraicChannelWithLeafletMesh : public AlgebraicMesh,
00330 public virtual ChannelWithLeafletMesh<ELEMENT>
00331 {
00332   public:
00333  
00334  ///\short Constructor: Pass pointer to GeomObject that represents the leaflet,
00335  /// the length of the domain to left and right of the leaflet, the
00336  /// height of the leaflet and the overall height of the channel,
00337  /// the number of element columns to the left and right of the leaflet,
00338  /// the number of rows of elements from the bottom of the channel to 
00339  /// the end of the leaflet, the number of rows of elements above the
00340  /// end of the leaflet. Timestepper defaults to Steady default 
00341  /// Timestepper defined in the Mesh base class
00342  AlgebraicChannelWithLeafletMesh(GeomObject* leaflet_pt, const double& lleft,
00343                                  const double& lright, const double& hleaflet,
00344                                  const double& htot,
00345                                  const unsigned& nleft, const unsigned& nright,
00346                                  const unsigned& ny1, const unsigned&  ny2,
00347                                  TimeStepper* time_stepper_pt=
00348                                  &Mesh::Default_TimeStepper) :
00349   ChannelWithLeafletMesh<ELEMENT>(leaflet_pt,lleft,lright,hleaflet,htot,
00350                                   nleft,nright,ny1,ny2,
00351                                   time_stepper_pt)
00352   { 
00353    // Store origin of leaflet for fast reference
00354    Vector<double> zeta(1);
00355    zeta[0]=0.0;
00356    Vector<double> r(2);
00357    this->Leaflet_pt->position(zeta,r);
00358    X_0 = r[0];
00359 
00360    // Store length of the leaflet for fast access (it's also available
00361    // through the domain, of course)
00362    Hleaflet = hleaflet;  
00363    
00364    // Add the geometric object to the list associated with this AlgebraicMesh
00365    AlgebraicMesh::add_geom_object_list_pt(leaflet_pt);
00366 
00367    //Setup algebraic node update operations
00368    setup_algebraic_node_update();
00369   }
00370  
00371 
00372  /// \short Destructor: empty
00373  virtual ~AlgebraicChannelWithLeafletMesh(){}
00374 
00375 
00376  /// \short Update the geometric references that are used 
00377  /// to update node after mesh adaptation.
00378  /// Empty -- no update of node update required without adaptivity
00379  void update_node_update(AlgebraicNode*& node_pt) {}
00380 
00381  /// \short Update nodal position at time level t (t=0: present; 
00382  /// t>0: previous)
00383  void algebraic_node_update(const unsigned& t, AlgebraicNode*& node_pt);
00384                    
00385 protected:
00386  
00387  /// Function to setup the algebraic node update
00388  void setup_algebraic_node_update();
00389 
00390  /// Update function for nodes in lower left region (I)
00391  void node_update_I(const unsigned& t, AlgebraicNode*& node_pt);
00392 
00393  /// Update function for nodes in lower right region (II)
00394  void node_update_II(const unsigned& t, AlgebraicNode*& node_pt);
00395 
00396  /// Update function for nodes in upper left region (III)
00397  void node_update_III(const unsigned& t, AlgebraicNode*& node_pt);
00398 
00399  /// Update function for nodes in upper right region (IV)
00400  void node_update_IV(const unsigned& t,AlgebraicNode*& node_pt);
00401 
00402  /// Helper function
00403  void slanted_bound_up(const unsigned& t,const Vector<double>& zeta,
00404                        Vector<double>& r);
00405  
00406  /// Origin of the wall (stored explicitly for reference in
00407  /// algebraic node update -- it's also stored independently in 
00408  /// domain....)
00409  double X_0; 
00410 
00411  /// Length of the leaflet (stored explicitly for reference in
00412  /// algebraic node update -- it's also stored independently in 
00413  /// domain....)
00414  double Hleaflet; 
00415 
00416 };
00417 
00418 ///////////////////////////////////////////////////////////////////////////
00419 ///////////////////////////////////////////////////////////////////////////
00420 ///////////////////////////////////////////////////////////////////////////
00421 
00422 
00423 
00424 //===================================================================
00425 /// Refineable version of algebraic ChannelWithLeafletMesh
00426 //===================================================================
00427 template<class ELEMENT>
00428 class RefineableAlgebraicChannelWithLeafletMesh : 
00429 public RefineableQuadMesh<ELEMENT>,
00430  public virtual AlgebraicChannelWithLeafletMesh<ELEMENT>
00431 {
00432 
00433 public:
00434  
00435  ///\short Constructor: Pass pointer to GeomObject that represents the leaflet,
00436  /// the length of the domain to left and right of the leaflet, the
00437  /// height of the leaflet and the overall height of the channel,
00438  /// the number of element columns to the left and right of the leaflet,
00439  /// the number of rows of elements from the bottom of the channel to 
00440  /// the end of the leaflet, the number of rows of elements above the
00441  /// end of the leaflet. Timestepper defaults to Steady default 
00442  /// Timestepper defined in the Mesh base class
00443  RefineableAlgebraicChannelWithLeafletMesh(
00444   GeomObject* leaflet_pt, 
00445   const double& lleft,
00446   const double& lright,
00447   const double& hleaflet,
00448   const double& htot,
00449   const unsigned& nleft, const unsigned& nright,
00450   const unsigned& ny1, const unsigned&  ny2,
00451   TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) :
00452   ChannelWithLeafletMesh<ELEMENT>(leaflet_pt,lleft,lright,hleaflet,htot,
00453                                   nleft,nright,ny1,ny2,
00454                                   time_stepper_pt),
00455   AlgebraicChannelWithLeafletMesh<ELEMENT>(leaflet_pt,lleft,lright,hleaflet,
00456                                            htot,nleft,nright,ny1,ny2,
00457                                            time_stepper_pt)
00458   {
00459    // Build quadtree forest
00460    this->setup_quadtree_forest();
00461   }
00462  
00463  /// \short Update the node update data for specified node following 
00464  /// any mesh adapation
00465  void update_node_update(AlgebraicNode*& node_pt);
00466 
00467 };
00468 
00469 
00470 }
00471 
00472 #endif
channel_with_leaflet_mesh.template.h
