constraint_elements.h

Go to the documentation of this file.

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
// Elements for adding constraints to a problem
// Base constraint element has:
//  Lagrange multipliers to enforce constraints - internal data
//  Indices of Lagrange multiplier dofs in internal data
//  External (constrained) degrees of freedom
//  Residual/Jacobian fill in
//  Constraint functions
 
// This represents the most general constraint element type
// - it can carry out procedures for automatically:
//  identifying constrained degrees of freedom
//  removing unenforcible constraints due to a lack of (unpinned) constrained degrees of freedom
//  fill in residual and jacobian entries of constrained degrees of freedom and lagrange multipliers for general constraint functions
// Provided the constraint functions are suitably formulated so that the rows/columns of the lagrange multiplier are not entirely zero
//  this will allow for general constraints to be applied to a problem, albeit in a rather inefficient way
 
// More specific constraint elements are implemented to allow for less flexible constraints but more explicit, readable, and efficient implementations
 
 
// Apply constraints to specific setups - i.e. node pairs, elements and local coordinate pairs
 
//  Node-Node constraint element
//   - a slightly more specific implementation with general constraint functions provided by function pointers but some additional member variables which
//     store the nodes being constrained and the local indices of the degrees of freedom of those nodes to assist in writing constraint functions
 
//  SolidNode-SolidNode constraint element - a more specific version of Node-Node mortaring element which only constrains solid degrees of freedom
 
//  Element-Element constraint element
//   - a slightly more specific implementation with general constraint functions provided by function pointers but some additional member variables which
//     store the elements being constrained, the local coordinates in those elements being constrained and the local indices of the degrees of freedom of the
//     elements/nodes of the elements to assist in writing constraint functions
 
 
// Apply equality constraint to specific setups - i.e. node pairs, elements and local coordinate pairs
 
// Mortaring elements
//  Node-Node mortaring element - a more specific implementation of the Node-Node constraint element which ensures all or some degrees of freedom of the two nodes are identical in value, including solid nodal
//   automatically determine if the requested constraints are enforcible
//   - if there are sufficient unpinned constrained degrees of freedom (must be at least one in each node for each constraint)
//   - if there are sufficient degrees of freedom for the specified degrees of freedom indices
//   - if the nodes are solid when mortaring of solid positions is specified
 
//  SolidNode-SolidNode mortaring element - a more specific version of Node-Node mortaring element which only ensures solid degrees of freedom are identical
 
 
 
 
//TODO need a way of identifying the mortaring space and building
// constraint elements. 
//TODO need face mortaring elements to apply to the mortaring boundary
 
 
// IMPORTANT advice when using constraint elements
// Since the feasibility of the application of constraints is dependent on the pinned status of the constrained degrees of freedom the creation of constraint elements should always happen after any steps in
//  which degrees of freedom in a problem are pinned.
 
#ifndef CONSTRAINT_ELEMENTS_HEADER
#define CONSTRAINT_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
//TODO why do I need to specify paths like this here but not elsewhere?
// oomph-lib headers
#include "../../../oomph-lib/src/generic/Vector.h"
#include "../../../oomph-lib/src/generic/nodes.h"
#include "../../../oomph-lib/src/generic/elements.h"
#include "../../../oomph-lib/src/generic/shape.h"
 
 
namespace oomph
{
 
 // Must call construct_lagrange_multipliers() in the constructor of any derived class
 class ConstraintElement : public virtual GeneralisedElement
 {
 public:
  // pass in all data which is to be included in the constraint calculations
  ConstraintElement(Vector<Data*>& data_pt) : Num_Constraints(0)
  {
   build(data_pt);
  }
  // Constructor but with a non-reference argument - used in node-node constraint elements
  // where a temporary vector of Data* is made
  ConstraintElement(Vector<Data*> data_pt) : Num_Constraints(0)
  {
   build(data_pt);
  }
 
  // Fill in residual and jacobian procedures
  // - Generic by finite differencing, override for more efficient implementation
  virtual void fill_in_contribution_to_residuals(Vector<double>& residuals);
 
  // Determine how many lagrange multipliers are needed and if the constraint is enforcible
  virtual void construct_lagrange_multipliers();
 
private:
  // Finish the element build
  void build(Vector<Data*>& data_pt)
  {
   // Add the data as external data
   Num_Constrained_Data = data_pt.size();
   Constrained_Data_Local_Index.resize(Num_Constrained_Data);
   for(unsigned i=0; i<Num_Constrained_Data; i++)
   {
    // Add the external data - we want to add it to the finite differencing
    Constrained_Data_Local_Index[i] = add_external_data(data_pt[i], true);
   }
  }
 
 protected:
  // Generic interface for evaluating constraint functions
  // - must re-size the error vector to the number of constraints being enforced
  virtual void evaluate_constraint_functions(Vector<double>& error);
 
  // Access functions for derived classes
  Vector<unsigned> get_constrained_data_local_index()
  {
   return Constrained_Data_Local_Index;
  }
  unsigned get_lagrange_multiplier_index()
  {
   return Lagrange_Multiplier_Index;
  }
  Data* lagrange_multiplier()
  {
   return internal_data_pt(Lagrange_Multiplier_Index);
  }
  const int lagrange_eqn(const unsigned& i)
  {
   // std::cout << "lagrange_eqn(" << i << ") ninternal_data " << ninternal_data() << std::endl;
  #ifdef RANGE_CHECKING
   if(i>=Num_Constraints)
   {
    throw OomphLibError(std::to_string(i) + " is greater than the number of constraints " + std::to_string(Num_Constraints),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  #endif
   return internal_local_eqn(Lagrange_Multiplier_Index, i);
  }
  Data* ext_data_pt(const unsigned& i)
  {
  #ifdef RANGE_CHECKING
   if(i>=Num_Constrained_Data)
   {
    throw OomphLibError(std::to_string(i) + " is greater than the number of constrained data " + std::to_string(Num_Constrained_Data),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  #endif
   return external_data_pt(Constrained_Data_Local_Index[i]);
  }
  const int ext_eqn(const unsigned& i, const unsigned& j)
  {
  #ifdef RANGE_CHECKING
   if(i>=Num_Constrained_Data)
   {
    throw OomphLibError(std::to_string(i) + " is greater than the number of constrained data " + std::to_string(Num_Constrained_Data),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  #endif
   return external_local_eqn(Constrained_Data_Local_Index[i], j);
  }
 
  // The index of the external data pt corresponding to the constrained data
  Vector<unsigned> Constrained_Data_Local_Index;
  // The number of constrained data
  unsigned Num_Constrained_Data;
  // The index of the lagrange multiplier internal data
  unsigned Lagrange_Multiplier_Index;
  // The number of constraints (and also the number of lagrange multipliers)
  unsigned Num_Constraints;
 };
} // End namespace 
 
#endif