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);
 
  448    CHKERR recordFreeEnergy_dC<FTensor::Tensor0<FTensor::PackPtr<double *, 1>>,
 
  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);
 
  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;
 
  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();
 
 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",
 
 1450      auto &list_of_operators = 
postProc.getOpPtrVector();
 
 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 .
#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
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 poisson_ratio
Poisson ratio.
double young_modulus
Young modulus.
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}}\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)
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.
Specialization for double precision scalar field values calculation.
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 parameter for U_t (see PETSc Time Solver documentation)
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.