Remodeling.hpp

Implementation of element for bone remodeling

/** \file Remodeling.hpp
* \ingroup bone_remodelling
* \example Remodeling.hpp
 
* \brief Implementation of element for bone remodeling
 
*/
 
/*
 * 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 __REMODELING_HPP__
#define __REMODELING_HPP__
 
namespace BoneRemodeling {
 
/**
 * \brief Implementation of bone remodeling finite element
 *
 * Implementation base on paper \cite waffenschmidt2012anisotropic
 * <http://biomechanics.stanford.edu/paper/IJSS12.pdf>
 *
 */
struct Remodeling {
 
  /**
   * \brief Volume finite element
   */
  struct Fe : public MoFEM::VolumeElementForcesAndSourcesCore {
    Fe(MoFEM::Interface &m_field)
        : MoFEM::VolumeElementForcesAndSourcesCore(m_field) {}
    /** \brief Set integrate rule
     */
    int getRule(int order) { return 2 * (order - 1) + 1; };
  };
 
  /**
   * \brief Not used at this stage. Could be used to do some calculations,
   * before assembly of local elements.
   */
  struct FePrePostProcessRhs : public FEMethod {
 
    MoFEMErrorCode preProcess() {
      MoFEMFunctionBegin;
      switch (ts_ctx) {
      case CTX_TSSETIFUNCTION: {
        snes_ctx = CTX_SNESSETFUNCTION;
        snes_f = ts_F;
        break;
      }
      default:
        break;
      }
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode postProcess() {
      MoFEMFunctionBegin;
      MoFEMFunctionReturn(0);
    }
  };
 
  /**
   * \brief Not used at this stage. Could be used to do some calculations,
   * before assembly of local elements.
   */
  struct FePrePostProcessLhs : public FEMethod {
 
    MoFEMErrorCode preProcess() {
      //
      MoFEMFunctionBegin;
      switch (ts_ctx) {
      case CTX_TSSETIJACOBIAN: {
        snes_ctx = CTX_SNESSETJACOBIAN;
        snes_B = ts_B;
        break;
      }
      default:
        break;
      }
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode postProcess() {
      MoFEMFunctionBegin;
      MoFEMFunctionReturn(0);
    }
  };
 
  /**
   * Data structure for storing material parameters and evaluated values at
   * integration points.
   *
   * @param  oRder  order of approximation
   * @param lambda Lame parameter  \f[ \lambda=\frac{E}{2(1+v)} \f]
   * @param mu Lame parameter \f[ \mu=\frac{E\cdot\nu}{(1-2\nu)(1+\nu)} \f]
   * @param c density evolution (growth) velocity [d/m^2] \f[ k_{p}^{\ast} \f]
   * in Ellen Kuhl's paper
   * @param m algorithmic exponent [-]
   * @param n porosity exponent [-]
   * @param rHo_ref reference density
   * @param pSi_ref reference free energy
   * @param R0 mass conduction coefficient
   */
  struct CommonData {
 
    Mat A;
    Vec F, D;
    // Approximation
    int oRder;
    BitRefLevel bitLevel;
 
    // DM
    DMType dm_name; ///< dm (problem) name
    DM dm;          ///< Discretization manager
    TS ts;          ///< Time solver
 
    boost::shared_ptr<Fe> feLhs; ///< FE to make left hand side
    boost::shared_ptr<Fe> feRhs; ///< FE to make right hand side
    boost::shared_ptr<FePrePostProcessRhs> preProcRhs;
    boost::shared_ptr<FePrePostProcessLhs> preProcLhs;
 
    boost::ptr_map<string, NeummanForcesSurface>
        neumannForces;                              ///< Forces on surface
    boost::ptr_map<string, NodalForce> nodalForces; ///< Nodal forces
    boost::ptr_map<string, EdgeForce> edgeForces;   ///< Forces on edges
 
    boost::shared_ptr<NonlinearElasticElement> elasticPtr;
    boost::shared_ptr<ElasticMaterials> elasticMaterialsPtr;
 
    // Material parameters
    double lambda; ///< Lame parameter
    double mu;     ///< Lame parameter
    double c;      ///< density evolution (growth) velocity [d/m^2]
    double m;      ///< algorithmic exponent [-]
    double n;      ///< porosity exponent [-]
 
    double rHo_ref; ///< reference density
    double rHo_max; ///< max density
    double rHo_min; ///< min density
    int b;          ///< b exponent for bell function
    double pSi_ref; ///< reference free energy
    double R0;      ///< mass conduction coefficient
 
    double
        cUrrent_psi; ///< current free energy for evaluating equilibrium state
    double
        cUrrent_mass; ///< current free energy for evaluating equilibrium state
    PetscBool with_adol_c;
    PetscBool is_atom_testing; ///< for atom tests
    PetscBool less_post_proc;  ///< reduce file size
    bool nOremodellingBlock;
    Range tEts_all;   // all blockset
    Range tEts_block; // blockset with no remodelling
 
    MoFEMErrorCode getParameters();
    inline double getCFromDensity(const double &rho);
    inline double getCFromDensityDiff(const double &rho);
    // aDouble values
    FTensor::Tensor2_symmetric<adouble, 3>
        aC; ///< right Cauchy-Green deformation tensor \f[ C=F^{T}\cdot F \f]
    adouble aPsi; ///< Elastic energy
                  ///< \f[\psi_{0}=\frac{\mu}{2}\left(\textrm{tr}(\mathbf{C})-3\right)-\mu\ln(J)+\frac{\lambda}{2}\ln^{2}(\ln
                  ///< J)^{2} \f]
 
    const int tAg;
    const int kEep;
 
    CommonData() : tAg(1), kEep(0) {}
 
    struct DataContainers {
      boost::shared_ptr<MatrixDouble>
          gradDispPtr; ///< Ptr to gradient of displacements matrix container
      boost::shared_ptr<VectorDouble>
          rhoPtr; ///< Ptr to density matrix container
      boost::shared_ptr<MatrixDouble> gradRhoPtr; ///< Gradient of density field
      VectorDouble vecPsi;                        ///< Elastic energy density
      MatrixDouble matS;                          ///< 2nd Piola stress
      MatrixDouble matP;                          ///< 1st Piola stress
      VectorDouble vecR;                          ///< Mass sorce
      MatrixDouble matMaterialTangent; ///< Tangent matrix (derivative for F and
                                       ///< 1st Piola stress)
      MatrixDouble
          matPushedMaterialTangent; ///< Pushed tangent matrix (derivative for F
                                    ///< and 1st Piola stress)
      MatrixDouble matGradientOfDeformation; ///< Gradient of deformation
      MatrixDouble matInvF;                          ///< inverse of deformation gradient
      VectorDouble vecDetF;                          ///< determinant of F
      VectorDouble vecRhoDt;                 ///< Time derivative of density
    };
    DataContainers data;
  };
 
  MoFEM::Interface &mField;
  CommonData &commonData;
 
  Remodeling(MoFEM::Interface &m_field, CommonData &common_data)
      : mField(m_field), commonData(common_data) {}
 
  /**
   * \brief Get parameters form line command or config file
 
   * Read command line and config file to setup material and model parameters.
 
   * @return Error code
   */
  MoFEMErrorCode getParameters();
 
  /**
   * \brief Set and add entities to approximation fields
   * @return Error code
   */
  MoFEMErrorCode addFields();
 
  /**
   * \brief Set and add finite elements
   * @return Error code
   */
  MoFEMErrorCode addElements();
 
  /**
   * \brief (Testing only) Set finite element to run mass transport problem only
   * @return Error code
   */
  MoFEMErrorCode addElementsTestingDensity();
 
  /**
   * \brief (Testing only) Set finite element to run elastic problem only
   * @return Error code
   */
  MoFEMErrorCode addElementsTestingElasticty();
 
  /**
   * \brief Finite elements to calculate tractions
   * @return Error code
   */
  MoFEMErrorCode addMomentumFluxes();
 
  /**
   * \brief Set problem and DM
   * @return Error code
   */
  MoFEMErrorCode buildDM();
 
  /**
   * \brief Solve problem set up in DM
   * @return    Error code
   */
  MoFEMErrorCode solveDM();
};
 
/**
* Calculate free energy
*
\f[\psi_{0}=\frac{\mu}{2}\left(\textrm{tr}(\mathbf{C})-3\right)-\mu\ln(J)+\frac{\lambda}{2}\ln^{2}(\ln
J) \f]
*
*/
template <class B1, class B2, class T>
inline MoFEMErrorCode freeEnergy(Remodeling::CommonData &common_data, T &C,
                                 B1 &psi) {
 
  MoFEMFunctionBeginHot;
  FTensor::Index<'I', 3> I;
  FTensor::Index<'J', 3> J;
  const double mu = common_data.mu;
  const double lambda = common_data.lambda;
 
  // Compressible Neo-Hookean
  B2 det;
  CHKERR determinantTensor3by3(C, det);
  det = sqrt(det);
  B2 traceC;
  traceC = C(I, I);
  B2 log_det = log(det);
  psi = 0.5 * mu * (traceC - 3.0) - mu * log_det +
        0.5 * lambda * log_det * log_det;
 
  // St. Venant–Kirchhoff Material
  // T E;
  // E(I,J) = C(I,J);
  // E(0,0) -= 1;
  // E(1,1) -= 1;
  // E(2,2) -= 1;
  // E(I,J) *= 0.5;
  // B2 traceE = E(I,I);
  // psi = 0.5*lambda*traceE*traceE+mu*E(I,J)*E(I,J);
 
  MoFEMFunctionReturnHot(0);
}
 
template <class B1, class T>
inline MoFEMErrorCode recordFreeEnergy_dC(Remodeling::CommonData &common_data,
                                          T &dC, B1 &psi) {
 
  MoFEMFunctionBeginHot;
  FTensor::Index<'I', 3> I;
  FTensor::Index<'J', 3> J;
  trace_on(common_data.tAg, common_data.kEep);
  common_data.aC(I, J) <<= dC(I, J); // Set independent variables
  CHKERR freeEnergy<adouble, adouble, FTensor::Tensor2_symmetric<adouble, 3>>(
      common_data, common_data.aC, common_data.aPsi);
  double psi_val;
  common_data.aPsi >>= psi_val; // Set dependent variables
  psi = psi_val;
  trace_off();
  MoFEMFunctionReturnHot(0);
}
/**
 * \brief Evaluate density derivative with respect to time
  in case of Backward Euler Method
  \f[
  \frac{\partial \rho}{\partial t} = \frac{\rho_{n+1} - \rho_{n}}{\Delta t}
  \f]
  The density within a finite element is approximated directly as
  \f[
  \rho = N r
  \f]
 */
struct OpGetRhoTimeDirevative
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Remodeling::CommonData &commonData;
 
  OpGetRhoTimeDirevative(Remodeling::CommonData &common_data);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
/**
 * \brief Evaluate physical equations at integration points.
 
 \f[
 F_{ij}=u_{i,j}+\delta_{ij} \\
 C_{ij}=F_{k,i}F_{k,j} \\
 P_{ij}=F_{i,k}S_{k,j} \\
 \psi_{0}=\frac{\mu}{2}\left[\textrm{tr}(C)-3\right]-\mu\ln(\sqrt{\det
 C})+\frac{\lambda}{2}\ln^{2}(\sqrt{\det C}) \\
 R=c\left[\left[\frac{\rho}{\rho_0}\right]^{-m}\Psi-\Psi_0\right] \\
 \f]
 
 */
struct OpCalculateStress
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Remodeling::CommonData &commonData;
  VectorDouble vecC;
 
  OpCalculateStress(Remodeling::CommonData &common_data);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
/**
 * \brief Used to post proc stresses, energy, source term.
 *
 * calculation of principal stresses. for e.g.
 *
 */
struct OpPostProcStress
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Remodeling::CommonData &commonData;
  moab::Interface &postProcMesh;
  std::vector<EntityHandle> &mapGaussPts;
 
  OpPostProcStress(moab::Interface &post_proc_mesh,
                   std::vector<EntityHandle> &map_gauss_pts,
                   Remodeling::CommonData &common_data);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
/**
 * Evaluate tangent material matrix
 
 \f[
 C_{IJ} = F_{iI}F_{iJ}\\
 P_{iI} = F_{iI}S_{IJ}(\mathbf{C})\\
 \delta \Psi =  \delta v_{iI}P_{iI} = \delta v_{iJ} F_{iI}S_{IJ}\\
 D_{iJkL} \delta u_{k,L} = S_{IJ}F_{iI,kL} \delta u_{k,L} +
 F_{iI}S_{IJ,KL}F_{kL}\delta u_{k,L}\\
 D_{iJkL} \delta u_{k,L} = S_{IJ}\delta_{ik}\delta_{IL} \delta u_{k,L} +
 F_{iI}S_{IJ,KL}F_{kL}\delta u_{k,L}\\
 D_{iJkL} \delta u_{k,L} = \left( S_{IJ}\delta_{ik}\delta_{IL} +
 F_{iI}S_{IJ,KL}F_{kL} \right)\delta u_{k,L}\\ \f]
 
 */
struct OpCalculateStressTangent
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Remodeling::CommonData &commonData;
 
  OpCalculateStressTangent(Remodeling::CommonData &common_data);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
struct OpCalculateStressTangentWithAdolc
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Remodeling::CommonData &commonData;
  VectorDouble vecC;
 
  OpCalculateStressTangentWithAdolc(Remodeling::CommonData &common_data);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
// TODO correct equations in documentation
/**
 * \brief Assemble residual for conservation of momentum (stresses)
 *
 \f[
  R^{\phi}=\intop_{V}\left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n}
  \mathbf{P}_{ij}\nabla N_{j}\,dV
 \f]
 */
struct OpAssmbleStressRhs
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Remodeling::CommonData &commonData;
  VectorDouble nF; ///< Vector of the right hand side (internal forces)
 
  OpAssmbleStressRhs(Remodeling::CommonData &common_data);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
/**
 * \brief Assemble residual for conservation of mass (density)
 \f[
 R^{\rho}=\intop_{V}N\frac{\partial\rho}{\partial t}-N\mathbf{R}-R_{0}\nabla
 N\nabla\rho\,dV \f]
 */
struct OpAssmbleRhoRhs
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Remodeling::CommonData &commonData;
  VectorDouble nF; ///< Vector of the right hand side (internal forces)
 
  OpAssmbleRhoRhs(Remodeling::CommonData &common_data);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};
 
/**
 * \brief Diagonal block of tangent matrix \f$K_{\rho\rho}\f$
  \f[
   K^{\rho \rho}=\intop_{V} \Delta t N_i N_j - R_0 \nabla N_i \nabla N_j -c
  \frac{n-m}{\rho_0}\left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n-m} \psi_0
 N_i N_j \,dV \f]
 */
struct OpAssmbleRhoLhs_dRho
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Remodeling::CommonData &commonData;
  MatrixDouble locK_rho_rho;
  MatrixDouble transLocK_rho_rho;
 
  OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data);
 
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data);
};
 
/**
 * \brief Off diagonal block of tangent matrix \f$K_{\rho u}\f$
  \f[
  K_{\rho u}=-\intop_{V} c \left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n-m}
 \nabla N_j P_{ij} N_i \,dV \f]
 */
struct OpAssmbleRhoLhs_dF
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Remodeling::CommonData &commonData;
  MatrixDouble locK_rho_F;
 
  OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data);
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data);
};
 
/**
 * \brief Off diagonal block of tangent matrix \f$K_{u u}\f$
  \f[
    K_{u u}=\left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n} \nabla N_j
 D_{ijkl}\nabla N_l S_{jl} \f]
 */
template <bool ADOLC>
struct OpAssmbleStressLhs_dF
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Remodeling::CommonData &commonData;
  MatrixDouble locK_P_F;
  MatrixDouble transLocK_P_F;
  FTensor::Tensor2<double, 3, 3> diffDiff;
 
  OpAssmbleStressLhs_dF(Remodeling::CommonData &common_data);
 
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data);
};
 
/**
 * \brief Off diagonal block of tangent matrix \f$K_{u \rho}\f$
/f[
 K_{u \rho}=\intop_{V} \left[\frac{n}{\rho_{0}}\right]
\left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n}  \nabla N_j P_{ij} N_i \,dV
 /f]
 */
struct OpAssmbleStressLhs_dRho
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Remodeling::CommonData &commonData;
  MatrixDouble locK_P_RHO;
 
  OpAssmbleStressLhs_dRho(Remodeling::CommonData &common_data);
 
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data);
};
 
/**
 * Element used to post-process results at each time step
 */
struct MonitorPostProc : public FEMethod {
 
  MoFEM::Interface &mField;
  Remodeling::CommonData &commonData;
  PostProcVolumeOnRefinedMesh postProc;
  PostProcVolumeOnRefinedMesh postProcElastic;
  // DensityMapFe densityMaps;
  PetscBool mass_postproc;
  PetscBool equilibrium_flg;
  double rate;
  bool iNit;
  int pRT;
 
  MonitorPostProc(MoFEM::Interface &m_field,
                  Remodeling::CommonData &common_data);
 
  MoFEMErrorCode preProcess();
  MoFEMErrorCode operator()();
  MoFEMErrorCode postProcess();
};
 
/**
 * Operator used to calculate potential energy at each time step, for
 * postprocessing
 */
struct OpMassAndEnergyCalculation
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  Vec energVec;
  Vec massVec;
  Remodeling::CommonData &commonData;
 
  OpMassAndEnergyCalculation(const string &field_name,
                             Remodeling::CommonData &common_data,
                             Vec energy_vec, Vec mass_vec)
      : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
            field_name, UserDataOperator::OPROW),
        energVec(energy_vec), massVec(mass_vec), commonData(common_data) {}
 
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        DataForcesAndSourcesCore::EntData &row_data);
};
 
} // namespace BoneRemodeling
 
#endif //__REMODELING_HPP__