MWLS.hpp

Go to the documentation of this file.

/** \file MWLS.hpp
 */
 
/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
#ifndef __MWLS_HPP__
#define __MWLS_HPP__
 
namespace FractureMechanics {
 
/**
 
Moving least weighted square approximation (MWLS)
 
\f[
u^h(x) = p_k(x)a_k(x)\\
J(x)=\frac{1}{2} \sum_i w(x-x_i)\left(p_j(x_i)a_j(x)-u_i\right)^2 = \\
J(x)=\frac{1}{2} \sum_i w(x-x_i)\left( p_j(x_i)a_j(x) p_k(x_i)a_k(x) -
2p_j(x_i)a_j(x)u_i + u_i^2 \right) \\ \frac{\partial J}{\partial a_j} =
\sum_i w(x_i-x)\left(p_j(x_i)p_k(x_i) a_k(x) -p_j(x_j)u_i\right)\\
A_{jk}(x) = \sum_i w(x_i-x)p_j(x_i)p_k(x_i)\\
B_{ji}(x) = w(x-x_i) p_j(x_i)\\
A_{jk}(x) a_k(x) = B_{ji}(x)u_i\\
a_k(x) = \left( A^{-1}_{kj}(x) B_{ji}(x) \right) u_i\\
u^h(x) = p_k(x) \left( A^{-1}_{kj}(x) B_{ji}(x) \right) u_i \\
u^h(x) = \left[ p_k(x) \left( A^{-1}_{kj}(x) B_{ji}(x) \right) \right] u_i \\
\phi_i = p_k(x) \left( A^{-1}_{kj}(x) B_{ji}(x) \right) \\
\phi_{i,l} =
p_{k,l}(x) \left( A^{-1}_{kj}(x) B_{ji}(x) \right)+
p_k(x) \left( A^{-1}_{kj,l}(x) B_{ji}(x) + A^{-1}_{kj}(x) B_{ji,l}(x)  \right)\\
\phi_{i,l} =
p_{k,l}(x) \left( A^{-1}_{kj}(x) B_{ji}(x) \right)+
p_k(x) \left( A^{-1}_{kj,l}(x) B_{ji}(x) + A^{-1}_{kj}(x) B_{ji,l}(x)  \right)\\
\left(\mathbf{A}\mathbf{A}^{-1}\right)_{,l} = \mathbf{0} \to
\mathbf{A}^{-1}_{,l} = -\mathbf{A}^{-1} \mathbf{A}_{,l} \mathbf{A}^{-1}\\
A_{jk,l} = \sum_i w_{,l}(x_i-x)p_j(x_i)p_k(x_i) \\
B_{ij,l} = w_{,l}(x_i-x)p_j(x_i)
\f]
 
*/
struct MWLSApprox {
 
  MoFEM::Interface &mField;
  moab::Core mwlsCore;
  moab::Interface &mwlsMoab;
  Vec F_lambda;
 
  boost::weak_ptr<DofEntity> arcLengthDof;
 
  MWLSApprox(MoFEM::Interface &m_field, Vec F_lambda = PETSC_NULL,
             boost::shared_ptr<DofEntity> arc_length_dof = nullptr);
  virtual ~MWLSApprox() = default;
 
  /**
   * @brief Get the Use Local Base At Material Reference Configuration object
   *
   * If true base is calculated at initail integration point coordinates and
   * base function evaluated at another point. Local derivatives are uses, since
   * weight function is at initial point.
   *
   * @return true
   * @return false
   */
  bool getUseGlobalBaseAtMaterialReferenceConfiguration() const {
    return (useGlobalBaseAtMaterialReferenceConfiguration == PETSC_TRUE);
  }
 
  /**
   * \brief get options from command line for MWLS
   * @return error code
   */
  MoFEMErrorCode getMWLSOptions();
 
  /**
   * Load mesh with data and do basic calculation
   * @param  file_name File name
   * @return           Error code
   */
  MoFEMErrorCode loadMWLSMesh(std::string file_name);
 
  /**
   * \brief get values form elements to nodes
   * @param  th tag with approximated data
   * @param  th tag with approximated data
   * @return    error code
   */
  MoFEMErrorCode getValuesToNodes(Tag th);
 
  /**
   * \brief Get node in influence domain
   * @param  material_coords material position
   * @return error code
   */
  MoFEMErrorCode getInfluenceNodes(double material_coords[3]);
 
  /**
  * \brief Calculate approximation function coefficients
  * @param  material_coords material position
  * @return          error code
 
  NOTE: This function remove form influenceNodes entities which
  are beyond influence radius.
 
  */
  MoFEMErrorCode calculateApproxCoeff(bool trim_influence_nodes,
                                      bool global_derivatives);
 
  /**
  * \brief evaluate approximation function at material point
  * @param  material_coords material position
  * @return          error code
 
  */
  MoFEMErrorCode evaluateApproxFun(double eval_point[3]);
 
  /**
  * \brief evaluate 1st derivative of approximation function at material point
  * @param  global_derivatives decide whether global derivative will be
  calculated
  * @return          error code
 
  */
  MoFEMErrorCode evaluateFirstDiffApproxFun(double eval_point[3],
                                            bool global_derivatives);
 
  /**
  * \brief evaluate 2nd derivative of approximation function at material point
  * @param  global_derivatives decide whether global derivative will be
  calculated
  * @return          error code
 
  */
  MoFEMErrorCode evaluateSecondDiffApproxFun(double eval_point[3],
                                             bool global_derivatives);
 
  /**
   * Get tag data on influence nodes
   * @param  th tag
   * @return    error code
   */
  MoFEMErrorCode getTagData(Tag th);
 
  /**
   * @brief Get the norm of the error at nod
   *
   * @param th
   * @param total_stress_at_node
   * @param total_stress_error_at_node
   * @param hydro_static_error_at_node
   * @param deviatoric_error_at_node
   * @param total_energy
   * @param total_energy_error
   * @param young_modulus
   * @param poisson_ratio
   * @return MoFEMErrorCode
   */
  MoFEMErrorCode getErrorAtNode(
      Tag th, double &total_stress_at_node, double &total_stress_error_at_node,
      double &hydro_static_error_at_node, double &deviatoric_error_at_node,
      double &total_energy, double &total_energy_error,
      const double young_modulus, const double poisson_ratio);
 
  typedef ublas::vector<double, ublas::bounded_array<double, 9>> VecVals;
  typedef ublas::matrix<double, ublas::row_major,
                        ublas::bounded_array<double, 27>>
      MatDiffVals;
  typedef ublas::matrix<double, ublas::row_major,
                        ublas::bounded_array<double, 81>>
      MatDiffDiffVals;
 
  /**
   * Get values at point
   * @return error code
   */
  const VecVals &getDataApprox();
 
  /**
   * Get derivatives of values at points
   * @return error code
   */
  const MatDiffVals &getDiffDataApprox();
 
  /**
   * Get second derivatives of values at points
   * @return error code
   */
  const MatDiffDiffVals &getDiffDiffDataApprox();
 
  // Stress at gauss point on finite element
  MatrixDouble stressAtGaussPts;
  MatrixDouble diffStressAtGaussPts;
 
  // Density at gauss point on finite element
  boost::shared_ptr<VectorDouble> rhoAtGaussPts;
  boost::shared_ptr<MatrixDouble> diffRhoAtGaussPts;
  boost::shared_ptr<MatrixDouble> diffDiffRhoAtGaussPts;
 
  MatrixDouble singularInitialDisplacement;
  MatrixDouble singularCurrentDisplacement;
  MatrixDouble mwlsMaterialPositions;
  MatrixDouble approxBasePoint;
  VectorDouble3 approxPointCoords;
  Range tetsInBlock;
  Range trisInBlock;
  Tag thMidNode;
 
  /**
   * @brief Store Coefficient of base functions at integration points
   *
   */
  std::map<EntityHandle, std::vector<MatrixDouble>> invABMap;
 
  /**
   * @brief Influence nodes on elements on at integration points
   *
   */
  std::map<EntityHandle, std::vector<std::vector<EntityHandle>>>
      influenceNodesMap;
 
  /**
   * @brief Store weight function radius at the nodes
   *
   */
  std::map<EntityHandle, std::vector<double>> dmNodesMap;
 
  template <class T> struct OpMWLSBase : public T::UserDataOperator {
    boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
    boost::shared_ptr<MatrixDouble> matGradPosAtPtsPtr;
    boost::weak_ptr<CrackFrontElement> feSingularPtr;
    boost::shared_ptr<MWLSApprox> mwlsApprox;
 
    OpMWLSBase(boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr,
               boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr,
               boost::shared_ptr<CrackFrontElement> fe_singular_ptr,
               boost::shared_ptr<MWLSApprox> mwls_approx, bool testing = false)
        : T::UserDataOperator("MESH_NODE_POSITIONS",
                              T::UserDataOperator::OPROW),
          matPosAtPtsPtr(mat_pos_at_pts_ptr),
          matGradPosAtPtsPtr(mat_grad_pos_at_pts_ptr),
          feSingularPtr(fe_singular_ptr), mwlsApprox(mwls_approx) {}
 
    virtual ~OpMWLSBase() = default;
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data);
 
  protected:
    /**
     * @brief Do specific work in operator.
     *
     * Is protected to prevent direct call. This can be only called by
     * doWork method.
     *
     * @param side
     * @param type
     * @param data
     * @return MoFEMErrorCode
     */
    virtual MoFEMErrorCode
    doMWLSWork(int side, EntityType type,
               DataForcesAndSourcesCore::EntData &data) = 0;
  };
 
  template <class T>
  struct OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl : OpMWLSBase<T> {
 
    const bool testing;
    OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl(
        boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr,
        boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr,
        boost::shared_ptr<CrackFrontElement> fe_singular_ptr,
        boost::shared_ptr<MWLSApprox> mwls_approx, bool testing = false)
        : OpMWLSBase<T>(mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
                        fe_singular_ptr, mwls_approx, testing),
          testing(testing) {}
 
  protected:
    MoFEMErrorCode doMWLSWork(int side, EntityType type,
                              DataForcesAndSourcesCore::EntData &data);
  };
 
  typedef OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl<
      VolumeElementForcesAndSourcesCore>
      OpMWLSCalculateBaseCoeffcientsAtGaussPts;
 
  typedef OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl<
      FaceElementForcesAndSourcesCore>
      OpMWLSCalculateBaseCoeffcientsForFaceAtGaussPts;
 
  struct OpMWLSRhoAtGaussUsingPrecalulatedCoeffs
      : public OpMWLSBase<VolumeElementForcesAndSourcesCore> {
 
    std::string rhoTagName;
    const bool calculateDerivative;
    const bool calculate2ndDerivative;
    const bool testing;
    OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
        boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr,
        boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr,
        boost::shared_ptr<CrackFrontElement> fe_singular_ptr,
        boost::shared_ptr<MWLSApprox> mwls_approx,
        const std::string rho_tag_name, const bool calculate_derivative,
        const bool calculate_2nd_derivative, bool testing = false)
        : OpMWLSBase<VolumeElementForcesAndSourcesCore>(
              mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_singular_ptr,
              mwls_approx, testing),
          rhoTagName(rho_tag_name), calculateDerivative(calculate_derivative),
          calculate2ndDerivative(calculate_2nd_derivative), testing(testing) {}
 
  protected:
    MoFEMErrorCode doMWLSWork(int side, EntityType type,
                              DataForcesAndSourcesCore::EntData &data);
  };
 
  /**
   * \brief Evaluate density at integration points
   */
  struct OpMWLSRhoAtGaussPts : OpMWLSBase<VolumeElementForcesAndSourcesCore> {
    std::string rhoTagName;
 
    const bool calculateDerivative;
    const bool calculate2ndDerivative;
    const bool testing;
    OpMWLSRhoAtGaussPts(boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr,
                        boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr,
                        boost::shared_ptr<CrackFrontElement> fe_singular_ptr,
                        boost::shared_ptr<MWLSApprox> mwls_approx,
                        const std::string rho_tag_name,
                        const bool calculate_derivative,
                        const bool calculate_2nd_derivative,
                        bool testing = false)
        : OpMWLSBase<VolumeElementForcesAndSourcesCore>(
              mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_singular_ptr,
              mwls_approx, testing),
          rhoTagName(rho_tag_name), calculateDerivative(calculate_derivative),
          calculate2ndDerivative(calculate_2nd_derivative), testing(testing) {}
 
  protected:
    MoFEMErrorCode doMWLSWork(int side, EntityType type,
                              DataForcesAndSourcesCore::EntData &data);
  };
 
  struct OpMWLSRhoPostProcess
      : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
    boost::shared_ptr<moab::Interface> postProcMeshPtr;
    boost::shared_ptr<std::vector<EntityHandle>> mapGaussPtsPtr;
    boost::shared_ptr<MWLSApprox> mwlsApproxPtr;
    PostProcVolumeOnRefinedMesh::CommonData &commonData;
 
    OpMWLSRhoPostProcess(
        boost::shared_ptr<moab::Interface> &post_proc_mesh,
        boost::shared_ptr<std::vector<EntityHandle>> &map_gauss_pts,
        boost::shared_ptr<MWLSApprox> mwls_approx,
        PostProcVolumeOnRefinedMesh::CommonData &common_data)
        : VolumeElementForcesAndSourcesCore::UserDataOperator(
              "SPATIAL_POSITION", UserDataOperator::OPROW),
          postProcMeshPtr(post_proc_mesh), mapGaussPtsPtr(map_gauss_pts),
          mwlsApproxPtr(mwls_approx), commonData(common_data) {}
 
    PetscErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data);
  };
 
  struct OpMWLSStressAtGaussUsingPrecalulatedCoeffs
      : public OpMWLSBase<VolumeElementForcesAndSourcesCore> {
 
    std::string stressTagName;
    const bool calculateDerivative;
    const bool testing;
 
    OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
        boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr,
        boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr,
        boost::shared_ptr<CrackFrontElement> fe_singular_ptr,
        boost::shared_ptr<MWLSApprox> mwls_approx,
        const std::string stress_tag_name, bool calculate_derivative,
        bool testing = false)
        : OpMWLSBase<VolumeElementForcesAndSourcesCore>(
              mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_singular_ptr,
              mwls_approx, testing),
          stressTagName(stress_tag_name),
          calculateDerivative(calculate_derivative), testing(testing) {}
 
  protected:
    MoFEMErrorCode doMWLSWork(int side, EntityType type,
                              DataForcesAndSourcesCore::EntData &data);
  };
  /**
   * \brief Evaluate stress at integration points
   */
  struct OpMWLSStressAtGaussPts
      : OpMWLSBase<VolumeElementForcesAndSourcesCore> {
    std::string stressTagName;
    const bool calculateDerivative;
    const bool testing;
    boost::shared_ptr<MatrixDouble> matGradPosAtPtsPtr;
    OpMWLSStressAtGaussPts(
        boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr,
        boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr,
        boost::shared_ptr<CrackFrontElement> fe_singular_ptr,
        boost::shared_ptr<MWLSApprox> mwls_approx,
        const std::string stress_tag_name, bool calculate_derivative,
        bool testing = false)
        : OpMWLSBase<VolumeElementForcesAndSourcesCore>(
              mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_singular_ptr,
              mwls_approx, testing),
          matGradPosAtPtsPtr(mat_pos_at_pts_ptr),
          stressTagName(stress_tag_name),
          calculateDerivative(calculate_derivative), testing(testing) {}
 
    virtual ~OpMWLSStressAtGaussPts() = default;
 
  protected:
    MoFEMErrorCode doMWLSWork(int side, EntityType type,
                              DataForcesAndSourcesCore::EntData &data);
  };
 
  struct OpMWLSStressPostProcess
      : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
    boost::shared_ptr<moab::Interface> postProcMeshPtr;
    boost::shared_ptr<std::vector<EntityHandle>> mapGaussPtsPtr;
    boost::shared_ptr<MWLSApprox> mwlsApproxPtr;
 
    OpMWLSStressPostProcess(
        boost::shared_ptr<moab::Interface> post_proc_mesh,
        boost::shared_ptr<std::vector<EntityHandle>> map_gauss_pts,
        boost::shared_ptr<MWLSApprox> mwls_approx)
        : VolumeElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
          postProcMeshPtr(post_proc_mesh), mapGaussPtsPtr(map_gauss_pts),
          mwlsApproxPtr(mwls_approx) {}
 
    PetscErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data);
  };
 
  /**
   * \brief Evaluate stress at integration points
   */
  struct OpMWLSStressAndErrorsAtGaussPts
      : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
    std::string stressTagName;
    boost::shared_ptr<MWLSApprox> mwlsApproxPtr;
    MoFEM::Interface &mField;
 
    OpMWLSStressAndErrorsAtGaussPts(boost::shared_ptr<MWLSApprox> mwls_approx,
                                    const std::string stress_tag_name,
                                    MoFEM::Interface &m_field)
        : VolumeElementForcesAndSourcesCore::UserDataOperator(
              "SPATIAL_POSITION", UserDataOperator::OPROW),
          mwlsApproxPtr(mwls_approx), stressTagName(stress_tag_name),
          mField(m_field) {}
 
  protected:
    MoFEMErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data);
  };
 
  /**
   * \brief Integrate force vector
   */
  struct OpMWLSSpatialStressRhs
      : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
    boost::shared_ptr<MatrixDouble> matGradPosAtPtsPtr;
    boost::shared_ptr<MWLSApprox> mwlsApprox;
    const bool testing;
    OpMWLSSpatialStressRhs(
        boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr,
        boost::shared_ptr<MWLSApprox> mwls_approx, const bool testing = false)
        : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
              "SPATIAL_POSITION", UserDataOperator::OPROW),
          matGradPosAtPtsPtr(mat_grad_pos_at_pts_ptr), mwlsApprox(mwls_approx),
          testing(testing) {}
    VectorDouble nF;
    MoFEMErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data);
  };
 
  /**
   * \brief Integrate force vector
   */
  struct OpMWLSMaterialStressRhs
      : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
    boost::shared_ptr<MatrixDouble> spaceGradPosAtPtsPtr;
    boost::shared_ptr<MatrixDouble> matGradPosAtPtsPtr;
    boost::shared_ptr<MWLSApprox> mwlsApprox;
    Range forcesOnlyOnEntitiesRow;
    VectorInt iNdices;
    Vec &F_lambda;
    OpMWLSMaterialStressRhs(
        boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr,
        boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr,
        boost::shared_ptr<MWLSApprox> mwls_approx,
        Range *forces_only_on_entities_row = NULL)
        : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
          spaceGradPosAtPtsPtr(space_grad_pos_at_pts_ptr),
          matGradPosAtPtsPtr(mat_grad_pos_at_pts_ptr), mwlsApprox(mwls_approx),
          F_lambda(mwls_approx->F_lambda) {
      if (forces_only_on_entities_row)
        forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
    }
    VectorDouble nF;
    MoFEMErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data);
  };
 
  struct OpMWLSSpatialStressLhs_DX
      : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
    boost::shared_ptr<MatrixDouble> matGradPosAtPtsPtr;
    boost::shared_ptr<MWLSApprox> mwlsApprox;
    VectorInt iNdices;
    OpMWLSSpatialStressLhs_DX(
        boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr,
        boost::shared_ptr<MWLSApprox> mwls_approx)
        : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
              "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
              UserDataOperator::OPROWCOL),
          matGradPosAtPtsPtr(mat_grad_pos_at_pts_ptr), mwlsApprox(mwls_approx) {
      sYmm = false;
    }
    MatrixDouble nA;
    MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type,
                          DataForcesAndSourcesCore::EntData &row_data,
                          DataForcesAndSourcesCore::EntData &col_data);
  };
 
  struct OpMWLSMaterialStressLhs_Dx
      : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
    boost::shared_ptr<MatrixDouble> spaceGradPosAtPtsPtr;
    boost::shared_ptr<MatrixDouble> matGradPosAtPtsPtr;
    boost::shared_ptr<MWLSApprox> mwlsApprox;
    Range forcesOnlyOnEntitiesRow;
    VectorInt iNdices;
    OpMWLSMaterialStressLhs_Dx(
        boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr,
        boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr,
        boost::shared_ptr<MWLSApprox> mwls_approx,
        Range *forces_only_on_entities_row = nullptr)
        : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", "SPATIAL_POSITION",
              UserDataOperator::OPROWCOL),
          spaceGradPosAtPtsPtr(space_grad_pos_at_pts_ptr),
          matGradPosAtPtsPtr(mat_grad_pos_at_pts_ptr), mwlsApprox(mwls_approx) {
      sYmm = false;
      if (forces_only_on_entities_row)
        forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
    }
    MatrixDouble nA;
    MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type,
                          DataForcesAndSourcesCore::EntData &row_data,
                          DataForcesAndSourcesCore::EntData &col_data);
  };
 
  struct OpMWLSMaterialStressLhs_DX
      : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
    boost::shared_ptr<MatrixDouble> spaceGradPosAtPtsPtr;
    boost::shared_ptr<MatrixDouble> matGradPosAtPtsPtr;
    boost::shared_ptr<MWLSApprox> mwlsApprox;
    Range forcesOnlyOnEntitiesRow;
    VectorInt iNdices;
    OpMWLSMaterialStressLhs_DX(
        boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr,
        boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr,
        boost::shared_ptr<MWLSApprox> mwls_approx,
        Range *forces_only_on_entities_row = nullptr)
        : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
              "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
              UserDataOperator::OPROWCOL),
          spaceGradPosAtPtsPtr(space_grad_pos_at_pts_ptr),
          matGradPosAtPtsPtr(mat_grad_pos_at_pts_ptr), mwlsApprox(mwls_approx) {
      sYmm = false;
      if (forces_only_on_entities_row)
        forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
    }
    MatrixDouble nA;
    MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type,
                          DataForcesAndSourcesCore::EntData &row_data,
                          DataForcesAndSourcesCore::EntData &col_data);
  };
 
  /// \brief This is used to set up analytical function for testing
  /// approximation
  boost::function<VectorDouble(const double, const double, const double)>
      functionTestHack;
 
  inline MatrixDouble &getInvAB() { return invAB; }
  inline std::vector<EntityHandle> &getInfluenceNodes() {
    return influenceNodes;
  }
  inline auto &getShortenDm() { return shortenDm; }
 
private:
  boost::shared_ptr<BVHTree> myTree;
 
  // 1 (1) x y z (4) xx yy zz xy xz yz (10)
  int nbBasePolynomials; ///< Number of base functions (1 - linear, 4 - linear,
                         ///< 10 - quadratic)
 
  double dmFactor;        ///< Relative value of weight function radius
  double shortenDm;       ///< Radius of weight function used to calculate base
                          ///< function.
  int maxElemsInDMFactor; ///< Maximal number of elements in the horizon is
                          ///< calculated from max_nb_elems = maxElemsInDMFactor
                          ///< *nbBasePolynomials. Value is set by command line
                          ///< -mwls_max_elems_factor
 
  int maxThreeDepth; ///< Maximal three depths
  int splitsPerDir;  ///< Split per direction in the tree
 
  PetscBool useGlobalBaseAtMaterialReferenceConfiguration;
  PetscBool useNodalData;
  PetscBool mwlsAnalyticalInternalStressTest;
 
  EntityHandle treeRoot;
  double maxEdgeL;
 
  Range levelNodes, levelEdges, levelTets;
  std::vector<EntityHandle> influenceNodes;
  std::map<EntityHandle, std::array<double, 3>> elemCentreMap;
 
  EntityHandle nearestInfluenceNode;
 
  ublas::symmetric_matrix<double, ublas::lower> A;
  ublas::symmetric_matrix<double, ublas::lower> testA;
  MatrixDouble invA;
  ublas::triangular_matrix<double, ublas::lower> L;
  MatrixDouble B;
  VectorDouble baseFun;
  VectorDouble approxFun;
  MatrixDouble invAB;
  VectorDouble baseFunInvA;
 
  ublas::symmetric_matrix<double, ublas::lower, ublas::row_major,
                          DoubleAllocator>
      diffA[3];
  ublas::symmetric_matrix<double, ublas::lower, ublas::row_major,
                          DoubleAllocator>
      diffInvA[3];
  MatrixDouble diffB[3];
  VectorDouble diffBaseFun[3];
  VectorDouble diffDiffBaseFun[9];
  MatrixDouble diffApproxFun;
  MatrixDouble diffDiffApproxFun;
 
  MatrixDouble influenceNodesData;
  VecVals outData;
 
  double outHydroStaticError;
 
  VecVals outDeviatoricDifference;
  double outDeviatoricError;
 
  MatDiffVals outDiffData;
  MatDiffDiffVals outDiffDiffData;
 
  template <class T = FTensor::Tensor1<double, 3>>
  MoFEMErrorCode calculateBase(const T &t_coords) {
    MoFEMFunctionBeginHot;
 
    baseFun.resize(nbBasePolynomials, false);
    baseFun.clear();
 
    baseFun[0] = 1;
    if (nbBasePolynomials > 1) {
      baseFun[1] = t_coords(0); // x
      baseFun[2] = t_coords(1); // y
      baseFun[3] = t_coords(2); // z
    }
 
    if (nbBasePolynomials > 4) {
      baseFun[4] = t_coords(0) * t_coords(0); // x^2
      baseFun[5] = t_coords(1) * t_coords(1); // y^2
      baseFun[6] = t_coords(2) * t_coords(2); // z^2
      baseFun[7] = t_coords(0) * t_coords(1); // xy
      baseFun[8] = t_coords(0) * t_coords(2); // xz
      baseFun[9] = t_coords(1) * t_coords(2); // yz
    }
 
    MoFEMFunctionReturnHot(0);
  }
 
  template <class T = FTensor::Tensor1<double, 3>>
  MoFEMErrorCode calculateDiffBase(const T &t_coords, const double dm) {
    MoFEMFunctionBeginHot;
 
    // Local base functions
    for (int d = 0; d != 3; ++d) {
      diffBaseFun[d].resize(nbBasePolynomials, false);
      diffBaseFun[d].clear();
    }
 
    for (int d = 0; d != 3; d++) {
      diffBaseFun[d][0] = 0;
    }
 
    if (nbBasePolynomials > 1) {
      // derivative dx
      diffBaseFun[0][1] = 1. / dm;
      // derivative dy
      diffBaseFun[1][2] = 1. / dm;
      // derivative dz
      diffBaseFun[2][3] = 1. / dm;
    }
 
    if (nbBasePolynomials > 4) {
      // dx derivative
      diffBaseFun[0][4] = 2 * t_coords(0) / dm;
      diffBaseFun[0][5] = 0;
      diffBaseFun[0][6] = 0;
      diffBaseFun[0][7] = t_coords(1) / dm;
      diffBaseFun[0][8] = t_coords(2) / dm;
      diffBaseFun[0][9] = 0;
      // dy derivative
      diffBaseFun[1][4] = 0;
      diffBaseFun[1][5] = 2 * t_coords(1) / dm;
      diffBaseFun[1][6] = 0;
      diffBaseFun[1][7] = t_coords(0) / dm;
      diffBaseFun[1][8] = 0;
      diffBaseFun[1][9] = t_coords(2) / dm;
      // dz derivative
      diffBaseFun[2][4] = 0;
      diffBaseFun[2][5] = 0;
      diffBaseFun[2][6] = 2 * t_coords(2) / dm;
      diffBaseFun[2][7] = 0;
      diffBaseFun[2][8] = t_coords(0) / dm;
      diffBaseFun[2][9] = t_coords(1) / dm;
    }
    MoFEMFunctionReturnHot(0);
  }
 
  template <class T = FTensor::Tensor1<double, 3>>
  MoFEMErrorCode calculateDiffDiffBase(const T &t_coords, const double dm) {
    MoFEMFunctionBeginHot;
 
    // Local base functions
    for (int d = 0; d != 9; ++d) {
      diffDiffBaseFun[d].resize(nbBasePolynomials, false);
      diffDiffBaseFun[d].clear();
    }
 
    // if (nbBasePolynomials > 1) {
    // }
 
    if (nbBasePolynomials > 4) {
      const double dm2 = dm * dm;
      // set the non-zero values
      // dxd? derivative
      diffDiffBaseFun[0][4] = 2. / dm2; // dxdx
      diffDiffBaseFun[1][7] = 1. / dm2; // dxdy
      diffDiffBaseFun[2][8] = 1. / dm2; // dxdz
      // dyd? derivative
      diffDiffBaseFun[4][5] = 2. / dm2; // dydy
      diffDiffBaseFun[3][7] = 1. / dm2; // dydx
      diffDiffBaseFun[5][9] = 1. / dm2; // dydz
      // dzd? derivative
      diffDiffBaseFun[8][6] = 2. / dm2; // dzdz
      diffDiffBaseFun[6][8] = 1. / dm2; // dzdx
      diffDiffBaseFun[7][9] = 1. / dm2; // dzdy
    }
    MoFEMFunctionReturnHot(0);
  }
 
  inline double evalWeight(const double r) {
    const double rr = r * r;
    const double rrr = rr * r;
    const double r4 = rrr * r;
    if (r < 1) {
      return 1 - 6 * rr + 8 * rrr - 3 * r4;
    } else {
      return 0;
    }
    // if(r<=0.5) {
    //   return 2./3-4*rr+4*rrr;
    // } else if(r<=1) {
    //   return 4./3-4*r+4*rr-(4./3)*rrr;
    // } else return 0;
  }
 
  inline double evalDiffWeight(const double r) {
    const double rr = r * r;
    const double rrr = rr * r;
    if (r < 1) {
      return -12 * r + 24 * rr - 12 * rrr;
    } else {
      return 0;
    }
    // if(r<=0.5) {
    //   return -8*r+12*rr;
    // } else if(r<=1) {
    //   return -4+8*r-4*rr;
    // } else return 0;
  }
 
  inline double evalDiffDiffWeight(const double r) {
    const double rr = r * r;
    if (r < 1) {
      return -36 * rr + 48 * r - 12;
    } else {
      return 0;
    }
  }
 
  std::vector<EntityHandle> treeTets;
  std::vector<EntityHandle> leafsTetsVecLeaf;
  std::vector<double> leafsTetsCentre;
  std::vector<EntityHandle> leafsVec;
  std::vector<double> leafsDist;
  std::vector<std::pair<double, EntityHandle>> distLeafs;
 
};
 
} // namespace FractureMechanics
 
#endif //__MWLS_HPP__