v0.14.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}\). More...
 
struct  OpAssmbleRhoLhs_dRho
 Diagonal block of tangent matrix \(K_{\rho\rho}\). More...
 
struct  OpAssmbleRhoRhs
 Assemble residual for conservation of mass (density) More...
 
struct  OpAssmbleRhs
 Assemble RHS vector f. More...
 
struct  OpAssmbleStressLhs_dF
 Off diagonal block of tangent matrix \(K_{u u}\). 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  OpCalculateLhs
 Assemble LHS matrix K. More...
 
struct  OpCalculateStress
 Evaluate physical equations at integration points. More...
 
struct  OpCalculateStressTangent
 
struct  OpCalculateStressTangentWithAdolc
 
struct  OpCalulatefRhoAtGaussPts
 Assemble local vector containing density data. More...
 
struct  OpGetRhoTimeDirevative
 Evaluate density derivative with respect to time in case of Backward Euler Method. More...
 
struct  OpMassAndEnergyCalculation
 
struct  OpMassCalculation
 Calculate mass before approximation. More...
 
struct  OpMassCalculationFromApprox
 Calculate mass after approximation. 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 
)
inline

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 
275  MoFEMFunctionBeginHot;
276  FTensor::Index<'I', 3> I;
277  FTensor::Index<'J', 3> J;
278  const double mu = common_data.mu;
279  const double lambda = common_data.lambda;
280 
281  // Compressible Neo-Hookean
282  B2 det;
283  CHKERR determinantTensor3by3(C, 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 
301  MoFEMFunctionReturnHot(0);
302 }

◆ recordFreeEnergy_dC()

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

Definition at line 305 of file Remodeling.hpp.

306  {
307 
308  MoFEMFunctionBeginHot;
309  FTensor::Index<'I', 3> I;
310  FTensor::Index<'J', 3> J;
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();
319  MoFEMFunctionReturnHot(0);
320 }

◆ 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 }

◆ 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 }

◆ 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 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
FTensor::Index
Definition: Index.hpp:23
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1511
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
Generated by Doxygen 1.8.17 and hosted at University of Glasgow