element_with_moving_nodes.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//====================================================================
//Functions for the ElementWithMovingNode class
#include "element_with_moving_nodes.h"
#include "geom_objects.h"
#include "algebraic_elements.h"

namespace oomph
{

 /////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////
 //Functions for the ElementWithMovingNodes class
 ///////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////

 //====================================================================
 /// Return a set of all geometric data associated with the element's node
 /// update function
 //======================================================================
 void ElementWithMovingNodes::assemble_set_of_all_geometric_data(
  std::set<Data*> &unique_geom_data_pt)
 {
  //First clear the set (just in case)
  unique_geom_data_pt.clear();
  
  //Get number of nodes
  const unsigned n_node = this->nnode();
  
  // Loop over all nodes
  for(unsigned n=0;n<n_node;n++)
   {
    //Cache pointer to the Node
    Node* const nod_pt = this->node_pt(n);
    
    //Is the node hanging
    const bool node_is_hanging = nod_pt->is_hanging();
    
    // Default number of master nodes
    unsigned nmaster=1;
    
    // Default: Node isn't hanging so it's its own master node
    Node* master_node_pt=nod_pt;
    
    // Cache the hanging point
    HangInfo* hang_info_pt = 0;
    
    // Find the number of master nodes if the node is hanging
    if (node_is_hanging) 
     {
      hang_info_pt = nod_pt->hanging_pt();
      nmaster = hang_info_pt->nmaster();
     }
    
    // Loop over all master nodes
    for (unsigned imaster=0;imaster<nmaster;imaster++)
     {
      //Get the master node
      if (node_is_hanging)
       {master_node_pt = hang_info_pt->master_node_pt(imaster);}
      

      //Find the number of data
      const unsigned n_geom_data = master_node_pt->ngeom_data();
      //If there are geometric data add them to the set
      if(n_geom_data > 0)
       {
        //Get vector of geometric data involved in the geometric 
        //change of this node
        Data** node_geom_data_pt = master_node_pt->all_geom_data_pt();
        
        for(unsigned i=0;i<n_geom_data;i++)
         {
          unique_geom_data_pt.insert(node_geom_data_pt[i]);
         }
       }
          
      // Find the number of geometric objects
      unsigned n_geom_obj = master_node_pt->ngeom_object();

      //If there are geometric objects, add them to the set
      if(n_geom_obj > 0)
       {
        // Get vector of geometric objects involved in the default 
        // update function for this (master) node.
        // Vector is constructed by copy operation.
        GeomObject** geom_object_pt = master_node_pt->all_geom_object_pt();
        
        //Loop over the geometric objects
        for(unsigned i=0;i<n_geom_obj;i++)
         {
          //Get the next geometric object
          GeomObject* geom_obj_pt=geom_object_pt[i];

          // Number of items of geometric data that the geometric
          // object depends on
          unsigned n_geom_data=geom_obj_pt->ngeom_data();
          
          // Loop over geometric data and add to set (use set to ensure 
          // that each one is only counted once)
          for (unsigned idata=0;idata<n_geom_data;idata++)
           {
            unique_geom_data_pt.insert(geom_obj_pt->geom_data_pt(idata));
           }         
         }
       } //End of add geom object loop
     }
   }
 }
 
 //=================================================================
 /// Construct the vector of (unique) geometric data
 //=================================================================
 void ElementWithMovingNodes::complete_setup_of_dependencies()
 {
  // This set will hold the pointers to all the unique (geometric) Data that 
  // affects the shape of this element
  std::set<Data*> unique_geom_data_pt;
  //Assemble that data
  this->assemble_set_of_all_geometric_data(unique_geom_data_pt);
  
  // Resize storage for the pointers to the Data items that are
  //involved in the element's node update operation.
  Geom_data_pt.clear();
  
  // Loop over all the unique remaining Data items involved in the
  // node update operations
  typedef std::set<Data*>::iterator IT;
  for (IT it=unique_geom_data_pt.begin();it!=unique_geom_data_pt.end();it++)
   {
    Geom_data_pt.push_back(*it);
   }
 }

 //==================================================================
 /// \short Function to describe the local dofs of the element[s]. The ostream 
 /// specifies the output stream to which the description 
 /// is written; the string stores the currently 
 /// assembled output that is ultimately written to the
 /// output stream by Data::describe_dofs(...); it is typically
 /// built up incrementally as we descend through the
 /// call hierarchy of this function when called from 
 /// Problem::describe_dofs(...)
 //==================================================================
 void ElementWithMovingNodes::
 describe_local_dofs(std::ostream& out,const std::string& current_string) const
 {
  // Call the standard finite element classification.
  FiniteElement::describe_local_dofs(out,current_string);
  
  //Set the number of data
  const unsigned n_geom_data = ngeom_data();

  //Loop over the node update data
  for(unsigned i=0;i<n_geom_data;i++)
   {
    // Pointer to geometric Data
    Data* data_pt=Geom_data_pt[i];
    
    std::stringstream conversion;
    conversion<<" of Geometric Data "<<i<<current_string;
    std::string in(conversion.str());
    data_pt->describe_dofs(out,in);
   }
 }


 //==================================================================
 /// Assign local equation numbers for the geometric data associated
 ///with the element.
 //==================================================================
 void ElementWithMovingNodes::assign_all_generic_local_eqn_numbers(
  const bool &store_local_dof_pt)
 {
  //Get local number of dofs so far
  unsigned local_eqn_number = this->ndof();
  
  //Set the number of data
  const unsigned n_geom_data = ngeom_data();
  
  // Reset number of geometric dofs
  Ngeom_dof=0;

  //If we have any geometric data
  if(n_geom_data > 0)
   {
    // Work out total number of values involved
    // Initialise from the first object
    unsigned n_total_values = Geom_data_pt[0]->nvalue();

    //Add the values from the other data
    for(unsigned i=1;i<n_geom_data;i++)
     {n_total_values += Geom_data_pt[i]->nvalue();}
    
    //If allocated delete the old storage
    if(Geometric_data_local_eqn)
     {
      delete[] Geometric_data_local_eqn[0];
      delete[] Geometric_data_local_eqn;
     }
    
    //If there are no values, we are done, null out the storage and
    //return
    if(n_total_values==0) {Geometric_data_local_eqn=0; return;}
    
    //Resize the storage for the geometric data local equation numbers
    //Firstly allocate the row pointers
    Geometric_data_local_eqn = new int*[n_geom_data];

    //Now allocate storage for all the equation numbers
    Geometric_data_local_eqn[0] = new int[n_total_values];

    //Initially all local equations are unclassified
    for(unsigned i=0;i<n_total_values;i++)
     {Geometric_data_local_eqn[0][i] = Data::Is_unclassified;}
    
    //Loop over the remaining rows and set their pointers
    for(unsigned i=1;i<n_geom_data;++i)
     {
      //Initially set the pointer to the i-th row to the pointer
      //to the i-1th row
      Geometric_data_local_eqn[i] = Geometric_data_local_eqn[i-1];

      //Now increase the row pointer by the number of values 
      //stored at the i-1th geometric data
      Geometric_data_local_eqn[i] += Geom_data_pt[i-1]->nvalue();
     }
        
    //A local queue to store the global equation numbers
    std::deque<unsigned long> global_eqn_number_queue;
    
    //Loop over the node update data
    for(unsigned i=0;i<n_geom_data;i++)
     {
      // Pointer to geometric Data
      Data* data_pt=Geom_data_pt[i];
      
      // Loop over values at this Data item
      unsigned n_value= data_pt->nvalue();
      for (unsigned j=0;j<n_value;j++)
       {
        // Get global equation number
        long eqn_number=data_pt->eqn_number(j);
        
        //If equation number positive 
        if (eqn_number>=0)
         {
          //Add the global equation number to our queue
          global_eqn_number_queue.push_back(eqn_number);
          //Add pointer to the dof to the queue if required
          if(store_local_dof_pt)
           {
            GeneralisedElement::Dof_pt_deque.push_back(
             data_pt->value_pt(j));
           }

          //Add to local value
          Geometric_data_local_eqn[i][j] = local_eqn_number;
          local_eqn_number++;

          // Bump up number of geometric dofs
          Ngeom_dof++;

         }
        else
         {
          //Set the local scheme to be pinned
          Geometric_data_local_eqn[i][j] = Data::Is_pinned;
         }
       }
     }
    
    //Now add our global equations numbers to the internal element storage
    this->add_global_eqn_numbers(global_eqn_number_queue,
                                 GeneralisedElement::Dof_pt_deque);
    //Clear the memory used in the deque
    if(store_local_dof_pt)
     {std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);}
   }
 }
 
 

 //==================================================================
 /// Calculate the node-update--related entries in the
 /// Jacobian. The vector passed
 /// in residuals has to contain the nonlinear residuals,
 /// evaluated for the current values of the unknowns, in
 /// case FDing is used to computed the derivatives.
 //==================================================================
 void ElementWithMovingNodes::fill_in_jacobian_from_geometric_data(
  Vector<double> &residuals, DenseMatrix<double> &jacobian)
 {
  if(!Bypass_fill_in_jacobian_from_geometric_data)
   {
    //Get number of Data items involved in node update operations
    const unsigned n_geometric_data = ngeom_data();
    
    //If there is nothing to be done, then leave
    if(n_geometric_data == 0) return;
    
    // Number of dofs
    const unsigned n_dof = this->ndof();
    
    // Number of nodes
    unsigned n_nod=this->nnode();
    
    // If there are no dofs, return
    if (n_nod==0) return;
    
    // Get nodal dimension from first node
    const unsigned dim_nod=node_pt(0)->ndim();
    
    // Number of shape controlling nodes for nonrefineable elements
    unsigned n_shape_controlling_node=nnode();
    
    // Are we dealing with a refineable element?
    RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(this);
    if (ref_el_pt!=0)
     {      
      // Adjust number of shape controlling nodes 
      n_shape_controlling_node=ref_el_pt->nshape_controlling_nodes();
     }
    
    // How are we going to evaluate the shape derivs?
    unsigned method=0;
    if (Method_for_shape_derivs==Shape_derivs_by_direct_fd)
     {
      method=0;
     }
    else if (Method_for_shape_derivs==Shape_derivs_by_chain_rule)
     {
      method=1;
     }
    else if (Method_for_shape_derivs==Shape_derivs_by_fastest_method)
     {
      // Direct FD-ing of residuals w.r.t. geometric dofs is likely to be 
      //faster if there are fewer geometric dofs than total nodal coordinates
      // (nodes x dim) in element:
      if (Ngeom_dof<(n_shape_controlling_node*dim_nod)) 
       {
        method=0;
       }
      else
       {
        method=1;
       }
      
     }
    
    // Choose method
    //===============
    switch(method)
     {
      
      // Direct FD:
      //-----------
     case 0:
      
     {
      //Create newres vector
      Vector<double> newres(n_dof);
      
      //Use the default finite difference step
      const double fd_step = GeneralisedElement::Default_fd_jacobian_step;
      
      //Integer storage for the local unknown
      int local_unknown=0;
      
      //Loop over the Data items that affect the node update operations
      for(unsigned i=0;i<n_geometric_data;i++)
       {
        //Loop over values
        unsigned n_value = Geom_data_pt[i]->nvalue();
        for(unsigned j=0;j<n_value;j++)
         {
          local_unknown = geometric_data_local_eqn(i,j);
          
          //If the value is free
          if(local_unknown >= 0)
           {
            //Get a pointer to the geometric data value
            double *value_pt = Geom_data_pt[i]->value_pt(j);
            
            //Save the old value
            double old_var = *value_pt;
            
            //Increment the variable
            *value_pt += fd_step;
            
            //Update the whole element (Bit inefficient)
            this->node_update();
            
            //Calculate the new residuals
            this->get_residuals(newres);
            
            //Now do finite differences
            for(unsigned m=0;m<n_dof;m++)
             {
              //Stick the entry into the Jacobian matrix
              jacobian(m,local_unknown) = (newres[m] - residuals[m])/fd_step;
             }
            
            //Reset the variable
            *value_pt = old_var;
            
            //We're relying on the total node update in the next loop
           }
         }
       }
      
      // Node update the element one final time to get things back to
      // the original state
      this->node_update();
     }
     
     break;
     
     // Chain rule
     //-----------
     case 1:
      
     {
      // Get derivatives of residuals w.r.t. all nodal coordinates
      RankThreeTensor<double> dresidual_dnodal_coordinates(
       n_dof,
       dim_nod,
       n_shape_controlling_node,
       0.0);
      
      // Use FD-version in base class?
      if (Evaluate_dresidual_dnodal_coordinates_by_fd)
       {
        if (ref_el_pt!=0)
         {
          ref_el_pt->RefineableElement::get_dresidual_dnodal_coordinates(
           dresidual_dnodal_coordinates);
         }
        else
         {
          FiniteElement::get_dresidual_dnodal_coordinates(
           dresidual_dnodal_coordinates);
         }
       }
      // Otherwise use the overloaded analytical version in derived
      // class (if it exists -- if it doesn't this just drops through
      // to the default implementation in FiniteElement).
      else
       {
        this->get_dresidual_dnodal_coordinates(dresidual_dnodal_coordinates);
       }
      
      // Get derivatives of nodal coordinates w.r.t. geometric dofs
      RankThreeTensor<double> dnodal_coordinates_dgeom_dofs(
       n_dof,
       dim_nod,
       n_shape_controlling_node,
       0.0);
      
      get_dnodal_coordinates_dgeom_dofs(dnodal_coordinates_dgeom_dofs);
      
      // Assemble Jacobian via chain rule
      for (unsigned l=0;l<n_dof;l++)
       {
        //Loop over the Data items that affect the node update operations
        for(unsigned i_data=0;i_data<n_geometric_data;i_data++)
         {
          //Loop over values
          unsigned n_value = Geom_data_pt[i_data]->nvalue();
          for(unsigned j_val=0;j_val<n_value;j_val++)
           {
            int k = geometric_data_local_eqn(i_data,j_val);
            
            //If the value is free
            if(k >= 0)
             {
              jacobian(l,k)=0.0;
              for (unsigned i=0;i<dim_nod;i++)
               {
                for (unsigned j=0;j<n_shape_controlling_node;j++)
                 {
                  jacobian(l,k)+=
                   dresidual_dnodal_coordinates(l,i,j)*
                   dnodal_coordinates_dgeom_dofs(k,i,j);
                 }
               }
             }
           }
         }
       }
     }
     
   break;
   
     default:
      
      std::ostringstream error_message;
      error_message << "Never get here: method " << method;
      throw OomphLibError(
       error_message.str(),
       OOMPH_CURRENT_FUNCTION,
       OOMPH_EXCEPTION_LOCATION);
      
     }
   }

 }
  
 
//======================================================================
/// \short Compute derivatives of the nodal coordinates w.r.t. 
/// to the geometric dofs. Default implementation by FD can be overwritten
/// for specific elements.
/// dnodal_coordinates_dgeom_dofs(l,i,j) = dX_{ij} / d s_l
//======================================================================
 void ElementWithMovingNodes::get_dnodal_coordinates_dgeom_dofs(
  RankThreeTensor<double>& dnodal_coordinates_dgeom_dofs)
 {
  //Get number of Data items involved in node update operations
  const unsigned n_geometric_data = ngeom_data();
  
  //If there is nothing to be done, then leave
  if(n_geometric_data == 0) {return;}
  
  // Number of nodes
  const unsigned n_nod=nnode();
  
  // If the element has no nodes (why??!!) return straightaway
  if (n_nod==0) return;
  
  // Get dimension from first node
  unsigned dim_nod=node_pt(0)->ndim();
  
  // Number of shape controlling nodes for nonrefineable elements
  unsigned n_shape_controlling_node=n_nod;
  
  // Are we dealing with a refineable element?
  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(this);
  if (ref_el_pt!=0)
   {      
    // Adjust number of shape controlling nodes 
    n_shape_controlling_node=ref_el_pt->nshape_controlling_nodes();
   }
  
  // Current and advanced nodal positions
  DenseMatrix<double> pos(dim_nod,n_shape_controlling_node);
  
  // Shape controlling nodes
  std::map<Node*,unsigned> local_shape_controlling_node_lookup;
  
  // Refineable element:
  if (ref_el_pt!=0)
   {
    local_shape_controlling_node_lookup=
     ref_el_pt->shape_controlling_node_lookup();
   }
  // Non-refineable element: the nodes themselves
  else
   {
    unsigned count=0;
    for (unsigned j=0;j<n_nod;j++)
     {
      local_shape_controlling_node_lookup[node_pt(j)]=count;
      count++;
     }
   }
  
  // Loop over all shape-controlling nodes to backup their original position
  for (std::map<Node*,unsigned>::iterator it=
        local_shape_controlling_node_lookup.begin();
       it!=local_shape_controlling_node_lookup.end();
       it++)
   {    
    // Get node
    Node* nod_pt=it->first;
    
    // Get its number
    unsigned node_number=it->second;
    
    // Backup
    for (unsigned i=0;i<dim_nod;i++)
     {
      pos(i,node_number)=nod_pt->position(i);
     }
   }
  
  
  //Integer storage for the local unknown
  int local_unknown=0;
  
  //Use the default finite difference step
  const double fd_step = GeneralisedElement::Default_fd_jacobian_step;
  
  //Loop over the Data items that affect the node update operations
  for(unsigned i=0;i<n_geometric_data;i++)
   {
    //Loop over values
    unsigned n_value = Geom_data_pt[i]->nvalue();
    for(unsigned j=0;j<n_value;j++)
     {
      local_unknown = geometric_data_local_eqn(i,j);
      
      //If the value is free
      if(local_unknown >= 0)
       {
        //Get a pointer to the geometric data value
        double *value_pt = Geom_data_pt[i]->value_pt(j);
        
        //Save the old value
        double old_var = *value_pt;
        
        //Increment the variable
        *value_pt += fd_step;
        
        //Update the whole element
        this->node_update();
        
        // Loop over all shape-controlling nodes
        for (std::map<Node*,unsigned>::iterator it=
              local_shape_controlling_node_lookup.begin();
             it!=local_shape_controlling_node_lookup.end();
             it++)
         {    
          // Get node
          Node* nod_pt=it->first;
          
          // Get its number
          unsigned node_number=it->second;
          
          // Get advanced position and FD
          for (unsigned ii=0;ii<dim_nod;ii++)
           {
            dnodal_coordinates_dgeom_dofs(local_unknown,ii,node_number)=
             (nod_pt->position(ii)-pos(ii,node_number))/fd_step;
           }
         }          
        
        //Reset the variable
        *value_pt = old_var;
        
        //We're relying on the total node update in the next loop       
       }
     }
   }
  // Node update the element one final time to get things back to
  // the original state
  this->node_update();
  
 }
 
}