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 //Include guards
00029 #ifndef OOMPH_QUARTER_CIRCLE_SECTOR_DOMAIN_HEADER
00030 #define OOMPH_QUARTER_CIRCLE_SECTOR_DOMAIN_HEADER
00031 
00032 // Generic oomph-lib includes
00033 #include "../generic/quadtree.h"
00034 #include "../generic/domain.h"
00035 #include "../generic/geom_objects.h"
00036 
00037 namespace oomph
00038 {
00039 
00040 
00041 
00042 
00043 //=================================================================
00044 /// \short Circular sector as domain. Domain is bounded by 
00045 /// curved boundary which is represented by a GeomObject. Domain is 
00046 /// parametrised by three macro elements as shown here:
00047 /// \image html DomainWithMacroElementSketchAndCoords.gif
00048 //=================================================================
00049 class QuarterCircleSectorDomain : public Domain 
00050 { 
00051 
00052 public: 
00053 
00054 
00055 
00056  /// \short Constructor: Pass boundary object and start and end coordinates
00057  /// and fraction along boundary object where outer ring is divided.
00058  QuarterCircleSectorDomain(GeomObject* boundary_geom_object_pt, 
00059                            const double& xi_lo,  
00060                            const double& fract_mid,                    
00061                            const double& xi_hi) :
00062   Xi_lo(xi_lo), Fract_mid(fract_mid), Xi_hi(xi_hi), 
00063   Wall_pt(boundary_geom_object_pt), BL_squash_fct_pt(&default_BL_squash_fct)
00064   {
00065    // There are three macro elements
00066    unsigned nmacro=3;
00067 
00068    // Resize
00069    Macro_element_pt.resize(nmacro);
00070 
00071    // Create macro elements
00072    for (unsigned i=0;i<nmacro;i++)
00073     {
00074      Macro_element_pt[i]=new QMacroElement<2>(this,i);
00075     }
00076   }
00077 
00078 
00079  /// Broken copy constructor
00080  QuarterCircleSectorDomain(const QuarterCircleSectorDomain&) 
00081   { 
00082    BrokenCopy::broken_copy("QuarterCircleSectorDomain");
00083   } 
00084  
00085  /// Broken assignment operator
00086  void operator=(const QuarterCircleSectorDomain&) 
00087   {
00088    BrokenCopy::broken_assign("QuarterCircleSectorDomain");
00089   }
00090 
00091  
00092  /// Destructor: Kill macro elements
00093  ~QuarterCircleSectorDomain()
00094   {
00095    for (unsigned i=0;i<3;i++)
00096     {
00097      delete Macro_element_pt[i];
00098     }
00099   }
00100  
00101  /// \short Typedef for function pointer for function that squashes
00102  /// the outer two macro elements towards 
00103  /// the wall by mapping the input value of the "radial" macro element
00104  /// coordinate to the return value
00105  typedef double (*BLSquashFctPt)(const double& s);
00106 
00107  
00108  /// \short Function pointer for function that squashes
00109  /// the outer two macro elements towards 
00110  /// the wall by mapping the input value of the "radial" macro element
00111  /// coordinate to the return value
00112  BLSquashFctPt& bl_squash_fct_pt()
00113   {
00114    return BL_squash_fct_pt;
00115   }
00116 
00117 
00118  /// \short Function that squashes the outer two macro elements towards 
00119  /// the wall by mapping the input value of the "radial" macro element
00120  /// coordinate to the return value
00121  double s_squashed(const double& s)
00122   {
00123    return BL_squash_fct_pt(s);
00124   }
00125 
00126  /// \short Vector representation of the  i_macro-th macro element
00127  /// boundary i_direct (N/S/W/E) at time level t 
00128  /// (t=0: present; t>0: previous):
00129  /// f(s). Note that the local coordinate \b s is a 1D
00130  /// Vector rather than a scalar -- this is unavoidable because
00131  /// this function implements the pure virtual function in the
00132  /// Domain base class.
00133  void macro_element_boundary(const unsigned& t,
00134                              const unsigned& i_macro,
00135                              const unsigned& i_direct,
00136                              const Vector<double>& s,
00137                              Vector<double>& f);
00138 
00139 
00140 private:
00141 
00142  /// Lower limit for the (1D) coordinates along the wall
00143  double Xi_lo;
00144 
00145  /// Fraction along wall where outer ring is to be divided
00146  double Fract_mid;
00147 
00148  /// Upper limit for the (1D) coordinates along the wall
00149  double Xi_hi;
00150 
00151  /// Pointer to geometric object that represents the curved wall
00152  GeomObject* Wall_pt;
00153 
00154  /// \short Function pointer for function that squashes
00155  /// the outer two macro elements towards 
00156  /// the wall by mapping the input value of the "radial" macro element
00157  /// coordinate to the return value
00158  BLSquashFctPt BL_squash_fct_pt;
00159 
00160  /// \short Default for function that squashes
00161  /// the outer two macro elements towards 
00162  /// the wall by mapping the input value of the "radial" macro element
00163  /// coordinate to the return value: Identity.
00164  static double default_BL_squash_fct(const double& s)
00165   {
00166    return s;
00167   }
00168  
00169  /// \short Boundary of top left macro element zeta \f$ \in [-1,1] \f$
00170  void r_top_left_N(const unsigned& t, const Vector<double>& zeta, 
00171                    Vector<double>& f);
00172 
00173  /// \short Boundary of top left macro element zeta \f$ \in [-1,1] \f$
00174  void r_top_left_W(const unsigned& t, const Vector<double>& zeta, 
00175                    Vector<double>& f);
00176 
00177  /// \short Boundary of top left macro element zeta \f$ \in [-1,1] \f$
00178  void r_top_left_S(const unsigned& t, const Vector<double>& zeta, 
00179                    Vector<double>& f);
00180 
00181  /// \short Boundary of top left macro element zeta \f$ \in [-1,1] \f$
00182  void r_top_left_E(const unsigned& t, const Vector<double>& zeta, 
00183                           Vector<double>& f);
00184 
00185  /// \short Boundary of bottom right macro element zeta \f$ \in [-1,1] \f$
00186  void r_bot_right_N(const unsigned& t, const Vector<double>& zeta, 
00187                     Vector<double>& f);
00188 
00189  /// \short Boundary of bottom right macro element zeta \f$ \in [-1,1] \f$
00190  void r_bot_right_W(const unsigned& t, const Vector<double>& zeta, 
00191                     Vector<double>& f);
00192 
00193  /// \short Boundary of bottom right macro element zeta \f$ \in [-1,1] \f$
00194  void r_bot_right_S(const unsigned& t, const Vector<double>& zeta, 
00195                     Vector<double>& f);
00196 
00197  /// \short Boundary of bottom right macro element zeta \f$ \in [-1,1] \f$
00198  void r_bot_right_E(const unsigned& t, const Vector<double>& zeta, 
00199                            Vector<double>& f);
00200 
00201  /// \short Boundary of central box macro element zeta \f$ \in [-1,1] \f$
00202  void r_centr_N(const unsigned& t, const Vector<double>& zeta, 
00203                 Vector<double>& f);
00204 
00205  /// \short Boundary of central box macro element zeta \f$ \in [-1,1] \f$
00206  void r_centr_E(const unsigned& t, const Vector<double>& zeta, 
00207                 Vector<double>& f);
00208 
00209  /// \short Boundary of central box macro element zeta \f$ \in [-1,1] \f$
00210  void r_centr_S(const unsigned& t, const Vector<double>& zeta, 
00211                 Vector<double>& f);
00212 
00213  /// \short Boundary of central box macro element zeta \f$ \in [-1,1] \f$
00214  void r_centr_W(const unsigned& t, const Vector<double>& zeta, 
00215                        Vector<double>& f);
00216 
00217 
00218 };
00219 
00220 
00221 /////////////////////////////////////////////////////////////////////////
00222 /////////////////////////////////////////////////////////////////////////
00223 /////////////////////////////////////////////////////////////////////////
00224 
00225 
00226 
00227 //=================================================================
00228 /// Vector representation of the  imacro-th macro element
00229 /// boundary idirect (N/S/W/E) at time level t (t=0: present; t>0: previous):
00230 /// f(s)
00231 //=================================================================
00232 void QuarterCircleSectorDomain::macro_element_boundary(
00233  const unsigned& t,
00234  const unsigned& imacro,
00235  const unsigned& idirect,
00236  const Vector<double>& s,
00237  Vector<double>& f)
00238 {
00239  
00240 
00241 #ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
00242    // Warn about time argument being moved to the front
00243    OomphLibWarning(
00244     "Order of function arguments has changed between versions 0.8 and 0.85",
00245     "QuarterCircleSectorDomain::macro_element_boundary(...)",
00246     OOMPH_EXCEPTION_LOCATION);
00247 #endif
00248 
00249  // Which macro element?
00250  // --------------------
00251  switch(imacro)
00252   {
00253 
00254    using namespace QuadTreeNames;
00255 
00256    // Macro element 0: Central box
00257   case 0:
00258    
00259    // Which direction?
00260    if (idirect==N)
00261     {
00262      r_centr_N(t,s,f);
00263     }
00264    else if (idirect==S)
00265     {
00266      r_centr_S(t,s,f);
00267     }
00268    else if (idirect==W)
00269     {
00270      r_centr_W(t,s,f);
00271     }
00272    else if (idirect==E)
00273     {
00274      r_centr_E(t,s,f);
00275     }
00276    else
00277     {
00278      std::ostringstream error_stream;
00279      error_stream << "idirect is " << idirect 
00280                   << " not one of N, S, E, W" <<  std::endl;
00281 
00282      throw OomphLibError(
00283       error_stream.str(),
00284       "QuarterCircleSectorDomain::macro_element_boundary()",
00285       OOMPH_EXCEPTION_LOCATION);
00286     }
00287    
00288    break;
00289 
00290    // Macro element 1: Bottom right
00291   case 1:
00292    
00293    // Which direction?
00294    if (idirect==N)
00295     {
00296      r_bot_right_N(t,s,f);
00297     }
00298    else if (idirect==S)
00299     {
00300      r_bot_right_S(t,s,f);
00301     }
00302    else if (idirect==W)
00303     {
00304      r_bot_right_W(t,s,f);
00305     }
00306    else if (idirect==E)
00307     {
00308      r_bot_right_E(t,s,f);
00309     }
00310    else
00311     {
00312      std::ostringstream error_stream;
00313      error_stream << "idirect is " << idirect 
00314                   << " not one of N, S, E, W" <<  std::endl;
00315      
00316      throw OomphLibError(
00317       error_stream.str(),
00318       "QuarterCircleSectorDomain::macro_element_boundary()",
00319       OOMPH_EXCEPTION_LOCATION);
00320     }
00321    
00322    break;   
00323    
00324    // Macro element 2:Top left
00325   case 2:
00326    
00327    // Which direction?
00328    if (idirect==N)
00329     {
00330      r_top_left_N(t,s,f);
00331     }
00332    else if (idirect==S)
00333     {
00334      r_top_left_S(t,s,f);
00335     }
00336    else if (idirect==W)
00337     {
00338      r_top_left_W(t,s,f);
00339     }
00340    else if (idirect==E)
00341     {
00342      r_top_left_E(t,s,f);
00343     }
00344    else
00345     {
00346      std::ostringstream error_stream;
00347      error_stream << "idirect is " << idirect 
00348                   << " not one of N, S, E, W" <<  std::endl;
00349 
00350      throw OomphLibError(
00351       error_stream.str(),
00352       "QuarterCircleSectorDomain::macro_element_boundary()",
00353       OOMPH_EXCEPTION_LOCATION);
00354     }
00355    
00356    break;
00357    
00358   default:
00359    
00360    // Error
00361    std::ostringstream error_stream;
00362    error_stream << "Wrong imacro " << imacro << std::endl;
00363 
00364    throw OomphLibError(error_stream.str(),
00365                        "QuarterCircleSectorDomain::macro_element_boundary()",
00366                        OOMPH_EXCEPTION_LOCATION);
00367   }
00368  
00369 }
00370 
00371 
00372 
00373 
00374 //=================================================================
00375 /// Northern edge of top left macro element \f$ s \in [-1,1] \f$
00376 //=================================================================
00377 void QuarterCircleSectorDomain::r_top_left_N(const unsigned& t, 
00378                                       const Vector<double>& s,
00379                                       Vector<double>& f)
00380  {                                   
00381   Vector<double> x(1);
00382 
00383   // Coordinate along wall
00384   x[0]=Xi_hi+(Fract_mid*(Xi_hi-Xi_lo)-Xi_hi)*0.5*(s[0]+1.0);
00385 
00386   Wall_pt->position(t,x,f);
00387  }                                                  
00388 
00389 
00390 
00391 
00392 //=================================================================
00393 /// Western edge of top left macro element \f$s \in [-1,1] \f$
00394 //=================================================================
00395  void QuarterCircleSectorDomain::r_top_left_W(const unsigned& t, 
00396                                       const Vector<double>& s, 
00397                                       Vector<double>& f)
00398  {                                   
00399   Vector<double> x(1);
00400       
00401   // Top left corner
00402   Vector<double> r_top(2);
00403   x[0]=Xi_hi;
00404      
00405   Wall_pt->position(t,x,r_top);
00406 
00407   f[0]=0.0;
00408   f[1]=0.5*r_top[1]*(1.0+s_squashed(0.5*(s[0]+1.0)));
00409  }                                                  
00410 
00411 
00412 
00413 //=================================================================
00414 /// Southern edge of top left macro element \f$ s \in [-1,1] \f$
00415 //=================================================================
00416  void QuarterCircleSectorDomain::r_top_left_S(const unsigned& t, 
00417                                       const Vector<double>& s, 
00418                                       Vector<double>& f)
00419  {                                   
00420   Vector<double> x(1);
00421 
00422   // Top left corner
00423   Vector<double> r_top(2);
00424   x[0]=Xi_hi;
00425   
00426   Wall_pt->position(t,x,r_top);
00427 
00428 
00429   // Bottom right corner
00430   Vector<double> r_bot(2);
00431   x[0]=0.0;
00432 
00433   Wall_pt->position(t,x,r_bot);
00434 
00435   f[0]=0.5*r_bot[0]*0.5*(s[0]+1.0);
00436   f[1]=0.5*r_top[1];
00437 
00438  }                                                  
00439 
00440 
00441 
00442 //=================================================================
00443 /// Eastern edge of top left macro element \f$ s \in [-1,1] \f$
00444 //=================================================================
00445 void QuarterCircleSectorDomain::r_top_left_E(const unsigned& t, 
00446                                              const Vector<double>& s, 
00447                                              Vector<double>& f)
00448 {                                   
00449  Vector<double> x(1);
00450 
00451   // Top left corner
00452   Vector<double> r_top(2);
00453   x[0]=Xi_hi;
00454         
00455   Wall_pt->position(t,x,r_top);
00456 
00457   // Bottom right corner
00458   Vector<double> r_bot(2);
00459   x[0]=Xi_lo;
00460 
00461   Wall_pt->position(t,x,r_bot);
00462                    
00463   // Halfway along wall
00464   Vector<double> r_half(2);
00465   x[0]=Xi_lo+Fract_mid*(Xi_hi-Xi_lo);
00466 
00467   Wall_pt->position(t,x,r_half);
00468 
00469   f[0]=0.5*(r_bot[0]+s_squashed(0.5*(s[0]+1.0))*(2.0*r_half[0]-r_bot[0]));
00470   f[1]=0.5*(r_top[1]+s_squashed(0.5*(s[0]+1.0))*(2.0*r_half[1]-r_top[1]));
00471 
00472  }                                                  
00473 
00474 
00475 
00476 
00477 
00478 
00479 //=================================================================
00480 /// Northern edge of bottom right macro element
00481 //=================================================================
00482  void QuarterCircleSectorDomain::r_bot_right_N(const unsigned& t, 
00483                                        const Vector<double>& s,
00484                                        Vector<double>& f)
00485  {                                   
00486   r_top_left_E(t,s,f); 
00487  }                                                  
00488 
00489 //=================================================================
00490 /// Western edge of bottom right macro element
00491 //=================================================================
00492 void QuarterCircleSectorDomain::r_bot_right_W(const unsigned& t, 
00493                                        const Vector<double>& s, 
00494                                        Vector<double>& f)
00495  {                                   
00496   Vector<double> x(1);
00497       
00498   // Top left corner
00499   Vector<double> r_top(2);
00500   x[0]=Xi_hi;
00501 
00502   Wall_pt->position(t,x,r_top);
00503 
00504   // Bottom right corner
00505   Vector<double> r_bot(2);
00506   x[0]=Xi_lo;  
00507        
00508   Wall_pt->position(t,x,r_bot);
00509 
00510   f[0]=0.5*r_bot[0];
00511   f[1]=0.5*r_top[1]*0.5*(s[0]+1.0);
00512  }                                                  
00513 
00514 //=================================================================
00515 /// Southern edge of bottom right macro element
00516 //=================================================================
00517  void QuarterCircleSectorDomain::r_bot_right_S(const unsigned& t, 
00518                                        const Vector<double>& s, 
00519                                        Vector<double>& f)
00520  {                                   
00521   Vector<double> x(1);
00522 
00523   // Bottom right corner
00524   Vector<double> r_bot(2);
00525   x[0]=Xi_lo;
00526   Wall_pt->position(t,x,r_bot);
00527 
00528 
00529   f[0]=0.5*r_bot[0]*(1.0+s_squashed(0.5*(s[0]+1.0)));
00530   f[1]=0.0;
00531 
00532  }                                                  
00533 
00534 //=================================================================
00535 /// Eastern edge of bottom right macro element
00536 //=================================================================
00537 void QuarterCircleSectorDomain::r_bot_right_E(const unsigned& t, 
00538                                               const Vector<double>& s,
00539                                               Vector<double>& f)
00540 {           
00541  Vector<double> x(1);
00542  
00543  // Coordinate along wall
00544  x[0]=Xi_lo+(Fract_mid*(Xi_hi-Xi_lo)-Xi_lo)*(s[0]+1.0)*0.5;         
00545  
00546  Wall_pt->position(t,x,f);
00547  
00548 }                                                  
00549 
00550 
00551 //=================================================================
00552 /// Northern edge of central box
00553 //=================================================================
00554  void QuarterCircleSectorDomain::r_centr_N(const unsigned& t, 
00555                                    const Vector<double>& s,
00556                                    Vector<double>& f)
00557  {                                   
00558   r_top_left_S(t,s,f); 
00559  }                                                  
00560 
00561 
00562 
00563 
00564 //=================================================================
00565 /// Eastern edge of central box
00566 //=================================================================
00567 void QuarterCircleSectorDomain::r_centr_E(const unsigned& t, 
00568                                           const Vector<double>& s,
00569                                    Vector<double>& f)
00570  {                                   
00571   r_bot_right_W(t,s,f); 
00572  }
00573 
00574 
00575 
00576 
00577 //=================================================================
00578 /// Southern edge of central box
00579 //=================================================================
00580 void QuarterCircleSectorDomain::r_centr_S(const unsigned& t, 
00581                                   const Vector<double>& s,
00582                                   Vector<double>& f)
00583 {                                   
00584  Vector<double> x(1);
00585  
00586  // Bottom right corner
00587  Vector<double> r_bot(2);
00588  x[0]=Xi_lo;
00589  Wall_pt->position(t,x,r_bot);
00590  
00591  f[0]=0.5*r_bot[0]*0.5*(s[0]+1.0);
00592  f[1]=0.0; 
00593 }                                                  
00594 
00595 
00596 
00597 
00598 //=================================================================
00599 /// Western  edge of central box
00600 //=================================================================
00601 void QuarterCircleSectorDomain::r_centr_W(const unsigned& t, 
00602                                           const Vector<double>& s,
00603                                           Vector<double>& f)
00604 {                                   
00605  
00606  Vector<double> x(1);
00607  
00608  // Top left corner
00609  Vector<double> r_top(2);
00610  x[0]=Xi_hi;
00611  Wall_pt->position(t,x,r_top);
00612  
00613  f[0]=0.0;
00614  f[1]=0.5*r_top[1]*0.5*(s[0]+1.0);
00615  
00616 }                                                  
00617 
00618 
00619 
00620 }
00621 
00622 #endif
quarter_circle_sector_domain.h
