one_d_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_ONE_D_MESH_TEMPLATE_CC
#define OOMPH_ONE_D_MESH_TEMPLATE_CC



// oomph-lib headers
#include "one_d_mesh.template.h"



namespace oomph
{


/// The generic mesh construction routine --- this contains all the hard
/// work and is called by all constructors
 template <class ELEMENT>
 void OneDMesh<ELEMENT>::build_mesh(TimeStepper* time_stepper_pt) 
 {
  // Set the length of the domain
  Length = Xmax - Xmin;
  
  // Set the number of boundaries -- there are 2 boundaries in a 1D mesh
  set_nboundary(2);
     
  // Allocate storage for the pointers to the elements
  Element_pt.resize(N);
     
  // Allocate memory for the first element
  Element_pt[0] = new ELEMENT;
     
  // Read out the number of nodes in the element (the member function
  // nnode_1d() is implemented in QElement)
  const unsigned n_node
   = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
     
  // We can now allocate storage for the pointers to the nodes in the mesh
  Node_pt.resize(1+(n_node-1)*N);
     
  // Initialise the node counter
  unsigned node_count = 0;

  // Initialise the minimum x coordinate in the mesh
  const double xinit = Xmin;
     
  // Calculate the length of the element
  const double el_length = Length/double(N);
     
  // Allocate storage for the local coordinate in the element
  Vector<double> s_fraction;
     
  // If the number of elements is 1, the first element is also the 
  // last element
  if(N==1)
   {
    // Set the first node
    // ------------------

    // Allocate memory for the node, using the element's own construct_node
    // function -- only the element knows what type of nodes it needs!
    Node_pt[node_count]
     = finite_element_pt(0)->construct_boundary_node(0,time_stepper_pt);
       
    // Set the position of the node
    node_pt(node_count)->x(0) = xinit;
       
    // Add the node to the boundary 0
    add_boundary_node(0,Node_pt[node_count]);
       
    // Increment the counter for the nodes
    node_count++;
       
    // Now build central nodes (ignore last one which needs special
    // ------------------------------------------------------------
    // treatment because it's on the boundary)
    // ---------------------------------------
    for(unsigned jnod=1;jnod<(n_node-1);jnod++)
     {
      // Allocate memory for nodes, as before
      Node_pt[node_count]
       = finite_element_pt(0)->construct_node(jnod,time_stepper_pt);
         
      // Get the local coordinate of the node
      finite_element_pt(0)->local_fraction_of_node(jnod,s_fraction);
         
      //Set the position of the node (linear mapping)
      node_pt(node_count)->x(0) = xinit + el_length*s_fraction[0];
         
      //Increment the node counter
      node_count++;
     }
       
    // New final node
    // --------------
       
    // Allocate memory for the node, as before
    Node_pt[node_count] = finite_element_pt(0)
     ->construct_boundary_node(n_node-1,time_stepper_pt);
       
    // Set the position of the node
    node_pt(node_count)->x(0) = xinit + Length;
       
    // Add the node to the boundary 1
    add_boundary_node(1,Node_pt[node_count]);
       
    // Increment the node counter
    node_count++;
   }

  // Otherwise, i.e. if there is more than one element, build all elements
  else
   {
    // -------------
    // FIRST ELEMENT 
    // -------------
       
    // Set the first node
    // ------------------
       
    // Allocate memory for the node, using the element's own construct_node
    // function -- only the element knows what type of nodes it needs!
    Node_pt[node_count]
     = finite_element_pt(0)->construct_boundary_node(0,time_stepper_pt);
       
    // Set the position of the node
    node_pt(node_count)->x(0) = xinit;
       
    // Add the node to the boundary 0
    add_boundary_node(0,Node_pt[node_count]);
       
    //Increment the counter for the nodes
    node_count++;

    // Now build the other nodes in the first element
    // ----------------------------------------------
       
    // Loop over the other nodes in the first element
    for(unsigned jnod=1;jnod<n_node;jnod++)
     {
      // Allocate memory for the nodes
      Node_pt[node_count]
       = finite_element_pt(0)->construct_node(jnod,time_stepper_pt);
         
      // Get the local coordinate of the node
      finite_element_pt(0)->local_fraction_of_node(jnod,s_fraction);
         
      // Set the position of the node (linear mapping)
      node_pt(node_count)->x(0) = xinit + el_length*s_fraction[0];
         
      // Increment the node counter
      node_count++;
     }
       
    // ----------------
    // CENTRAL ELEMENTS
    // ----------------

    // Loop over central elements in mesh
    for(unsigned e=1;e<(N-1);e++)
     {
      // Allocate memory for the new element
      Element_pt[e] = new ELEMENT;
         
      // The first node is the same as the last node in the neighbouring
      // element on the left
      finite_element_pt(e)->node_pt(0)
       = finite_element_pt(e-1)->node_pt((n_node-1));
         
      // Loop over the remaining nodes in the element
      for(unsigned jnod=1;jnod<n_node;jnod++)
       {
        // Allocate memory for the nodes, as before
        Node_pt[node_count]
         = finite_element_pt(e)->construct_node(jnod,time_stepper_pt);
           
        // Get the local coordinate of the nodes
        finite_element_pt(e)->local_fraction_of_node(jnod,s_fraction);
           
        // Set the position of the node (linear mapping)
        node_pt(node_count)->x(0) = xinit + el_length*(e + s_fraction[0]);
           
        // Increment the node counter
        node_count++;
       }
     } // End of loop over central elements
       
       
       // FINAL ELEMENT
       //--------------
       
       // Allocate memory for element
    Element_pt[N-1] = new ELEMENT;
       
    // The first node is the same as the last node in the neighbouring
    // element on the left
    finite_element_pt(N-1)->node_pt(0)
     = finite_element_pt(N-2)->node_pt(n_node-1);
       
    // New central nodes (ignore last one which needs special treatment
    // because it's on the boundary)
    for(unsigned jnod=1;jnod<(n_node-1);jnod++)
     {
      // Allocate memory for nodes, as before
      Node_pt[node_count]
       = finite_element_pt(N-1)->construct_node(jnod,time_stepper_pt);
         
      // Get the local coordinate of the node
      finite_element_pt(N-1)->local_fraction_of_node(jnod,s_fraction);
         
      // Set the position of the node
      node_pt(node_count)->x(0) = xinit + el_length*(N-1 + s_fraction[0]);
         
      // Increment the node counter
      node_count++;
     }
       
    // New final node
    // --------------
       
    // Allocate memory for the node, as before
    Node_pt[node_count] = finite_element_pt(N-1)
     ->construct_boundary_node(n_node-1,time_stepper_pt);
       
    // Set the position of the node
    node_pt(node_count)->x(0) = xinit + Length;
       
    // Add the node to the boundary 1
    add_boundary_node(1,Node_pt[node_count]);
       
    // Increment the node counter
    node_count++;
   }
 }

} // End of namespace

#endif