v0.9.0
Classes | Functions
BoneRemodeling Namespace Reference

Classes

struct  CommonData
 Data structure for storing global and local matrix K and vector f. More...
 
struct  DataFromMetaIO
 Load data from MetaImage file, translate grayscale values to densities. More...
 
struct  DensityMapFe
 
struct  MonitorPostProc
 
struct  OpAssmbleRhoLhs_dF
 Off diagonal block of tangent matrix \(K_{\rho u}\)

\[ K_{\rho u}=-\intop_{V} c \left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n-m} \nabla N_j P_{ij} N_i \,dV \]

. More...

 
struct  OpAssmbleRhoLhs_dRho
 Diagonal block of tangent matrix \(K_{\rho\rho}\)

\[ 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 \]

. More...

 
struct  OpAssmbleRhoRhs
 Assemble residual for conservation of mass (density)

\[ R^{\rho}=\intop_{V}N\frac{\partial\rho}{\partial t}-N\mathbf{R}-R_{0}\nabla N\nabla\rho\,dV \]

. More...

 
struct  OpAssmbleRhs
 Assemble RHS vector f.

\[ \mathbf{f}_{e}=\int_{V_{e}}\phi q^{\ast}\,dV_{e} \]

. More...

 
struct  OpAssmbleStressLhs_dF
 Off diagonal block of tangent matrix \(K_{u u}\)

\[ K_{u u}=\left[\frac{\rho_{0}}{\rho_{0}^{\ast}}\right]^{n} \nabla N_j D_{ijkl}\nabla N_l S_{jl} \]

. More...

 
struct  OpAssmbleStressLhs_dRho
 Off diagonal block of tangent matrix \(K_{u \rho}\) /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]. More...
 
struct  OpAssmbleStressRhs
 Assemble residual for conservation of momentum (stresses) More...
 
struct  OpCalculateStress
 Evaluate physical equations at integration points. More...
 
struct  OpCalculateStressTangent
 
struct  OpCalculateStressTangentWithAdolc
 
struct  OpCalulatefRhoAtGaussPts
 Assemble local vector containing density data. More...
 
struct  OpCalulateLhs
 Assemble LHS matrix K.

\[ K=\int_{V}\phi^{T}\phi\,dV \]

. More...

 
struct  OpGetRhoTimeDirevative
 Evaluate density derivative with respect to time in case of Backward Euler Method

\[ \frac{\partial \rho}{\partial t} = \frac{\rho_{n+1} - \rho_{n}}{\Delta t} \]

The density within a finite element is approximated directly as

\[ \rho = N r \]

. More...

 
struct  OpMassAndEnergyCalculation
 
struct  OpMassCalculation
 Calculate mass before approximation.

\[ M=\int_{V}\rho\,dV \]

. More...

 
struct  OpMassCalculationFromApprox
 Calculate mass after approximation.

\[ M=\int_{V}\rho N\,dV \]

. More...

 
struct  OpPostProcStress
 Used to post proc stresses, energy, source term. More...
 
struct  OpVolumeCalculation
 Calculate volume of the model. More...
 
struct  Remodeling
 Implementation of bone remodeling finite element. More...
 
struct  SurfaceKDTree
 Create KDTree on surface of the mesh and calculate distance. More...
 

Functions

template<class B1 , class B2 , class T >
MoFEMErrorCode freeEnergy (Remodeling::CommonData &common_data, T &C, B1 &psi)
 
template<class B1 , class T >
MoFEMErrorCode recordFreeEnergy_dC (Remodeling::CommonData &common_data, T &dC, B1 &psi)
 
PetscErrorCode sphereSurfaceIntegration28 (MatrixDouble &gauss_pts, VectorDouble &weights)
 28 integration points on sphere More...
 
PetscErrorCode sphereSurfaceIntegration21 (MatrixDouble &gauss_pts, VectorDouble &weights)
 21 integration points on sphere More...
 
PetscErrorCode sphereSurfaceIntegration61 (MatrixDouble &gauss_pts, VectorDouble &weights)
 61 integration points on sphere More...
 

Function Documentation

◆ freeEnergy()

template<class B1 , class B2 , class T >
MoFEMErrorCode BoneRemodeling::freeEnergy ( Remodeling::CommonData common_data,
T C,
B1 &  psi 
)

Calculate free energy

\[\psi_{0}=\frac{\mu}{2}\left(\textrm{tr}(\mathbf{C})-3\right)-\mu\ln(J)+\frac{\lambda}{2}\ln^{2}(\ln J) \]

Examples
Remodeling.hpp.

Definition at line 272 of file Remodeling.hpp.

273  {
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 }
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:415
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define CHKERR
Inline error check.
Definition: definitions.h:596

◆ recordFreeEnergy_dC()

template<class B1 , class T >
MoFEMErrorCode BoneRemodeling::recordFreeEnergy_dC ( Remodeling::CommonData common_data,
T dC,
B1 &  psi 
)
Examples
Remodeling.hpp.

Definition at line 305 of file Remodeling.hpp.

306  {
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 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define CHKERR
Inline error check.
Definition: definitions.h:596

◆ sphereSurfaceIntegration21()

PetscErrorCode RemodelingFE::sphereSurfaceIntegration21 ( MatrixDouble gauss_pts,
VectorDouble weights 
)

21 integration points on sphere

Parameters
gauss_ptsreferenced matrix consisting normals at integration points
weightsWeights at gauss points
Returns
Error code

Definition at line 69 of file SphereSurfaceIntegration.cpp.

71  {
72  PetscFunctionBegin;
73  weights.resize(21,false);
74  gauss_pts.resize(21,3,false);
75  for(int i = 0; i != 3; i++ ) {
76  weights[ i ] = 0.02652141274;
77  }
78  for(int i = 3; i != 9; i++ ) {
79  weights[ i ] = 0.01993014153;
80  }
81  for(int i = 9; i != 21; i++ ) {
82  weights[ i ] = 0.02507124272;
83  }
84  gauss_pts(0,0) = gauss_pts(1,1) = gauss_pts(2,2) = 1.;
85  gauss_pts(0,1) = gauss_pts(0,2) = gauss_pts(1,0) = gauss_pts(1,2) = gauss_pts(2,0) = gauss_pts(2,1) = 0.;
86  gauss_pts(3,0) = gauss_pts(3,1) = gauss_pts(4,0) = gauss_pts(5,0) = gauss_pts(5,2) = gauss_pts(6,0) = gauss_pts(7,1) = gauss_pts(7,2) = gauss_pts(8,1) = 0.7071067812;
87  gauss_pts(4,1) = gauss_pts(6,2) = gauss_pts(8,2) = -0.7071067812;
88  gauss_pts(3,2) = gauss_pts(4,2) = gauss_pts(5,1) = gauss_pts(6,1) = gauss_pts(7,0) = gauss_pts(8,0) = 0.;
89  double help [ 3 ] [ 3 ] =
90  { { 0.3879072746, 0.3879072746, 0.8360956240 },
91  { 0.3879072746, 0.8360956240, 0.3879072746 },
92  { 0.8360956240, 0.3879072746, 0.3879072746 } };
93 
94  for ( i = 0; i < 3; i++ ) {
95  i4 = i * 4;
96  gauss_pts(9 + i4,0) = help [ i ] [ 0 ];
97  gauss_pts(9 + i4,1) = help [ i ] [ 1 ];
98  gauss_pts(9 + i4,2) = help [ i ] [ 2 ];
99 
100  gauss_pts(10 + i4,0) = help [ i ] [ 0 ];
101  gauss_pts(10 + i4,1) = help [ i ] [ 1 ];
102  gauss_pts(10 + i4,2) = -help [ i ] [ 2 ];
103 
104  gauss_pts(11 + i4,0) = help [ i ] [ 0 ];
105  gauss_pts(11 + i4,1) = -help [ i ] [ 1 ];
106  gauss_pts(11 + i4,2) = help [ i ] [ 2 ];
107 
108  gauss_pts(12 + i4,0) = help [ i ] [ 0 ];
109  gauss_pts(12 + i4,1) = -help [ i ] [ 1 ];
110  gauss_pts(12 + i4,2) = -help [ i ] [ 2 ];
111  }
112  PetscFunctionReturn(0);
113 }
static char help[]

◆ sphereSurfaceIntegration28()

PetscErrorCode RemodelingFE::sphereSurfaceIntegration28 ( MatrixDouble gauss_pts,
VectorDouble weights 
)

28 integration points on sphere

Parameters
gauss_ptsreferenced matrix consisting normals at integration points
weightsWeights at gauss points
Returns
Error code

Definition at line 30 of file SphereSurfaceIntegration.cpp.

32  {
33  PetscFunctionBegin;
34  gauss_pts.resize(28,3,false);
35  weights.resize(28,false);
36  for(int i = 0; i != 4; i++ ) {
37  weights[ i ] = 0.0160714276E0;
38  weights[ i + 4 ] = 0.0204744730E0;
39  weights[ i + 8 ] = 0.0204744730E0;
40  weights[ i + 12 ] = 0.0204744730E0;
41  weights[ i + 16 ] = 0.0158350505E0;
42  weights[ i + 20 ] = 0.0158350505E0;
43  weights[ i + 24 ] = 0.0158350505E0;
44  }
45  double help [ 7 ] [ 3 ] = {
46  { .577350259E0, .577350259E0, .577350259E0 }, { .935113132E0, .250562787E0, .250562787E0 },
47  { .250562787E0, .935113132E0, .250562787E0 }, { .250562787E0, .250562787E0, .935113132E0 },
48  { .186156720E0, .694746614E0, .694746614E0 }, { .694746614E0, .186156720E0, .694746614E0 },
49  { .694746614E0, .694746614E0, .186156720E0 }
50  };
51  for(int i = 0; i != 7; i++ ) {
52  i4 = i * 4;
53  gauss_pts(i4,0) = help [ i ] [ 0 ];
54  gauss_pts(i4,1) = help [ i ] [ 1 ];
55  gauss_pts(i4,2) = help [ i ] [ 2 ];
56  gauss_pts(i4 + 1,0) = help [ i ] [ 0 ];
57  gauss_pts(i4 + 1,1) = help [ i ] [ 1 ];
58  gauss_pts(i4 + 1,2) = -help [ i ] [ 2 ];
59  gauss_pts(i4 + 2,0) = help [ i ] [ 0 ];
60  gauss_pts(i4 + 2,1) = -help [ i ] [ 1 ];
61  gauss_pts(i4 + 2,2) = help [ i ] [ 2 ];
62  gauss_pts(i4 + 3,0) = help [ i ] [ 0 ];
63  gauss_pts(i4 + 3,1) = -help [ i ] [ 1 ];
64  gauss_pts(i4 + 3,2) = -help [ i ] [ 2 ];
65  }
66  PetscFunctionReturn(0);
67 }
static char help[]

◆ sphereSurfaceIntegration61()

PetscErrorCode RemodelingFE::sphereSurfaceIntegration61 ( MatrixDouble gauss_pts,
VectorDouble weights 
)

61 integration points on sphere

Parameters
gauss_ptsreferenced matrix consisting normals at integration points
weightsWeights at gauss points
Returns
Error code

Definition at line 115 of file SphereSurfaceIntegration.cpp.

117  {
118  PetscFunctionBegin;
119  double help [ 61 ] [ 4 ] = {
120  { 1.000000000000, 0.000000000000, 0.000000000000, 0.00795844204678 },
121  { 0.745355992500, 0.0, 0.666666666667, 0.00795844204678 },
122  { 0.745355992500, -0.577350269190, -0.333333333333, 0.00795844204678 },
123  { 0.745355992500, 0.577350269190, -0.333333333333, 0.00795844204678 },
124  { 0.333333333333, 0.577350269190, 0.745355992500, 0.00795844204678 },
125  { 0.333333333333, -0.577350269190, 0.745355992500, 0.00795844204678 },
126  { 0.333333333333, -0.934172358963, 0.127322003750, 0.00795844204678 },
127  { 0.333333333333, -0.356822089773, -0.872677996250, 0.00795844204678 },
128  { 0.333333333333, 0.356822089773, -0.872677996250, 0.00795844204678 },
129  { 0.333333333333, 0.934172358963, 0.127322003750, 0.00795844204678 },
130  { 0.794654472292, -0.525731112119, 0.303530999103, 0.0105155242892 },
131  { 0.794654472292, 0.0, -0.607061998207, 0.0105155242892 },
132  { 0.794654472292, 0.525731112119, 0.303530999103, 0.0105155242892 },
133  { 0.187592474085, 0.0, 0.982246946377, 0.0105155242892 },
134  { 0.187592474085, -0.850650808352, -0.491123473188, 0.0105155242892 },
135  { 0.187592474085, 0.850650808352, -0.491123473188, 0.0105155242892 },
136  { 0.934172358963, 0.0, 0.356822089773, 0.0100119364272 },
137  { 0.934172358963, -0.309016994375, -0.178411044887, 0.0100119364272 },
138  { 0.934172358963, 0.309016994375, -0.178411044887, 0.0100119364272 },
139  { 0.577350269190, 0.309016994375, 0.755761314076, 0.0100119364272 },
140  { 0.577350269190, -0.309016994375, 0.755761314076, 0.0100119364272 },
141  { 0.577350269190, -0.809016994375, -0.110264089708, 0.0100119364272 },
142  { 0.577350269190, -0.5, -0.645497224368, 0.0100119364272 },
143  { 0.577350269190, 0.5, -0.645497224368, 0.0100119364272 },
144  { 0.577350269190, 0.809016994375, -0.110264089708, 0.0100119364272 },
145  { 0.356822089773, -0.809016994375, 0.467086179481, 0.0100119364272 },
146  { 0.356822089773, 0.0, -0.934172358963, 0.0100119364272 },
147  { 0.356822089773, 0.809016994375, 0.467086179481, 0.0100119364272 },
148  { 0.0, 0.5, 0.866025403784, 0.0100119364272 },
149  { 0.0, -1.0, 0.0, 0.0100119364272 },
150  { 0.0, 0.5, -0.866025403784, 0.0100119364272 },
151  { 0.947273580412, -0.277496978165, 0.160212955043, 0.00690477957966 },
152  { 0.812864676392, -0.277496978165, 0.512100034157, 0.00690477957966 },
153  { 0.595386501297, -0.582240127941, 0.553634669695, 0.00690477957966 },
154  { 0.595386501297, -0.770581752342, 0.227417407053, 0.00690477957966 },
155  { 0.812864676392, -0.582240127941, -0.015730584514, 0.00690477957966 },
156  { 0.492438766306, -0.753742692223, -0.435173546254, 0.00690477957966 },
157  { 0.274960591212, -0.942084316623, -0.192025554687, 0.00690477957966 },
158  { -0.076926487903, -0.942084316623, -0.326434458707, 0.00690477957966 },
159  { -0.076926487903, -0.753742692223, -0.652651721349, 0.00690477957966 },
160  { 0.274960591212, -0.637341166847, -0.719856173359, 0.00690477957966 },
161  { 0.947273580412, 0.0, -0.320425910085, 0.00690477957966 },
162  { 0.812864676392, -0.304743149777, -0.496369449643, 0.00690477957966 },
163  { 0.595386501297, -0.188341624401, -0.781052076747, 0.00690477957966 },
164  { 0.595386501297, 0.188341624401, -0.781052076747, 0.00690477957966 },
165  { 0.812864676392, 0.304743149777, -0.496369449643, 0.00690477957966 },
166  { 0.492438766306, 0.753742692223, -0.435173546254, 0.00690477957966 },
167  { 0.274960591212, 0.637341166847, -0.719856173359, 0.00690477957966 },
168  { -0.076926487903, 0.753742692223, -0.652651721349, 0.00690477957966 },
169  { -0.076926487903, 0.942084316623, -0.326434458707, 0.00690477957966 },
170  { 0.274960591212, 0.942084316623, -0.192025554687, 0.00690477957966 },
171  { 0.947273580412, 0.277496978165, 0.160212955043, 0.00690477957966 },
172  { 0.812864676392, 0.582240127941, -0.015730584514, 0.00690477957966 },
173  { 0.595386501297, 0.770581752342, 0.227417407053, 0.00690477957966 },
174  { 0.595386501297, 0.582240127941, 0.553634669695, 0.00690477957966 },
175  { 0.812864676392, 0.277496978165, 0.512100034157, 0.00690477957966 },
176  { 0.492438766306, 0.0, 0.870347092509, 0.00690477957966 },
177  { 0.274960591212, 0.304743149777, 0.911881728046, 0.00690477957966 },
178  { -0.076926487903, 0.188341624401, 0.979086180056, 0.00690477957966 },
179  { -0.076926487903, -0.188341624401, 0.979086180056, 0.00690477957966 },
180  { 0.274960591212, -0.304743149777, 0.911881728046, 0.00690477957966 }
181  };
182  for ( i = 0; i < numberOfMicroplanes; i++ ) {
183  gauss_pts(i,0 ) = help [ i ] [ 0 ];
184  gauss_pts(i,1 ) = help [ i ] [ 1 ];
185  gauss_pts(i,2 ) = help [ i ] [ 2 ];
186  weights[i] = help [ i ] [ 3 ];
187  }
188  PetscFunctionReturn(0);
189 }
static char help[]