brick_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 "map_matrix.h"
#include "brick_mesh.h"

namespace oomph
{

//====================================================================
/// Helper namespace for generation of brick from tet mesh
//====================================================================
namespace BrickFromTetMeshHelper
{

 /// Tolerance for mismatch during setup of boundary coordinates
 double Face_position_tolerance=1.0e-12;

}



///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////


//================================================================
/// Setup lookup schemes which establish which elements are located
/// next to which boundaries (Doc to outfile if it's open).
//================================================================
void BrickMeshBase::setup_boundary_element_info(std::ostream &outfile)
{
 //Martina's fixed and commented version of the code.
 
 // Define variable doc and set the initial value to False.
 bool doc=false;
    
 // If file declared in outfile exists,set doc=true to enable writing to stream
 if (outfile) doc=true;
 
 // Number of boundaries. Gives the value assigned by the function nboundary()
 unsigned nbound=nboundary();

 if(doc)
   {
    outfile << "The number of boundaries is " << nbound << "\n";
   }

 // Wipe/allocate storage for arrays
 // Command clear will reomve all elements of vectors
 // Command resize will adapt pointer to correct boundary size. Internal data
 Boundary_element_pt.clear();
 Face_index_at_boundary.clear();
 Boundary_element_pt.resize(nbound);
 Face_index_at_boundary.resize(nbound);


 // Temporary vector of vectors of pointers to elements on the boundaries: 
 // vector_of_boundary_element_pt[i] is the vector of pointers to all
 // elements that have nodes on boundary i. 
 // This is not a set to ensure UNIQUE ordering on every processor
 // There are a number of elements with nodes on each boundary.
 // For each boundary we store the appropriate elements in a vector.
 // There is therefore a vector for each boundary, which we store in
 // another vector, i.e. a matrix (vector of vectors).
 // Then vector_of_boundary_element_pt[0] is a vector of pointers to elements with
 // nodes on boundary 0 and vector_of_boundary_element_pt[0][0] is the first
 // element with nodes on boundary 0.
 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
 // Calling pointer to FiniteElement and pointer to  Vector<int> defining
 // the matrix to identify each element on boundary
    MapMatrixMixed<unsigned,FiniteElement*,Vector<int>* >
    boundary_identifier;
 

 // Temporary container to store pointers to temporary vectors
 // so they can be deleted. Creating storage to store these
 // temporary vectors of vectors of pointers previously defined
 Vector<Vector<int>*> tmp_vect_pt; 
 
 // Loop over elements
 //-------------------
 unsigned nel=nelement();
 for (unsigned e=0;e<nel;e++)
  {
   // Get pointer to element
   // and put it in local storage fe_pt.
    FiniteElement* fe_pt=finite_element_pt(e);
  
    //print out values of all elements to doc
      if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;

   // Loop over the element's nodes and find out which boundaries they're on
   // ----------------------------------------------------------------------
   // Return number of nodes along one edge of the current element defined by fe_pt
   // hence, loop over nodes of each element in 3D-dimension
 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++)
      {
       for (unsigned i2=0;i2<nnode_1d;i2++)
        {
         // Local node number, which is an assumed ordering defined in our underlying
         // finite element scheme.
         unsigned j=i0+i1*nnode_1d+i2*nnode_1d*nnode_1d;
         
        // Get pointer to the set of boundaries that this node lives on
        // for each node reset boundaries_pt to 0,
        // create pointer to each node in each element
        // and give boundaries_pt a value
         std::set<unsigned>* boundaries_pt=0;
         fe_pt->node_pt(j)->get_boundaries_pt(boundaries_pt);
         
     // If the node lives on some boundaries....
         // If not equal to 0, node is on boundary
         // Loop through values of the current node stored in boundaries_pt
      if (boundaries_pt!=0)
          {
              
            // Loop over the entries in the set "it" (name we use to refer to the iterator)
           for (std::set<unsigned>::iterator it=boundaries_pt->begin();
                it!=boundaries_pt->end();++it)
            {

             // What's the boundary?
            // Define boundary_id to have the value pointed to by boundaries_pt
             unsigned boundary_id=*it;

             // Add pointer to finite element to vector for the appropriate 
             // boundary 

             // Does the pointer already exist 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 faces of the element
             // coincide with the boundary since each face of the element must 
             // have all four corner nodes on the boundary.
             if (boundary_identifier(boundary_id,fe_pt)==0)
              {
                  
                // Here we make our vector of integers and keep track of the pointer to it
                // The vector stores information about local node position
                // for each boundary element
               Vector<int>* tmp_pt=new Vector<int>;
                  
               // Add to the local scheme that will allow us to delete it later
               // at the very end of the function
               tmp_vect_pt.push_back(tmp_pt);
                  
               // Add the pointer to the storage scheme defined by boundary_id
               // and pointer to finite element
               boundary_identifier(boundary_id,fe_pt)=tmp_pt;
              }
            
             // Are we at a corner node? (local for each boundary element)
             // i0,i1,i2 represent the current 3D coordinates- which local corner node.
             if (((i0==0)||(i0==nnode_1d-1))&&((i1==0)||(i1==nnode_1d-1))
                 &&((i2==0)||(i2==nnode_1d-1)))
              {
               // Create index to represent position relative to s_0
               // s_0,s_1,s_2 are local 3D directions
               (*boundary_identifier(boundary_id,fe_pt)).
                push_back(1*(2*i0/(nnode_1d-1)-1));               
               
               // Create index to represent position relative to s_1
               (*boundary_identifier(boundary_id,fe_pt)).
                push_back(2*(2*i1/(nnode_1d-1)-1));
 
               // Create index to represent position relative to s_2
               (*boundary_identifier(boundary_id,fe_pt)).
                push_back(3*(2*i2/(nnode_1d-1)-1));
              }
            }
          }
        }
      }     
    }  
  }
   

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

 // Loop over boundaries
 //---------------------
 for (unsigned i=0;i<nbound;i++)
  {
   // Loop over elements on given boundary
   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 for each element
     FiniteElement* fe_pt=(*it);
     
     // Initialise count for boundary identities (-3,-2,-1,1,2,3)
     std::map<int,unsigned> count;
     
     // Loop over coordinates in 3D dimension
     for (int ii=0;ii<3;ii++)
      {
       // Loop over upper/lower end of coordinates
       // count -/+ for each direction separately
       for (int sign=-1;sign<3;sign+=2)
        {
         count[(ii+1)*sign]=0;
            
         // Initialise map of counts to 0 before loop
         // count boundary indicator for each element
         // and its nodes
        }
      }
     
     // Loop over boundary indicators for this element/boundary
     // count nodes for any one fixed boundary determined by
     // boundary indicator
     unsigned n_indicators=(*boundary_identifier(i,fe_pt)).size();
     for (unsigned j=0;j<n_indicators;j++)
      {
       count[(*boundary_identifier(i,fe_pt))[j] ]++;
      }
     
     // Determine the correct boundary indicator by checking that it 
     // occurs four times (since four corner nodes of the element's boundary
     // need to be located on the boundary domain)
     int indicator=-10;
     
     
     // Loop over coordinates
     for (int ii=0;ii<3;ii++)
      {
       // Loop over upper/lower end of coordinates
       for (int sign=-1;sign<3;sign+=2)
        {
         // If an index occurs four times then a face is on the boundary
         //But we can have multiple faces on the same boundary, so add all the ones that we find!
         if (count[(ii+1)*sign]==4)
          {
           indicator=(ii+1)*sign;
           
           // Copy into the member data structures
           Boundary_element_pt[i].push_back(*it);
           Face_index_at_boundary[i].push_back(indicator);
          }
        }
      }
    }
  }
 
 // 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 along boundary is "
               <<  Face_index_at_boundary[i][e] << std::endl;
      }
    }
  }
 

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


 // Cleanup temporary vectors
 unsigned n=tmp_vect_pt.size();
 for (unsigned i=0;i<n;i++)
  {
   delete tmp_vect_pt[i];
  }

 return;
 
 
 // Doc?
 /*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 of pointers to elements on the boundaries: 
 // vector_of_boundary_element_pt[i] is the vector of pointers to all
 // elements that have nodes on boundary i. 
 // This is not a set to ensure UNIQUE ordering on every processor
 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;

 // Tmp container to store pointers to tmp vectors so they can be deleted
 Vector<Vector<int>*> tmp_vect_pt; 

 // 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;

   // 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++)
      {
       for (unsigned i2=0;i2<nnode_1d;i2++)
        {
         // Local node number
         unsigned j=i0+i1*nnode_1d+i2*nnode_1d*nnode_1d;
         
         // Get pointer to set 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)
          {
           for (std::set<unsigned>::iterator it=boundaries_pt->begin();
                it!=boundaries_pt->end();++it)
            {

             // What's the boundary?
             unsigned boundary_id=*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 (brick!) element must 
             // have exactly four corner nodes on the boundary.
             if (boundary_identifier(boundary_id,fe_pt)==0)
              {
               Vector<int>* tmp_pt=new Vector<int>;
               tmp_vect_pt.push_back(tmp_pt); 
               boundary_identifier(boundary_id,fe_pt)=tmp_pt;
              }
            
             // Are we at a corner node?
             if (((i0==0)||(i0==nnode_1d-1))&&((i1==0)||(i1==nnode_1d-1))
                 &&((i2==0)||(i2==nnode_1d-1)))
              {
               // Create index to represent position relative to s_0
               (*boundary_identifier(boundary_id,fe_pt)).
                push_back(1*(2*i0/(nnode_1d-1)-1));               
               
               // Create index to represent position relative to s_1
               (*boundary_identifier(boundary_id,fe_pt)).
                push_back(2*(2*i1/(nnode_1d-1)-1));
 
               // Create index to represent position relative to s_2
               (*boundary_identifier(boundary_id,fe_pt)).
                push_back(3*(2*i2/(nnode_1d-1)-1));
              }
            }
          }
        }
      }     
    }  
  }
   

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

 // Loop over boundaries
 //---------------------
 for (unsigned i=0;i<nbound;i++)
  {
   // Loop over elements on given boundary
   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 (-3,-2,-1,1,2,3)
     std::map<int,unsigned> count;
     
     // Loop over coordinates
     for (int ii=0;ii<3;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] ]++;
      }
     
     // Determine the correct boundary indicator by checking that it 
     // occurs four times (since four 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 (int ii=0;ii<3;ii++)
      {
       // Loop over upper/lower end of coordinates
       for (int sign=-1;sign<3;sign+=2)
        {
         if (count[(ii+1)*sign]==4)
          {
           // Check that we haven't found multiple boundaries
           if (found)
            {
             throw OomphLibError(
              "Trouble: Multiple boundary identifiers!\n",
              OOMPH_CURRENT_FUNCTION,
              OOMPH_EXCEPTION_LOCATION);
            }
           found=true;
           indicator=(ii+1)*sign;
          }
        }
      }
     
     // Check if we've found one boundary
     if (found)
      {
       
       // Store element
       Boundary_element_pt[i].push_back(*it);
       
       // Now convert boundary indicator into information required
       // for FaceElements
       switch (indicator)
        {
        case -3:
         
         // s_2 is fixed at -1.0:
         Face_index_at_boundary[i].push_back(-3);
         break;
         
        case -2:
         
         // s_1 is fixed at -1.0:
         Face_index_at_boundary[i].push_back(-2);
         break;
         
        case -1:
         
         // s_0 is fixed at -1.0:
         Face_index_at_boundary[i].push_back(-1);
         break;
         
         
        case 1:
         
         // s_0 is fixed at 1.0:
         Face_index_at_boundary[i].push_back(1);
         break;
         
        case 2:
         
         // s_1 is fixed at 1.0:
         Face_index_at_boundary[i].push_back(2);
         break;
         
        case 3:
         
         // s_2 is fixed at 1.0:
         Face_index_at_boundary[i].push_back(3);
         break;
         
         
        default:
         
         throw OomphLibError("Never get here",
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
        }
       
      }
    }
  }
 
 // 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 along boundary is "
               <<  Face_index_at_boundary[i][e] << std::endl;
      }
    }
  }
 

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


 // Cleanup temporary vectors
 unsigned n=tmp_vect_pt.size();
 for (unsigned i=0;i<n;i++)
  {
   delete tmp_vect_pt[i];
   }*/

}


}