face_mesh_project.h

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 guard to prevent multiple inclusions of the header
#ifndef OOMPH_FACE_MESH_PROJECT_HEADER
#define OOMPH_FACE_MESH_PROJECT_HEADER

namespace oomph
{


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


 //======================================================================
 /// \short Class that makes the finite element specified as template argument
 /// projectable -- on the assumption that all fields are interpolated
 /// by isoparametric Lagrange interpolation between the nodes.
 //======================================================================
 template <class ELEMENT>
 class GenericLagrangeInterpolatedProjectableElement :
  public virtual ProjectableElement<ELEMENT>
 {
  
 public:

  /// Constructor
  GenericLagrangeInterpolatedProjectableElement()
   {
    Boundary_id=UINT_MAX;
   }

 /// \short Nodal value of boundary coordinate
 double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i)
  const
  {
   //Vector in which to hold the intrinsic coordinate
   Vector<double> zeta(this->dim());
   
#ifdef PARANOID
   if (Boundary_id==UINT_MAX)
    {
     std::ostringstream error_message;
     error_message 
      << "Boundary_id is (still) UINT_MAX -- please set\n"
      << "the actual value with set_boundary_id(...)\n";
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   
   //Get the k-th generalised boundary coordinate at node n
   this->node_pt(n)->get_coordinates_on_boundary(
    Boundary_id,k,zeta);
   
   //Return the individual coordinate
   return zeta[i];
  }

  /// Boundary id 
  void set_boundary_id(const unsigned& boundary_id)
   {
    Boundary_id=boundary_id;
   }

  /// Boundary id 
  unsigned boundary_id() const
   {
#ifdef PARANOID
    if (Boundary_id==UINT_MAX)
     {
      std::ostringstream error_message;
      error_message 
       << "Boundary_id is (still) UINT_MAX -- please set\n"
       << "the actual value with set_boundary_id(...)\n";
      throw OomphLibError(error_message.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
#endif
    return Boundary_id;
   }

  /// \short Specify the values associated with field fld. 
  /// The information is returned in a vector of pairs which comprise 
  /// the Data object and the value within it, that correspond to field fld. 
  Vector<std::pair<Data*,unsigned> > data_values_of_field(const unsigned& fld)
   { 
    // Create the vector
    unsigned nnod=this->nnode();
    Vector<std::pair<Data*,unsigned> > data_values(nnod);
    
    // Loop over all nodes
    for (unsigned j=0;j<nnod;j++)
     {
      // Add the data value associated field: The node itself
      data_values[j]=std::make_pair(this->node_pt(j),fld);
     }
    
    // Return the vector
    return data_values;
   }
  
  /// \short Number of fields to be projected.
  unsigned nfields_for_projection()
   {
    return this->node_pt(0)->nvalue();
   }
  
  /// \short Number of history values to be stored for fld-th field
  /// (includes current value!). Extract from first node but assume it's
  /// the same for all.
  unsigned nhistory_values_for_projection(const unsigned &fld)
   {
    return this->node_pt(0)->ntstorage();   
   }
  
  ///\short Number of positional history values (Note: count includes 
  /// current value!). Extract from first node but assume it's
  /// the same for all.
  unsigned nhistory_values_for_coordinate_projection()
   {
    return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
   }
  
  
  /// \short Return Jacobian of mapping and shape functions of field fld
  /// at local coordinate s. 
  double jacobian_and_shape_of_field(const unsigned &fld, 
                                     const Vector<double> &s, 
                                     Shape &psi)
   {
    this->shape(s,psi);
    return this->J_eulerian(s);
   }
  
  
  
  /// \short Return interpolated field fld at local coordinate s, at time level
  /// t (t=0: present; t>0: history values)
  double get_field(const unsigned &t, 
                   const unsigned &fld,
                   const Vector<double>& s)
   {
    //Local shape function
    unsigned n_node=this->nnode();
    Shape psi(n_node);
    
    //Find values of shape function
    this->shape(s,psi);
    
    //Initialise value of u
    double interpolated_u = 0.0;
    
    //Sum over the local nodes
    for(unsigned l=0;l<n_node;l++) 
     {
      interpolated_u += this->nodal_value(t,l,fld)*psi[l];
     }
    return interpolated_u;     
   }
  
  /// \short Return number of values in field fld
  unsigned nvalue_of_field(const unsigned &fld)
   {
    return this->nnode();
   }
  
  
  /// \short Return local equation number of value j in field fld. Assumed to be
  /// the local nodal equation.
  int local_equation(const unsigned &fld,
                     const unsigned &j)
   {
    return this->nodal_local_eqn(j,fld);     
   }
  

 private:

  /// Boundary id
  unsigned Boundary_id;

 };
 
 

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

 
 //======================================================================
 /// Class that makes backup (via a deep copy) of a mesh, keeping alive
 /// enough information to allow the solution that is currently stored
 /// on the mesh to be projected onto another mesh sometime in the 
 /// future (when the original mesh may already have been deleted).
 /// This is mainly useful for the projection of additional 
 /// nodal values (such as Lagrange multipliers) created by FaceElements.
 /// ASSUMPTION: All fields in the element are represented by isoparametric
 ///             Lagrange interpolation between the nodal values.
 ///             Any fields that do not fall into this category will not
 ///             be copied across correctly and if you're unlucky 
 ///             the code may die...).
 //======================================================================
  template <class GEOMETRIC_ELEMENT> 
 class BackupMeshForProjection : public virtual Mesh
 {

 public:

  /// \short Constructor: Pass existing mesh and the boundary ID (need to find
  /// the boundary coordinates that are used for the projection.
  /// Optional final argument specifies the ID of the field (i.e. the
  /// index of the relevant nodal value!) to be projected.
  /// If omitted, we project all of them. 
   BackupMeshForProjection(Mesh* mesh_pt, const unsigned& boundary_id,
                           const unsigned& id_of_field_to_be_projected=UINT_MAX)
    : Boundary_id(boundary_id), 
   ID_of_field_to_be_projected(id_of_field_to_be_projected)
   {
    // Find unique nodes (via elements because Node_pt vector in
    // original mesh may not have been filled (typical for most
    // face element meshes)
    unsigned nel=mesh_pt->nelement();
    Element_pt.reserve(nel);
    Node_pt.reserve(mesh_pt->nnode());
    for (unsigned e=0;e<nel;e++)
     {
      FiniteElement* el_pt=mesh_pt->finite_element_pt(e);
      if (el_pt!=0)
       {
        // Make new element
#ifdef PARANOID
        if (dynamic_cast<GEOMETRIC_ELEMENT*>(mesh_pt->element_pt(e))==0)
         {
          std::ostringstream error_message;
          error_message 
           << "Element is of wrong type " << typeid(el_pt).name() 
           << " doesn't match template parameter!" << std::endl;
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);

          if (el_pt->ninternal_data()!=0)
           {
            std::ostringstream error_message;
            error_message 
             << "Internal data will NOT be projected across!\n"
             << "If you want this functionality you'll have to \n"
             << "implement it yourself" << std::endl;
            OomphLibWarning(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
           }
         }
#endif
        
        // Make a new element
        GenericLagrangeInterpolatedProjectableElement<GEOMETRIC_ELEMENT>* 
         new_el_pt=new GenericLagrangeInterpolatedProjectableElement<GEOMETRIC_ELEMENT>;

        // Set boundary ID
        new_el_pt->set_boundary_id(Boundary_id);

        // Set nodal dimension
        unsigned nodal_dim=el_pt->node_pt(0)->ndim();
        new_el_pt->set_nodal_dimension(nodal_dim);

        // Add it to mesh
        add_element_pt(new_el_pt);
        
        // Create new nodes if needed
        unsigned nnod=el_pt->nnode();
        for (unsigned j=0;j<nnod;j++)
         {
          Node* old_node_pt=el_pt->node_pt(j);
          if (New_node_pt[old_node_pt]==0)
           {
            Node* new_nod_pt=0;
            

#ifdef PARANOID
            // Check boundary node-ness
            if (!old_node_pt->is_on_boundary())
             {
              std::ostringstream error_message;
              error_message 
               << "Node isn't on a boundary!"
               << std::endl;
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
             }
#endif
            

            // How many values are we projecting? Default: All
            unsigned nval=old_node_pt->nvalue();
            unsigned first_index_in_old_node=0;
            if (ID_of_field_to_be_projected!=UINT_MAX)
             {
              nval=dynamic_cast<BoundaryNodeBase*>
               (old_node_pt)->nvalue_assigned_by_face_element
               (ID_of_field_to_be_projected);
              first_index_in_old_node=dynamic_cast<BoundaryNodeBase*>
               (old_node_pt)->index_of_first_value_assigned_by_face_element
               (ID_of_field_to_be_projected);
             }
            
            // Build new node
            new_nod_pt=new BoundaryNode<Node>(old_node_pt->time_stepper_pt(),
                                              old_node_pt->ndim(), 
                                              old_node_pt->nposition_type(),
                                              nval);
            
            // Copy data across
            unsigned n_time = old_node_pt->ntstorage();
            for(unsigned t=0;t<n_time;t++)
             {
              for(unsigned i=0;i<nval;i++) 
               {
                new_nod_pt->set_value
                 (t,i,old_node_pt->value(t,first_index_in_old_node+i));
               }
             }
            
            // Copy nodal positions
            unsigned n_dim=old_node_pt->ndim();
            for(unsigned i=0;i<n_dim;i++)
             {
              new_nod_pt->x(i)=old_node_pt->x(i);
             }

            // Add to boundary
            new_nod_pt->add_to_boundary(Boundary_id);
            
            // Transfer boundary coordinates
#ifdef PARANOID
            if (!old_node_pt->is_on_boundary(Boundary_id))
             {
              std::ostringstream error_message;
              error_message 
               << "Boundary ID specified as " << Boundary_id 
               << " but node isn't actually on that boundary!"
               << std::endl;
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
             }
#endif
            
            // Get vector of coordinates on mesh boundary from old node
            unsigned n=old_node_pt->ncoordinates_on_boundary(Boundary_id);
            Vector<double> zeta(n);
            old_node_pt->get_coordinates_on_boundary(Boundary_id,zeta);
            
            // Set for new node
            new_nod_pt->set_coordinates_on_boundary(Boundary_id,zeta);
            
            // Add node
            Node_pt.push_back(new_nod_pt);
          
            // Setup association
            New_node_pt[old_node_pt]=new_nod_pt;
           }
          
          // Set node pointer from new element
          new_el_pt->node_pt(j)=New_node_pt[old_node_pt];
         }
       }
     }
   }
  
  /// \short Project the solution that was present in the original mesh
  /// and from which this backup mesh was created onto the mesh
  /// pointed to by new_mesh_pt. Note that elements in the new mesh do 
  /// not have to be projectable. The original mesh may by now have
  /// been deleted.
  void project_onto_new_mesh(Mesh* new_mesh_pt)
  {
   // Make copy of new mesh that we can project onto
   BackupMeshForProjection<GEOMETRIC_ELEMENT>* projectable_new_mesh_pt=
    new BackupMeshForProjection<GEOMETRIC_ELEMENT>(new_mesh_pt,Boundary_id,
                                         ID_of_field_to_be_projected);
   
   // Create projection problem
   ProjectionProblem<GenericLagrangeInterpolatedProjectableElement
    <GEOMETRIC_ELEMENT> >* 
    proj_problem_pt=new 
    ProjectionProblem<GenericLagrangeInterpolatedProjectableElement
    <GEOMETRIC_ELEMENT> >;
   
   // Set the mesh we want to project onto
   proj_problem_pt->mesh_pt()=projectable_new_mesh_pt;
   
   // Project from projectable copy of original mesh -- this one!
   bool dont_project_positions=true;
   proj_problem_pt->project(this,dont_project_positions);
   
   // Copy nodal values onto the corresponding nodes in the new mesh
   projectable_new_mesh_pt->copy_onto_original_mesh();
   
   // Kill!
   delete proj_problem_pt;
   delete projectable_new_mesh_pt;
  }
  
  
  /// \short Copy nodal values back onto original mesh from which this 
  /// mesh was built. This obviously only makes sense if the original
  /// mesh still exists!
  void copy_onto_original_mesh()
   {

    for (std::map<Node*,Node*>::iterator it=New_node_pt.begin();
         it!=New_node_pt.end();it++)
     {      
      // Get old node (in the previously existing mesh)
      Node* old_node_pt=(*it).first;
      
      // ...and corresponding new one (in the mesh where we did the projection
      Node* new_node_pt=(*it).second; 
      
      // How many values are we moving across?
      unsigned nval=old_node_pt->nvalue();
      unsigned first_index_in_old_node=0;
      if (ID_of_field_to_be_projected!=UINT_MAX)
       {
        nval=dynamic_cast<BoundaryNodeBase*>(old_node_pt)->
         nvalue_assigned_by_face_element
         (ID_of_field_to_be_projected);
        first_index_in_old_node=dynamic_cast<BoundaryNodeBase*>(old_node_pt)->
         index_of_first_value_assigned_by_face_element
         (ID_of_field_to_be_projected);
       }
      
      // Copy data across (include offset in orig mesh)
      unsigned n_time = old_node_pt->ntstorage();
      for(unsigned t=0;t<n_time;t++)
       {
        for(unsigned i=0;i<nval;i++) 
         {
          old_node_pt->set_value
           (t,first_index_in_old_node+i,new_node_pt->value(t,i));
         }
       }

     }
   }


 private:

  /// Map returning new node, labeled by node point in original mesh
  std::map<Node*,Node*> New_node_pt;

  /// Boundary id 
  unsigned Boundary_id;

  /// \short ID of field to be projected (UINT_MAX if there isn't one, in which 
  /// case we're doing all of them.
  unsigned ID_of_field_to_be_projected;

 };


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

}

#endif