v0.9.1
Remodeling.hpp
Go to the documentation of this file.
1 /** \file Remodeling.hpp
2 * \ingroup bone_remodelling
3 * \example Remodeling.hpp
4 
5 * \brief Implementation of element for bone remodeling
6 
7 */
8 
9 /*
10  * This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #ifndef __REMODELING_HPP__
25 #define __REMODELING_HPP__
26 
27 namespace BoneRemodeling {
28 
29 /**
30  * \brief Implementation of bone remodeling finite element
31  *
32  * Implementation base on paper \cite waffenschmidt2012anisotropic
33  * <http://biomechanics.stanford.edu/paper/IJSS12.pdf>
34  *
35  */
36 struct Remodeling {
37 
38  /**
39  * \brief Volume finite element
40  */
42  Fe(MoFEM::Interface &m_field)
44  /** \brief Set integrate rule
45  */
46  int getRule(int order) { return 2 * (order - 1) + 1; };
47  };
48 
49  /**
50  * \brief Not used at this stage. Could be used to do some calculations,
51  * before assembly of local elements.
52  */
53  struct FePrePostProcessRhs : public FEMethod {
54 
57  switch (ts_ctx) {
58  case CTX_TSSETIFUNCTION: {
59  snes_ctx = CTX_SNESSETFUNCTION;
60  snes_f = ts_F;
61  break;
62  }
63  default:
64  break;
65  }
67  }
68 
72  }
73  };
74 
75  /**
76  * \brief Not used at this stage. Could be used to do some calculations,
77  * before assembly of local elements.
78  */
79  struct FePrePostProcessLhs : public FEMethod {
80 
82  //
84  switch (ts_ctx) {
85  case CTX_TSSETIJACOBIAN: {
86  snes_ctx = CTX_SNESSETJACOBIAN;
87  snes_B = ts_B;
88  break;
89  }
90  default:
91  break;
92  }
94  }
95 
99  }
100  };
101 
102  /**
103  * Data structure for storing material parameters and evaluated values at
104  * integration points.
105  *
106  * @param oRder order of approximation
107  * @param lambda Lame parameter \f[ \lambda=\frac{E}{2(1+v)} \f]
108  * @param mu Lame parameter \f[ \mu=\frac{E\cdot\nu}{(1-2\nu)(1+\nu)} \f]
109  * @param c density evolution (growth) velocity [d/m^2] \f[ k_{p}^{\ast} \f]
110  * in Ellen Kuhl's paper
111  * @param m algorithmic exponent [-]
112  * @param n porosity exponent [-]
113  * @param rHo_ref reference density
114  * @param pSi_ref reference free energy
115  * @param R0 mass conduction coefficient
116  */
117  struct CommonData {
118 
119  Mat A;
120  Vec F, D;
121  // Approximation
122  int oRder;
124 
125  // DM
126  DMType dm_name; ///< dm (problem) name
127  DM dm; ///< Discretization manager
128  TS ts; ///< Time solver
129 
130  boost::shared_ptr<Fe> feLhs; ///< FE to make left hand side
131  boost::shared_ptr<Fe> feRhs; ///< FE to make right hand side
132  boost::shared_ptr<FePrePostProcessRhs> preProcRhs;
133  boost::shared_ptr<FePrePostProcessLhs> preProcLhs;
134 
135  boost::ptr_map<string, NeummanForcesSurface>
136  neumannForces; ///< Forces on surface
137  boost::ptr_map<string, NodalForce> nodalForces; ///< Nodal forces
138  boost::ptr_map<string, EdgeForce> edgeForces; ///< Forces on edges
139 
140  boost::shared_ptr<NonlinearElasticElement> elasticPtr;
141  boost::shared_ptr<ElasticMaterials> elasticMaterialsPtr;
142 
143  // Material parameters
144  double lambda; ///< Lame parameter
145  double mu; ///< Lame parameter
146  double c; ///< density evolution (growth) velocity [d/m^2]
147  double m; ///< algorithmic exponent [-]
148  double n; ///< porosity exponent [-]
149 
150  double rHo_ref; ///< reference density
151  double rHo_max; ///< max density
152  double rHo_min; ///< min density
153  int b; ///< b exponent for bell function
154  double pSi_ref; ///< reference free energy
155  double R0; ///< mass conduction coefficient
156 
157  double
158  cUrrent_psi; ///< current free energy for evaluating equilibrium state
159  double
160  cUrrent_mass; ///< current free energy for evaluating equilibrium state
161  PetscBool with_adol_c;
162  PetscBool is_atom_testing; ///< for atom tests
163  PetscBool less_post_proc; ///< reduce file size
165  Range tEts_all; // all blockset
166  Range tEts_block; // blockset with no remodelling
167 
169  inline double getCFromDensity(const double &rho);
170  inline double getCFromDensityDiff(const double &rho);
171  // aDouble values
173  aC; ///< right Cauchy-Green deformation tensor \f[ C=F^{T}\cdot F \f]
174  adouble aPsi; ///< Elastic energy
175  ///< \f[\psi_{0}=\frac{\mu}{2}\left(\textrm{tr}(\mathbf{C})-3\right)-\mu\ln(J)+\frac{\lambda}{2}\ln^{2}(\ln
176  ///< J)^{2} \f]
177 
178  const int tAg;
179  const int kEep;
180 
181  CommonData() : tAg(1), kEep(0) {}
182 
183  struct DataContainers {
184  boost::shared_ptr<MatrixDouble>
185  gradDispPtr; ///< Ptr to gradient of displacements matrix container
186  boost::shared_ptr<VectorDouble>
187  rhoPtr; ///< Ptr to density matrix container
188  boost::shared_ptr<MatrixDouble> gradRhoPtr; ///< Gradient of density field
189  VectorDouble vecPsi; ///< Elastic energy density
190  MatrixDouble matS; ///< 2nd Piola stress
191  MatrixDouble matP; ///< 1st Piola stress
192  VectorDouble vecR; ///< Mass sorce
193  MatrixDouble matMaterialTangent; ///< Tangent matrix (derivative for F and
194  ///< 1st Piola stress)
196  matPushedMaterialTangent; ///< Pushed tangent matrix (derivative for F
197  ///< and 1st Piola stress)
198  MatrixDouble matGradientOfDeformation; ///< Gradient of deformation
199  MatrixDouble matInvF; ///< inverse of deformation gradient
200  VectorDouble vecDetF; ///< determinant of F
201  VectorDouble vecRhoDt; ///< Time derivative of density
202  };
204  };
205 
208 
209  Remodeling(MoFEM::Interface &m_field, CommonData &common_data)
210  : mField(m_field), commonData(common_data) {}
211 
212  /**
213  * \brief Get parameters form line command or config file
214 
215  * Read command line and config file to setup material and model parameters.
216 
217  * @return Error code
218  */
220 
221  /**
222  * \brief Set and add entities to approximation fields
223  * @return Error code
224  */
226 
227  /**
228  * \brief Set and add finite elements
229  * @return Error code
230  */
232 
233  /**
234  * \brief (Testing only) Set finite element to run mass transport problem only
235  * @return Error code
236  */
238 
239  /**
240  * \brief (Testing only) Set finite element to run elastic problem only
241  * @return Error code
242  */
244 
245  /**
246  * \brief Finite elements to calculate tractions
247  * @return Error code
248  */
250 
251  /**
252  * \brief Set problem and DM
253  * @return Error code
254  */
256 
257  /**
258  * \brief Solve problem set up in DM
259  * @return Error code
260  */
262 };
263 
264 /**
265 * Calculate free energy
266 *
267 \f[\psi_{0}=\frac{\mu}{2}\left(\textrm{tr}(\mathbf{C})-3\right)-\mu\ln(J)+\frac{\lambda}{2}\ln^{2}(\ln
268 J) \f]
269 *
270 */
271 template <class B1, class B2, class T>
273  B1 &psi) {
274 
278  const double mu = common_data.mu;
279  const double lambda = common_data.lambda;
280 
281  // Compressible Neo-Hookean
282  B2 det;
284  det = sqrt(det);
285  B2 traceC;
286  traceC = C(I, I);
287  B2 log_det = log(det);
288  psi = 0.5 * mu * (traceC - 3.0) - mu * log_det +
289  0.5 * lambda * log_det * log_det;
290 
291  // St. Venant–Kirchhoff Material
292  // T E;
293  // E(I,J) = C(I,J);
294  // E(0,0) -= 1;
295  // E(1,1) -= 1;
296  // E(2,2) -= 1;
297  // E(I,J) *= 0.5;
298  // B2 traceE = E(I,I);
299  // psi = 0.5*lambda*traceE*traceE+mu*E(I,J)*E(I,J);
300 
302 }
303 
304 template <class B1, class T>
306  T &dC, B1 &psi) {
307 
311  trace_on(common_data.tAg, common_data.kEep);
312  common_data.aC(I, J) <<= dC(I, J); // Set independent variables
313  CHKERR freeEnergy<adouble, adouble, FTensor::Tensor2_symmetric<adouble, 3>>(
314  common_data, common_data.aC, common_data.aPsi);
315  double psi_val;
316  common_data.aPsi >>= psi_val; // Set dependent variables
317  psi = psi_val;
318  trace_off();
320 }
321 /**
322  * \brief Evaluate density derivative with respect to time
323  in case of Backward Euler Method
324  \f[
325  \frac{\partial \rho}{\partial t} = \frac{\rho_{n+1} - \rho_{n}}{\Delta t}
326  \f]
327  The density within a finite element is approximated directly as
328  \f[
329  \rho = N r
330  \f]
331  */
334 
336 
338 
339  MoFEMErrorCode doWork(int side, EntityType type,
341 };
342 
343 /**
344  * \brief Evaluate physical equations at integration points.
345 
346  \f[
347  F_{ij}=u_{i,j}+\delta_{ij} \\
348  C_{ij}=F_{k,i}F_{k,j} \\
349  P_{ij}=F_{i,k}S_{k,j} \\
350  \psi_{0}=\frac{\mu}{2}\left[\textrm{tr}(C)-3\right]-\mu\ln(\sqrt{\det
351  C})+\frac{\lambda}{2}\ln^{2}(\sqrt{\det C}) \\
352  R=c\left[\left[\frac{\rho}{\rho_0}\right]^{-m}\Psi-\Psi_0\right] \\
353  \f]
354 
355  */
358 
361 
363 
364  MoFEMErrorCode doWork(int side, EntityType type,
366 };
367 
368 /**
369  * \brief Used to post proc stresses, energy, source term.
370  *
371  * calculation of principal stresses. for e.g.
372  *
373  */
376 
379  std::vector<EntityHandle> &mapGaussPts;
380 
381  OpPostProcStress(moab::Interface &post_proc_mesh,
382  std::vector<EntityHandle> &map_gauss_pts,
383  Remodeling::CommonData &common_data);
384 
385  MoFEMErrorCode doWork(int side, EntityType type,
387 };
388 
389 /**
390  * Evaluate tangent material matrix
391 
392  \f[
393  C_{IJ} = F_{iI}F_{iJ}\\
394  P_{iI} = F_{iI}S_{IJ}(\mathbf{C})\\
395  \delta \Psi = \delta v_{iI}P_{iI} = \delta v_{iJ} F_{iI}S_{IJ}\\
396  D_{iJkL} \delta u_{k,L} = S_{IJ}F_{iI,kL} \delta u_{k,L} +
397  F_{iI}S_{IJ,KL}F_{kL}\delta u_{k,L}\\
398  D_{iJkL} \delta u_{k,L} = S_{IJ}\delta_{ik}\delta_{IL} \delta u_{k,L} +
399  F_{iI}S_{IJ,KL}F_{kL}\delta u_{k,L}\\
400  D_{iJkL} \delta u_{k,L} = \left( S_{IJ}\delta_{ik}\delta_{IL} +
401  F_{iI}S_{IJ,KL}F_{kL} \right)\delta u_{k,L}\\ \f]
402 
403  */
406 
408 
410 
411  MoFEMErrorCode doWork(int side, EntityType type,
413 };
414 
417 
420 
422 
423  MoFEMErrorCode doWork(int side, EntityType type,
425 };
426 // TODO correct equations in documentation
427 /**
428  * \brief Assemble residual for conservation of momentum (stresses)
429  *
430  \f[
431  R^{\phi}=\intop_{V}\left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n}
432  \mathbf{P}_{ij}\nabla N_{j}\,dV
433  \f]
434  */
437 
439  VectorDouble nF; ///< Vector of the right hand side (internal forces)
440 
442 
443  MoFEMErrorCode doWork(int side, EntityType type,
445 };
446 
447 /**
448  * \brief Assemble residual for conservation of mass (density)
449  \f[
450  R^{\rho}=\intop_{V}N\frac{\partial\rho}{\partial t}-N\mathbf{R}-R_{0}\nabla
451  N\nabla\rho\,dV \f]
452  */
455 
457  VectorDouble nF; ///< Vector of the right hand side (internal forces)
458 
460 
461  MoFEMErrorCode doWork(int side, EntityType type,
463 };
464 
465 /**
466  * \brief Diagonal block of tangent matrix \f$K_{\rho\rho}\f$
467  \f[
468  K^{\rho \rho}=\intop_{V} \Delta t N_i N_j - R_0 \nabla N_i \nabla N_j -c
469  \frac{n-m}{\rho_0}\left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n-m} \psi_0
470  N_i N_j \,dV \f]
471  */
474 
478 
480 
481  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
482  EntityType col_type,
485 };
486 
487 /**
488  * \brief Off diagonal block of tangent matrix \f$K_{\rho u}\f$
489  \f[
490  K_{\rho u}=-\intop_{V} c \left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n-m}
491  \nabla N_j P_{ij} N_i \,dV \f]
492  */
495 
498 
500  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
501  EntityType col_type,
504 };
505 
506 /**
507  * \brief Off diagonal block of tangent matrix \f$K_{u u}\f$
508  \f[
509  K_{u u}=\left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n} \nabla N_j
510  D_{ijkl}\nabla N_l S_{jl} \f]
511  */
512 template <bool ADOLC>
515 
520 
522 
523  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
524  EntityType col_type,
527 };
528 
529 /**
530  * \brief Off diagonal block of tangent matrix \f$K_{u \rho}\f$
531 /f[
532  K_{u \rho}=\intop_{V} \left[\frac{n}{\rho_{0}}\right]
533 \left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n} \nabla N_j P_{ij} N_i \,dV
534  /f]
535  */
538 
541 
543 
544  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
545  EntityType col_type,
548 };
549 
550 /**
551  * Element used to post-process results at each time step
552  */
553 struct MonitorPostProc : public FEMethod {
554 
559  // DensityMapFe densityMaps;
560  PetscBool mass_postproc;
561  PetscBool equilibrium_flg;
562  double rate;
563  bool iNit;
564  int pRT;
565 
567  Remodeling::CommonData &common_data);
568 
572 };
573 
574 /**
575  * Operator used to calculate potential energy at each time step, for
576  * postprocessing
577  */
580 
581  Vec energVec;
582  Vec massVec;
584 
585  OpMassAndEnergyCalculation(const string &field_name,
586  Remodeling::CommonData &common_data,
587  Vec energy_vec, Vec mass_vec)
589  field_name, UserDataOperator::OPROW),
590  energVec(energy_vec), massVec(mass_vec), commonData(common_data) {}
591 
592  MoFEMErrorCode doWork(int row_side, EntityType row_type,
594 };
595 
596 } // namespace BoneRemodeling
597 
598 #endif //__REMODELING_HPP__
PetscBool equilibrium_flg
Definition: Remodeling.hpp:561
Evaluate density derivative with respect to time in case of Backward Euler Method The density within...
Definition: Remodeling.hpp:332
moab::Interface & postProcMesh
Definition: Remodeling.hpp:378
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:377
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:820
MoFEMErrorCode getParameters()
Get parameters form line command or config file.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
boost::shared_ptr< VectorDouble > rhoPtr
Ptr to density matrix container.
Definition: Remodeling.hpp:187
Deprecated interface functions.
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:456
OpCalculateStress(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:280
MoFEMErrorCode freeEnergy(Remodeling::CommonData &common_data, T &C, B1 &psi)
Definition: Remodeling.hpp:272
PetscBool less_post_proc
reduce file size
Definition: Remodeling.hpp:163
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:482
OpCalculateStressTangent(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:717
double cUrrent_psi
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:158
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:892
Off diagonal block of tangent matrix .
Definition: Remodeling.hpp:493
MoFEMErrorCode addElementsTestingDensity()
(Testing only) Set finite element to run mass transport problem only
MoFEMErrorCode preProcess()
Definition: Remodeling.hpp:81
Assemble residual for conservation of mass (density) .
Definition: Remodeling.hpp:453
boost::shared_ptr< MatrixDouble > gradRhoPtr
Gradient of density field.
Definition: Remodeling.hpp:188
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:475
int getRule(int order)
Set integrate rule.
Definition: Remodeling.hpp:46
Definition: Remodeling.hpp:553
Diagonal block of tangent matrix .
Definition: Remodeling.hpp:472
VectorDouble vecRhoDt
Time derivative of density.
Definition: Remodeling.hpp:201
Off diagonal block of tangent matrix .
Definition: Remodeling.hpp:513
double pSi_ref
reference free energy
Definition: Remodeling.hpp:154
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
MoFEMErrorCode addMomentumFluxes()
Finite elements to calculate tractions.
bool iNit
Definition: Remodeling.hpp:563
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
double cUrrent_mass
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:160
MoFEMErrorCode postProcess()
Definition: Remodeling.hpp:69
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:732
double rate
Definition: Remodeling.hpp:562
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Remodeling.cpp:977
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpMassAndEnergyCalculation(const string &field_name, Remodeling::CommonData &common_data, Vec energy_vec, Vec mass_vec)
Definition: Remodeling.hpp:585
PetscBool mass_postproc
Definition: Remodeling.hpp:560
FTensor::Index< 'I', 3 > I
Definition: PlasticOps.hpp:70
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:550
OpCalculateStressTangentWithAdolc(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:382
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:516
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
MoFEMErrorCode solveDM()
Solve problem set up in DM.
MoFEMErrorCode addElements()
Set and add finite elements.
VectorDouble nF
Vector of the right hand side (internal forces)
Definition: Remodeling.hpp:457
double getCFromDensity(const double &rho)
Definition: Remodeling.cpp:218
MonitorPostProc(MoFEM::Interface &m_field, Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Off diagonal block of tangent matrix /f[ K_{u \rho}=\intop_{V} \left[\frac{n}{\rho_{0}}\right] \left...
Definition: Remodeling.hpp:536
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpGetRhoTimeDirevative(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:232
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:407
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:71
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:335
boost::shared_ptr< ElasticMaterials > elasticMaterialsPtr
Definition: Remodeling.hpp:141
MoFEMErrorCode preProcess()
MoFEMErrorCode preProcess()
Definition: Remodeling.hpp:55
MatrixDouble matGradientOfDeformation
Gradient of deformation.
Definition: Remodeling.hpp:198
OpAssmbleStressLhs_dRho(Remodeling::CommonData &common_data)
FTensor::Tensor2_symmetric< adouble, 3 > aC
right Cauchy-Green deformation tensor
Definition: Remodeling.hpp:173
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:539
OpAssmbleStressLhs_dF(Remodeling::CommonData &common_data)
Assemble residual for conservation of momentum (stresses)
Definition: Remodeling.hpp:435
boost::shared_ptr< MatrixDouble > gradDispPtr
Ptr to gradient of displacements matrix container.
Definition: Remodeling.hpp:185
MoFEMErrorCode operator()()
MoFEM::Interface & mField
Definition: Remodeling.hpp:555
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:239
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:359
MoFEMErrorCode addElementsTestingElasticty()
(Testing only) Set finite element to run elastic problem only
MoFEM::Interface & mField
Definition: Remodeling.hpp:206
Evaluate physical equations at integration points.
Definition: Remodeling.hpp:356
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
FTensor::Tensor2< double, 3, 3 > diffDiff
Definition: Remodeling.hpp:519
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::vector< EntityHandle > & mapGaussPts
Definition: Remodeling.hpp:379
double m
algorithmic exponent [-]
Definition: Remodeling.hpp:147
MoFEMErrorCode postProcess()
Definition: Remodeling.hpp:96
VectorDouble vecPsi
Elastic energy density.
Definition: Remodeling.hpp:189
double getCFromDensityDiff(const double &rho)
Definition: Remodeling.cpp:225
Post processing.
MoFEMErrorCode addFields()
Set and add entities to approximation fields.
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:496
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
Definition: Remodeling.hpp:79
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
Definition: Remodeling.hpp:136
int pRT
Definition: Remodeling.hpp:564
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
Definition: Remodeling.hpp:137
PetscBool is_atom_testing
for atom tests
Definition: Remodeling.hpp:162
#define CHKERR
Inline error check.
Definition: definitions.h:602
PostProcVolumeOnRefinedMesh postProc
Definition: Remodeling.hpp:557
MoFEMErrorCode recordFreeEnergy_dC(Remodeling::CommonData &common_data, T &dC, B1 &psi)
Definition: Remodeling.hpp:305
OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:969
constexpr int order
MoFEMErrorCode buildDM()
Set problem and DM.
OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data)
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
Fe(MoFEM::Interface &m_field)
Definition: Remodeling.hpp:42
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
boost::shared_ptr< FePrePostProcessLhs > preProcLhs
Definition: Remodeling.hpp:133
double n
porosity exponent [-]
Definition: Remodeling.hpp:148
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:438
DM dm
Discretization manager.
Definition: Remodeling.hpp:127
OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:886
Implementation of bone remodeling finite element.
Definition: Remodeling.hpp:36
double c
density evolution (growth) velocity [d/m^2]
Definition: Remodeling.hpp:146
Used to post proc stresses, energy, source term.
Definition: Remodeling.hpp:374
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
Definition: Remodeling.hpp:130
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
MoFEMErrorCode postProcess()
Volume finite element.
Definition: Remodeling.hpp:41
PostProcVolumeOnRefinedMesh postProcElastic
Definition: Remodeling.hpp:558
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:556
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Definition: Remodeling.hpp:140
VectorDouble nF
Vector of the right hand side (internal forces)
Definition: Remodeling.hpp:439
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
Definition: Remodeling.hpp:138
Remodeling(MoFEM::Interface &m_field, CommonData &common_data)
Definition: Remodeling.hpp:209
MatrixDouble matInvF
inverse of deformation gradient
Definition: Remodeling.hpp:199
int b
b exponent for bell function
Definition: Remodeling.hpp:153
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
Definition: Remodeling.hpp:53
boost::shared_ptr< FePrePostProcessRhs > preProcRhs
Definition: Remodeling.hpp:132
double R0
mass conduction coefficient
Definition: Remodeling.hpp:155
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:396
OpAssmbleStressRhs(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:814
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
Definition: Remodeling.hpp:131