eighth_sphere_mesh.template.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//====================================================================
#ifndef OOMPH_EIGHTH_SPHERE_MESH_TEMPLATE_CC
#define OOMPH_EIGHTH_SPHERE_MESH_TEMPLATE_CC

#include "eighth_sphere_mesh.template.h"


namespace oomph
{

//======================================================================
/// Constructor for the eighth of a sphere mesh: Pass timestepper;
/// defaults to static default timestepper.
//======================================================================
template<class ELEMENT>
EighthSphereMesh<ELEMENT>::EighthSphereMesh(const double& radius,
                                            TimeStepper* time_stepper_pt) :
 Radius(radius)
{ 
 
 // Mesh can only be built with 3D Qelements.
 MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(3);
 
 //Set the number of boundaries
 this->set_nboundary(4);

 // Provide storage for the four elements
 this->Element_pt.resize(4);

 //Set the domain pointer: Pass the radius of the sphere
 Domain_pt = new EighthSphereDomain(Radius);


 Vector<double> s(3), s_fraction(3);
 Vector<double> r(3);

 // Create first element
 //--------------------
 this->Element_pt[0]= new ELEMENT;

 // Give it nodes
 
 // How many 1D nodes are there?
 unsigned nnode1d = this->finite_element_pt(0)->nnode_1d();
 
 // Loop over nodes
 for (unsigned i = 0; i < nnode1d; i++)
  {
   for (unsigned j = 0; j < nnode1d; j++)
    {
     for (unsigned k = 0; k < nnode1d; k++)
      {
       unsigned jnod = k * nnode1d * nnode1d + j * nnode1d + i;
       
       //If we're on a boundary, make a boundary node
       if((i==0) || (j==0) || (k==0))
        {
         // Allocate the node according to the element's wishes
         this->Node_pt.push_back(
          this->finite_element_pt(0)->
          construct_boundary_node(jnod,time_stepper_pt));         
        }
       //Otherwise it's a normal node
       else
        {
         // Allocate the node according to the element's wishes
         this->Node_pt.push_back(this->finite_element_pt(0)->
                                 construct_node(jnod,time_stepper_pt));  
        }

       // Work out the node's coordinates in the finite element's local
       // coordinate system
       this->finite_element_pt(0)->local_fraction_of_node(jnod,s_fraction);
         

       //Get the position of the node from macro element mapping
       s[0]=-1.0+2.0*s_fraction[0];  
       s[1]=-1.0+2.0*s_fraction[1];  
       s[2]=-1.0+2.0*s_fraction[2];
       Domain_pt->macro_element_pt(0)->macro_map(s,r);

       // Assign coordinates
       this->finite_element_pt(0)->node_pt(jnod)->x(0) = r[0];
       this->finite_element_pt(0)->node_pt(jnod)->x(1) = r[1];
       this->finite_element_pt(0)->node_pt(jnod)->x(2) = r[2];
       
       // Add the node to the appropriate boundary
       if(i==0)
        add_boundary_node(0,this->finite_element_pt(0)->node_pt(jnod));
       if(j==0)
        add_boundary_node(1,this->finite_element_pt(0)->node_pt(jnod));
       if(k==0)
        add_boundary_node(2,this->finite_element_pt(0)->node_pt(jnod));
         
      }
    }
  }

 
 // Create a second element
 //------------------------ 
 this->Element_pt[1]= new ELEMENT;;

 // Give it nodes
 
 // Loop over nodes 
 for (unsigned i = 0; i < nnode1d; i++)
  {
   for (unsigned j = 0; j < nnode1d; j++)
    {
     for (unsigned k = 0; k < nnode1d; k++)
      {
       unsigned jnod = k*nnode1d*nnode1d+j*nnode1d+i;
       
       // If node has not been yet created, create it
       if(i>0)
        {
         //If the node is on a boundary, make a boundary node
         if((i==nnode1d-1) || (j==0) || (k==0))
          {
           this->
            Node_pt.push_back(this->
                              finite_element_pt(1)->
                              construct_boundary_node(jnod,time_stepper_pt));
          }
         //Otherwise make a normal node
         else
          {
           this->Node_pt.push_back(this->finite_element_pt(1)->
                                   construct_node(jnod,time_stepper_pt));
          }

         // Work out the node's coordinates in the finite element's local
         // coordinate system
         this->finite_element_pt(1)->local_fraction_of_node(jnod,s_fraction);

         //Get the position of the node from macro element mapping
         s[0]=-1.0+2.0*s_fraction[0];
         s[1]=-1.0+2.0*s_fraction[1];
         s[2]=-1.0+2.0*s_fraction[2];
         Domain_pt->macro_element_pt(1)->macro_map(s,r);
          
         // Assign coordinate
         this->finite_element_pt(1)->node_pt(jnod)->x(0) = r[0];
         this->finite_element_pt(1)->node_pt(jnod)->x(1) = r[1];
         this->finite_element_pt(1)->node_pt(jnod)->x(2) = r[2];
         
         // Add the node to the approriate boundary
         if(j==0)
          add_boundary_node(1,this->finite_element_pt(1)->node_pt(jnod));
         if(k==0)
          add_boundary_node(2,this->finite_element_pt(1)->node_pt(jnod));
         if(i==nnode1d-1)
          add_boundary_node(3,this->finite_element_pt(1)->node_pt(jnod));
        }
              
       // ...else use the node already created
       else
        {
         this->finite_element_pt(1)->node_pt(jnod)=
          this->finite_element_pt(0)->node_pt(jnod+nnode1d-1);
        }
      }
    }
  }
 
   
 
 
 // Create a third element
 //------------------------
 this->Element_pt[2]=  new ELEMENT;

 // Give it nodes
 
 // Loop over nodes
 for (unsigned i = 0; i < nnode1d; i++)
  {
   for (unsigned j = 0; j < nnode1d; j++)
    {
     for (unsigned k = 0; k < nnode1d; k++)
      {
       unsigned jnod=k*nnode1d*nnode1d+j*nnode1d+i;
         
       // If the node has not been yet created, create it
       if((i<nnode1d-1) && (j>0))
        {
         //If it's on a boundary, make a boundary node
         if((i==0) || (j==nnode1d-1) || (k==0))
          {
           // Allocate the node according to the element's wishes
           this->
            Node_pt.push_back(this->
                              finite_element_pt(2)->
                              construct_boundary_node(jnod,time_stepper_pt));
          }
         //Otherwise allocate a normal node
         else
          {
           // Allocate the node according to the element's wishes
           this->Node_pt.push_back(this->finite_element_pt(2)->
                                   construct_node(jnod,time_stepper_pt));
          }

         // Work out the node's coordinates in the finite element's local
         // coordinate system
         this->finite_element_pt(2)->local_fraction_of_node(jnod,s_fraction);
         
         //Get the position of the node from macro element mapping
         s[0]=-1.0+2.0*s_fraction[0];
         s[1]=-1.0+2.0*s_fraction[1];
         s[2]=-1.0+2.0*s_fraction[2];
         Domain_pt->macro_element_pt(2)->macro_map(s,r);
          
         // Assign coordinates
         this->finite_element_pt(2)->node_pt(jnod)->x(0) = r[0];
         this->finite_element_pt(2)->node_pt(jnod)->x(1) = r[1];
         this->finite_element_pt(2)->node_pt(jnod)->x(2) = r[2];
         
         // Add the node to the appropriate boundary
         if(i==0)
          add_boundary_node(0,this->finite_element_pt(2)->node_pt(jnod));
         if(k==0)
          add_boundary_node(2,this->finite_element_pt(2)->node_pt(jnod));
         if(j==nnode1d-1)
          add_boundary_node(3,this->finite_element_pt(2)->node_pt(jnod));
        }
       
       // ...else use the nodes already created
       else
        {
         // If the node belongs to the element 0
         if(j==0)
          this->finite_element_pt(2)->node_pt(jnod)=
           this->finite_element_pt(0)->node_pt(jnod+nnode1d*(nnode1d-1));
         
           // ...else it belongs to the element 1
         else
          if(i==nnode1d-1)
           this->finite_element_pt(2)->node_pt(jnod)=
            this->finite_element_pt(1)->node_pt(nnode1d*nnode1d*k+j+i*nnode1d);
        }
      }
    }
  }
   
 
 //Create the fourth element
 //-------------------------
 this->Element_pt[3]= new ELEMENT;

 // Give it nodes

 // Loop over nodes
 for (unsigned i = 0; i < nnode1d; i++)
  {
   for (unsigned j = 0; j < nnode1d; j++)
    {
     for (unsigned k = 0; k < nnode1d; k++)
      {
       unsigned jnod=k*nnode1d*nnode1d+j*nnode1d+i;           

       // If the node has not been yet created, create it
       if((k>0) && (i<nnode1d-1) && (j<nnode1d-1))
        {
         //If it's on a boundary, allocate a boundary node
         if((i==0) || (j==0) || (k==nnode1d-1))
          {
           // Allocate the node according to the element's wishes
           this->
            Node_pt.push_back(this->
                              finite_element_pt(3)->
                              construct_boundary_node(jnod,time_stepper_pt));
          }
         //Otherwise allocate a normal node
         else
          {
           // Allocate the node according to the element's wishes
           this->Node_pt.push_back(this->finite_element_pt(3)->
                                   construct_node(jnod,time_stepper_pt));
          }

         // Work out the node's coordinates in the finite element's local
         // coordinate system
         this->finite_element_pt(3)->local_fraction_of_node(jnod,s_fraction);
         
         //Get the position of the node from macro element mapping
         s[0]=-1.0+2.0*s_fraction[0];
         s[1]=-1.0+2.0*s_fraction[1];
         s[2]=-1.0+2.0*s_fraction[2];
         Domain_pt->macro_element_pt(3)->macro_map(s,r);
          
         // Assign coordinates
         this->finite_element_pt(3)->node_pt(jnod)->x(0) = r[0];
         this->finite_element_pt(3)->node_pt(jnod)->x(1) = r[1];
         this->finite_element_pt(3)->node_pt(jnod)->x(2) = r[2];
         
         // Add the node to the appropriate boundary
         if(i==0)
          add_boundary_node(0,this->finite_element_pt(3)->node_pt(jnod));
         if(j==0)
          add_boundary_node(1,this->finite_element_pt(3)->node_pt(jnod));
         if(k==nnode1d-1)
          add_boundary_node(3,this->finite_element_pt(3)->node_pt(jnod));
        }

       // ...otherwise the node was already created: use it.
       else
        {
         // if k=0 then the node belongs to element 0
         if(k==0)
          {
           this->finite_element_pt(3)->node_pt(jnod)= 
            this->finite_element_pt(0)->
            node_pt(jnod+nnode1d*nnode1d*(nnode1d-1));
          }
         else
          // else if i==nnode1d-1 the node already exists in element 1
          if(i==nnode1d-1)
           {
            this->finite_element_pt(3)->node_pt(jnod)=
             this->finite_element_pt(1)->
             node_pt(k+i*nnode1d*nnode1d+j*nnode1d);
           }
          else
           // else, the node exists in element 2
           {
            this->finite_element_pt(3)->node_pt(jnod)=
             this->finite_element_pt(2)->
             node_pt(i+k*nnode1d+j*nnode1d*nnode1d);
           }
        }
      }
    }
  }

 // Setup boundary element lookup schemes
 setup_boundary_element_info();

}

}
#endif