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 "map_matrix.h"
00029 #include "brick_mesh.h"
00030 
00031 namespace oomph
00032 {
00033 
00034 //================================================================
00035 /// Setup lookup schemes which establish which elements are located
00036 /// next to which boundaries (Doc to outfile if it's open).
00037 //================================================================
00038 void BrickMeshBase::setup_boundary_element_info(std::ostream &outfile)
00039 {
00040 
00041  // Doc?
00042  bool doc=false;
00043  if (outfile) doc=true;
00044  
00045  // Number of boundaries
00046  unsigned nbound=nboundary();
00047  
00048  // Wipe/allocate storage for arrays
00049  Boundary_element_pt.clear();
00050  Face_index_at_boundary.clear();
00051  Boundary_element_pt.resize(nbound);
00052  Face_index_at_boundary.resize(nbound);
00053  
00054  // Temporary vector of sets of pointers to elements on the boundaries: 
00055  // set_of_boundary_element_pt[i] is the set of pointers to all
00056  // elements that have nodes on boundary i. 
00057  Vector<std::set<FiniteElement*> > set_of_boundary_element_pt;
00058  set_of_boundary_element_pt.resize(nbound);
00059  
00060  // Matrix map for working out the fixed local coord for elements on boundary
00061  MapMatrixMixed<unsigned,FiniteElement*,Vector<int>* > 
00062   boundary_identifier;
00063 
00064  // Tmp container to store pointers to tmp vectors so they can be deleted
00065  Vector<Vector<int>*> tmp_vect_pt; 
00066 
00067  // Loop over elements
00068  //-------------------
00069  unsigned nel=nelement();
00070  for (unsigned e=0;e<nel;e++)
00071   {
00072    // Get pointer to element
00073    FiniteElement* fe_pt=finite_element_pt(e);
00074 
00075    if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
00076 
00077    // Loop over the element's nodes and find out which boundaries they're on
00078    // ----------------------------------------------------------------------
00079    unsigned nnode_1d=fe_pt->nnode_1d();
00080 
00081    // Loop over nodes in order
00082    for (unsigned i0=0;i0<nnode_1d;i0++)
00083     {
00084      for (unsigned i1=0;i1<nnode_1d;i1++)
00085       {
00086        for (unsigned i2=0;i2<nnode_1d;i2++)
00087         {
00088          // Local node number
00089          unsigned j=i0+i1*nnode_1d+i2*nnode_1d*nnode_1d;
00090          
00091          // Get pointer to set of boundaries that this
00092          // node lives on
00093          std::set<unsigned>* boundaries_pt=0;
00094          fe_pt->node_pt(j)->get_boundaries_pt(boundaries_pt);
00095          
00096          // If the node lives on some boundaries....
00097          if (boundaries_pt!=0)
00098           {
00099            for (std::set<unsigned>::iterator it=boundaries_pt->begin();
00100                 it!=boundaries_pt->end();++it)
00101             {
00102 
00103              // What's the boundary?
00104              unsigned boundary_id=*it;
00105 
00106              // Add pointer to finite element to set for the appropriate 
00107              // boundary -- storage in set makes sure we don't count elements
00108              // multiple times
00109              set_of_boundary_element_pt[boundary_id].insert(fe_pt);
00110              
00111              // For the current element/boundary combination, create
00112              // a vector that stores an indicator which element boundaries
00113              // the node is located (boundary_identifier=-/+1 for nodes
00114              // on the left/right boundary; boundary_identifier=-/+2 for nodes
00115              // on the lower/upper boundary. We determine these indices
00116              // for all corner nodes of the element and add them to a vector
00117              // to a vector. This allows us to decide which face of the element
00118              // coincides with the boundary since the (brick!) element must 
00119              // have exactly four corner nodes on the boundary.
00120              if (boundary_identifier(boundary_id,fe_pt)==0)
00121               {
00122                Vector<int>* tmp_pt=new Vector<int>;
00123                tmp_vect_pt.push_back(tmp_pt); 
00124                boundary_identifier(boundary_id,fe_pt)=tmp_pt;
00125               }
00126             
00127              // Are we at a corner node?
00128              if (((i0==0)||(i0==nnode_1d-1))&&((i1==0)||(i1==nnode_1d-1))
00129                  &&((i2==0)||(i2==nnode_1d-1)))
00130               {
00131                // Create index to represent position relative to s_0
00132                (*boundary_identifier(boundary_id,fe_pt)).
00133                 push_back(1*(2*i0/(nnode_1d-1)-1));               
00134                
00135                // Create index to represent position relative to s_1
00136                (*boundary_identifier(boundary_id,fe_pt)).
00137                 push_back(2*(2*i1/(nnode_1d-1)-1));
00138  
00139                // Create index to represent position relative to s_2
00140                (*boundary_identifier(boundary_id,fe_pt)).
00141                 push_back(3*(2*i2/(nnode_1d-1)-1));
00142               }
00143             }
00144           }
00145         }
00146       }     
00147     }  
00148   }
00149    
00150 
00151  // Now copy everything across into permanent arrays
00152  //-------------------------------------------------
00153 
00154 
00155  // Loop over boundaries
00156  //---------------------
00157  for (unsigned i=0;i<nbound;i++)
00158   {
00159 
00160    // Number of elements on this boundary (currently stored in a set)
00161    unsigned nel=set_of_boundary_element_pt[i].size();
00162    
00163    // Loop over elements on given boundary
00164    typedef std::set<FiniteElement*>::iterator IT;
00165    for (IT it=set_of_boundary_element_pt[i].begin();
00166         it!=set_of_boundary_element_pt[i].end();
00167         it++)
00168     {
00169      // Push back into permanent array
00170      Boundary_element_pt[i].push_back(*it);
00171     }
00172    
00173    // Allocate storage for the face identifiers
00174    Face_index_at_boundary[i].resize(nel);
00175    
00176    // Loop over elements on this boundary
00177    //-------------------------------------
00178    for (unsigned e=0;e<nel;e++)
00179     {
00180      // Recover pointer to element
00181      FiniteElement* fe_pt=Boundary_element_pt[i][e];
00182      
00183      // Initialise count for boundary identiers (-3,-2,-1,1,2,3)
00184      std::map<int,unsigned> count;
00185      
00186      // Loop over coordinates
00187      for (int ii=0;ii<3;ii++)
00188       {
00189        // Loop over upper/lower end of coordinates
00190        for (int sign=-1;sign<3;sign+=2)
00191         {
00192          count[(ii+1)*sign]=0;
00193         }
00194       }
00195      
00196      // Loop over boundary indicators for this element/boundary
00197      unsigned n_indicators=(*boundary_identifier(i,fe_pt)).size();
00198      for (unsigned j=0;j<n_indicators;j++)
00199       {
00200        count[(*boundary_identifier(i,fe_pt))[j] ]++;
00201       }
00202 
00203      // Determine the correct boundary indicator by checking that it 
00204      // occurs four times (since four corner nodes of the element's boundary
00205      // need to be located on the domain boundary
00206      int indicator=-10;
00207 
00208      //Check that we're finding exactly one boundary indicator
00209      bool found=false;
00210 
00211      // Loop over coordinates
00212      for (int ii=0;ii<3;ii++)
00213       {
00214        // Loop over upper/lower end of coordinates
00215        for (int sign=-1;sign<3;sign+=2)
00216         {
00217          if (count[(ii+1)*sign]==4)
00218           {
00219            // Check that we haven't found multiple boundaries
00220            if (found)
00221             {
00222              throw OomphLibError(
00223               "Trouble: Multiple boundary identifiers!\n",
00224               "BrickMeshBase::setup_boundary_element_info()",
00225               OOMPH_EXCEPTION_LOCATION);
00226             }
00227            found=true;
00228            indicator=(ii+1)*sign;
00229           }
00230         }
00231       }
00232      
00233      // Check if we've found one boundary
00234      if (!found)
00235       {
00236        std::ostringstream error_stream;
00237        error_stream
00238         << "Failed to find a boundary identifier.\n\n" 
00239         << "NOTE: This function is not implemented in full generality\n "
00240         << "      and this error may be thrown erroneously:\n\n"
00241         << "Cases where fewer than four vertex nodes of a brick element \n"
00242         << "that has some nodes located on a mesh boundary are located \n"
00243         << "on the same boundary MAY in fact arise legally, for instance \n"
00244         << "if only an edge of an element is located on the mesh boundary.\n"
00245         << "However, there is currently a problem, either with our code\n"
00246         << "or with the intel compiler, in which this error also gets \n"
00247         << "thrown when there ARE four vertex nodes on the same domain boundary!\n"
00248         << "This problem seems to occur only with 3d brick meshes in which\n"
00249         << "Nodes were initially created as \"ordinary\" Nodes and then\n"
00250         << "converted to BoundaryNodes, using the \n"
00251         << "Mesh::convert_to_boundary_node(...) function. \n\n"
00252         << "There are currently two options:\n"
00253         << " (1) If you have developed a new mesh in which this error is\n"
00254         << "     thrown erroneously, rewrite this function: Elements should\n"
00255         << "     only be added to the boundary element vector IF it has been\n"
00256         << "     established that four of their vertex nodes are located on a mesh\n"
00257         << "     boundary. We will implement this rewrite ourselves \n"
00258         << "     once we have tracked down the problem (or established that \n"
00259         << "     it is caused by the intel compiler -- wishful thinking?)\n"
00260         << " (2) If you use one of oomph-lib's existing 3D meshes then \n"
00261         << "     the error should not occur. In this case a workaround\n"
00262         << "     is to avoid the conversion to boundary nodes, by\n"
00263         << "     compiling the meshes (or indeed the entire library)\n "
00264         << "     with the C++ compiler flag CONVERT_BOUNDARY_NODE_IS_BROKEN.\n"
00265         << "     With this flag all Nodes in the affected meshes are created\n"
00266         << "     as BoundaryNodes and the problem is avoided -- at the\n"
00267         << "     expense of a slight additional memory overhead.\n";
00268        throw OomphLibError(error_stream.str(),
00269                            "BrickMeshBase::setup_boundary_element_info()",
00270                            OOMPH_EXCEPTION_LOCATION);
00271       }
00272      
00273 
00274      // Now convert boundary indicator into information required
00275      // for FaceElements
00276      switch (indicator)
00277       {
00278       case -3:
00279    
00280        // s_2 is fixed at -1.0:
00281        Face_index_at_boundary[i][e] = -3;
00282        break;
00283       
00284       case -2:
00285    
00286        // s_1 is fixed at -1.0:
00287        Face_index_at_boundary[i][e] = -2;
00288        break;
00289        
00290       case -1:
00291        
00292        // s_0 is fixed at -1.0:
00293        Face_index_at_boundary[i][e] = -1;
00294        break;
00295        
00296        
00297       case 1:
00298        
00299        // s_0 is fixed at 1.0:
00300        Face_index_at_boundary[i][e] = 1;
00301        break;
00302        
00303       case 2:
00304        
00305        // s_1 is fixed at 1.0:
00306        Face_index_at_boundary[i][e] = 2;
00307        break;
00308     
00309       case 3:
00310        
00311        // s_2 is fixed at 1.0:
00312        Face_index_at_boundary[i][e] = 3;
00313        break;
00314     
00315    
00316       default:
00317 
00318        throw OomphLibError("Never get here",
00319                            "BrickMeshBase::setup_boundary_element_info()",
00320                            OOMPH_EXCEPTION_LOCATION);
00321       }
00322      
00323     }
00324   }
00325  
00326 
00327  // Doc?
00328  //-----
00329  if (doc)
00330   {
00331    // Loop over boundaries
00332    for (unsigned i=0;i<nbound;i++)
00333     {
00334      unsigned nel=Boundary_element_pt[i].size();
00335      outfile << "Boundary: " << i
00336              << " is adjacent to " << nel
00337              << " elements" << std::endl;
00338      
00339      // Loop over elements on given boundary
00340      for (unsigned e=0;e<nel;e++)
00341       {
00342        FiniteElement* fe_pt=Boundary_element_pt[i][e];
00343        outfile << "Boundary element:" <<  fe_pt
00344                << " Face index along boundary is "
00345                <<  Face_index_at_boundary[i][e] << std::endl;
00346       }
00347     }
00348   }
00349  
00350 
00351  // Lookup scheme has now been setup yet
00352  Lookup_for_elements_next_boundary_is_setup=true;
00353 
00354 
00355  // Cleanup temporary vectors
00356  unsigned n=tmp_vect_pt.size();
00357  for (unsigned i=0;i<n;i++)
00358   {
00359    delete tmp_vect_pt[i];
00360   }
00361 
00362 }
00363 
00364 
00365 }
brick_mesh.cc
