v0.5.86
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
*/
Fe(MoFEM::Interface &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 {
PetscErrorCode preProcess() {
PetscFunctionBegin;
switch (ts_ctx) {
break;
}
default:
break;
}
PetscFunctionReturn(0);
}
PetscErrorCode postProcess() {
PetscFunctionBegin;
PetscFunctionReturn(0);
}
};
/**
* \brief Not used at this stage. Could be used to do some calculations,
* before assembly of local elements.
*/
struct FePrePostProcessLhs: public FEMethod {
PetscErrorCode preProcess() {
//
PetscFunctionBegin;
switch (ts_ctx) {
break;
}
default:
break;
}
PetscFunctionReturn(0);
}
PetscErrorCode postProcess() {
PetscFunctionBegin;
PetscFunctionReturn(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 pSi_ref; ///< reference free energy
double R0; ///< mass conduction coefficient
double cUrrent_psi; ///< current free energy for evauluating equilibrium state
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
PetscErrorCode getParameters();
// 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
VectorDouble vecRhoDt; ///< Time derivative of density
};
DataContainers data;
};
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
*/
PetscErrorCode getParameters();
/**
* \brief Set and add entities to approximation fields
* @return Error code
*/
PetscErrorCode addFields();
/**
* \brief Set and add finite elements
* @return Error code
*/
PetscErrorCode addElements();
/**
* \brief (Testing only) Set finite element to run mass transport problem only
* @return Error code
*/
PetscErrorCode addElementsTestingDensity();
/**
* \brief (Testing only) Set finite element to run elastic problem only
* @return Error code
*/
PetscErrorCode addElementsTestingElasticty();
/**
* \brief Finite elements to calculate tractions
* @return Error code
*/
PetscErrorCode addMomentumFluxes();
/**
* \brief Set problem and DM
* @return Error code
*/
PetscErrorCode buildDM();
/**
* \brief Solve problem set up in DM
* @return Error code
*/
PetscErrorCode 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 PetscErrorCode freeEnergy(
Remodeling::CommonData &common_data,T &C,B1 &psi
) {
PetscFunctionBegin;
const double mu = common_data.mu;
const double lambda = common_data.lambda;
// Compressible Neo-Hookean
B2 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);
PetscFunctionReturn(0);
}
template<class B1,class T>
inline PetscErrorCode recordFreeEnergy_dC(Remodeling::CommonData &common_data,T &dC,B1 &psi) {
PetscFunctionBegin;
trace_on(common_data.tAg,common_data.kEep);
common_data.aC(I,J) <<= dC(I,J); // Set independent variables
ierr = 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();
PetscFunctionReturn(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);
PetscErrorCode 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]
*/
Remodeling::CommonData &commonData;
OpCalculateStress(Remodeling::CommonData &common_data);
PetscErrorCode doWork(
int side,EntityType type,DataForcesAndSourcesCore::EntData &data
);
};
/**
* \brief Used to post proc stresses, energy, source term.
*
* calculation of principal stresses. for e.g.
*
*/
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
);
PetscErrorCode 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 OpCalulateStressTangent: public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
Remodeling::CommonData &commonData;
OpCalulateStressTangent(Remodeling::CommonData &common_data);
PetscErrorCode 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]
*/
Remodeling::CommonData &commonData;
VectorDouble nF; ///< Vector of the right hand side (internal forces)
OpAssmbleStressRhs(Remodeling::CommonData &common_data);
PetscErrorCode 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]
*/
Remodeling::CommonData &commonData;
VectorDouble nF; ///< Vector of the right hand side (internal forces)
OpAssmbleRhoRhs(Remodeling::CommonData &common_data);
PetscErrorCode 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]
*/
Remodeling::CommonData &commonData;
MatrixDouble locK_rho_rho;
MatrixDouble transLocK_rho_rho;
OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data);
PetscErrorCode 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]
*/
Remodeling::CommonData &commonData;
MatrixDouble locK_rho_F;
OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data);
PetscErrorCode 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]
*/
struct OpAssmbleStressLhs_dF: public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
Remodeling::CommonData &commonData;
MatrixDouble locK_P_F;
MatrixDouble transLocK_P_F;
OpAssmbleStressLhs_dF(Remodeling::CommonData &common_data);
PetscErrorCode 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);
PetscErrorCode 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 {
Remodeling::CommonData &commonData;
PostProcVolumeOnRefinedMesh postProcElastic;
//DensityMapFe densityMaps;
PetscBool mass_postproc;
PetscBool equilibrium_flg;
double rate;
bool iNit;
int pRT;
MoFEM::Interface &m_field,Remodeling::CommonData &common_data
);
PetscErrorCode preProcess();
PetscErrorCode operator()();
PetscErrorCode 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
):
energVec(energy_vec),
massVec(mass_vec),
commonData(common_data){
}
PetscErrorCode doWork(int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
);
};
}
#endif //__REMODELING_HPP__