quad_mesh.cc

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
#include <algorithm>
#include "map_matrix.h"
#include "quad_mesh.h"


namespace oomph
{


//================================================================
/// Setup lookup schemes which establish which elements are located
/// next to which boundaries (Doc to outfile if it's open).
//================================================================
void QuadMeshBase::setup_boundary_element_info(std::ostream &outfile)
{

 bool doc=false;
 if (outfile) doc=true;
 
 // Number of boundaries
 unsigned nbound=nboundary();
 
 // Wipe/allocate storage for arrays
 Boundary_element_pt.clear();
 Face_index_at_boundary.clear();
 Boundary_element_pt.resize(nbound);
 Face_index_at_boundary.resize(nbound);
 
 // Temporary vector of vectors to pointers to elements on the boundaries: 
 // This is not a set to ensure UNIQUE ordering
 Vector<Vector<FiniteElement*> > vector_of_boundary_element_pt;
 vector_of_boundary_element_pt.resize(nbound);
 
 // Matrix map for working out the fixed local coord for elements on boundary
 MapMatrixMixed<unsigned,FiniteElement*,Vector<int>* > 
  boundary_identifier;
 

 // Loop over elements
 //-------------------
 unsigned nel=nelement();
 for (unsigned e=0;e<nel;e++)
  {
   // Get pointer to element
   FiniteElement* fe_pt=finite_element_pt(e);
   
   if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
   
   // Only include 2D elements! Some meshes contain interface elements too.
   if (fe_pt->dim()==2)
    {
     // Loop over the element's nodes and find out which boundaries they're on
     // ----------------------------------------------------------------------
     unsigned nnode_1d=fe_pt->nnode_1d();
     
     // Loop over nodes in order
     for (unsigned i0=0;i0<nnode_1d;i0++)
      {
       for (unsigned i1=0;i1<nnode_1d;i1++)
        {
         // Local node number
         unsigned j=i0+i1*nnode_1d;
         
         // Get pointer to vector of boundaries that this
         // node lives on
         std::set<unsigned>* boundaries_pt=0;
         fe_pt->node_pt(j)->get_boundaries_pt(boundaries_pt);

         // If the node lives on some boundaries....
         if (boundaries_pt!=0)
          {
           // Loop over boundaries
           //unsigned nbound=(*boundaries_pt).size();
           for (std::set<unsigned>::iterator it=boundaries_pt->begin();
                it != boundaries_pt->end();++it)
            {
             // Add pointer to finite element to vector for the appropriate 
             // boundary

             // Does the pointer already exits in the vector
             Vector<FiniteElement*>::iterator b_el_it =
              std::find(vector_of_boundary_element_pt[*it].begin(),
                        vector_of_boundary_element_pt[*it].end(),
                        fe_pt);
             //Only insert if we have not found it (i.e. got to the end)
             if(b_el_it == vector_of_boundary_element_pt[*it].end())
              {
               vector_of_boundary_element_pt[*it].push_back(fe_pt);
              }

             // For the current element/boundary combination, create
             // a vector that stores an indicator which element boundaries
             // the node is located (boundary_identifier=-/+1 for nodes
             // on the left/right boundary; boundary_identifier=-/+2 for nodes
             // on the lower/upper boundary. We determine these indices
             // for all corner nodes of the element and add them to a vector
             // to a vector. This allows us to decide which face of the element
             // coincides with the boundary since the (quad!) element must 
             // have exactly two corner nodes on the boundary.
             if (boundary_identifier(*it,fe_pt)==0)
              {
               boundary_identifier(*it,fe_pt)=new Vector<int>;
              }
             
             // Are we at a corner node?
             if (((i0==0)||(i0==nnode_1d-1))&&((i1==0)||(i1==nnode_1d-1)))
              {
               // Create index to represent position relative to s_0
               (*boundary_identifier(*it,fe_pt)).
                push_back(1*(2*i0/(nnode_1d-1)-1));               
               
               // Create index to represent position relative to s_1
               (*boundary_identifier(*it,fe_pt)).
                push_back(2*(2*i1/(nnode_1d-1)-1));
              } 

            }
          }
         //else
         // {
         //  oomph_info << "...does not live on any boundaries " << std::endl;
         // }
  
        }
      }
    }
    //else
    //{
     //oomph_info << "Element " << e << " does not qualify" << std::endl;
    //}
  }

 // Now copy everything across into permanent arrays
 //-------------------------------------------------

 // Note: vector_of_boundary_element_pt contains all elements
 // that have (at least) one corner node on a boundary -- can't copy 
 // them across into Boundary_element_pt
 // yet because some of them might have only one node on the
 // the boundary in which case they don't qualify as 
 // boundary elements!

 // Loop over boundaries
 //---------------------
 for (unsigned i=0;i<nbound;i++)
  {
   // Number of elements on this boundary
   // nel is unused, so I've commented it out - RWhite.
//   unsigned nel=vector_of_boundary_element_pt[i].size();
    
   // Allocate storage for the face identifiers
   //Face_index_at_boundary[i].resize(nel);
   
   // Loop over elements that have at least one corner node on this boundary
   //-----------------------------------------------------------------------
   //unsigned e_count=0;
   typedef Vector<FiniteElement*>::iterator IT;
   for (IT it=vector_of_boundary_element_pt[i].begin();
        it!=vector_of_boundary_element_pt[i].end();
        it++)
    {    
     // Recover pointer to element
     FiniteElement* fe_pt=*it;

     // Initialise count for boundary identiers (-2,-1,1,2)
     std::map<int,unsigned> count;
     
     // Loop over coordinates
     for (unsigned ii=0;ii<2;ii++)
      {
       // Loop over upper/lower end of coordinates
       for (int sign=-1;sign<3;sign+=2)
        {
         count[(ii+1)*sign]=0;
        }
      }
     
     // Loop over boundary indicators for this element/boundary
     unsigned n_indicators=(*boundary_identifier(i,fe_pt)).size();
     for (unsigned j=0;j<n_indicators;j++)
      {
       count[(*boundary_identifier(i,fe_pt))[j] ]++;
      }
     delete boundary_identifier(i,fe_pt);

     // Determine the correct boundary indicator by checking that it 
     // occurs twice (since two corner nodes of the element's boundary
     // need to be located on the domain boundary
     int indicator=-10;

     //Check that we're finding exactly one boundary indicator
     //bool found=false;

     // Loop over coordinates
     for (unsigned ii=0;ii<2;ii++)
      {
       // Loop over upper/lower end of coordinates
       for (int sign=-1;sign<3;sign+=2)
        {
         //If an index occurs twice then that face is on the boundary
         //But we can have multiple faces on the same boundary id, so just add all the ones that we
         //have
         if (count[(ii+1)*sign]==2)
          {
           // Check that we haven't found multiple boundaries
           /*if (found)
            {
             std::string error_message=
              "Trouble: Multiple boundary identifiers!\n";
             error_message += 
              "Elements should only have at most 2 corner ";
             error_message += 
              "nodes on any one boundary.\n";

             throw OomphLibError(
              error_message,
              OOMPH_CURRENT_FUNCTION,
              OOMPH_EXCEPTION_LOCATION);
            }
            found=true;*/
           indicator=(ii+1)*sign;

           //Copy into the data structure
           Boundary_element_pt[i].push_back(*it);
           Face_index_at_boundary[i].push_back(indicator);
          }
        }
      }
     
     // Element has exactly two corner nodes on boundary
     /*if (found)
      {
       // Add to permanent storage
       Boundary_element_pt[i].push_back(fe_pt);

       // Now convert boundary indicator into information required
       // for FaceElements
       switch (indicator)
        {
         //South face
        case -2:
         
         // s_1 is fixed at -1.0:
         Face_index_at_boundary[i][e_count] = -2;
         break;
         
         //West face
        case -1:
         
         // s_0 is fixed at -1.0:
         Face_index_at_boundary[i][e_count] = -1;
         break;
         
         //East face
        case 1:
         
         // s_0 is fixed at 1.0:
         Face_index_at_boundary[i][e_count] = 1;
         break;
         
         //North face
        case 2:
         
         // s_1 is fixed at 1.0:
         Face_index_at_boundary[i][e_count] = 2;
         break;
         
        default:

         throw OomphLibError("Never get here",
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
                             }

       // Increment counter
       e_count++;
       }*/
     
    }
  }
 


 // Doc?
 //-----
 if (doc)
  {
   // Loop over boundaries
   for (unsigned i=0;i<nbound;i++)
    {
     unsigned nel=Boundary_element_pt[i].size();
     outfile << "Boundary: " << i
             << " is adjacent to " << nel
             << " elements" << std::endl;
     
     // Loop over elements on given boundary
     for (unsigned e=0;e<nel;e++)
      {
       FiniteElement* fe_pt=Boundary_element_pt[i][e];
       outfile << "Boundary element:" <<  fe_pt
               << " Face index on boundary is " 
               <<  Face_index_at_boundary[i][e] << std::endl;
      }
    }
  }
 

 // Lookup scheme has now been setup
 Lookup_for_elements_next_boundary_is_setup=true;


}

}