text: // Copyright (c) 2019    Patrick Diehl
  //
  // Distributed under the GNU GENERAL PUBLIC LICENSE, Version 3.0.
  // (See accompanying file LICENSE.txt)
  
  #ifndef MATERIAL_PD_ELASTIC_STATEMATERIAL_H
  #define MATERIAL_PD_ELASTIC_STATEMATERIAL_H
  
  #include "baseMaterial.h"
  #include "util/point.h"
  #include "util/matrixBlaze.h"
  #include "geometry/neighbor.h"
  #include "fe/mesh.h"
  #include "geometry/volumeCorrection.h"
  
  #include <vector>
  
  #include <hpx/include/parallel_algorithm.hpp>
  
  // forward declaration
  namespace inp {
  struct MaterialDeck;
  struct OutputDeck;
  }
  
  // forward declaration
  namespace data {
  class DataManager;
  }
  
  // forward declaration of class
  namespace fe {
  class Mesh;
  } // namespace fe
  
  namespace geometry {
  class VolumeCorrection;
  } // namespace geometry
  
  namespace util {
  class StateBasedHelperFunctions;
  } // namespace material
  
  namespace material {
  
  namespace pd {
  
  /*! @brief A class implementing the state-based elastic material model
   *
   * This class implements the the state-based elastic material model
   * as described by [Silling](https://doi.org/10.1007/s10659-007-9125-1)
   * and implementation details are described by [Littlewood](https://www.sandia.gov/~djlittl/docs/PeridynamicSoftwareRoadmap.pdf).
   *
   * @see https://doi.org/10.1007/s10659-007-9125-1
   */
  class ElasticState : public BaseMaterial {
  
   public:
    /*!
     * @brief Constructor
     * @param deck Pointer to the input deck
     * @param DataManager Pointer to the data manager object
     */
    ElasticState(inp::MaterialDeck *deck, data::DataManager *dataManager);
  
    /*!
     * @brief Returns energy and force state between node i and node j
     * @param i Id of node i
     * @param j Local id in the neighbor list of node i
     * @return Value Pair of energy and force
     */
  
    std::pair<util::Point3, double> getBondEF(size_t i, size_t j);
  
    /*!
     * @brief Returns critical bond strain between node i and node j
     *
     * @param i Id of node i
     * @param j Id of node j
     * @return Value Critical strain
     */
    double getSc(size_t i, size_t j);
  
    /*!
     * @brief Computes the stress tensor for one node
     * @param i Id of node i
     * @return The stress tensor for node i
     */
    util::Matrix33 getStress(size_t i);
  
    /*!
     * @brief Computes the strain vector for one node
     * @param i ID of node i
     * @return The strain tensor for node I
     */
    util::Matrix33 getStrain(size_t i);
  
    /*!
     * @brief return the factor for two-dimensional problems
     */
    double getFactor2D();
  
    void update();
  
    /*!
     * @brief Get direction of bond force
     * @return vector Unit vector along the bond force
     */
    util::Point3
    getBondForceDirection(const util::Point3 &dx,
                          const util::Point3 &du) const override {
      return dx / dx.length();
    }
  
    /*!
     * @brief Returns the value of influence function
     *
     * @param r Reference (initial) bond length
     * @return value Influence function at r
     */
    double getInfFn(const double &r) const {return 1.; };
  
   private:
    /*!
     * @brief Computes elastic state-based material parameters from elastic constants
     *
     * Either Young's modulus E or bulk modulus K, and Poisson ratio \f$ \nu \f$,
     * are needed.
     *
     * @param deck Input material deck
     * @param M Moment of influence function
     */
    void computeParameters(inp::MaterialDeck *deck, size_t dim);
  
    /**
     * @name Helper functions to compute the stress and strain
     */
    /**@{*/
  
    /*!
     * @brief Computes the deformation gradient for node i
     * @param i Id of node
     * @return The deformation gradient for node i
     *
     */
    util::Matrix33 deformation_gradient(size_t i);
  
    /*!
     * @brief Computes the x vector state between node i and node j
     * @param i Id of node
     * @param j Id of node
     * @return The x vector state
     *
     */
    util::Point3 X_vector_state(size_t i, size_t j);
  
    /*!
     * @brief Computes the y vector state between node i and node j
     * @param i Id of node
     * @param j Id of node
     * @return The y vector state
     *
     */
    util::Point3 Y_vector_state(size_t i, size_t j);
  
    /*!
     * @brief Computes the K shape tensor for node i
     * @param i Id of node
     * @return The  K shape tensor for node i
     *
     */
    util::Matrix33 K_shape_tensor(size_t i);
  
    /*!
     * @brief Computes the K modulus tensor
     * @param i Node i
     * @param j Neighbor j
     * @param k Neighbor k
     * @param m Index for the volume correction
     * @return The K modulus tensor
     *
     */
    util::Matrix33 K_modulus_tensor(size_t i, size_t j, size_t k, size_t m);
  
    /*!
     * @brief Provide the image of x under the Dirac Delta Function
     * @param deck The input deck
     * @param problem The corresponding problem
     * @param x Node x
     * @param i Node i
     * @param j Neighbor j
     * @param k Neighbor k
     * @param m Index for the volume correction
     * @return 1 if x is a null-vector, otherwise 0
     */
    double dirac_delta(util::Point3 x, size_t i, size_t j, size_t m);
  
    /*@}*/
  
  
  
    /*! @brief Correction factor for using plain stress */
    double d_factor2D;
  
    /*! @brief Dimension of the problem */
    size_t dim;
  
    double horizon;
  
    /*! @brief Compute strain energy */
    bool strainEnergy;
  
  
    /**
     * @name Pointers for the function parameters
     */
    /**@{*/
  
    /*! @brief Pointer to the material deck */
    const inp::MaterialDeck *d_deck;
  
    /*! @brief Pointer to the nodes */
    //const std::vector<util::Point3> *d_nodes;
  
    /*! @brief Pointer to the weighted volumes of the nodes */
    //const std::vector<double> *d_weightedVolume;
  
    /*! @brief Pointer to the volume correction of the nodes */
    //const std::vector<std::vector<double>> *d_volumeCorrection;
  
    /*! @brief Pointer to the volumes of the nodes */
    //const std::vector<double> * d_volumes;
  
    /*! @brief Pointer to the neighbors of the nodes */
    //geometry::Neighbor *d_neighbors;
  
    data::DataManager *d_dataManager_p;
  
    /*@}*/
  
  
  };
  
  } // namespace pd
  
  } // namespace material
  
  #endif // MATERIAL_PD_ELASTIC_STATEMATERIAL_H
  