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_QUARTER_TUBE_MESH_TEMPLATE_CC
00029 #define OOMPH_QUARTER_TUBE_MESH_TEMPLATE_CC
00030 
00031 #include "quarter_tube_mesh.template.h"
00032 
00033 
00034 namespace oomph
00035 {
00036 
00037 
00038 //====================================================================
00039 /// \short Constructor for deformable quarter tube mesh class. 
00040 /// The domain is specified by the GeomObject that 
00041 /// identifies boundary 3. Pass pointer to geometric object that
00042 /// specifies the wall, start and end coordinates on the 
00043 /// geometric object, and the fraction along
00044 /// which the dividing line is to be placed, and the timestepper.
00045 /// Timestepper defaults to Static dummy timestepper.
00046 //====================================================================
00047 template <class ELEMENT>
00048 QuarterTubeMesh<ELEMENT>::QuarterTubeMesh(GeomObject* wall_pt,
00049                                           const Vector<double>& xi_lo, 
00050                                           const double& fract_mid,
00051                                           const Vector<double>& xi_hi,
00052                                           const unsigned& nlayer,
00053                                           TimeStepper* time_stepper_pt) :  
00054  Wall_pt(wall_pt), Xi_lo(xi_lo), Fract_mid(fract_mid), Xi_hi(xi_hi)
00055 {
00056  
00057  // Build macro element-based domain
00058  Domain_pt = new QuarterTubeDomain(wall_pt,xi_lo,fract_mid,xi_hi,nlayer);
00059 
00060  //Set the number of boundaries
00061  set_nboundary(5);
00062 
00063  //We have only bothered to parametrise boundary 3
00064  Boundary_coordinate_exists[3] = true;
00065 
00066  // Allocate the store for the elements
00067  unsigned nelem=3*nlayer;
00068  Element_pt.resize(nelem);
00069 
00070  // Create  dummy element so we can determine the number of nodes 
00071  ELEMENT* dummy_el_pt=new ELEMENT;
00072 
00073  // Read out the number of linear points in the element
00074  unsigned n_p = dummy_el_pt->nnode_1d();
00075 
00076  // Kill the element
00077  delete dummy_el_pt;
00078 
00079  // Can now allocate the store for the nodes 
00080  unsigned nnodes_total=
00081   (n_p*n_p+(n_p-1)*n_p+(n_p-1)*(n_p-1))*(1+nlayer*(n_p-1));
00082  Node_pt.resize(nnodes_total);
00083 
00084 
00085  Vector<double> s(3);
00086  Vector<double> r(3);
00087 
00088  //Storage for the intrinsic boundary coordinate
00089  Vector<double> zeta(2);
00090 
00091  // Loop over elements and create all nodes
00092  for (unsigned ielem=0;ielem<nelem;ielem++)
00093   {
00094 
00095    // Create element
00096    Element_pt[ielem] = new ELEMENT;
00097 
00098    // Loop over rows in z/s_2-direction
00099    for (unsigned i2=0;i2<n_p;i2++)
00100     {
00101 
00102      // Loop over rows in y/s_1-direction
00103      for (unsigned i1=0;i1<n_p;i1++)
00104       {
00105        
00106        // Loop over rows in x/s_0-direction
00107        for (unsigned i0=0;i0<n_p;i0++)
00108         {
00109 
00110          // Local node number
00111          unsigned jnod_local=i0+i1*n_p+i2*n_p*n_p;
00112 
00113          // Create the node 
00114          Node* node_pt=finite_element_pt(ielem)->
00115           construct_node(jnod_local,time_stepper_pt);
00116 
00117          //Set the position of the node from macro element mapping
00118          s[0]=-1.0+2.0*double(i0)/double(n_p-1);  
00119          s[1]=-1.0+2.0*double(i1)/double(n_p-1);  
00120          s[2]=-1.0+2.0*double(i2)/double(n_p-1); 
00121          Domain_pt->macro_element_pt(ielem)->macro_map(s,r);
00122 
00123          node_pt->x(0) = r[0];
00124          node_pt->x(1) = r[1];
00125          node_pt->x(2) = r[2];
00126 
00127         }
00128       }
00129     }
00130 
00131   }
00132 
00133  // Initialise number of global nodes
00134  unsigned node_count=0;
00135 
00136  // Tolerance for node killing:
00137  double node_kill_tol=1.0e-12;
00138 
00139  // Check for error in node killing
00140  bool stopit=false;
00141 
00142  // Loop over elements
00143  for (unsigned ielem=0;ielem<nelem;ielem++)
00144   {
00145 
00146    // Which layer?
00147    unsigned ilayer=unsigned(ielem/3);
00148 
00149    // Which macro element?
00150    switch(ielem%3)
00151     {
00152      
00153      // Macro element 0: Central box
00154      //-----------------------------
00155     case 0:
00156    
00157 
00158      // Loop over rows in z/s_2-direction
00159      for (unsigned i2=0;i2<n_p;i2++)
00160       {
00161        
00162        // Loop over rows in y/s_1-direction
00163        for (unsigned i1=0;i1<n_p;i1++)
00164         {
00165          
00166          // Loop over rows in x/s_0-direction
00167          for (unsigned i0=0;i0<n_p;i0++)
00168           {
00169            
00170            // Local node number
00171            unsigned jnod_local=i0+i1*n_p+i2*n_p*n_p;
00172            
00173            // Has the node been killed?
00174            bool killed=false;
00175 
00176            // First layer of all nodes in s_2 direction gets killed
00177            // and re-directed to nodes in previous element layer
00178            if ((i2==0)&&(ilayer>0))
00179             {
00180              
00181              // Neighbour element
00182              unsigned ielem_neigh=ielem-3;
00183              
00184              // Node in neighbour element
00185              unsigned i0_neigh=i0;
00186              unsigned i1_neigh=i1;
00187              unsigned i2_neigh=n_p-1;
00188              unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
00189              
00190              // Check:
00191              for (unsigned i=0;i<3;i++)
00192               {
00193                 double error=std::abs(
00194                      finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
00195                      finite_element_pt(ielem_neigh)->
00196                      node_pt(jnod_local_neigh)->x(i));
00197                  if (error>node_kill_tol)
00198                  {
00199                   oomph_info << "Error in node killing for i " 
00200                        << i << " " << error << std::endl;
00201                   stopit=true;
00202                  }
00203               }
00204              
00205              // Kill node
00206              delete finite_element_pt(ielem)->node_pt(jnod_local);  
00207              killed=true;
00208              
00209              // Set pointer to neighbour:
00210              finite_element_pt(ielem)->node_pt(jnod_local)=
00211               finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);  
00212              
00213             }
00214     
00215            // No duplicate node: Copy across to mesh
00216            if (!killed)
00217             {
00218              Node_pt[node_count]=finite_element_pt(ielem)->node_pt(jnod_local);
00219 
00220              // Set boundaries:
00221 
00222              // Back: Boundary 0
00223              if ((i2==0)&&(ilayer==0))
00224               {
00225                this->convert_to_boundary_node(Node_pt[node_count]);
00226                add_boundary_node(0,Node_pt[node_count]);
00227               }
00228 
00229              // Front: Boundary 4
00230              if ((i2==n_p-1)&&(ilayer==nlayer-1))
00231               {
00232                this->convert_to_boundary_node(Node_pt[node_count]);
00233                add_boundary_node(4,Node_pt[node_count]);
00234               }
00235 
00236              // Left symmetry plane: Boundary 1
00237              if (i0==0)
00238               {
00239                this->convert_to_boundary_node(Node_pt[node_count]);
00240                add_boundary_node(1,Node_pt[node_count]);
00241               }
00242 
00243              // Bottom symmetry plane: Boundary 2
00244              if (i1==0)
00245               {
00246                this->convert_to_boundary_node(Node_pt[node_count]);
00247                add_boundary_node(2,Node_pt[node_count]);
00248               }
00249 
00250              // Increment node counter
00251              node_count++;
00252             }
00253            
00254 
00255           }
00256         }
00257       }
00258      
00259      
00260      break;
00261      
00262      // Macro element 1: Lower right box
00263      //---------------------------------
00264     case 1: 
00265 
00266 
00267      // Loop over rows in z/s_2-direction
00268      for (unsigned i2=0;i2<n_p;i2++)
00269       {
00270        
00271        // Loop over rows in y/s_1-direction
00272        for (unsigned i1=0;i1<n_p;i1++)
00273         {
00274          
00275          // Loop over rows in x/s_0-direction
00276          for (unsigned i0=0;i0<n_p;i0++)
00277           {
00278            
00279            // Local node number
00280            unsigned jnod_local=i0+i1*n_p+i2*n_p*n_p;
00281            
00282            // Has the node been killed?
00283            bool killed=false;
00284 
00285            // First layer of all nodes in s_2 direction gets killed
00286            // and re-directed to nodes in previous element layer
00287            if ((i2==0)&&(ilayer>0))
00288             {
00289              
00290              // Neighbour element
00291              unsigned ielem_neigh=ielem-3;
00292              
00293              // Node in neighbour element
00294              unsigned i0_neigh=i0;
00295              unsigned i1_neigh=i1;
00296              unsigned i2_neigh=n_p-1;
00297              unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
00298              
00299 
00300              // Check:
00301              for (unsigned i=0;i<3;i++)
00302               {
00303                 double error=std::abs(
00304                      finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
00305                      finite_element_pt(ielem_neigh)->
00306                      node_pt(jnod_local_neigh)->x(i));
00307                 if (error>node_kill_tol)
00308                  {
00309                   oomph_info << "Error in node killing for i " 
00310                        << i << " " << error << std::endl;
00311                   stopit=true;
00312                  }
00313               }
00314              
00315              // Kill node
00316              delete finite_element_pt(ielem)->node_pt(jnod_local);  
00317              killed=true;
00318              
00319              // Set pointer to neighbour:
00320              finite_element_pt(ielem)->node_pt(jnod_local)=
00321               finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);  
00322              
00323             }
00324            // Not in first layer:
00325            else
00326             {
00327              
00328              // Duplicate node: kill and set pointer to central element
00329              if (i0==0)
00330               {
00331                
00332                // Neighbour element
00333                unsigned ielem_neigh=ielem-1;
00334                
00335                // Node in neighbour element
00336                unsigned i0_neigh=n_p-1;
00337                unsigned i1_neigh=i1;
00338                unsigned i2_neigh=i2;
00339                unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
00340                
00341 
00342                // Check:
00343                for (unsigned i=0;i<3;i++)
00344                 {
00345                  double error=std::abs(
00346                   finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
00347                   finite_element_pt(ielem_neigh)->
00348                   node_pt(jnod_local_neigh)->x(i));
00349                if (error>node_kill_tol)
00350                   {
00351                    oomph_info << "Error in node killing for i " 
00352                         << i << " " << error << std::endl;
00353                    stopit=true;
00354                   }
00355                 }
00356                
00357                // Kill node
00358                delete finite_element_pt(ielem)->node_pt(jnod_local);  
00359                killed=true;
00360                
00361                // Set pointer to neighbour:
00362                finite_element_pt(ielem)->node_pt(jnod_local)=
00363                 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);  
00364                
00365               }
00366             }
00367 
00368            // No duplicate node: Copy across to mesh
00369            if (!killed)
00370             {
00371              Node_pt[node_count]=finite_element_pt(ielem)->node_pt(jnod_local);
00372 
00373              // Set boundaries:
00374              
00375              // Back: Boundary 0
00376              if ((i2==0)&&(ilayer==0))
00377               {
00378                this->convert_to_boundary_node(Node_pt[node_count]);
00379                add_boundary_node(0,Node_pt[node_count]);
00380               }
00381 
00382              // Front: Boundary 4
00383              if ((i2==n_p-1)&&(ilayer==nlayer-1))
00384               {
00385                this->convert_to_boundary_node(Node_pt[node_count]);
00386                add_boundary_node(4,Node_pt[node_count]);
00387               }
00388 
00389              // Bottom symmetry plane: Boundary 2
00390              if (i1==0)
00391               {
00392                this->convert_to_boundary_node(Node_pt[node_count]);
00393                add_boundary_node(2,Node_pt[node_count]);
00394               }
00395 
00396              // Tube wall: Boundary 3
00397              if (i0==n_p-1)
00398               {
00399                this->convert_to_boundary_node(Node_pt[node_count]);
00400                add_boundary_node(3,Node_pt[node_count]);
00401 
00402 
00403                // Get axial boundary coordinate
00404                zeta[0]=Xi_lo[0]+
00405                 (double(ilayer)+double(i2)/double(n_p-1))*
00406                 (Xi_hi[0]-Xi_lo[0])/double(nlayer);
00407 
00408                // Get azimuthal boundary coordinate
00409                zeta[1]=Xi_lo[1]+
00410                 double(i1)/double(n_p-1)*0.5*(Xi_hi[1]-Xi_lo[1]);
00411 
00412                Node_pt[node_count]->set_coordinates_on_boundary(3,zeta);
00413                 
00414               }
00415 
00416              // Increment node counter
00417              node_count++;
00418             }
00419            
00420           }
00421         }
00422       }
00423 
00424      break;
00425 
00426 
00427      // Macro element 2: Top left box
00428      //--------------------------------
00429     case 2: 
00430 
00431      // Loop over rows in z/s_2-direction
00432      for (unsigned i2=0;i2<n_p;i2++)
00433       {
00434        
00435        // Loop over rows in y/s_1-direction
00436        for (unsigned i1=0;i1<n_p;i1++)
00437         {
00438          
00439          // Loop over rows in x/s_0-direction
00440          for (unsigned i0=0;i0<n_p;i0++)
00441           {
00442            
00443            // Local node number
00444            unsigned jnod_local=i0+i1*n_p+i2*n_p*n_p;
00445            
00446            // Has the node been killed?
00447            bool killed=false;
00448            
00449            // First layer of all nodes in s_2 direction gets killed
00450            // and re-directed to nodes in previous element layer
00451            if ((i2==0)&&(ilayer>0))
00452             {
00453              
00454              // Neighbour element
00455              unsigned ielem_neigh=ielem-3;
00456              
00457              // Node in neighbour element
00458              unsigned i0_neigh=i0;
00459              unsigned i1_neigh=i1;
00460              unsigned i2_neigh=n_p-1;
00461              unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
00462              
00463              // Check:
00464              for (unsigned i=0;i<3;i++)
00465               {
00466                double error=std::abs(
00467                 finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
00468                 finite_element_pt(ielem_neigh)->
00469                 node_pt(jnod_local_neigh)->x(i));
00470                if (error>node_kill_tol)
00471                 {
00472                  oomph_info << "Error in node killing for i " 
00473                       << i << " " << error << std::endl;
00474                  stopit=true;
00475                 }
00476               }
00477              
00478              // Kill node
00479              delete finite_element_pt(ielem)->node_pt(jnod_local);  
00480              killed=true;
00481              
00482              // Set pointer to neighbour:
00483              finite_element_pt(ielem)->node_pt(jnod_local)=
00484               finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);  
00485              
00486             }
00487            // Not in first layer:
00488            else
00489             {
00490              
00491              // Duplicate node: kill and set pointer to node in bottom right
00492              // element
00493              if (i0==n_p-1)
00494               {
00495                
00496                // Neighbour element
00497                unsigned ielem_neigh=ielem-1;
00498                
00499                // Node in neighbour element
00500                unsigned i0_neigh=i1;
00501                unsigned i1_neigh=n_p-1;
00502                unsigned i2_neigh=i2;
00503                unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
00504                
00505                // Check:
00506                for (unsigned i=0;i<3;i++)
00507                 {
00508                  double error=std::abs(
00509                   finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
00510                   finite_element_pt(ielem_neigh)->
00511                   node_pt(jnod_local_neigh)->x(i));
00512                  if (error>node_kill_tol)
00513                   {
00514                    oomph_info << "Error in node killing for i " 
00515                         << i << " " << error << std::endl;
00516                    stopit=true;
00517                   }
00518                 }
00519                
00520                // Kill node
00521                delete finite_element_pt(ielem)->node_pt(jnod_local);  
00522                killed=true;
00523 
00524                // Set pointer to neighbour:
00525                finite_element_pt(ielem)->node_pt(jnod_local)=
00526                 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);  
00527              
00528               }
00529 
00530 
00531              // Duplicate node: kill and set pointer to central element
00532              if ((i1==0)&&(i0!=n_p-1))
00533               {
00534                
00535                // Neighbour element
00536                unsigned ielem_neigh=ielem-2;
00537                
00538                // Node in neighbour element
00539                unsigned i0_neigh=i0;
00540                unsigned i1_neigh=n_p-1;
00541                unsigned i2_neigh=i2;
00542                unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p+i2_neigh*n_p*n_p;
00543                
00544                // Check:
00545                for (unsigned i=0;i<3;i++)
00546                 {
00547                  double error=std::abs(
00548                   finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
00549                   finite_element_pt(ielem_neigh)->
00550                   node_pt(jnod_local_neigh)->x(i));
00551                  if (error>node_kill_tol) 
00552                   {
00553                    oomph_info << "Error in node killing for i " 
00554                         << i << " " << error << std::endl;
00555                    stopit=true;
00556                   }
00557                 }
00558 
00559                // Kill node
00560                delete finite_element_pt(ielem)->node_pt(jnod_local);  
00561                killed=true;
00562 
00563                // Set pointer to neighbour:
00564                finite_element_pt(ielem)->node_pt(jnod_local)=
00565                 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);  
00566              
00567               }
00568              
00569              // No duplicate node: Copy across to mesh
00570              if (!killed)
00571               {
00572                Node_pt[node_count]=finite_element_pt(ielem)->
00573                 node_pt(jnod_local);
00574 
00575                // Set boundaries:
00576                
00577                // Back: Boundary 0
00578                if ((i2==0)&&(ilayer==0))
00579                 {
00580                  this->convert_to_boundary_node(Node_pt[node_count]);
00581                  add_boundary_node(0,Node_pt[node_count]);
00582                 }
00583 
00584                // Front: Boundary 4
00585                if ((i2==n_p-1)&&(ilayer==nlayer-1))
00586                 {
00587                  this->convert_to_boundary_node(Node_pt[node_count]);
00588                  add_boundary_node(4,Node_pt[node_count]);
00589                 }
00590                
00591                // Left symmetry plane: Boundary 1
00592                if (i0==0)
00593                 {
00594                  this->convert_to_boundary_node(Node_pt[node_count]);
00595                  add_boundary_node(1,Node_pt[node_count]);
00596                 }
00597 
00598 
00599                // Tube wall: Boundary 3
00600                if (i1==n_p-1)
00601                 {
00602                  this->convert_to_boundary_node(Node_pt[node_count]);
00603                  add_boundary_node(3,Node_pt[node_count]);
00604 
00605 
00606                  // Get axial boundary coordinate
00607                  zeta[0]=Xi_lo[0]+
00608                   (double(ilayer)+double(i2)/double(n_p-1))*
00609                   (Xi_hi[0]-Xi_lo[0])/double(nlayer);
00610                  
00611                  // Get azimuthal boundary coordinate
00612                  zeta[1]=Xi_hi[1]-
00613                   double(i0)/double(n_p-1)*0.5*(Xi_hi[1]-Xi_lo[1]);
00614                  
00615                  Node_pt[node_count]->set_coordinates_on_boundary(3,zeta);
00616 
00617                 }
00618 
00619                // Increment node counter
00620                node_count++;
00621 
00622               }
00623              
00624             }
00625           }
00626         }
00627       }
00628 
00629      break;
00630     }
00631   }
00632  
00633  // Terminate if there's been an error
00634  if (stopit)
00635   {
00636    std::ostringstream error_stream;
00637    error_stream << "Error in killing nodes\n"
00638                 << "The most probable cause is that the domain is not\n"
00639                 << "compatible with the mesh.\n"
00640                 << "For the QuarterTubeMesh, the domain must be\n"
00641                 << "topologically consistent with a quarter tube with a\n"
00642                 << "non-curved centreline.\n";
00643    throw OomphLibError(error_stream.str(),
00644                        "QuatrerTubeMesh::QuarterTubeMesh()",
00645                        OOMPH_EXCEPTION_LOCATION);
00646   }
00647 
00648  // Setup boundary element lookup schemes
00649  setup_boundary_element_info();
00650 
00651 }
00652 
00653 ///////////////////////////////////////////////////////////////////////
00654 ///////////////////////////////////////////////////////////////////////
00655 // Algebraic-mesh-version of RefineableQuarterTubeMesh
00656 ///////////////////////////////////////////////////////////////////////
00657 ///////////////////////////////////////////////////////////////////////
00658 
00659 
00660 //======================================================================
00661 /// Setup algebraic node update data, based on 3 regions, each
00662 /// undergoing a different node update strategy. These regions are
00663 /// defined by the three MacroElements in each of the nlayer slices
00664 /// of the QuarterTubeDomain used to build the mesh.
00665 /// The Mesh is suspended from the `wall' GeomObject pointed to
00666 /// by wall_pt. The lower right edge of the mesh is located at the
00667 /// wall's coordinate xi[1]==xi_lo[1], the upper left edge at
00668 /// xi[1]=xi_hi[1], i.e. a view looking down the tube length.
00669 /// The dividing line between the two outer regions is located
00670 /// at the fraction fract_mid between these two coordinates.
00671 /// Node updating strategy:
00672 /// - the starting cross sectional shape along the tube length is
00673 ///   assumed to be uniform
00674 /// - the cross sectional shape of the central region remains
00675 ///   rectangular; the position of its top right corner is located
00676 ///   at a fixed fraction of the starting width and height of the
00677 ///   domain measured at xi[1]==xi_lo[1] and xi[1]==xi_hi[1]
00678 ///   respectively. Nodes in this region are located at fixed
00679 ///   horizontal and vertical fractions of the region.
00680 /// - Nodes in the two outer regions (bottom right and top left)
00681 ///   are located on straight lines running from the edges of the
00682 ///   central region to the outer wall.
00683 //======================================================================
00684 template <class ELEMENT>
00685 void AlgebraicRefineableQuarterTubeMesh<ELEMENT>::
00686 setup_algebraic_node_update()
00687 {
00688 
00689 #ifdef PARANOID
00690  /// Pointer to first algebraic element in central region
00691  AlgebraicElementBase* algebraic_element_pt=
00692   dynamic_cast<AlgebraicElementBase*>(Mesh::element_pt(0));
00693 
00694  if (algebraic_element_pt==0)
00695   {
00696    std::ostringstream error_message;
00697    error_message <<
00698     "Element in AlgebraicRefineableQuarterTubeMesh must be\n ";
00699    error_message << "derived from AlgebraicElementBase\n";
00700    error_message<< "but it is of type:  "
00701                 << typeid(Mesh::element_pt(0)).name() << std::endl;
00702    std::string function_name =
00703     "AlgebraicRefineableQuarterTubeMesh::";
00704    function_name += "setup_algebraic_node_update()";
00705    throw OomphLibError(error_message.str(),function_name,
00706                        OOMPH_EXCEPTION_LOCATION);
00707   }
00708 #endif
00709 
00710  // Find number of nodes in an element from the zeroth element
00711  unsigned nnodes_3d = Mesh::finite_element_pt(0)->nnode();
00712 
00713  // also find number of nodes in 1d line and 2d slice
00714  unsigned nnodes_1d = Mesh::finite_element_pt(0)->nnode_1d();
00715  unsigned nnodes_2d = nnodes_1d*nnodes_1d;
00716 
00717  // find node number of a top-left and a bottom-right node in an element
00718  // (orientation: looking down tube)
00719  unsigned tl_node = nnodes_2d-nnodes_1d;
00720  unsigned br_node = nnodes_1d-1;
00721 
00722  // find x & y distances to top-right node in element 0 - this is the same
00723  // node as the top-left node of element 1
00724  double x_c_element= Mesh::finite_element_pt(1)->node_pt(tl_node)->x(0);
00725  double y_c_element= Mesh::finite_element_pt(1)->node_pt(tl_node)->x(1);
00726 
00727  // Get x-distance to bottom-right edge of wall, i.e. coord of node
00728  // at bottom-right of bottom-right of element 1
00729  double x_wall= Mesh::finite_element_pt(1)->node_pt(br_node)->x(0);
00730 
00731  // Get y-distance to top-left edge of wall, i.e. coord of node
00732  // at top-left of element 2
00733  double y_wall= Mesh::finite_element_pt(2)->node_pt(tl_node)->x(1);
00734 
00735  // Establish fractional widths in central region
00736  Lambda_x = Centre_box_size*x_c_element/x_wall;
00737  Lambda_y = Centre_box_size*y_c_element/y_wall;
00738 
00739  // how many elements are there?
00740  unsigned nelements = Mesh::nelement();
00741 
00742  // loop through the elements
00743  for (unsigned e=0; e<nelements; e++)
00744   {
00745    // get pointer to element
00746    FiniteElement* el_pt=Mesh::finite_element_pt(e);
00747 
00748    // set region id
00749    unsigned region_id = e%3;
00750 
00751    // find the first node for which to set up node update info - must
00752    // bypass the first nnodes_2d nodes after the first 3 elements
00753    unsigned nstart=nnodes_2d;
00754    if (e<3)
00755     {
00756      nstart=0;
00757     }
00758 
00759    // loop through the nodes,
00760    for (unsigned n=nstart; n<nnodes_3d; n++)
00761     {
00762      // find z coordinate of node
00763      // NOTE: to implement axial spacing replace z by z_spaced where
00764      // z_spaced=axial_spacing_fct(z) when finding the GeomObjects
00765      // and local coords below 
00766      // BUT store z as the third reference value since this is the value
00767      // required by update_node_update()
00768      double z = el_pt->node_pt(n)->x(2);
00769 
00770      // Find wall GeomObject and the GeomObject local coordinates
00771      // at bottom-right edge of wall
00772      Vector<double> xi(2);
00773      xi[0] = z;
00774      xi[1] = QuarterTubeMesh<ELEMENT>::Xi_lo[1];
00775      GeomObject* obj_br_pt;
00776      Vector<double> s_br(2);
00777      this->Wall_pt->locate_zeta(xi,obj_br_pt,s_br);
00778 
00779      // Find wall GeomObject/local coordinate
00780      // at top-left edge of wall
00781      xi[1] = QuarterTubeMesh<ELEMENT>::Xi_hi[1];
00782      GeomObject* obj_tl_pt;
00783      Vector<double> s_tl(2);
00784      this->Wall_pt->locate_zeta(xi,obj_tl_pt,s_tl);
00785 
00786      // Element 0: central region
00787      //--------------------------
00788      if (region_id==0)
00789       {
00790        // Nodal coordinates in x and y direction
00791        double x=el_pt->node_pt(n)->x(0);
00792        double y=el_pt->node_pt(n)->x(1);
00793 
00794        // The update function requires two geometric objects
00795        Vector<GeomObject*> geom_object_pt(2);
00796 
00797        // Wall GeomObject at bottom right end of wall mesh:
00798        geom_object_pt[0] = obj_br_pt;
00799 
00800        // Wall GeomObject at top left end of wall mesh:
00801        geom_object_pt[1] = obj_tl_pt;
00802 
00803        // The update function requires seven parameters:
00804        Vector<double> ref_value(7);
00805 
00806        // First reference value: fractional x-position inside region
00807        ref_value[0]=x/x_c_element;
00808 
00809        // Second cfractional y-position inside region
00810        ref_value[1]=y/y_c_element;
00811 
00812        // Third reference value: initial z coordinate of node
00813        ref_value[2]=z;
00814 
00815        // Fourth and fifth reference values:
00816        // local coordinate in GeomObject at bottom-right of wall.
00817        // Note: must recompute this reference for new nodes
00818        // since values interpolated from parent nodes will
00819        // be wrong (this is true for all local coordinates
00820        // within GeomObjects)
00821        ref_value[3]=s_br[0];
00822        ref_value[4]=s_br[1];
00823 
00824        // Sixth and seventh reference values:
00825        // local coordinate in GeomObject at top-left of wall.
00826        ref_value[5]=s_tl[0];
00827        ref_value[6]=s_tl[1];
00828 
00829        // Setup algebraic update for node: Pass update information
00830        static_cast<AlgebraicNode*>(el_pt->node_pt(n))->
00831         add_node_update_info(
00832          Central_region,    // enumerated ID
00833          this,              // mesh
00834          geom_object_pt,    // vector of geom objects
00835          ref_value);        // vector of ref. values
00836       }
00837 
00838      // Element 1: Bottom right region
00839      //-------------------------------
00840      else if (region_id==1)
00841       {
00842        // Fractional distance between nodes
00843        double fract = 1.0/double(nnodes_1d-1);
00844 
00845        // Fraction in the s_0-direction
00846        double rho_0 = fract * double(n%nnodes_1d);
00847 
00848        // Fraction in the s_1-direction
00849        double rho_1 = fract * double((n/nnodes_1d)%nnodes_1d);
00850 
00851        // Find coord of reference point on wall
00852        xi[1]=QuarterTubeMesh<ELEMENT>::Xi_lo[1]
00853         +rho_1*QuarterTubeMesh<ELEMENT>::Fract_mid*(
00854          QuarterTubeMesh<ELEMENT>::Xi_hi[1] -
00855          QuarterTubeMesh<ELEMENT>::Xi_lo[1]);
00856 
00857        // Identify wall GeomObject and local coordinate of
00858        // reference point in GeomObject
00859        GeomObject* obj_wall_pt;
00860        Vector<double> s_wall(2);
00861        this->Wall_pt->locate_zeta(xi,obj_wall_pt,s_wall);
00862 
00863        // The update function requires three geometric objects
00864        Vector<GeomObject*> geom_object_pt(3);
00865 
00866        // Wall element at bottom-right end of wall mesh:
00867        geom_object_pt[0] = obj_br_pt;
00868 
00869        // Wall element at top left end of wall mesh:
00870        geom_object_pt[1] = obj_tl_pt;
00871 
00872        // Wall element at that contians reference point:
00873        geom_object_pt[2] = obj_wall_pt;
00874 
00875        // The update function requires nine parameters:
00876        Vector<double> ref_value(9);
00877 
00878        // First reference value: fractional s0-position inside region
00879        ref_value[0]=rho_0;
00880 
00881        // Second reference value: fractional s1-position inside region
00882        ref_value[1]=rho_1;
00883 
00884        // Thrid reference value: initial z coord of node
00885        ref_value[2]=z;
00886 
00887        // Fourth and fifth reference values:
00888        //   local coordinates at bottom-right of wall.
00889        ref_value[3]=s_br[0];
00890        ref_value[4]=s_br[1];
00891 
00892        // Sixth and seventh reference values:
00893        //   local coordinates at top-left of wall.
00894        ref_value[5]=s_tl[0];
00895        ref_value[6]=s_tl[1];
00896 
00897        // Eighth and ninth reference values:
00898        //   local coordinate of reference point.
00899        ref_value[7]=s_wall[0];
00900        ref_value[8]=s_wall[1];
00901 
00902        // Setup algebraic update for node: Pass update information
00903        static_cast<AlgebraicNode*>(el_pt->node_pt(n))->
00904         add_node_update_info(
00905          Lower_right_region,     // enumerated ID
00906          this,                   // mesh
00907          geom_object_pt,         // vector of geom objects
00908          ref_value);             // vector of ref. vals
00909       }
00910 
00911      // Element 2: Top left region
00912      //---------------------------
00913      else if (region_id==2)
00914       {
00915        // Fractional distance between nodes
00916        double fract = 1.0/double(nnodes_1d-1);
00917 
00918        // Fraction in the s_0-direction
00919        double rho_0 = fract * double(n%nnodes_1d);
00920 
00921        // Fraction in the s_1-direction
00922        double rho_1 = fract * double((n/nnodes_1d)%nnodes_1d);
00923 
00924        // Find coord of reference point on wall
00925        xi[1]=QuarterTubeMesh<ELEMENT>::Xi_hi[1]
00926         +rho_0*(1.0-QuarterTubeMesh<ELEMENT>::Fract_mid)*
00927         (QuarterTubeMesh<ELEMENT>::Xi_lo[1]-
00928          QuarterTubeMesh<ELEMENT>::Xi_hi[1]);
00929 
00930        // Identify GeomObject and local coordinate at
00931        // reference point on wall
00932        GeomObject* obj_wall_pt;
00933        Vector<double> s_wall(2);
00934        this->Wall_pt->locate_zeta(xi,obj_wall_pt,s_wall);
00935 
00936        // The update function requires three geometric objects
00937        Vector<GeomObject*> geom_object_pt(3);
00938 
00939        // Wall element at bottom-right of wall mesh:
00940        geom_object_pt[0] = obj_br_pt;
00941 
00942        // Wall element at top-left of wall mesh:
00943        geom_object_pt[1] = obj_tl_pt;
00944 
00945        // Wall element contianing reference point:
00946        geom_object_pt[2] = obj_wall_pt;
00947 
00948        // The update function requires nine parameters:
00949        Vector<double> ref_value(9);
00950 
00951        // First reference value: fractional s0-position inside region
00952        ref_value[0]=rho_0;
00953 
00954        // Second reference value: fractional s1-position inside region
00955        ref_value[1]=rho_1;
00956 
00957        // Third reference value: initial z coord
00958        ref_value[2]=z;
00959 
00960        // Fourth and fifth reference values: 
00961        //   local coordinates in bottom-right GeomObject
00962        ref_value[3]=s_br[0];
00963        ref_value[4]=s_br[1];
00964 
00965        // Sixth and seventh reference values:
00966        //   local coordinates in top-left GeomObject
00967        ref_value[5]=s_tl[0];
00968        ref_value[6]=s_tl[1];
00969 
00970        // Eighth and ninth reference values: 
00971        //   local coordinates at reference point
00972        ref_value[7]=s_wall[0];
00973        ref_value[8]=s_wall[1];
00974 
00975        // Setup algebraic update for node: Pass update information
00976        static_cast<AlgebraicNode*>(el_pt->node_pt(n))->
00977         add_node_update_info(
00978          Upper_left_region,      // Enumerated ID
00979          this,                   // mesh
00980          geom_object_pt,         // vector of geom objects
00981          ref_value);             // vector of ref. vals
00982 
00983       }
00984     }
00985   }
00986 }
00987 
00988 //======================================================================
00989 /// \short Algebraic update function: Update in central region according
00990 /// to wall shape at time level t (t=0: present; t>0: previous)
00991 //======================================================================
00992 template <class ELEMENT>
00993 void AlgebraicRefineableQuarterTubeMesh<ELEMENT>::
00994 node_update_central_region(const unsigned& t, AlgebraicNode*& node_pt)
00995 {
00996  
00997 #ifdef PARANOID
00998  // We're updating the nodal positions (!) at time level t
00999  // and determine them by evaluating the wall GeomObject's
01000  // position at that time level. I believe this only makes sense
01001  // if the t-th history value in the positional timestepper
01002  // actually represents previous values (rather than some
01003  // generalised quantity). Hence if this function is called with
01004  // t>nprev_values(), we issue a warning and terminate the execution.
01005  // It *might* be possible that the code still works correctly
01006  // even if this condition is violated (e.g. if the GeomObject's
01007  // position() function returns the appropriate "generalised"
01008  // position value that is required by the timestepping scheme but it's
01009  // probably worth flagging this up and forcing the user to manually switch
01010  // off this warning if he/she is 100% sure that this is kosher.
01011  if (t>node_pt->position_time_stepper_pt()->nprev_values())
01012   {
01013    std::string error_message =
01014     "Trying to update the nodal position at a time level\n";
01015    error_message += "beyond the number of previous values in the nodes'\n";
01016    error_message += "position timestepper. This seems highly suspect!\n";
01017    error_message += "If you're sure the code behaves correctly\n";
01018    error_message += "in your application, remove this warning \n";
01019    error_message += "or recompile with PARNOID switched off.\n";
01020    
01021    std::string function_name =
01022     "AlgebraicRefineableQuarterTubeMesh::";
01023    function_name += "node_update_central_region()",
01024     throw OomphLibError(error_message,function_name,
01025                         OOMPH_EXCEPTION_LOCATION);
01026   }
01027 #endif
01028  
01029  // Extract references for update in central region by copy construction
01030  Vector<double> ref_value(node_pt->vector_ref_value(Central_region));
01031  
01032  // Extract geometric objects for update in central region by copy construction
01033  Vector<GeomObject*> geom_object_pt(node_pt->
01034                                     vector_geom_object_pt(Central_region));
01035  
01036  // First reference value: fractional x-position of node inside region
01037  double rho_x=ref_value[0];
01038  
01039  // Second reference value: fractional y-position of node inside region
01040  double rho_y=ref_value[1];
01041  
01042  // Wall position in bottom right corner:
01043  
01044  // Pointer to GeomObject
01045  GeomObject* obj_br_pt = geom_object_pt[0];
01046 
01047  // Local coordinate at bottom-right reference point:
01048  Vector<double> s_br(2);
01049  s_br[0]=ref_value[3];
01050  s_br[1]=ref_value[4];
01051 
01052  // Get wall position
01053  unsigned n_dim = obj_br_pt->ndim();
01054  Vector<double> r_br(n_dim);
01055  obj_br_pt->position(t,s_br,r_br);
01056 
01057  // Wall position in top left corner:
01058  
01059  // Pointer to GeomObject
01060  GeomObject* obj_tl_pt=geom_object_pt[1];
01061  
01062  // Local coordinate:
01063  Vector<double> s_tl(2);
01064  s_tl[0]=ref_value[5];
01065  s_tl[1]=ref_value[6];
01066 
01067  // Get wall position
01068  Vector<double> r_tl(n_dim);
01069  obj_tl_pt->position(t,s_tl,r_tl);
01070  
01071  // Assign new nodal coordinate 
01072  node_pt->x(t,0)=r_br[0]*Lambda_x*rho_x;
01073  node_pt->x(t,1)=r_tl[1]*Lambda_y*rho_y;
01074  node_pt->x(t,2)=(r_br[2]+r_tl[2])*0.5;
01075 }
01076 
01077 
01078 //====================================================================
01079 /// Algebraic update function: Update in lower-right region according
01080 /// to wall shape at time level t (t=0: present; t>0: previous)
01081 //====================================================================
01082 template <class ELEMENT>
01083 void AlgebraicRefineableQuarterTubeMesh<ELEMENT>::
01084 node_update_lower_right_region(const unsigned& t, AlgebraicNode*& node_pt)
01085 {
01086 
01087 #ifdef PARANOID
01088  // We're updating the nodal positions (!) at time level t
01089  // and determine them by evaluating the wall GeomObject's
01090  // position at that gime level. I believe this only makes sense
01091  // if the t-th history value in the positional timestepper
01092  // actually represents previous values (rather than some
01093  // generalised quantity). Hence if this function is called with
01094  // t>nprev_values(), we issue a warning and terminate the execution.
01095  // It *might* be possible that the code still works correctly
01096  // even if this condition is violated (e.g. if the GeomObject's
01097  // position() function returns the appropriate "generalised"
01098  // position value that is required by the timestepping scheme but it's
01099  // probably worth flagging this up and forcing the user to manually switch
01100  // off this warning if he/she is 100% sure that this is kosher.
01101  if (t>node_pt->position_time_stepper_pt()->nprev_values())
01102   {
01103    std::string error_message =
01104     "Trying to update the nodal position at a time level";
01105    error_message += "beyond the number of previous values in the nodes'";
01106    error_message += "position timestepper. This seems highly suspect!";
01107    error_message += "If you're sure the code behaves correctly";
01108    error_message += "in your application, remove this warning ";
01109    error_message += "or recompile with PARNOID switched off.";
01110 
01111    std::string function_name =
01112     "AlgebraicRefineableQuarterTubeSectorMesh::";
01113    function_name += "node_update_lower_right_region()",
01114     throw OomphLibError(error_message,function_name,
01115                         OOMPH_EXCEPTION_LOCATION);
01116   }
01117 #endif
01118 
01119  // Extract references for update in lower-right region by copy construction
01120  Vector<double> ref_value(node_pt->vector_ref_value(Lower_right_region));
01121 
01122  // Extract geometric objects for update in lower-right region 
01123  // by copy construction
01124  Vector<GeomObject*> 
01125   geom_object_pt(node_pt->vector_geom_object_pt(Lower_right_region));
01126 
01127  // First reference value: fractional s0-position of node inside region
01128  double rho_0=ref_value[0];
01129 
01130  // Use s_squashed to modify fractional s0 position
01131  rho_0 = this->Domain_pt->s_squashed(rho_0);
01132 
01133  // Second reference value: fractional s1-position of node inside region
01134  double rho_1=ref_value[1];
01135 
01136  // Wall position in bottom right corner:
01137 
01138  // Pointer to GeomObject
01139  GeomObject* obj_br_pt=geom_object_pt[0];
01140 
01141  // Local coordinate
01142  Vector<double> s_br(2);
01143  s_br[0]=ref_value[3];
01144  s_br[1]=ref_value[4];
01145 
01146  // Eulerian dimension
01147  unsigned n_dim=obj_br_pt->ndim();
01148 
01149  // Get wall position
01150  Vector<double> r_br(n_dim);
01151  obj_br_pt->position(t,s_br,r_br);
01152 
01153  // Wall position in top left corner:
01154 
01155  // Pointer to GeomObject
01156  GeomObject* obj_tl_pt=geom_object_pt[1];
01157 
01158  // Local coordinate
01159  Vector<double> s_tl(2);
01160  s_tl[0]=ref_value[5];
01161  s_tl[1]=ref_value[6];
01162 
01163  Vector<double> r_tl(n_dim);
01164  obj_tl_pt->position(t,s_tl,r_tl);
01165 
01166  // Wall position at reference point:
01167 
01168  // Pointer to GeomObject
01169  GeomObject* obj_wall_pt=geom_object_pt[2];
01170 
01171  // Local coordinate
01172  Vector<double> s_wall(2);
01173  s_wall[0]=ref_value[7];
01174  s_wall[1]=ref_value[8];
01175 
01176  Vector<double> r_wall(n_dim);
01177  obj_wall_pt->position(t,s_wall,r_wall);
01178 
01179  // Position Vector to corner of the central region
01180  Vector<double> r_corner(n_dim);
01181  r_corner[0]=Lambda_x*r_br[0];
01182  r_corner[1]=Lambda_y*r_tl[1];
01183  r_corner[2]=(r_br[2]+r_tl[2])*0.5;
01184 
01185  // Position Vector to left edge of region
01186  Vector<double> r_left(n_dim);
01187  r_left[0]=r_corner[0];
01188  r_left[1]=rho_1*r_corner[1];
01189  r_left[2]=r_corner[2];
01190 
01191  // Assign new nodal coordinate
01192  node_pt->x(t,0)=r_left[0]+rho_0*(r_wall[0]-r_left[0]);
01193  node_pt->x(t,1)=r_left[1]+rho_0*(r_wall[1]-r_left[1]);
01194  node_pt->x(t,2)=r_left[2]+rho_0*(r_wall[2]-r_left[2]);
01195 
01196 }
01197 //====================================================================
01198 /// Algebraic update function: Update in upper left region according
01199 /// to wall shape at time level t (t=0: present; t>0: previous)
01200 //====================================================================
01201 template <class ELEMENT>
01202 void AlgebraicRefineableQuarterTubeMesh<ELEMENT>::
01203 node_update_upper_left_region(const unsigned& t, AlgebraicNode*& node_pt)
01204 {
01205 
01206 #ifdef PARANOID
01207  // We're updating the nodal positions (!) at time level t
01208  // and determine them by evaluating the wall GeomObject's
01209  // position at that gime level. I believe this only makes sense
01210  // if the t-th history value in the positional timestepper
01211  // actually represents previous values (rather than some
01212  // generalised quantity). Hence if this function is called with
01213  // t>nprev_values(), we issue a warning and terminate the execution.
01214  // It *might* be possible that the code still works correctly
01215  // even if this condition is violated (e.g. if the GeomObject's
01216  // position() function returns the appropriate "generalised"
01217  // position value that is required by the timestepping scheme but it's
01218  // probably worth flagging this up and forcing the user to manually switch
01219  // off this warning if he/she is 100% sure that this is kosher.
01220  if (t>node_pt->position_time_stepper_pt()->nprev_values())
01221   {
01222    std::string error_message =
01223     "Trying to update the nodal position at a time level";
01224    error_message += "beyond the number of previous values in the nodes'";
01225    error_message += "position timestepper. This seems highly suspect!";
01226    error_message += "If you're sure the code behaves correctly";
01227    error_message += "in your application, remove this warning ";
01228    error_message += "or recompile with PARNOID switched off.";
01229 
01230    std::string function_name =
01231     "AlgebraicRefineableQuarterTubeMesh::";
01232    function_name += "node_update_upper_left_region()";
01233 
01234    throw OomphLibError(error_message,function_name,
01235                        OOMPH_EXCEPTION_LOCATION);
01236   }
01237 #endif
01238 
01239  // Extract references for update in upper-left region by copy construction
01240  Vector<double> ref_value(node_pt->vector_ref_value(Upper_left_region));
01241 
01242  // Extract geometric objects for update in upper-left region 
01243  // by copy construction
01244  Vector<GeomObject*> 
01245   geom_object_pt(node_pt->vector_geom_object_pt(Upper_left_region));
01246 
01247  // First reference value: fractional s0-position of node inside region
01248  double rho_0=ref_value[0];
01249 
01250  // Second  reference value: fractional s1-position of node inside region
01251  double rho_1=ref_value[1];
01252 
01253  // Use s_squashed to modify fractional s1 position
01254  rho_1 = this->Domain_pt->s_squashed(rho_1);
01255 
01256  // Wall position in bottom right corner:
01257 
01258  // Pointer to GeomObject
01259  GeomObject* obj_br_pt=geom_object_pt[0];
01260 
01261  // Eulerian dimension
01262  unsigned n_dim=obj_br_pt->ndim();
01263 
01264  // Local coordinate
01265  Vector<double> s_br(2);
01266  s_br[0]=ref_value[3];
01267  s_br[1]=ref_value[4];
01268 
01269  // Get wall position
01270  Vector<double> r_br(n_dim);
01271  obj_br_pt->position(t,s_br,r_br);
01272 
01273  // Wall position in top left corner:
01274 
01275  // Pointer to GeomObject
01276  GeomObject* obj_tl_pt=node_pt->geom_object_pt(1);
01277 
01278  // Local coordinate
01279  Vector<double> s_tl(2);
01280  s_tl[0]=ref_value[5];
01281  s_tl[1]=ref_value[6];
01282 
01283  Vector<double> r_tl(n_dim);
01284  obj_tl_pt->position(t,s_tl,r_tl);
01285 
01286  // Wall position at reference point:
01287 
01288  // Pointer to GeomObject
01289  GeomObject* obj_wall_pt=node_pt->geom_object_pt(2);
01290 
01291  // Local coordinate
01292  Vector<double> s_wall(2);
01293  s_wall[0]=ref_value[7];
01294  s_wall[1]=ref_value[8];
01295 
01296  Vector<double> r_wall(n_dim);
01297  obj_wall_pt->position(t,s_wall,r_wall);
01298 
01299  // Position vector to corner of the central region
01300  Vector<double> r_corner(n_dim);
01301  r_corner[0]=Lambda_x*r_br[0];
01302  r_corner[1]=Lambda_y*r_tl[1];
01303  r_corner[2]=(r_br[2]+r_tl[2])*0.5;
01304 
01305  // Position Vector to top face of central region
01306  Vector<double> r_top(n_dim);
01307  r_top[0]=rho_0*r_corner[0];
01308  r_top[1]=r_corner[1];
01309  r_top[2]=r_corner[2];
01310 
01311  // Assign new nodal coordinate
01312  node_pt->x(t,0)=r_top[0]+rho_1*(r_wall[0]-r_top[0]);
01313  node_pt->x(t,1)=r_top[1]+rho_1*(r_wall[1]-r_top[1]);
01314  node_pt->x(t,2)=r_top[2]+rho_1*(r_wall[2]-r_top[2]);
01315 
01316 }
01317 
01318 
01319 //======================================================================
01320 /// Update algebraic update function for nodes in region defined by
01321 /// region_id.
01322 //======================================================================
01323 template <class ELEMENT>
01324 void AlgebraicRefineableQuarterTubeMesh<ELEMENT>::
01325 update_node_update_in_region(AlgebraicNode*& node_pt, int& region_id)
01326 {
01327  // Extract references by copy construction
01328  Vector<double> ref_value(node_pt->vector_ref_value(region_id));
01329 
01330  // Extract geometric objects for update by copy construction
01331  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt(region_id));
01332 
01333  // Now remove the update info to allow overwriting below
01334  node_pt->kill_node_update_info(region_id);
01335 
01336  // Fixed reference value: Start coordinate on wall
01337  double xi_lo = QuarterTubeMesh<ELEMENT>::Xi_lo[1];
01338 
01339  // Fixed reference value: Fractional position of dividing line
01340  double fract_mid = QuarterTubeMesh<ELEMENT>::Fract_mid;
01341 
01342  // Fixed reference value: End coordinate on wall
01343  double xi_hi = QuarterTubeMesh<ELEMENT>::Xi_hi[1];
01344 
01345  // get initial z-coordinate of node
01346  // NOTE: use modified z found using axial_spacing_fct() 
01347  // to implement axial spacing
01348  double z = ref_value[2];
01349 
01350  // Update reference to wall point in bottom right corner
01351  //------------------------------------------------------
01352 
01353  // Wall position in bottom right corner:
01354  Vector<double> xi(2);
01355  xi[0]=z;
01356  xi[1]=xi_lo;
01357 
01358  // Find corresponding wall element/local coordinate
01359  GeomObject* obj_br_pt;
01360  Vector<double> s_br(2);
01361  this->Wall_pt->locate_zeta(xi,obj_br_pt,s_br);
01362 
01363  // Wall element at bottom right end of wall mesh:
01364  geom_object_pt[0] = obj_br_pt;
01365 
01366  // Local coordinate in this wall element.
01367  ref_value[3]=s_br[0];
01368  ref_value[4]=s_br[1];
01369 
01370 
01371  // Update reference to wall point in upper left corner
01372  //----------------------------------------------------
01373 
01374  // Wall position in top left corner
01375  xi[1]=xi_hi;
01376 
01377  // Find corresponding wall element/local coordinate
01378  GeomObject* obj_tl_pt;
01379  Vector<double> s_tl(2);
01380  this->Wall_pt->locate_zeta(xi,obj_tl_pt,s_tl);
01381 
01382  // Wall element at top left end of wall mesh:
01383  geom_object_pt[1] = obj_tl_pt;
01384 
01385  // Local coordinate in this wall element.
01386  ref_value[5]=s_tl[0];
01387  ref_value[6]=s_tl[1];
01388 
01389  if (region_id!=Central_region)
01390   {
01391    // Update reference to reference point on wall
01392    //--------------------------------------------
01393 
01394    // Reference point on wall
01395    if (region_id==Lower_right_region)
01396     {
01397      // Second reference value: fractional s1-position of node inside region
01398      double rho_1=ref_value[1];
01399 
01400      // position of reference point
01401      xi[1]=xi_lo+rho_1*fract_mid*(xi_hi-xi_lo);
01402     }
01403    else // case of Upper_left region
01404     {
01405      // First reference value: fractional s0-position of node inside region
01406      double rho_0=ref_value[0];
01407 
01408      // position of reference point
01409      xi[1]=xi_hi-rho_0*(1.0-fract_mid)*(xi_hi-xi_lo);
01410 
01411     }
01412 
01413    // Identify wall element number and local coordinate of
01414    // reference point on wall
01415    GeomObject* obj_wall_pt;
01416    Vector<double> s_wall(2);
01417    this->Wall_pt->locate_zeta(xi,obj_wall_pt,s_wall);
01418 
01419    // Wall element at that contians reference point:
01420    geom_object_pt[2] = obj_wall_pt;
01421 
01422    // Local coordinate in this wall element.
01423    ref_value[7]=s_wall[0];
01424    ref_value[8]=s_wall[1];
01425   }
01426 
01427  // Setup algebraic update for node: Pass update information
01428  node_pt->add_node_update_info(
01429   region_id,              // Enumerated ID
01430   this,                   // mesh
01431   geom_object_pt,         // vector of geom objects
01432   ref_value);             // vector of ref. vals
01433 
01434 }
01435 
01436 
01437 }
01438 #endif
quarter_tube_mesh.template.cc
