node_element_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.
 
#ifndef NODE_ELEMENT_CONSTRAINT_ELEMENTS_HEADER
#define NODE_ELEMENT_CONSTRAINT_ELEMENTS_HEADER
 
#include "constraint_elements.h"
 
namespace oomph
{
 // It's difficult to write a general interface to interpolate values in an element.
 // So we only implement solid dof only mortaring element.
 // Mortars the position of a node to a local coordinate in an element.
 class NodeElementSolidOnlyMortaringElement : public virtual ConstraintElement
 {
 public:
  NodeElementSolidOnlyMortaringElement(SolidNode* node_pt, FiniteElement* elem_pt, Vector<double>& s)
  : ConstraintElement(NodeElementSolidOnlyMortaringElement::node_and_vector_to_data_pt(node_pt, elem_pt)),
    Local_Coord(s),
    Node_pt(node_pt),
    Element_pt(elem_pt),
    Dim(node_pt->ndim())
  {
   if(Dim != elem_pt->dim())
   {
    throw OomphLibError("Mortaring element and node do not have the same number of dimensions " + std::to_string(Dim) + " " + std::to_string(elem_pt->dim()),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
   }
   // Build the lagrange multipliers
   Num_Constraints = Dim;
   Lagrange_Multiplier_Index = add_internal_data(new Data(Dim), false);
  }
  // We know exactly how to fill in the residual and jacobian
  void fill_in_contribution_to_residuals(Vector<double>& residuals) override
  {
   fill_in_contribution_to_jacobian_mortared_nodes(residuals, GeneralisedElement::Dummy_matrix, 0);
  }
  void fill_in_contribution_to_jacobian(Vector<double>& residuals, DenseMatrix<double>& jacobian) override
  {
   fill_in_contribution_to_jacobian_mortared_nodes(residuals, jacobian, 1);
  }
 
  void position_in_element(Vector<double>& x)
  {
   // The number of nodes in the element
   const unsigned n_node = Element_pt->nnode();
   
   // The shape function
   Shape psi(n_node);
   Element_pt->shape(Local_Coord, psi);
 
   // Get the interpolated element coordinate
   for(unsigned i=0; i<Dim; i++)
   {
    x[i] = 0.0;
    for(unsigned l=0; l<n_node; l++)
    {
     // Get the coordinate of the node
     x[i] += psi[l] * ext_data_pt(1+l)->value(i);
    }
   }
  }
 
  // Used for when pruning redundant mortaring elements from the mesh
  SolidNode* solid_node_pt()
  {
   return Node_pt;
  }
 
 private:
  // Fill in rsidual and jacobian generic
  void fill_in_contribution_to_jacobian_mortared_nodes(Vector<double>& residuals, DenseMatrix<double>& jacobian, const bool& flag)
  {
   // The number of nodes in the element
   const unsigned n_node = Element_pt->nnode();
   
   // The shape function
   Shape psi(n_node);
   Element_pt->shape(Local_Coord, psi);
 
   // Loop over the dimensions
   for(unsigned i=0; i<Dim; i++)
   {
    // eqn numbers
    const int lambda_eqn = lagrange_eqn(i);
    const int node_eqn = ext_eqn(0,i);
    Vector<int> element_node_eqn(n_node);
 
    // Get the lagrange multiplier value
    const double lambda = lagrange_multiplier()->value(i);
 
    // The node coordinate
    const double x_node = ext_data_pt(0)->value(i);
 
    // Get the interpolated element coordinate
    double x_elem = 0.0;
    for(unsigned l=0;l<n_node;l++)
    {
     // Get the eqn number
     element_node_eqn[l] = ext_eqn(1+l, i);
     // Get the coordinate of the node
     x_elem += psi[l] * ext_data_pt(1+l)->value(i);
    }
 
    // Fill in the residuals and jacobian entries
 
    // Lambda residual
    if(lambda_eqn>=0)
    {
     residuals[lambda_eqn] += (x_node - x_elem);
     if(flag)
     {
      if(node_eqn>=0)
      {
       jacobian(lambda_eqn, node_eqn) += 1.0;
      }
      for(unsigned l=0;l<n_node;l++)
      {
       if(element_node_eqn[l]>=0)
       {
        jacobian(lambda_eqn, element_node_eqn[l]) -= psi[l];
       }
      }
     }
    }
 
    // Node residual
    if(node_eqn>=0)
    {
     residuals[node_eqn] += lambda;
     if(flag)
     {
      if(lambda_eqn>=0)
      {
       jacobian(node_eqn, lambda_eqn) += 1.0;
      }
     }
    }
 
    // Element node residual
    for(unsigned l=0;l<n_node;l++)
    {
     if(element_node_eqn[l]>=0)
     {
      residuals[element_node_eqn[l]] -= lambda*psi[l];
      if(flag)
      {
       if(lambda_eqn>=0)
       {
        jacobian(element_node_eqn[l], lambda_eqn) -= psi[l];
       }
      }
     }
    }
   }
  }
 
 
  // Get the vector of data from the node and the finite element
  static Vector<Data*> node_and_vector_to_data_pt(SolidNode* node_pt, FiniteElement* elem_pt)
  {
   const unsigned n_node_element = elem_pt->nnode();
   Vector<Data*> data_pt(1 + n_node_element);
   data_pt[0] = node_pt->variable_position_pt();
   for(unsigned i=0; i<n_node_element; i++)
   {
    data_pt[i+1] = dynamic_cast<SolidNode*>(elem_pt->node_pt(i))->variable_position_pt();
   }
   return data_pt;
  }
 
 private:
  Vector<double> Local_Coord;
  SolidNode* Node_pt;
  FiniteElement* Element_pt;
  unsigned Dim;
 };
} // End namespace
#endif