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 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);
218 double mid_rho = 0.5 * (rHo_max + rHo_min);
219 double var_h = (
rho - mid_rho) / (rHo_max - mid_rho);
221 return 1 / (1 + (b != 0) * pow(var_h, 2 * b));
225 double mid_rho = 0.5 * (rHo_max + rHo_min);
226 double var_h = (
rho - mid_rho) / (rHo_max - mid_rho);
227 return (b != 0) * (-2) * b * pow(var_h, 2 * b - 1) /
228 ((rHo_max - mid_rho) * pow(pow(var_h, 2 * b) + 1, 2));
239 DataForcesAndSourcesCore::EntData &data) {
243 const int nb_gauss_pts = data.getN().size1();
246 const int nb_base_functions = data.getN().size2();
248 auto base_function = data.getFTensor0N();
250 if (type == MBVERTEX) {
254 int nb_dofs = data.getIndices().size();
261 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
263 for (; bb < nb_dofs; bb++) {
264 const double field_data = array[data.getLocalIndices()[bb]];
265 drho_dt += field_data * base_function;
269 for (; bb != nb_base_functions; bb++)
282 commonData(common_data) {
294 DataForcesAndSourcesCore::EntData &data) {
298 if (type != MBVERTEX)
303 "Gradient at integration points of displacement is not calculated");
310 "Density at integration points is not calculated");
325 matS.resize(nb_gauss_pts, 6,
false);
330 matF.resize(9, nb_gauss_pts,
false);
344 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
347 F(
i,
I) = grad(
i,
I);
355 const double log_det = log(det_f);
356 psi =
F(
i,
J) *
F(
i,
J) - 3 - 2 * log_det;
358 psi += 0.5 *
lambda * log_det * log_det;
360 const double coef =
lambda * log_det -
mu;
361 P(
i,
J) =
mu *
F(
i,
J) + coef * invF(
J,
i);
362 S(
i,
I) = invF(
i,
J) ^ P(
J,
I);
365 R =
c * kp * (pow(
rho / rho_ref,
n -
m) * psi - psi_ref);
385 commonData(common_data) {
396 int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
400 if (type != MBVERTEX)
406 vecC.resize(6,
false);
409 dS_dC.resize(6, 6 * nb_gauss_pts,
false);
414 matS.resize(nb_gauss_pts, 6,
false);
417 dP_dF.resize(81, nb_gauss_pts,
false);
419 matF.resize(9, nb_gauss_pts,
false);
434 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
437 F(
i,
I) = grad(
i,
I);
457 "ADOL-C function evaluation with error");
468 const int is = 6 * gg;
469 double *hessian_ptr[6] = {&dS_dC(0, is), &dS_dC(1, is), &dS_dC(2, is),
470 &dS_dC(3, is), &dS_dC(4, is), &dS_dC(5, is)};
476 "ADOL-C function evaluation with error");
493 for (
int ii = 0; ii < 6; ii++) {
494 for (
int jj = 0; jj < ii; jj++) {
495 dS_dC(jj, is + ii) = dS_dC(ii, is + jj);
550 std::vector<EntityHandle> &map_gauss_pts,
554 commonData(common_data), postProcMesh(post_proc_mesh),
555 mapGaussPts(map_gauss_pts) {
567 DataForcesAndSourcesCore::EntData &data) {
571 if (type != MBVERTEX)
574 bzero(def_VAL, 9 *
sizeof(
double));
578 th_piola2, MB_TAG_CREAT | MB_TAG_SPARSE,
583 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
587 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
591 Tag th_prin_stress_vect1, th_prin_stress_vect2, th_prin_stress_vect3;
597 th_piola1, MB_TAG_CREAT | MB_TAG_SPARSE,
601 principal, MB_TAG_CREAT | MB_TAG_SPARSE,
605 th_prin_stress_vect1,
606 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
608 th_prin_stress_vect2,
609 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
611 th_prin_stress_vect3,
612 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
615 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
619 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
633 VectorDouble3 eigen_values(3);
636 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
645 for (
int ii = 0; ii != 3; ii++) {
646 for (
int jj = 0; jj != 3; jj++) {
647 stress(ii, jj) = S(ii, jj);
655 for (
int ii = 0; ii != 3; ii++) {
656 for (
int jj = 0; jj != 3; jj++) {
657 stress(ii, jj) = P(ii, jj);
664 const double t2[] = {grad_rho(0), grad_rho(1), grad_rho(2)};
668 noalias(eigen_vectors) = stress;
669 VectorDouble3 prin_stress_vect1(3);
670 VectorDouble3 prin_stress_vect2(3);
671 VectorDouble3 prin_stress_vect3(3);
674 int n = 3, lda = 3, info, lwork = -1;
676 info =
lapack_dsyev(
'V',
'U',
n, &(eigen_vectors.data()[0]), lda,
677 &(eigen_values.data()[0]), &wkopt, lwork);
679 SETERRQ1(PETSC_COMM_SELF, 1,
680 "something is wrong with lapack_dsyev info = %d", info);
683 info =
lapack_dsyev(
'V',
'U',
n, &(eigen_vectors.data()[0]), lda,
684 &(eigen_values.data()[0]), work, lwork);
686 SETERRQ1(PETSC_COMM_SELF, 1,
687 "something is wrong with lapack_dsyev info = %d", info);
689 for (
int ii = 0; ii < 3; ii++) {
690 prin_stress_vect1[ii] = eigen_vectors(0, ii);
691 prin_stress_vect2[ii] = eigen_vectors(1, ii);
692 prin_stress_vect3[ii] = eigen_vectors(2, ii);
698 1, &*prin_stress_vect1.begin());
700 1, &*prin_stress_vect2.begin());
702 1, &*prin_stress_vect3.begin());
720 commonData(common_data) {
732 DataForcesAndSourcesCore::EntData &data) {
736 if (type != MBVERTEX)
748 dP_dF.resize(81, nb_gauss_pts,
false);
750 matF.resize(9, nb_gauss_pts,
false);
766 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
769 F(
i,
I) = grad(
i,
I);
777 const double log_det = log(det_f);
778 psi =
F(
i,
J) *
F(
i,
J) - 3 - 2 * log_det;
780 psi += 0.5 *
lambda * log_det * log_det;
782 const double var =
lambda * log_det -
mu;
783 P(
i,
J) =
mu *
F(
i,
J) + var * invF(
J,
i);
785 const double coef =
mu -
lambda * log_det;
789 invF(
J,
i) * (invF(
L,
k) *
lambda) + invF(
L,
i) * (invF(
J,
k) * coef);
791 D2(0, 0, 0, 0) +=
mu;
792 D2(0, 1, 0, 1) +=
mu;
793 D2(0, 2, 0, 2) +=
mu;
794 D2(1, 0, 1, 0) +=
mu;
795 D2(1, 1, 1, 1) +=
mu;
796 D2(1, 2, 1, 2) +=
mu;
797 D2(2, 0, 2, 0) +=
mu;
798 D2(2, 1, 2, 1) +=
mu;
799 D2(2, 2, 2, 2) +=
mu;
816 commonData(common_data) {}
820 DataForcesAndSourcesCore::EntData &data) {
824 const int nb_dofs = data.getIndices().size();
827 nF.resize(nb_dofs,
false);
830 const int nb_gauss_pts = data.getN().size1();
831 const int nb_base_functions = data.getN().size2();
832 auto diff_base_functions = data.getFTensor1DiffN<3>();
848 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
852 double a = w * pow(
rho / rho_ref,
n);
855 for (; bb != nb_dofs / 3; bb++) {
856 double *ptr = &
nF[3 * bb];
858 t1(
i) +=
a * P(
i,
I) * diff_base_functions(
I);
861 ++diff_base_functions;
865 for (; bb != nb_base_functions; bb++) {
866 ++diff_base_functions;
877 nb_dofs, &*data.getIndices().begin(), &*
nF.begin(), ADD_VALUES);
885 commonData(common_data) {}
889 DataForcesAndSourcesCore::EntData &data) {
893 const int nb_dofs = data.getIndices().size();
896 nF.resize(nb_dofs,
false);
899 const int nb_gauss_pts = data.getN().size1();
900 const int nb_base_functions = data.getN().size2();
901 auto base_functions = data.getFTensor0N();
902 auto diff_base_functions = data.getFTensor1DiffN<3>();
925 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
931 for (; bb != nb_dofs; bb++) {
932 nF[bb] += w * base_functions * drho_dt;
933 nF[bb] -= w * base_functions *
R;
934 nF[bb] -= w * R0 * diff_base_functions(
I) * grad_rho(
I);
936 ++diff_base_functions;
940 for (; bb != nb_base_functions; bb++) {
942 ++diff_base_functions;
955 nb_dofs, &*data.getIndices().begin(), &*
nF.begin(), ADD_VALUES);
963 commonData(common_data) {
970 DataForcesAndSourcesCore::EntData &row_data,
971 DataForcesAndSourcesCore::EntData &col_data) {
975 const int row_nb_dofs = row_data.getIndices().size();
978 const int col_nb_dofs = col_data.getIndices().size();
986 const int row_nb_gauss_pts = row_data.getN().size1();
987 if (!row_nb_gauss_pts)
989 const int row_nb_base_functions = row_data.getN().size2();
1003 auto row_base_functions = row_data.getFTensor0N();
1004 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1006 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
1013 const double a =
c * kp * ((
n -
m) /
rho) * pow(
rho / rho_ref,
n -
m) * psi;
1014 const double a_diff =
1015 c * kp_diff * (pow(
rho / rho_ref,
n -
m) * psi - psi_ref);
1018 for (; row_bb != row_nb_dofs; row_bb++) {
1019 auto col_base_functions = col_data.getFTensor0N(gg, 0);
1020 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1021 for (
int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1023 ts_a * w * row_base_functions * col_base_functions;
1025 w * R0 * row_diff_base_functions(
I) * col_diff_base_functions(
I);
1027 a * w * row_base_functions * col_base_functions;
1029 a_diff * w * row_base_functions * col_base_functions;
1031 ++col_base_functions;
1032 ++col_diff_base_functions;
1034 ++row_base_functions;
1035 ++row_diff_base_functions;
1037 for (; row_bb != row_nb_base_functions; row_bb++) {
1038 ++row_base_functions;
1039 ++row_diff_base_functions;
1048 &*row_data.getIndices().begin(), col_nb_dofs,
1053 if (row_side != col_side || row_type != col_type) {
1057 &*col_data.getIndices().begin(), row_nb_dofs,
1058 &*row_data.getIndices().begin(),
1068 commonData(common_data) {
1074 EntityType col_type,
1075 DataForcesAndSourcesCore::EntData &row_data,
1076 DataForcesAndSourcesCore::EntData &col_data) {
1080 const int row_nb_dofs = row_data.getIndices().size();
1083 const int col_nb_dofs = col_data.getIndices().size();
1088 locK_rho_F.resize(row_nb_dofs, col_nb_dofs,
false);
1091 const int row_nb_gauss_pts = row_data.getN().size1();
1092 const int col_nb_base_functions = col_data.getN().size2();
1106 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>();
1108 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
1113 const double a = w *
c * kp * pow(
rho / rho_ref,
n -
m);
1116 for (; col_bb != col_nb_dofs / 3; col_bb++) {
1117 t1(
i) =
a * P(
i,
I) * col_diff_base_functions(
I);
1118 auto row_base_functions = row_data.getFTensor0N(gg, 0);
1119 for (
int row_bb = 0; row_bb != row_nb_dofs; row_bb++) {
1123 k(
i) -= row_base_functions * t1(
i);
1124 ++row_base_functions;
1126 ++col_diff_base_functions;
1128 for (; col_bb != col_nb_base_functions; col_bb++) {
1129 ++col_diff_base_functions;
1137 &*row_data.getIndices().begin(), col_nb_dofs,
1138 &*col_data.getIndices().begin(), &
locK_rho_F(0, 0),
1143template <
bool ADOLC>
1148 commonData(common_data) {
1151template <
bool ADOLC>
1153 int row_side,
int col_side, EntityType row_type, EntityType col_type,
1154 DataForcesAndSourcesCore::EntData &row_data,
1155 DataForcesAndSourcesCore::EntData &col_data) {
1159 const int row_nb_dofs = row_data.getIndices().size();
1162 const int col_nb_dofs = col_data.getIndices().size();
1166 const bool diagonal_block = (row_type == col_type) && (row_side == col_side);
1169 locK_P_F.resize(row_nb_dofs, col_nb_dofs,
false);
1172 const int row_nb_gauss_pts = row_data.getN().size1();
1173 const int row_nb_base_functions = row_data.getN().size2();
1175 MatrixDouble &matS = commonData.data.matS;
1176 MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
1182 const double n = commonData.n;
1183 const double rho_ref = commonData.rHo_ref;
1193 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1201 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
1204 double w = getVolume() * getGaussPts()(3, gg);
1205 const double a = w * pow(
rho / rho_ref,
n);
1208 for (; row_bb != row_nb_dofs / 3; row_bb++) {
1210 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1212 const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
1214 for (; col_bb != final_bb; col_bb++) {
1217 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1218 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1219 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1220 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1221 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1222 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1223 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1224 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1225 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1228 a * row_diff_base_functions(
J) * col_diff_base_functions(
L);
1233 t_assemble(
i,
k) += diffDiff(
J,
L) * D2(
i,
J,
k,
L);
1238 t0 = diffDiff(
J,
L) * S(
J,
L);
1239 t_assemble(0, 0) += t0;
1240 t_assemble(1, 1) += t0;
1241 t_assemble(2, 2) += t0;
1245 ++col_diff_base_functions;
1248 ++row_diff_base_functions;
1250 for (; row_bb != row_nb_base_functions; row_bb++) {
1251 ++row_diff_base_functions;
1262 if (diagonal_block) {
1263 for (
int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
1265 for (; col_bb != row_bb + 1; col_bb++) {
1267 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1268 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1269 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1270 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1271 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1272 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1273 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1274 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1275 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1277 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 0),
1278 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 1),
1279 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 2),
1280 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 0),
1281 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 1),
1282 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 2),
1283 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 0),
1284 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 1),
1285 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 2));
1286 t_off_side(
i,
k) = t_assemble(
k,
i);
1291 const int *row_ind = &*row_data.getIndices().begin();
1292 const int *col_ind = &*col_data.getIndices().begin();
1294 col_ind, &locK_P_F(0, 0), ADD_VALUES);
1296 if (row_type != col_type || row_side != col_side) {
1297 transLocK_P_F.resize(col_nb_dofs, row_nb_dofs,
false);
1298 noalias(transLocK_P_F) = trans(locK_P_F);
1300 row_ind, &transLocK_P_F(0, 0), ADD_VALUES);
1310 commonData(common_data) {
1316 EntityType col_type,
1317 DataForcesAndSourcesCore::EntData &row_data,
1318 DataForcesAndSourcesCore::EntData &col_data) {
1322 const int row_nb_dofs = row_data.getIndices().size();
1325 const int col_nb_dofs = col_data.getIndices().size();
1330 locK_P_RHO.resize(row_nb_dofs, col_nb_dofs,
false);
1344 const int row_nb_gauss_pts = row_data.getN().size1();
1345 const int row_nb_base_functions = row_data.getN().size2();
1346 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1348 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
1352 double a = w * (
n /
rho) * pow(
rho / rho_ref,
n);
1355 for (; row_bb != row_nb_dofs / 3; row_bb++) {
1356 t1(
i) =
a * row_diff_base_functions(
I) * P(
i,
I);
1357 auto base_functions = col_data.getFTensor0N(gg, 0);
1358 for (
int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1365 k(
i) += t1(
i) * base_functions;
1368 ++row_diff_base_functions;
1370 for (; row_bb != row_nb_base_functions; row_bb++) {
1371 ++row_diff_base_functions;
1379 &*row_data.getIndices().begin(), col_nb_dofs,
1380 &*col_data.getIndices().begin(), &
locK_P_RHO(0, 0),
1388 : mField(m_field), commonData(common_data), postProc(m_field),
1389 postProcElastic(m_field),
1408 PetscBool analysis_mesh_flg = PETSC_FALSE;
1411 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Bone remodeling post-process",
1416 "is used for testing, calculates mass and energy at each time step",
"",
1419 CHKERR PetscOptionsBool(
"-analysis_mesh",
1420 "saves analysis mesh at each time step",
"",
1421 analysis_mesh_flg, &analysis_mesh_flg, PETSC_NULL);
1423 CHKERR PetscOptionsScalar(
1424 "-equilibrium_stop_rate",
1425 "is used to stop calculations when equilibium state is achieved",
"",
1431 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
1432 "Bone remodeling post-process",
"none");
1435 CHKERR PetscOptionsInt(
"-my_output_prt",
1436 "frequncy how often results are dumped on hard disk",
1452 list_of_operators.push_back(
1468 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
1482 int mass_vec_ghost[] = {0};
1484 1, 1, mass_vec_ghost, &mass_vec);
1486 CHKERR VecZeroEntries(mass_vec);
1487 CHKERR VecDuplicate(mass_vec, &energ_vec);
1505 CHKERR VecAssemblyBegin(mass_vec);
1506 CHKERR VecAssemblyEnd(mass_vec);
1508 CHKERR VecSum(mass_vec, &mass_sum);
1510 CHKERR VecAssemblyBegin(energ_vec);
1511 CHKERR VecAssemblyEnd(energ_vec);
1513 CHKERR VecSum(energ_vec, &energ_sum);
1515 CHKERR PetscPrintf(PETSC_COMM_WORLD,
1516 "Time: %g Mass: %6.5g Elastic energy: %6.5g \n", ts_t,
1517 mass_sum, energ_sum);
1518 CHKERR VecDestroy(&mass_vec);
1519 CHKERR VecDestroy(&energ_vec);
1522 double equilibrium_rate =
1524 double equilibrium_mass_rate =
1526 if (equilibrium_rate <
rate) {
1529 "Energy equilibrium state is achieved! Difference = %0.6g %%. \n",
1530 equilibrium_rate * 100);
1532 if (equilibrium_mass_rate <
rate) {
1535 "Mass equilibrium state is achieved! Difference = %0.6g %%. \n",
1536 equilibrium_mass_rate * 100);
1544#if PETSC_VERSION_LE(3, 8, 0)
1545 CHKERR TSGetTimeStepNumber(ts, &step);
1547 CHKERR TSGetStepNumber(ts, &step);
1550 if ((step) %
pRT == 0) {
1553 sss <<
"out_" << step <<
".h5m";
1556 "PARALLEL=WRITE_PART");
1557 if (analysis_mesh_flg) {
1559 ttt <<
"analysis_mesh_" << step <<
".h5m";
1561 "PARALLEL=WRITE_PART");
1568 ss <<
"out_elastic_" << step <<
".h5m";
1570 "PARALLEL=WRITE_PART");
1579 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Bone remodeling",
"none");
1582 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"", 2,
1589 string name = it->getName();
1590 if (name.compare(0, 14,
"NO_REMODELLING") == 0) {
1624 bool add_rho_field =
false;
1632 add_rho_field =
true;
1672 if (add_rho_field) {
1684 "MESH_NODE_POSITIONS");
1702 "MESH_NODE_POSITIONS");
1714 "MESH_NODE_POSITIONS");
1729 "DISPLACEMENTS",
"MESH_NODE_POSITIONS",
false,
true);
1747 list_of_operators_rhs.push_back(
1749 list_of_operators_rhs.push_back(
1764 list_of_operators_lhs.push_back(
1766 list_of_operators_rhs.push_back(
1772 list_of_operators_lhs.push_back(
1774 list_of_operators_lhs.push_back(
1778 list_of_operators_lhs.push_back(
1792 CHKERR MetaNeummanForces::addNeumannBCElements(
mField,
"DISPLACEMENTS");
1797 CHKERR MetaNeummanForces::setMomentumFluxOperators(
1801 PETSC_NULL,
"DISPLACEMENTS");
1806 for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
1809 mit->second->methodsOp.push_back(
1812 for (boost::ptr_map<string, NodalForce>::iterator mit =
1815 mit->second->methodsOp.push_back(
1818 for (boost::ptr_map<string, EdgeForce>::iterator mit =
1821 mit->second->methodsOp.push_back(
1869 CHKERR TSSetIFunction(ts,
F, PETSC_NULL, PETSC_NULL);
1870 CHKERR TSSetIJacobian(ts, A, A, PETSC_NULL, PETSC_NULL);
1872#if PETSC_VERSION_GE(3, 8, 0)
1873 CHKERR TSSetMaxTime(ts, ftime);
1875 CHKERR TSSetFromOptions(ts);
1879 CHKERR TSGetSNES(ts, &snes);
1881 CHKERR SNESGetKSP(snes, &ksp);
1883 CHKERR KSPGetPC(ksp, &pc);
1884 PetscBool is_pcfs = PETSC_FALSE;
1885 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1889 if (is_pcfs == PETSC_TRUE) {
1894 problem_ptr->
getName(),
ROW,
"DISPLACEMENTS", 0, 3, &is_disp);
1896 problem_ptr->
getName(),
ROW,
"RHO", 0, 1, &is_rho);
1899 CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
1900 CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
1901 CHKERR ISDestroy(&is_disp);
1902 CHKERR ISDestroy(&is_rho);
1914#if PETSC_VERSION_GE(3, 7, 0)
1915 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1918 CHKERR TSGetTime(ts, &ftime);
1919 PetscInt steps, snesfails, rejects, nonlinits, linits;
1920#if PETSC_VERSION_LE(3, 8, 0)
1921 CHKERR TSGetTimeStepNumber(ts, &steps);
1923 CHKERR TSGetStepNumber(ts, &steps);
1926 CHKERR TSGetSNESFailures(ts, &snesfails);
1927 CHKERR TSGetStepRejections(ts, &rejects);
1928 CHKERR TSGetSNESIterations(ts, &nonlinits);
1929 CHKERR TSGetKSPIterations(ts, &linits);
1930 PetscPrintf(PETSC_COMM_WORLD,
1931 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
1933 steps, rejects, snesfails, ftime, nonlinits, linits);
1944 int row_side, EntityType row_type,
1945 DataForcesAndSourcesCore::EntData &row_data) {
1948 if (row_type != MBVERTEX)
1951 const int nb_gauss_pts = row_data.getN().size1();
1958 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1961 energy += vol * psi;
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#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)
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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
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 addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2 or H1 field gradient.
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)
const double n
refractive index of diffusive medium
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', DIM1 > J
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
MoFEMErrorCode recordFreeEnergy_dC(Remodeling::CommonData &common_data, T &dC, B1 &psi)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
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.
double poisson_ratio
Poisson ratio.
FTensor::Index< 'm', 3 > m
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}...
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)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpCalculateStressTangentWithAdolc(Remodeling::CommonData &common_data)
Remodeling::CommonData & commonData
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Remodeling::CommonData & commonData
OpCalculateStressTangent(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)
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 field rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar field 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 reference 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.