25#include <boost/program_options.hpp>
27namespace po = boost::program_options;
35#error "MoFEM need to be compiled with ADOL-C"
38#define TENSOR1_VEC_PTR(VEC) &VEC[0], &VEC[1], &VEC[2]
40#define SYMMETRIC_TENSOR4_MAT_PTR(MAT) \
41 &MAT(0, 0), &MAT(0, 1), &MAT(0, 2), &MAT(0, 3), &MAT(0, 4), &MAT(0, 5), \
42 &MAT(1, 0), &MAT(1, 1), &MAT(1, 2), &MAT(1, 3), &MAT(1, 4), &MAT(1, 5), \
43 &MAT(2, 0), &MAT(2, 1), &MAT(2, 2), &MAT(2, 3), &MAT(2, 4), &MAT(2, 5), \
44 &MAT(3, 0), &MAT(3, 1), &MAT(3, 2), &MAT(3, 3), &MAT(3, 4), &MAT(3, 5), \
45 &MAT(4, 0), &MAT(4, 1), &MAT(4, 2), &MAT(4, 3), &MAT(4, 4), &MAT(4, 5), \
46 &MAT(5, 0), &MAT(5, 1), &MAT(5, 2), &MAT(5, 3), &MAT(5, 4), &MAT(5, 5)
48#define TENSOR4_MAT_PTR(MAT) &MAT(0, 0), MAT.size2()
50#define TENSOR2_MAT_PTR(MAT) \
51 &MAT(0, 0), &MAT(1, 0), &MAT(2, 0), &MAT(3, 0), &MAT(4, 0), &MAT(5, 0), \
52 &MAT(6, 0), &MAT(7, 0), &MAT(8, 0)
54#define SYMMETRIC_TENSOR2_MAT_PTR(MAT) \
55 &MAT(0, 0), &MAT(0, 1), &MAT(0, 2), &MAT(0, 3), &MAT(0, 4), &MAT(0, 5)
57#define SYMMETRIC_TENSOR2_VEC_PTR(VEC) \
58 &VEC[0], &VEC[1], &VEC[2], &VEC[3], &VEC[4], &VEC[5]
64 PetscBool flg_config = PETSC_FALSE;
67 char my_config_file_name[255];
87 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Bone remodeling parameters",
91 CHKERR PetscOptionsString(
"-my_config",
"configuration file name",
"",
92 "my_config.in", my_config_file_name, 255,
96 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Config file: %s loaded.\n",
99 ifstream ini_file(my_config_file_name);
100 if (!ini_file.is_open()) {
101 SETERRQ(PETSC_COMM_SELF, 1,
"*** -my_config does not exist ***");
103 po::variables_map vm;
104 po::options_description config_file_options;
107 config_file_options.add_options()(
108 "young_modulus", po::value<double>(&
young_modulus)->default_value(1));
109 config_file_options.add_options()(
110 "poisson_ratio", po::value<double>(&
poisson_ratio)->default_value(1));
111 config_file_options.add_options()(
112 "c", po::value<double>(&
c)->default_value(1));
113 config_file_options.add_options()(
114 "m", po::value<double>(&
m)->default_value(1));
115 config_file_options.add_options()(
116 "n", po::value<double>(&
n)->default_value(1));
117 config_file_options.add_options()(
118 "rHo_ref", po::value<double>(&
rHo_ref)->default_value(1));
119 config_file_options.add_options()(
120 "pSi_ref", po::value<double>(&
pSi_ref)->default_value(1));
121 config_file_options.add_options()(
122 "R0", po::value<double>(&
R0)->default_value(1));
124 po::parsed_options parsed =
125 parse_config_file(ini_file, config_file_options,
true);
129 }
catch (
const std::exception &ex) {
130 std::ostringstream ss;
131 ss << ex.what() << std::endl;
136 CHKERR PetscOptionsScalar(
"-young_modulus",
"get young modulus",
"",
139 CHKERR PetscOptionsScalar(
"-poisson_ratio",
"get poisson_ratio",
"",
142 CHKERR PetscOptionsScalar(
"-c",
"density evolution (growth) velocity [d/m^2]",
143 "",
c, &
c, PETSC_NULL);
145 CHKERR PetscOptionsScalar(
"-m",
"algorithmic exponent",
"",
m, &
m,
148 CHKERR PetscOptionsScalar(
"-n",
"porosity exponent",
"",
n, &
n, PETSC_NULL);
150 CHKERR PetscOptionsInt(
"-b",
"bell function exponent",
"",
b, &
b, PETSC_NULL);
152 CHKERR PetscOptionsScalar(
"-rho_ref",
"reference bone density",
"",
rHo_ref,
155 CHKERR PetscOptionsScalar(
"-rho_max",
"reference bone density",
"",
rHo_max,
157 CHKERR PetscOptionsScalar(
"-rho_min",
"reference bone density",
"",
rHo_min,
160 CHKERR PetscOptionsScalar(
"-psi_ref",
"reference energy density",
"",
pSi_ref,
163 CHKERR PetscOptionsScalar(
"-r0",
"mass source",
"",
R0, &
R0, PETSC_NULL);
165 CHKERR PetscOptionsBool(
"-my_is_atom_test",
166 "is used with testing, exit with error when diverged",
169 CHKERR PetscOptionsBool(
"-less_post_proc",
170 "is used to reduce output file size",
"",
175 "calculate the material tangent with automatic differentiation",
"",
181 CHKERR PetscPrintf(PETSC_COMM_WORLD,
182 "Young's modulus E[Pa]: %4.3g\n",
184 CHKERR PetscPrintf(PETSC_COMM_WORLD,
185 "Poisson ratio nu[-] %4.3g\n",
187 CHKERR PetscPrintf(PETSC_COMM_WORLD,
188 "Lame coefficient lambda[Pa]: %4.3g\n",
lambda);
189 CHKERR PetscPrintf(PETSC_COMM_WORLD,
190 "Lame coefficient mu[Pa]: %4.3g\n",
mu);
191 CHKERR PetscPrintf(PETSC_COMM_WORLD,
192 "Density growth velocity [d/m2] %4.3g\n",
c);
193 CHKERR PetscPrintf(PETSC_COMM_WORLD,
194 "Algorithmic exponent m[-]: %4.3g\n",
m);
195 CHKERR PetscPrintf(PETSC_COMM_WORLD,
196 "Porosity exponent n[-]: %4.3g\n",
n);
197 CHKERR PetscPrintf(PETSC_COMM_WORLD,
198 "Reference density rHo_ref[kg/m3]: %4.3g\n",
rHo_ref);
199 CHKERR PetscPrintf(PETSC_COMM_WORLD,
200 "Reference energy pSi_ref[Pa]: %4.3g\n",
pSi_ref);
201 CHKERR PetscPrintf(PETSC_COMM_WORLD,
202 "Mass conduction R0[kg/m3s]: %4.3g\n",
R0);
203 CHKERR PetscPrintf(PETSC_COMM_WORLD,
204 "Bell function coefficent b[-]: %4.3g\n",
b);
206 CHKERR PetscPrintf(PETSC_COMM_WORLD,
207 "Max density rHo_max[kg/m3]: %4.3g\n",
rHo_max);
208 CHKERR PetscPrintf(PETSC_COMM_WORLD,
209 "Min density rHo_min[kg/m3]: %4.3g\n",
rHo_min);
212 ierr = PetscOptionsEnd();
219 double mid_rho = 0.5 * (rHo_max + rHo_min);
220 double var_h = (
rho - mid_rho) / (rHo_max - mid_rho);
222 return 1 / (1 + (b != 0) * pow(var_h, 2 * b));
226 double mid_rho = 0.5 * (rHo_max + rHo_min);
227 double var_h = (
rho - mid_rho) / (rHo_max - mid_rho);
228 return (b != 0) * (-2) * b * pow(var_h, 2 * b - 1) /
229 ((rHo_max - mid_rho) * pow(pow(var_h, 2 * b) + 1, 2));
240 DataForcesAndSourcesCore::EntData &data) {
244 const int nb_gauss_pts = data.getN().size1();
247 const int nb_base_functions = data.getN().size2();
249 auto base_function = data.getFTensor0N();
251 if (type == MBVERTEX) {
255 int nb_dofs = data.getIndices().size();
262 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
264 for (; bb < nb_dofs; bb++) {
265 const double field_data = array[data.getLocalIndices()[bb]];
266 drho_dt += field_data * base_function;
270 for (; bb != nb_base_functions; bb++)
283 commonData(common_data) {
295 DataForcesAndSourcesCore::EntData &data) {
299 if (type != MBVERTEX)
304 "Gradient at integration points of displacement is not calculated");
311 "Density at integration points is not calculated");
326 matS.resize(nb_gauss_pts, 6,
false);
331 matF.resize(9, nb_gauss_pts,
false);
345 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
348 F(
i,
I) = grad(
i,
I);
356 const double log_det = log(det_f);
357 psi =
F(
i,
J) *
F(
i,
J) - 3 - 2 * log_det;
359 psi += 0.5 *
lambda * log_det * log_det;
361 const double coef =
lambda * log_det -
mu;
362 P(
i,
J) =
mu *
F(
i,
J) + coef * invF(
J,
i);
363 S(
i,
I) = invF(
i,
J) ^ P(
J,
I);
366 R =
c * kp * (pow(
rho / rho_ref,
n -
m) * psi - psi_ref);
386 commonData(common_data) {
397 int side,
EntityType type, DataForcesAndSourcesCore::EntData &data) {
401 if (type != MBVERTEX)
407 vecC.resize(6,
false);
410 dS_dC.resize(6, 6 * nb_gauss_pts,
false);
415 matS.resize(nb_gauss_pts, 6,
false);
418 dP_dF.resize(81, nb_gauss_pts,
false);
420 matF.resize(9, nb_gauss_pts,
false);
435 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
438 F(
i,
I) = grad(
i,
I);
449 CHKERR recordFreeEnergy_dC<FTensor::Tensor0<FTensor::PackPtr<double *, 1>>,
458 "ADOL-C function evaluation with error");
469 const int is = 6 * gg;
470 double *hessian_ptr[6] = {&dS_dC(0, is), &dS_dC(1, is), &dS_dC(2, is),
471 &dS_dC(3, is), &dS_dC(4, is), &dS_dC(5, is)};
477 "ADOL-C function evaluation with error");
494 for (
int ii = 0; ii < 6; ii++) {
495 for (
int jj = 0; jj < ii; jj++) {
496 dS_dC(jj, is + ii) = dS_dC(ii, is + jj);
551 std::vector<EntityHandle> &map_gauss_pts,
555 commonData(common_data), postProcMesh(post_proc_mesh),
556 mapGaussPts(map_gauss_pts) {
568 DataForcesAndSourcesCore::EntData &data) {
572 if (type != MBVERTEX)
575 bzero(def_VAL, 9 *
sizeof(
double));
579 th_piola2, MB_TAG_CREAT | MB_TAG_SPARSE,
584 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
588 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
592 Tag th_prin_stress_vect1, th_prin_stress_vect2, th_prin_stress_vect3;
598 th_piola1, MB_TAG_CREAT | MB_TAG_SPARSE,
602 principal, MB_TAG_CREAT | MB_TAG_SPARSE,
606 th_prin_stress_vect1,
607 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
609 th_prin_stress_vect2,
610 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
612 th_prin_stress_vect3,
613 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
616 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
620 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
637 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
646 for (
int ii = 0; ii != 3; ii++) {
647 for (
int jj = 0; jj != 3; jj++) {
648 stress(ii, jj) = S(ii, jj);
656 for (
int ii = 0; ii != 3; ii++) {
657 for (
int jj = 0; jj != 3; jj++) {
658 stress(ii, jj) = P(ii, jj);
665 const double t2[] = {grad_rho(0), grad_rho(1), grad_rho(2)};
669 noalias(eigen_vectors) = stress;
675 int n = 3, lda = 3, info, lwork = -1;
677 info =
lapack_dsyev(
'V',
'U',
n, &(eigen_vectors.data()[0]), lda,
678 &(eigen_values.data()[0]), &wkopt, lwork);
680 SETERRQ1(PETSC_COMM_SELF, 1,
681 "something is wrong with lapack_dsyev info = %d", info);
684 info =
lapack_dsyev(
'V',
'U',
n, &(eigen_vectors.data()[0]), lda,
685 &(eigen_values.data()[0]), work, lwork);
687 SETERRQ1(PETSC_COMM_SELF, 1,
688 "something is wrong with lapack_dsyev info = %d", info);
690 for (
int ii = 0; ii < 3; ii++) {
691 prin_stress_vect1[ii] = eigen_vectors(0, ii);
692 prin_stress_vect2[ii] = eigen_vectors(1, ii);
693 prin_stress_vect3[ii] = eigen_vectors(2, ii);
699 1, &*prin_stress_vect1.begin());
701 1, &*prin_stress_vect2.begin());
703 1, &*prin_stress_vect3.begin());
721 commonData(common_data) {
733 DataForcesAndSourcesCore::EntData &data) {
737 if (type != MBVERTEX)
749 dP_dF.resize(81, nb_gauss_pts,
false);
751 matF.resize(9, nb_gauss_pts,
false);
767 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
770 F(
i,
I) = grad(
i,
I);
778 const double log_det = log(det_f);
779 psi =
F(
i,
J) *
F(
i,
J) - 3 - 2 * log_det;
781 psi += 0.5 *
lambda * log_det * log_det;
783 const double var =
lambda * log_det -
mu;
784 P(
i,
J) =
mu *
F(
i,
J) + var * invF(
J,
i);
786 const double coef =
mu -
lambda * log_det;
790 invF(
J,
i) * (invF(
L,
k) *
lambda) + invF(
L,
i) * (invF(
J,
k) * coef);
792 D2(0, 0, 0, 0) +=
mu;
793 D2(0, 1, 0, 1) +=
mu;
794 D2(0, 2, 0, 2) +=
mu;
795 D2(1, 0, 1, 0) +=
mu;
796 D2(1, 1, 1, 1) +=
mu;
797 D2(1, 2, 1, 2) +=
mu;
798 D2(2, 0, 2, 0) +=
mu;
799 D2(2, 1, 2, 1) +=
mu;
800 D2(2, 2, 2, 2) +=
mu;
817 commonData(common_data) {}
821 DataForcesAndSourcesCore::EntData &data) {
825 const int nb_dofs = data.getIndices().size();
828 nF.resize(nb_dofs,
false);
831 const int nb_gauss_pts = data.getN().size1();
832 const int nb_base_functions = data.getN().size2();
833 auto diff_base_functions = data.getFTensor1DiffN<3>();
849 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
853 double a = w * pow(
rho / rho_ref,
n);
856 for (; bb != nb_dofs / 3; bb++) {
857 double *ptr = &
nF[3 * bb];
859 t1(
i) +=
a * P(
i,
I) * diff_base_functions(
I);
862 ++diff_base_functions;
866 for (; bb != nb_base_functions; bb++) {
867 ++diff_base_functions;
878 nb_dofs, &*data.getIndices().begin(), &*
nF.begin(), ADD_VALUES);
886 commonData(common_data) {}
890 DataForcesAndSourcesCore::EntData &data) {
894 const int nb_dofs = data.getIndices().size();
897 nF.resize(nb_dofs,
false);
900 const int nb_gauss_pts = data.getN().size1();
901 const int nb_base_functions = data.getN().size2();
902 auto base_functions = data.getFTensor0N();
903 auto diff_base_functions = data.getFTensor1DiffN<3>();
926 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
932 for (; bb != nb_dofs; bb++) {
933 nF[bb] += w * base_functions * drho_dt;
934 nF[bb] -= w * base_functions *
R;
935 nF[bb] -= w * R0 * diff_base_functions(
I) * grad_rho(
I);
937 ++diff_base_functions;
941 for (; bb != nb_base_functions; bb++) {
943 ++diff_base_functions;
956 nb_dofs, &*data.getIndices().begin(), &*
nF.begin(), ADD_VALUES);
964 commonData(common_data) {
971 DataForcesAndSourcesCore::EntData &row_data,
972 DataForcesAndSourcesCore::EntData &col_data) {
976 const int row_nb_dofs = row_data.getIndices().size();
979 const int col_nb_dofs = col_data.getIndices().size();
987 const int row_nb_gauss_pts = row_data.getN().size1();
988 if (!row_nb_gauss_pts)
990 const int row_nb_base_functions = row_data.getN().size2();
1004 auto row_base_functions = row_data.getFTensor0N();
1005 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1007 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
1014 const double a =
c * kp * ((
n -
m) /
rho) * pow(
rho / rho_ref,
n -
m) * psi;
1015 const double a_diff =
1016 c * kp_diff * (pow(
rho / rho_ref,
n -
m) * psi - psi_ref);
1019 for (; row_bb != row_nb_dofs; row_bb++) {
1020 auto col_base_functions = col_data.getFTensor0N(gg, 0);
1021 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1022 for (
int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1024 ts_a * w * row_base_functions * col_base_functions;
1026 w * R0 * row_diff_base_functions(
I) * col_diff_base_functions(
I);
1028 a * w * row_base_functions * col_base_functions;
1030 a_diff * w * row_base_functions * col_base_functions;
1032 ++col_base_functions;
1033 ++col_diff_base_functions;
1035 ++row_base_functions;
1036 ++row_diff_base_functions;
1038 for (; row_bb != row_nb_base_functions; row_bb++) {
1039 ++row_base_functions;
1040 ++row_diff_base_functions;
1049 &*row_data.getIndices().begin(), col_nb_dofs,
1054 if (row_side != col_side || row_type != col_type) {
1058 &*col_data.getIndices().begin(), row_nb_dofs,
1059 &*row_data.getIndices().begin(),
1069 commonData(common_data) {
1076 DataForcesAndSourcesCore::EntData &row_data,
1077 DataForcesAndSourcesCore::EntData &col_data) {
1081 const int row_nb_dofs = row_data.getIndices().size();
1084 const int col_nb_dofs = col_data.getIndices().size();
1089 locK_rho_F.resize(row_nb_dofs, col_nb_dofs,
false);
1092 const int row_nb_gauss_pts = row_data.getN().size1();
1093 const int col_nb_base_functions = col_data.getN().size2();
1107 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>();
1109 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
1114 const double a = w *
c * kp * pow(
rho / rho_ref,
n -
m);
1117 for (; col_bb != col_nb_dofs / 3; col_bb++) {
1118 t1(
i) =
a * P(
i,
I) * col_diff_base_functions(
I);
1119 auto row_base_functions = row_data.getFTensor0N(gg, 0);
1120 for (
int row_bb = 0; row_bb != row_nb_dofs; row_bb++) {
1124 k(
i) -= row_base_functions * t1(
i);
1125 ++row_base_functions;
1127 ++col_diff_base_functions;
1129 for (; col_bb != col_nb_base_functions; col_bb++) {
1130 ++col_diff_base_functions;
1138 &*row_data.getIndices().begin(), col_nb_dofs,
1139 &*col_data.getIndices().begin(), &
locK_rho_F(0, 0),
1144template <
bool ADOLC>
1149 commonData(common_data) {
1152template <
bool ADOLC>
1155 DataForcesAndSourcesCore::EntData &row_data,
1156 DataForcesAndSourcesCore::EntData &col_data) {
1160 const int row_nb_dofs = row_data.getIndices().size();
1163 const int col_nb_dofs = col_data.getIndices().size();
1167 const bool diagonal_block = (row_type == col_type) && (row_side == col_side);
1170 locK_P_F.resize(row_nb_dofs, col_nb_dofs,
false);
1173 const int row_nb_gauss_pts = row_data.getN().size1();
1174 const int row_nb_base_functions = row_data.getN().size2();
1177 MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
1183 const double n = commonData.n;
1184 const double rho_ref = commonData.rHo_ref;
1194 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1202 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
1205 double w = getVolume() * getGaussPts()(3, gg);
1206 const double a = w * pow(
rho / rho_ref,
n);
1209 for (; row_bb != row_nb_dofs / 3; row_bb++) {
1211 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1213 const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
1215 for (; col_bb != final_bb; col_bb++) {
1218 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1219 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1220 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1221 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1222 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1223 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1224 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1225 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1226 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1229 a * row_diff_base_functions(
J) * col_diff_base_functions(
L);
1234 t_assemble(
i,
k) += diffDiff(
J,
L) * D2(
i,
J,
k,
L);
1239 t0 = diffDiff(
J,
L) * S(
J,
L);
1240 t_assemble(0, 0) += t0;
1241 t_assemble(1, 1) += t0;
1242 t_assemble(2, 2) += t0;
1246 ++col_diff_base_functions;
1249 ++row_diff_base_functions;
1251 for (; row_bb != row_nb_base_functions; row_bb++) {
1252 ++row_diff_base_functions;
1263 if (diagonal_block) {
1264 for (
int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
1266 for (; col_bb != row_bb + 1; col_bb++) {
1268 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1269 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1270 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1271 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1272 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1273 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1274 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1275 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1276 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1278 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 0),
1279 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 1),
1280 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 2),
1281 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 0),
1282 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 1),
1283 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 2),
1284 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 0),
1285 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 1),
1286 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 2));
1287 t_off_side(
i,
k) = t_assemble(
k,
i);
1292 const int *row_ind = &*row_data.getIndices().begin();
1293 const int *col_ind = &*col_data.getIndices().begin();
1295 col_ind, &locK_P_F(0, 0), ADD_VALUES);
1297 if (row_type != col_type || row_side != col_side) {
1298 transLocK_P_F.resize(col_nb_dofs, row_nb_dofs,
false);
1299 noalias(transLocK_P_F) = trans(locK_P_F);
1301 row_ind, &transLocK_P_F(0, 0), ADD_VALUES);
1311 commonData(common_data) {
1318 DataForcesAndSourcesCore::EntData &row_data,
1319 DataForcesAndSourcesCore::EntData &col_data) {
1323 const int row_nb_dofs = row_data.getIndices().size();
1326 const int col_nb_dofs = col_data.getIndices().size();
1331 locK_P_RHO.resize(row_nb_dofs, col_nb_dofs,
false);
1345 const int row_nb_gauss_pts = row_data.getN().size1();
1346 const int row_nb_base_functions = row_data.getN().size2();
1347 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1349 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
1353 double a = w * (
n /
rho) * pow(
rho / rho_ref,
n);
1356 for (; row_bb != row_nb_dofs / 3; row_bb++) {
1357 t1(
i) =
a * row_diff_base_functions(
I) * P(
i,
I);
1358 auto base_functions = col_data.getFTensor0N(gg, 0);
1359 for (
int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1366 k(
i) += t1(
i) * base_functions;
1369 ++row_diff_base_functions;
1371 for (; row_bb != row_nb_base_functions; row_bb++) {
1372 ++row_diff_base_functions;
1380 &*row_data.getIndices().begin(), col_nb_dofs,
1381 &*col_data.getIndices().begin(), &
locK_P_RHO(0, 0),
1389 : mField(m_field), commonData(common_data), postProc(m_field),
1390 postProcElastic(m_field),
1409 PetscBool analysis_mesh_flg = PETSC_FALSE;
1412 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Bone remodeling post-process",
1417 "is used for testing, calculates mass and energy at each time step",
"",
1420 CHKERR PetscOptionsBool(
"-analysis_mesh",
1421 "saves analysis mesh at each time step",
"",
1422 analysis_mesh_flg, &analysis_mesh_flg, PETSC_NULL);
1424 CHKERR PetscOptionsScalar(
1425 "-equilibrium_stop_rate",
1426 "is used to stop calculations when equilibium state is achieved",
"",
1428 ierr = PetscOptionsEnd();
1433 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
1434 "Bone remodeling post-process",
"none");
1437 CHKERR PetscOptionsInt(
"-my_output_prt",
1438 "frequncy how often results are dumped on hard disk",
1440 ierr = PetscOptionsEnd();
1453 auto &list_of_operators =
postProc.getOpPtrVector();
1455 list_of_operators.push_back(
1471 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
1485 int mass_vec_ghost[] = {0};
1487 1, 1, mass_vec_ghost, &mass_vec);
1489 CHKERR VecZeroEntries(mass_vec);
1490 CHKERR VecDuplicate(mass_vec, &energ_vec);
1508 CHKERR VecAssemblyBegin(mass_vec);
1509 CHKERR VecAssemblyEnd(mass_vec);
1511 CHKERR VecSum(mass_vec, &mass_sum);
1513 CHKERR VecAssemblyBegin(energ_vec);
1514 CHKERR VecAssemblyEnd(energ_vec);
1516 CHKERR VecSum(energ_vec, &energ_sum);
1518 CHKERR PetscPrintf(PETSC_COMM_WORLD,
1519 "Time: %g Mass: %6.5g Elastic energy: %6.5g \n", ts_t,
1520 mass_sum, energ_sum);
1521 CHKERR VecDestroy(&mass_vec);
1522 CHKERR VecDestroy(&energ_vec);
1525 double equilibrium_rate =
1527 double equilibrium_mass_rate =
1529 if (equilibrium_rate <
rate) {
1532 "Energy equilibrium state is achieved! Difference = %0.6g %%. \n",
1533 equilibrium_rate * 100);
1535 if (equilibrium_mass_rate <
rate) {
1538 "Mass equilibrium state is achieved! Difference = %0.6g %%. \n",
1539 equilibrium_mass_rate * 100);
1547#if PETSC_VERSION_LE(3, 8, 0)
1548 CHKERR TSGetTimeStepNumber(ts, &step);
1550 CHKERR TSGetStepNumber(ts, &step);
1553 if ((step) %
pRT == 0) {
1556 sss <<
"out_" << step <<
".h5m";
1559 "PARALLEL=WRITE_PART");
1560 if (analysis_mesh_flg) {
1562 ttt <<
"analysis_mesh_" << step <<
".h5m";
1564 "PARALLEL=WRITE_PART");
1571 ss <<
"out_elastic_" << step <<
".h5m";
1573 "PARALLEL=WRITE_PART");
1582 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Bone remodeling",
"none");
1585 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"", 2,
1588 ierr = PetscOptionsEnd();
1593 string name = it->getName();
1594 if (name.compare(0, 14,
"NO_REMODELLING") == 0) {
1628 bool add_rho_field =
false;
1636 add_rho_field =
true;
1676 if (add_rho_field) {
1688 "MESH_NODE_POSITIONS");
1706 "MESH_NODE_POSITIONS");
1718 "MESH_NODE_POSITIONS");
1733 "DISPLACEMENTS",
"MESH_NODE_POSITIONS",
false,
true);
1751 list_of_operators_rhs.push_back(
1753 list_of_operators_rhs.push_back(
1768 list_of_operators_lhs.push_back(
1770 list_of_operators_rhs.push_back(
1776 list_of_operators_lhs.push_back(
1778 list_of_operators_lhs.push_back(
1782 list_of_operators_lhs.push_back(
1796 CHKERR MetaNeummanForces::addNeumannBCElements(
mField,
"DISPLACEMENTS");
1801 CHKERR MetaNeummanForces::setMomentumFluxOperators(
1805 PETSC_NULL,
"DISPLACEMENTS");
1810 for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
1813 mit->second->methodsOp.push_back(
1816 for (boost::ptr_map<string, NodalForce>::iterator mit =
1819 mit->second->methodsOp.push_back(
1822 for (boost::ptr_map<string, EdgeForce>::iterator mit =
1825 mit->second->methodsOp.push_back(
1873 CHKERR TSSetIFunction(ts,
F, PETSC_NULL, PETSC_NULL);
1874 CHKERR TSSetIJacobian(ts,
A,
A, PETSC_NULL, PETSC_NULL);
1876#if PETSC_VERSION_GE(3, 8, 0)
1877 CHKERR TSSetMaxTime(ts, ftime);
1879 CHKERR TSSetFromOptions(ts);
1883 CHKERR TSGetSNES(ts, &snes);
1885 CHKERR SNESGetKSP(snes, &ksp);
1887 CHKERR KSPGetPC(ksp, &pc);
1888 PetscBool is_pcfs = PETSC_FALSE;
1889 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1893 if (is_pcfs == PETSC_TRUE) {
1898 problem_ptr->
getName(),
ROW,
"DISPLACEMENTS", 0, 3, &is_disp);
1900 problem_ptr->
getName(),
ROW,
"RHO", 0, 1, &is_rho);
1903 CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
1904 CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
1905 CHKERR ISDestroy(&is_disp);
1906 CHKERR ISDestroy(&is_rho);
1918#if PETSC_VERSION_GE(3, 7, 0)
1919 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1922 CHKERR TSGetTime(ts, &ftime);
1923 PetscInt steps, snesfails, rejects, nonlinits, linits;
1924#if PETSC_VERSION_LE(3, 8, 0)
1925 CHKERR TSGetTimeStepNumber(ts, &steps);
1927 CHKERR TSGetStepNumber(ts, &steps);
1930 CHKERR TSGetSNESFailures(ts, &snesfails);
1931 CHKERR TSGetStepRejections(ts, &rejects);
1932 CHKERR TSGetSNESIterations(ts, &nonlinits);
1933 CHKERR TSGetKSPIterations(ts, &linits);
1934 PetscPrintf(PETSC_COMM_WORLD,
1935 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
1937 steps, rejects, snesfails, ftime, nonlinits, linits);
1949 DataForcesAndSourcesCore::EntData &row_data) {
1952 if (row_type != MBVERTEX)
1955 const int nb_gauss_pts = row_data.getN().size1();
1962 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1965 energy += vol * psi;
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT)
#define TENSOR2_MAT_PTR(MAT)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT)
@ MOFEM_STD_EXCEPTION_THROW
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define TENSOR4_MAT_PTR(MAT)
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
constexpr IntegrationType I
double young_modulus
Young modulus.
Remodeling::CommonData & commonData
MoFEMErrorCode postProcess()
MoFEMErrorCode operator()()
MoFEMErrorCode preProcess()
MoFEM::Interface & mField
PetscBool equilibrium_flg
PostProcVolumeOnRefinedMesh postProcElastic
PostProcVolumeOnRefinedMesh postProc
MonitorPostProc(MoFEM::Interface &m_field, Remodeling::CommonData &common_data)
Off diagonal block of tangent matrix .
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data)
Remodeling::CommonData & commonData
Diagonal block of tangent matrix .
MatrixDouble transLocK_rho_rho
OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
Remodeling::CommonData & commonData
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
MatrixDouble locK_rho_rho
Assemble residual for conservation of mass (density)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
VectorDouble nF
Vector of the right hand side (internal forces)
Remodeling::CommonData & commonData
Off diagonal block of tangent matrix .
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssmbleStressLhs_dF(Remodeling::CommonData &common_data)
Off diagonal block of tangent matrix /f[ K_{u \rho}=\intop_{V} \left[\frac{n}{\rho_{0}}\right] \left...
Remodeling::CommonData & commonData
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssmbleStressLhs_dRho(Remodeling::CommonData &common_data)
Assemble residual for conservation of momentum (stresses)
Remodeling::CommonData & commonData
VectorDouble nF
Vector of the right hand side (internal forces)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpAssmbleStressRhs(Remodeling::CommonData &common_data)
Evaluate physical equations at integration points.
Remodeling::CommonData & commonData
OpCalculateStress(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Remodeling::CommonData & commonData
OpCalculateStressTangent(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpCalculateStressTangentWithAdolc(Remodeling::CommonData &common_data)
Remodeling::CommonData & commonData
Evaluate density derivative with respect to time in case of Backward Euler Method.
OpGetRhoTimeDirevative(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Remodeling::CommonData & commonData
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Remodeling::CommonData & commonData
Used to post proc stresses, energy, source term.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, Remodeling::CommonData &common_data)
Remodeling::CommonData & commonData
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
VectorDouble vecDetF
determinant of F
VectorDouble vecPsi
Elastic energy density.
VectorDouble vecR
Mass sorce.
boost::shared_ptr< MatrixDouble > gradDispPtr
Ptr to gradient of displacements matrix container.
boost::shared_ptr< VectorDouble > rhoPtr
Ptr to density matrix container.
MatrixDouble matPushedMaterialTangent
MatrixDouble matP
1st Piola stress
MatrixDouble matS
2nd Piola stress
boost::shared_ptr< MatrixDouble > gradRhoPtr
Gradient of density field.
MatrixDouble matMaterialTangent
MatrixDouble matInvF
inverse of deformation gradient
VectorDouble vecRhoDt
Time derivative of density.
MatrixDouble matGradientOfDeformation
Gradient of deformation.
double getCFromDensity(const double &rho)
double pSi_ref
reference free energy
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
double n
porosity exponent [-]
int b
b exponent for bell function
double getCFromDensityDiff(const double &rho)
double m
algorithmic exponent [-]
DMType dm_name
dm (problem) name
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
PetscBool is_atom_testing
for atom tests
double c
density evolution (growth) velocity [d/m^2]
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
MoFEMErrorCode getParameters()
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
double rHo_max
max density
PetscBool less_post_proc
reduce file size
DM dm
Discretization manager.
boost::shared_ptr< ElasticMaterials > elasticMaterialsPtr
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
double rHo_min
min density
double rHo_ref
reference density
double lambda
Lame parameter.
double R0
mass conduction coefficient
double cUrrent_psi
current free energy for evaluating equilibrium state
boost::shared_ptr< NonlinearElasticElement > elasticPtr
double cUrrent_mass
current free energy for evaluating equilibrium state
MoFEM::Interface & mField
MoFEMErrorCode addElements()
Set and add finite elements.
MoFEMErrorCode solveDM()
Solve problem set up in DM.
MoFEMErrorCode buildDM()
Set problem and DM.
MoFEMErrorCode getParameters()
Get parameters form line command or config file.
MoFEMErrorCode addFields()
Set and add entities to approximation fields.
MoFEMErrorCode addMomentumFluxes()
Finite elements to calculate tractions.
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual int get_comm_rank() const =0
bool & doPrisms
\deprectaed
bool sYmm
If true assume that matrix is symmetric structure.
bool & doVertices
\deprectaed If false skip vertices
bool & doEdges
\deprectaed If false skip edges
bool & doQuads
\deprectaed
Deprecated interface functions.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Section manager is used to create indexes and sections.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
keeps basic data about problem
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
Interface for Time Stepping (TS) solver.
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
double getVolume() const
element volume (linear geometry)
Volume finite element base.
structure grouping operators and data used for calculation of nonlinear elastic element
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Force scale operator for reading two columns.