59                         {
   60 
   62 
   63  PetscBool flg_config = PETSC_FALSE;
   64  is_atom_testing = PETSC_FALSE;
   65  with_adol_c = PETSC_FALSE;
   66  char my_config_file_name[255];
   70                    
   73  rHo_ref = 1000;   
   74  pSi_ref = 2.75e4; 
   75                    
   76  rHo_max = rHo_ref + 0.5 * rHo_ref;
   77  rHo_min = rHo_ref - 0.5 * rHo_ref;
   78  b = 0;
   79  R0 = 0.0; 
   80  cUrrent_psi = 0.0;
   81  cUrrent_mass = 0.0;
   82  nOremodellingBlock = false;
   83  less_post_proc = PETSC_FALSE;
   84 
   86  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling parameters",
   87                           "none");
   88 
   89  
   90  CHKERR PetscOptionsString(
"-my_config", 
"configuration file name", 
"",
 
   91                            "my_config.in", my_config_file_name, 255,
   92                            &flg_config);
   93 
   94  if (flg_config) {
   95    CHKERR PetscPrintf(PETSC_COMM_WORLD, 
"Config file: %s loaded.\n",
 
   96                       my_config_file_name);
   97    try {
   98      ifstream ini_file(my_config_file_name);
   99      if (!ini_file.is_open()) {
  100        SETERRQ(PETSC_COMM_SELF, 1, "*** -my_config does not exist ***");
  101      }
  102      po::variables_map vm;
  103      po::options_description config_file_options;
  104      
  105 
  106      config_file_options.add_options()(
  107          "young_modulus", po::value<double>(&
young_modulus)->default_value(1));
 
  108      config_file_options.add_options()(
  109          "poisson_ratio", po::value<double>(&
poisson_ratio)->default_value(1));
 
  110      config_file_options.add_options()(
  111          "c", po::value<double>(&
c)->default_value(1));
 
  112      config_file_options.add_options()(
  113          "m", po::value<double>(&
m)->default_value(1));
 
  114      config_file_options.add_options()(
  115          "n", po::value<double>(&
n)->default_value(1));
 
  116      config_file_options.add_options()(
  117          "rHo_ref", po::value<double>(&rHo_ref)->default_value(1));
  118      config_file_options.add_options()(
  119          "pSi_ref", po::value<double>(&pSi_ref)->default_value(1));
  120      config_file_options.add_options()(
  121          "R0", po::value<double>(&R0)->default_value(1));
  122 
  123      po::parsed_options parsed =
  124          parse_config_file(ini_file, config_file_options, true);
  125      store(parsed, vm);
  126      po::notify(vm);
  127 
  128    } catch (const std::exception &ex) {
  129      std::ostringstream ss;
  130      ss << ex.what() << std::endl;
  132    }
  133  }
  134 
  135  CHKERR PetscOptionsScalar(
"-young_modulus", 
"get young modulus", 
"",
 
  137 
  138  CHKERR PetscOptionsScalar(
"-poisson_ratio", 
"get poisson_ratio", 
"",
 
  140 
  141  CHKERR PetscOptionsScalar(
"-c", 
"density evolution (growth) velocity [d/m^2]",
 
  142                            "", 
c, &
c, PETSC_NULL);
 
  143 
  144  CHKERR PetscOptionsScalar(
"-m", 
"algorithmic exponent", 
"", 
m, &
m,
 
  145                            PETSC_NULL);
  146 
  147  CHKERR PetscOptionsScalar(
"-n", 
"porosity exponent", 
"", 
n, &
n, PETSC_NULL);
 
  148 
  149  CHKERR PetscOptionsInt(
"-b", 
"bell function exponent", 
"", b, &b, PETSC_NULL);
 
  150 
  151  CHKERR PetscOptionsScalar(
"-rho_ref", 
"reference bone density", 
"", rHo_ref,
 
  152                            &rHo_ref, PETSC_NULL);
  153 
  154  CHKERR PetscOptionsScalar(
"-rho_max", 
"reference bone density", 
"", rHo_max,
 
  155                            &rHo_max, PETSC_NULL);
  156  CHKERR PetscOptionsScalar(
"-rho_min", 
"reference bone density", 
"", rHo_min,
 
  157                            &rHo_min, PETSC_NULL);
  158 
  159  CHKERR PetscOptionsScalar(
"-psi_ref", 
"reference energy density", 
"", pSi_ref,
 
  160                            &pSi_ref, PETSC_NULL);
  161 
  162  CHKERR PetscOptionsScalar(
"-r0", 
"mass source", 
"", R0, &R0, PETSC_NULL);
 
  163 
  164  CHKERR PetscOptionsBool(
"-my_is_atom_test",
 
  165                          "is used with testing, exit with error when diverged",
  166                          "", is_atom_testing, &is_atom_testing, PETSC_NULL);
  167 
  168  CHKERR PetscOptionsBool(
"-less_post_proc",
 
  169                          "is used to reduce output file size", "",
  170                          less_post_proc, &less_post_proc, PETSC_NULL);
  171 
  173      "-with_adolc",
  174      "calculate the material tangent with automatic differentiation", "",
  175      with_adol_c, &with_adol_c, PETSC_NULL);
  176 
  179 
  180  CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  181                     "Young's modulus E[Pa]:                %4.3g\n",
  183  CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  184                     "Poisson ratio nu[-]                   %4.3g\n",
  186  CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  187                     "Lame coefficient lambda[Pa]:          %4.3g\n", 
lambda);
 
  188  CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  189                     "Lame coefficient mu[Pa]:              %4.3g\n", 
mu);
 
  190  CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  191                     "Density growth velocity [d/m2]        %4.3g\n", 
c);
 
  192  CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  193                     "Algorithmic exponent m[-]:            %4.3g\n", 
m);
 
  194  CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  195                     "Porosity exponent n[-]:               %4.3g\n", 
n);
 
  196  CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  197                     "Reference density rHo_ref[kg/m3]:     %4.3g\n", rHo_ref);
  198  CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  199                     "Reference energy pSi_ref[Pa]:         %4.3g\n", pSi_ref);
  200  CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  201                     "Mass conduction R0[kg/m3s]:           %4.3g\n", R0);
  202  CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  203                     "Bell function coefficent b[-]:        %4.3g\n", b);
  204  if (b != 0) {
  205    CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  206                       "Max density rHo_max[kg/m3]:         %4.3g\n", rHo_max);
  207    CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  208                       "Min density rHo_min[kg/m3]:         %4.3g\n", rHo_min);
  209  }
  210 
  211  PetscOptionsEnd();
  212 
  214}
  215 
  216inline double Remodeling::CommonData::getCFromDensity(
const double &
rho) {
 
  217  double mid_rho = 0.5 * (rHo_max + rHo_min);
  218  double var_h = (
rho - mid_rho) / (rHo_max - mid_rho);
 
  219 
  220  return 1 / (1 + (b != 0) * pow(var_h, 2 * b));
  221}
  222 
  223inline double Remodeling::CommonData::getCFromDensityDiff(
const double &
rho) {
 
  224  double mid_rho = 0.5 * (rHo_max + rHo_min);
  225  double var_h = (
rho - mid_rho) / (rHo_max - mid_rho);
 
  226  return (b != 0) * (-2) * b * pow(var_h, 2 * b - 1) /
  227         ((rHo_max - mid_rho) * pow(pow(var_h, 2 * b) + 1, 2));
  228}
  229 
  230OpGetRhoTimeDirevative::OpGetRhoTimeDirevative(
  231    Remodeling::CommonData &common_data)
  234      commonData(common_data) {}
  235 
  237OpGetRhoTimeDirevative::doWork(int side, EntityType type,
  238                               DataForcesAndSourcesCore::EntData &data) {
  239 
  241 
  242  const int nb_gauss_pts = data.getN().size1();
  243  if (!nb_gauss_pts)
  245  const int nb_base_functions = data.getN().size2();
  246 
  247  auto base_function = data.getFTensor0N();
  248  commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
  249  if (type == MBVERTEX) {
  250    commonData.data.vecRhoDt.clear();
  251  }
  252 
  253  int nb_dofs = data.getIndices().size();
  254  if (nb_dofs == 0)
  256  double *array;
  257  CHKERR VecGetArray(getFEMethod()->ts_u_t, &array);
 
  259 
  260  for (int gg = 0; gg < nb_gauss_pts; gg++) {
  261    int bb = 0;
  262    for (; bb < nb_dofs; bb++) {
  263      const double field_data = array[data.getLocalIndices()[bb]];
  264      drho_dt += field_data * base_function;
  265      ++base_function;
  266    }
  267    
  268    for (; bb != nb_base_functions; bb++)
  269      ++base_function;
  270    ++drho_dt;
  271  }
  272 
  273  CHKERR VecRestoreArray(getFEMethod()->ts_u_t, &array);
 
  274 
  276}
  277 
  281      commonData(common_data) {
  282  
  283  doVertices = true;
  284  doEdges = false;
  285  doQuads = false;
  286  doTris = false;
  287  doTets = false;
  288  doPrisms = false;
  289}
  290 
  293                          DataForcesAndSourcesCore::EntData &data) {
  294 
  296 
  297  if (type != MBVERTEX)
  299 
  300  if (!commonData.data.gradDispPtr) {
  302            "Gradient at integration points of displacement is not calculated");
  303  }
  304  auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
  305 
  306  const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
  307  if (!commonData.data.rhoPtr) {
  309            "Density at integration points is not calculated");
  310  }
  311 
  313 
  314  commonData.data.vecPsi.resize(nb_gauss_pts, false);
  316  commonData.data.vecR.resize(nb_gauss_pts, false);
  318 
  319  commonData.data.matInvF.resize(9, nb_gauss_pts, false);
  320  auto invF = getFTensor2FromMat<3, 3>(commonData.data.matInvF);
  321  commonData.data.vecDetF.resize(nb_gauss_pts, false);
  324  matS.resize(nb_gauss_pts, 6, false);
  326  commonData.data.matP.resize(9, nb_gauss_pts, false);
  327  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
 
  328  MatrixDouble &matF = commonData.data.matGradientOfDeformation;
 
  329  matF.resize(9, nb_gauss_pts, false);
  331 
  332  const double c = commonData.c;
 
  333  const double n = commonData.n;
 
  334  const double m = commonData.m;
 
  335  const double rho_ref = commonData.rHo_ref;
  336  const double psi_ref = commonData.pSi_ref;
  337  const double lambda = commonData.lambda;
 
  338  const double mu = commonData.mu;
 
  342 
  343  for (int gg = 0; gg != nb_gauss_pts; gg++) {
  344 
  345    
  346    F(
i, 
I) = grad(
i, 
I);
 
  350 
  353 
  354    const double log_det = log(det_f);
  355    psi = 
F(
i, 
J) * 
F(
i, 
J) - 3 - 2 * log_det;
 
  357    psi += 0.5 * 
lambda * log_det * log_det;
 
  358 
  359    const double coef = 
lambda * log_det - 
mu;
 
  361    S(
i, 
I) = invF(
i, 
J) ^ 
P(
J, 
I);
 
  362 
  363    const double kp = commonData.getCFromDensity(
rho);
 
  364    R = 
c * kp * (pow(
rho / rho_ref, 
n - 
m) * psi - psi_ref);
 
  365 
  366    ++det_f;
  367    ++invF;
  368    ++grad;
  370    ++psi;
  371    ++S;
  375  }
  376 
  378}
  379 
  380OpCalculateStressTangentWithAdolc::OpCalculateStressTangentWithAdolc(
  381    Remodeling::CommonData &common_data)
  384      commonData(common_data) {
  385  
  386  doVertices = true;
  387  doEdges = false;
  388  doQuads = false;
  389  doTris = false;
  390  doTets = false;
  391  doPrisms = false;
  392}
  393 
  395    int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
  396 
  398 
  399  if (type != MBVERTEX)
  401 
  402  auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
  403  const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
  404 
  405  vecC.resize(6, false);
  407  MatrixDouble &dS_dC = commonData.data.matMaterialTangent;
 
  408  dS_dC.resize(6, 6 * nb_gauss_pts, false);
  409  dS_dC.clear();
  411 
  413  matS.resize(nb_gauss_pts, 6, false);
  415  MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
 
  416  dP_dF.resize(81, nb_gauss_pts, false);
  417  MatrixDouble &matF = commonData.data.matGradientOfDeformation;
 
  418  matF.resize(9, nb_gauss_pts, false);
  421  commonData.data.matP.resize(9, nb_gauss_pts, false);
  422  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
 
  423  commonData.data.vecPsi.resize(nb_gauss_pts, false);
  425 
  432 
  433  for (int gg = 0; gg != nb_gauss_pts; gg++) {
  434 
  435    
  436    F(
i, 
I) = grad(
i, 
I);
 
  440 
  441    
  442    
  443    
  445 
  446    
  447    CHKERR recordFreeEnergy_dC<FTensor::Tensor0<FTensor::PackPtr<double *, 1>>,
 
  449        commonData, C, psi);
  450 
  451    int r_g = ::gradient(commonData.tAg, 6, &vecC[0], &matS(gg, 0));
  452    if (r_g < 0) {
  453      
  454      
  456              "ADOL-C function evaluation with error");
  457    }
  458 
  459    
  460    S(0, 0) *= 2;
  461    S(1, 1) *= 2;
  462    S(2, 2) *= 2;
  463 
  464    
  466 
  467    const int is = 6 * gg;
  468    double *hessian_ptr[6] = {&dS_dC(0, is), &dS_dC(1, is), &dS_dC(2, is),
  469                              &dS_dC(3, is), &dS_dC(4, is), &dS_dC(5, is)};
  470    int r_h = ::hessian(commonData.tAg, 6, &vecC[0], hessian_ptr);
  471    if (r_h < 0) {
  472      
  473      
  475              "ADOL-C function evaluation with error");
  476    }
  477 
  478    
  479    
  480    
  481    
  482    
  483    
  484    
  485    
  486    
  487    
  488    
  489    
  490    
  491    
  492    for (int ii = 0; ii < 6; ii++) {
  493      for (int jj = 0; jj < ii; jj++) {
  494        dS_dC(jj, is + ii) = dS_dC(ii, is + jj);
  495      }
  496    }
  497 
  498    
  499    
  500    
  501    D1(0, 0, 0, 0) *= 4;
  502    D1(1, 1, 1, 1) *= 4;
  503    D1(2, 2, 2, 2) *= 4;
  504 
  505    D1(0, 0, 1, 1) *= 4;
  506    D1(1, 1, 0, 0) *= 4;
  507    D1(1, 1, 2, 2) *= 4;
  508    D1(2, 2, 1, 1) *= 4;
  509    D1(0, 0, 2, 2) *= 4;
  510    D1(2, 2, 0, 0) *= 4;
  511 
  512    D1(0, 0, 0, 1) *= 2;
  513    D1(0, 0, 0, 2) *= 2;
  514    D1(0, 0, 1, 2) *= 2;
  515    D1(0, 1, 0, 0) *= 2;
  516    D1(0, 2, 0, 0) *= 2;
  517    D1(1, 2, 0, 0) *= 2;
  518 
  519    D1(1, 1, 0, 1) *= 2;
  520    D1(1, 1, 0, 2) *= 2;
  521    D1(1, 1, 1, 2) *= 2;
  522    D1(0, 1, 1, 1) *= 2;
  523    D1(0, 2, 1, 1) *= 2;
  524    D1(1, 2, 1, 1) *= 2;
  525 
  526    D1(2, 2, 0, 1) *= 2;
  527    D1(2, 2, 0, 2) *= 2;
  528    D1(2, 2, 1, 2) *= 2;
  529    D1(0, 1, 2, 2) *= 2;
  530    D1(0, 2, 2, 2) *= 2;
  531    D1(1, 2, 2, 2) *= 2;
  532 
  534    
  535 
  536    ++grad;
  537    ++S;
  539    ++D1;
  540    ++D2;
  542    ++psi;
  543  }
  544 
  546}
  547 
  549                                   std::vector<EntityHandle> &map_gauss_pts,
  550                                   Remodeling::CommonData &common_data)
  553      commonData(common_data), postProcMesh(post_proc_mesh),
  554      mapGaussPts(map_gauss_pts) {
  555  
  556  doVertices = true;
  557  doEdges = false;
  558  doQuads = false;
  559  doTris = false;
  560  doTets = false;
  561  doPrisms = false;
  562}
  563 
  566                         DataForcesAndSourcesCore::EntData &data) {
  567 
  569 
  570  if (type != MBVERTEX)
  572  double def_VAL[9];
  573  bzero(def_VAL, 9 * sizeof(double));
  574 
  577                                     th_piola2, MB_TAG_CREAT | MB_TAG_SPARSE,
  578                                     def_VAL);
  579 
  582                                     MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
  583 
  586                                     MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
  587 
  590  Tag th_prin_stress_vect1, th_prin_stress_vect2, th_prin_stress_vect3;
 
  594 
  596                                       th_piola1, MB_TAG_CREAT | MB_TAG_SPARSE,
  597                                       def_VAL);
  598 
  600                                       principal, MB_TAG_CREAT | MB_TAG_SPARSE,
  601                                       def_VAL);
  602 
  604                                       th_prin_stress_vect1,
  605                                       MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
  607                                       th_prin_stress_vect2,
  608                                       MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
  610                                       th_prin_stress_vect3,
  611                                       MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
  612 
  614                                       MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
  615 
  617                                       th_grad_rho,
  618                                       MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
  619  }
  620 
  624  auto P = getFTensor2FromMat<3, 3>(
commonData.data.matP);
 
  627 
  629  auto grad_rho = getFTensor1FromMat<3>(*
commonData.data.gradRhoPtr);
 
  630 
  634 
  635  for (int gg = 0; gg != nb_gauss_pts; gg++) {
  636 
  637    double val = psi;
  641 
  642    
  643    
  644    for (int ii = 0; ii != 3; ii++) {
  645      for (int jj = 0; jj != 3; jj++) {
  646        stress(ii, jj) = S(ii, jj);
  647      }
  648    }
  650                                     &stress(0, 0));
  651 
  653 
  654      for (int ii = 0; ii != 3; ii++) {
  655        for (int jj = 0; jj != 3; jj++) {
  656          stress(ii, jj) = 
P(ii, jj);
 
  657        }
  658      }
  660                                       &stress(0, 0));
  663      const double t2[] = {grad_rho(0), grad_rho(1), grad_rho(2)};
  665 
  666      
  667      noalias(eigen_vectors) = stress;
  671      
  672      
  673      int n = 3, lda = 3, info, lwork = -1;
 
  674      double wkopt;
  675      info = 
lapack_dsyev(
'V', 
'U', 
n, &(eigen_vectors.data()[0]), lda,
 
  676                          &(eigen_values.data()[0]), &wkopt, lwork);
  677      if (info != 0)
  678        SETERRQ1(PETSC_COMM_SELF, 1,
  679                 "something is wrong with lapack_dsyev info = %d", info);
  681      double work[lwork];
  682      info = 
lapack_dsyev(
'V', 
'U', 
n, &(eigen_vectors.data()[0]), lda,
 
  683                          &(eigen_values.data()[0]), work, lwork);
  684      if (info != 0)
  685        SETERRQ1(PETSC_COMM_SELF, 1,
  686                 "something is wrong with lapack_dsyev info = %d", info);
  687 
  688      for (int ii = 0; ii < 3; ii++) {
  689        prin_stress_vect1[ii] = eigen_vectors(0, ii);
  690        prin_stress_vect2[ii] = eigen_vectors(1, ii);
  691        prin_stress_vect3[ii] = eigen_vectors(2, ii);
  692      }
  693 
  695                                       &eigen_values[0]);
  697                                       1, &*prin_stress_vect1.begin());
  699                                       1, &*prin_stress_vect2.begin());
  701                                       1, &*prin_stress_vect3.begin());
  702    }
  703 
  704    ++S;
  706    ++psi;
  709    ++grad_rho;
  710  }
  711 
  713}
  714 
  715OpCalculateStressTangent::OpCalculateStressTangent(
  716    Remodeling::CommonData &common_data)
  719      commonData(common_data) {
  720  
  721  doVertices = true;
  722  doEdges = false;
  723  doQuads = false;
  724  doTris = false;
  725  doTets = false;
  726  doPrisms = false;
  727}
  728 
  730OpCalculateStressTangent::doWork(int side, EntityType type,
  731                                DataForcesAndSourcesCore::EntData &data) {
  732 
  734 
  735  if (type != MBVERTEX)
  737 
  738  auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
  739  const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
  740 
  741  commonData.data.matInvF.resize(9, nb_gauss_pts, false);
  742  auto invF = getFTensor2FromMat<3, 3>(commonData.data.matInvF);
  743  commonData.data.vecDetF.resize(nb_gauss_pts, false);
  745 
  746  MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
 
  747  dP_dF.resize(81, nb_gauss_pts, false);
  748  MatrixDouble &matF = commonData.data.matGradientOfDeformation;
 
  749  matF.resize(9, nb_gauss_pts, false);
  752  commonData.data.matP.resize(9, nb_gauss_pts, false);
  753  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
 
  754  commonData.data.vecPsi.resize(nb_gauss_pts, false);
  762  const double lambda = commonData.lambda;
 
  763  const double mu = commonData.mu;
 
  764 
  765  for (int gg = 0; gg != nb_gauss_pts; gg++) {
  766 
  767    
  768    F(
i, 
I) = grad(
i, 
I);
 
  772 
  775 
  776    const double log_det = log(det_f);
  777    psi = 
F(
i, 
J) * 
F(
i, 
J) - 3 - 2 * log_det;
 
  779    psi += 0.5 * 
lambda * log_det * log_det;
 
  780 
  781    const double var = 
lambda * log_det - 
mu;
 
  783 
  784    const double coef = 
mu - 
lambda * log_det;
 
  785    
  786 
  788        invF(
J, 
i) * (invF(
L, 
k) * 
lambda) + invF(
L, 
i) * (invF(
J, 
k) * coef);
 
  789 
  790    D2(0, 0, 0, 0) += 
mu;
 
  791    D2(0, 1, 0, 1) += 
mu;
 
  792    D2(0, 2, 0, 2) += 
mu;
 
  793    D2(1, 0, 1, 0) += 
mu;
 
  794    D2(1, 1, 1, 1) += 
mu;
 
  795    D2(1, 2, 1, 2) += 
mu;
 
  796    D2(2, 0, 2, 0) += 
mu;
 
  797    D2(2, 1, 2, 1) += 
mu;
 
  798    D2(2, 2, 2, 2) += 
mu;
 
  799 
  800    ++det_f;
  801    ++invF;
  802    ++grad;
  804    ++D2;
  806    ++psi;
  807  }
  808 
  810}
  811 
  812OpAssmbleStressRhs::OpAssmbleStressRhs(Remodeling::CommonData &common_data)
  815      commonData(common_data) {}
  816 
  818OpAssmbleStressRhs::doWork(int side, EntityType type,
  819                           DataForcesAndSourcesCore::EntData &data) {
  820 
  822 
  823  const int nb_dofs = data.getIndices().size();
  824  if (!nb_dofs)
  826  nF.resize(nb_dofs, false);
  827  nF.clear();
  828 
  829  const int nb_gauss_pts = data.getN().size1();
  830  const int nb_base_functions = data.getN().size2();
  831  auto diff_base_functions = data.getFTensor1DiffN<3>();
  832  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
 
  834 
  835  
  836  
  837  
  838  
  839 
  840  const double n = commonData.n;
 
  841  const double rho_ref = commonData.rHo_ref;
  842 
  846 
  847  for (int gg = 0; gg != nb_gauss_pts; gg++) {
  848 
  849    
  850    double w = getVolume() * getGaussPts()(3, gg);
 
  851    double a = 
w * pow(
rho / rho_ref, 
n);
 
  852 
  853    int bb = 0;
  854    for (; bb != nb_dofs / 3; bb++) {
  855      double *ptr = &nF[3 * bb];
  857      t1(
i) += 
a * 
P(
i, 
I) * diff_base_functions(
I);
 
  858      
  859      
  860      ++diff_base_functions;
  861    }
  862    
  863    
  864    for (; bb != nb_base_functions; bb++) {
  865      ++diff_base_functions;
  866    }
  868    
  869    
  871  }
  872 
  873  
  875      getFEMethod()->ts_F, 
  876      nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
  877 
  879}
  880 
  881OpAssmbleRhoRhs::OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
  884      commonData(common_data) {}
  885 
  887OpAssmbleRhoRhs::doWork(int side, EntityType type,
  888                        DataForcesAndSourcesCore::EntData &data) {
  889 
  891 
  892  const int nb_dofs = data.getIndices().size();
  893  if (!nb_dofs)
  895  nF.resize(nb_dofs, false);
  896  nF.clear();
  897 
  898  const int nb_gauss_pts = data.getN().size1();
  899  const int nb_base_functions = data.getN().size2();
  900  auto base_functions = data.getFTensor0N();
  901  auto diff_base_functions = data.getFTensor1DiffN<3>();
  902  if (!commonData.data.rhoPtr) {
  904  }
  906  if (!commonData.data.gradRhoPtr) {
  908  }
  909  auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
  911 
  912  if (commonData.data.vecRhoDt.size() != nb_gauss_pts) {
  913    commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
  914    commonData.data.vecRhoDt.clear();
  915    
  916    
  917  }
  919 
  920  const double R0 = commonData.R0;
  921 
  923 
  924  for (int gg = 0; gg != nb_gauss_pts; gg++) {
  925 
  926    
  927    double w = getVolume() * getGaussPts()(3, gg);
 
  928 
  929    int bb = 0;
  930    for (; bb != nb_dofs; bb++) {
  931      nF[bb] += 
w * base_functions * drho_dt; 
 
  932      nF[bb] -= 
w * base_functions * 
R;       
 
  933      nF[bb] -= 
w * R0 * diff_base_functions(
I) * grad_rho(
I); 
 
  934      ++base_functions;
  935      ++diff_base_functions;
  936    }
  937    
  938    
  939    for (; bb != nb_base_functions; bb++) {
  940      ++base_functions;
  941      ++diff_base_functions;
  942    }
  943 
  944    
  946    ++grad_rho;
  947    ++drho_dt;
  949  }
  950 
  951  
  953      getFEMethod()->ts_F, 
  954      nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
  955 
  957}
  958 
  959OpAssmbleRhoLhs_dRho::OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
  962      commonData(common_data) {
  963  sYmm = true;
  964}
  965 
  967OpAssmbleRhoLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
  968                             EntityType col_type,
  969                             DataForcesAndSourcesCore::EntData &row_data,
  970                             DataForcesAndSourcesCore::EntData &col_data) {
  971 
  973 
  974  const int row_nb_dofs = row_data.getIndices().size();
  975  if (!row_nb_dofs)
  977  const int col_nb_dofs = col_data.getIndices().size();
  978  if (!col_nb_dofs)
  980 
  981  
  982  locK_rho_rho.resize(row_nb_dofs, col_nb_dofs, false);
  983  locK_rho_rho.clear();
  984 
  985  const int row_nb_gauss_pts = row_data.getN().size1();
  986  if (!row_nb_gauss_pts)
  988  const int row_nb_base_functions = row_data.getN().size2();
  989 
  990  const double ts_a = getFEMethod()->ts_a;
  991  const double R0 = commonData.R0;
  992  const double c = commonData.c;
 
  993  const double n = commonData.n;
 
  994  const double m = commonData.m;
 
  995  const double rho_ref = commonData.rHo_ref;
  996  const double psi_ref = commonData.pSi_ref;
  997 
  999 
 1002  auto row_base_functions = row_data.getFTensor0N();
 1003  auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
 1004 
 1005  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
 1006 
 1007    
 1008    double w = getVolume() * getGaussPts()(3, gg);
 
 1009 
 1010    const double kp = commonData.getCFromDensity(
rho);
 
 1011    const double kp_diff = commonData.getCFromDensityDiff(
rho);
 
 1012    const double a = 
c * kp * ((
n - 
m) / 
rho) * pow(
rho / rho_ref, 
n - 
m) * psi;
 
 1013    const double a_diff =
 1014        c * kp_diff * (pow(
rho / rho_ref, 
n - 
m) * psi - psi_ref);
 
 1015 
 1016    int row_bb = 0;
 1017    for (; row_bb != row_nb_dofs; row_bb++) {
 1018      auto col_base_functions = col_data.getFTensor0N(gg, 0);
 1019      auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
 1020      for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
 1021        locK_rho_rho(row_bb, col_bb) +=
 1022            ts_a * 
w * row_base_functions * col_base_functions;
 
 1023        locK_rho_rho(row_bb, col_bb) -=
 1024            w * R0 * row_diff_base_functions(
I) * col_diff_base_functions(
I);
 
 1025        locK_rho_rho(row_bb, col_bb) -=
 1026            a * 
w * row_base_functions * col_base_functions;
 
 1027        locK_rho_rho(row_bb, col_bb) -=
 1028            a_diff * 
w * row_base_functions * col_base_functions;
 
 1029 
 1030        ++col_base_functions;
 1031        ++col_diff_base_functions;
 1032      }
 1033      ++row_base_functions;
 1034      ++row_diff_base_functions;
 1035    }
 1036    for (; row_bb != row_nb_base_functions; row_bb++) {
 1037      ++row_base_functions;
 1038      ++row_diff_base_functions;
 1039    }
 1040    ++psi;
 1042  }
 1043 
 1044  
 1045 
 1047                      &*row_data.getIndices().begin(), col_nb_dofs,
 1048                      &*col_data.getIndices().begin(), &locK_rho_rho(0, 0),
 1049                      ADD_VALUES);
 1050 
 1051  
 1052  if (row_side != col_side || row_type != col_type) {
 1053    transLocK_rho_rho.resize(col_nb_dofs, row_nb_dofs, false);
 1054    noalias(transLocK_rho_rho) = trans(locK_rho_rho);
 1056                        &*col_data.getIndices().begin(), row_nb_dofs,
 1057                        &*row_data.getIndices().begin(),
 1058                        &transLocK_rho_rho(0, 0), ADD_VALUES);
 1059  }
 1060 
 1062}
 1063 
 1064OpAssmbleRhoLhs_dF::OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data)
 1067      commonData(common_data) {
 1068  sYmm = false;
 1069}
 1070 
 1072OpAssmbleRhoLhs_dF::doWork(int row_side, int col_side, EntityType row_type,
 1073                           EntityType col_type,
 1074                           DataForcesAndSourcesCore::EntData &row_data,
 1075                           DataForcesAndSourcesCore::EntData &col_data) {
 1076 
 1078 
 1079  const int row_nb_dofs = row_data.getIndices().size();
 1080  if (!row_nb_dofs)
 1082  const int col_nb_dofs = col_data.getIndices().size();
 1083  if (!col_nb_dofs)
 1085 
 1086  
 1087  locK_rho_F.resize(row_nb_dofs, col_nb_dofs, false);
 1088  locK_rho_F.clear();
 1089 
 1090  const int row_nb_gauss_pts = row_data.getN().size1();
 1091  const int col_nb_base_functions = col_data.getN().size2();
 1092 
 1093  const double c = commonData.c;
 
 1094  const double n = commonData.n;
 
 1095  const double m = commonData.m;
 
 1096  const double rho_ref = commonData.rHo_ref;
 1097 
 1100 
 1102  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
 
 1103 
 1105  auto col_diff_base_functions = col_data.getFTensor1DiffN<3>();
 1106 
 1107  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
 1108 
 1109    
 1110    double w = getVolume() * getGaussPts()(3, gg);
 
 1111    const double kp = commonData.getCFromDensity(
rho);
 
 1112    const double a = 
w * 
c * kp * pow(
rho / rho_ref, 
n - 
m);
 
 1113 
 1114    int col_bb = 0;
 1115    for (; col_bb != col_nb_dofs / 3; col_bb++) {
 1116      t1(
i) = 
a * 
P(
i, 
I) * col_diff_base_functions(
I);
 
 1117      auto row_base_functions = row_data.getFTensor0N(gg, 0);
 1118      for (int row_bb = 0; row_bb != row_nb_dofs; row_bb++) {
 1120                                        &locK_rho_F(row_bb, 3 * col_bb + 1),
 1121                                        &locK_rho_F(row_bb, 3 * col_bb + 2));
 1122        k(
i) -= row_base_functions * t1(
i);
 
 1123        ++row_base_functions;
 1124      }
 1125      ++col_diff_base_functions;
 1126    }
 1127    for (; col_bb != col_nb_base_functions; col_bb++) {
 1128      ++col_diff_base_functions;
 1129    }
 1130 
 1133  }
 1134 
 1136                      &*row_data.getIndices().begin(), col_nb_dofs,
 1137                      &*col_data.getIndices().begin(), &locK_rho_F(0, 0),
 1138                      ADD_VALUES);
 1139 
 1141}
 1142template <bool ADOLC>
 1143OpAssmbleStressLhs_dF<ADOLC>::OpAssmbleStressLhs_dF(
 1144    Remodeling::CommonData &common_data)
 1147      commonData(common_data) {
 1148  sYmm = true;
 1149}
 1150template <bool ADOLC>
 1152    int row_side, int col_side, EntityType row_type, EntityType col_type,
 1153    DataForcesAndSourcesCore::EntData &row_data,
 1154    DataForcesAndSourcesCore::EntData &col_data) {
 1155 
 1157 
 1158  const int row_nb_dofs = row_data.getIndices().size();
 1159  if (!row_nb_dofs)
 1161  const int col_nb_dofs = col_data.getIndices().size();
 1162  if (!col_nb_dofs)
 1164 
 1165  const bool diagonal_block = (row_type == col_type) && (row_side == col_side);
 1166 
 1167  
 1168  locK_P_F.resize(row_nb_dofs, col_nb_dofs, false);
 1169  locK_P_F.clear();
 1170 
 1171  const int row_nb_gauss_pts = row_data.getN().size1();
 1172  const int row_nb_base_functions = row_data.getN().size2();
 1173 
 1175  MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
 
 1176  
 1177  
 1178 
 1179  double t0;
 1180 
 1181  const double n = commonData.n;
 
 1182  const double rho_ref = commonData.rHo_ref;
 1183 
 1191 
 1192  auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
 1196 
 1197  
 1198  
 1199 
 1200  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
 1201 
 1202    
 1203    double w = getVolume() * getGaussPts()(3, gg);
 
 1204    const double a = 
w * pow(
rho / rho_ref, 
n);
 
 1205 
 1206    int row_bb = 0;
 1207    for (; row_bb != row_nb_dofs / 3; row_bb++) {
 1208 
 1209      auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
 1210 
 1211      const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
 1212      int col_bb = 0;
 1213      for (; col_bb != final_bb; col_bb++) {
 1214 
 1216            &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
 1217            &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
 1218            &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
 1219            &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
 1220            &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
 1221            &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
 1222            &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
 1223            &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
 1224            &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
 1225 
 1227            a * row_diff_base_functions(
J) * col_diff_base_functions(
L);
 
 1228 
 1229        
 1230        
 1231        
 1232        t_assemble(
i, 
k) += diffDiff(
J, 
L) * D2(
i, 
J, 
k, 
L);
 
 1233 
 1234        if (ADOLC) {
 1235          
 1236          
 1237          t0 = diffDiff(
J, 
L) * S(
J, 
L);
 
 1238          t_assemble(0, 0) += t0;
 1239          t_assemble(1, 1) += t0;
 1240          t_assemble(2, 2) += t0;
 1241        }
 1242 
 1243        
 1244        ++col_diff_base_functions;
 1245      }
 1246 
 1247      ++row_diff_base_functions;
 1248    }
 1249    for (; row_bb != row_nb_base_functions; row_bb++) {
 1250      ++row_diff_base_functions;
 1251    }
 1252 
 1253    
 1254    ++D2;
 1255    ++S;
 1256    
 1258  }
 1259 
 1260  
 1261  if (diagonal_block) {
 1262    for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
 1263      int col_bb = 0;
 1264      for (; col_bb != row_bb + 1; col_bb++) {
 1266            &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
 1267            &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
 1268            &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
 1269            &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
 1270            &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
 1271            &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
 1272            &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
 1273            &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
 1274            &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
 1276            &locK_P_F(3 * col_bb + 0, 3 * row_bb + 0),
 1277            &locK_P_F(3 * col_bb + 0, 3 * row_bb + 1),
 1278            &locK_P_F(3 * col_bb + 0, 3 * row_bb + 2),
 1279            &locK_P_F(3 * col_bb + 1, 3 * row_bb + 0),
 1280            &locK_P_F(3 * col_bb + 1, 3 * row_bb + 1),
 1281            &locK_P_F(3 * col_bb + 1, 3 * row_bb + 2),
 1282            &locK_P_F(3 * col_bb + 2, 3 * row_bb + 0),
 1283            &locK_P_F(3 * col_bb + 2, 3 * row_bb + 1),
 1284            &locK_P_F(3 * col_bb + 2, 3 * row_bb + 2));
 1285        t_off_side(
i, 
k) = t_assemble(
k, 
i);
 
 1286      }
 1287    }
 1288  }
 1289 
 1290  const int *row_ind = &*row_data.getIndices().begin();
 1291  const int *col_ind = &*col_data.getIndices().begin();
 1293                      col_ind, &locK_P_F(0, 0), ADD_VALUES);
 1294 
 1295  if (row_type != col_type || row_side != col_side) {
 1296    transLocK_P_F.resize(col_nb_dofs, row_nb_dofs, false);
 1297    noalias(transLocK_P_F) = trans(locK_P_F);
 1299                        row_ind, &transLocK_P_F(0, 0), ADD_VALUES);
 1300  }
 1301 
 1303}
 1304 
 1305OpAssmbleStressLhs_dRho::OpAssmbleStressLhs_dRho(
 1306    Remodeling::CommonData &common_data)
 1309      commonData(common_data) {
 1310  sYmm = false; 
 1311}
 1312 
 1314OpAssmbleStressLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
 1315                                EntityType col_type,
 1316                                DataForcesAndSourcesCore::EntData &row_data,
 1317                                DataForcesAndSourcesCore::EntData &col_data) {
 1318 
 1320 
 1321  const int row_nb_dofs = row_data.getIndices().size();
 1322  if (!row_nb_dofs)
 1324  const int col_nb_dofs = col_data.getIndices().size();
 1325  if (!col_nb_dofs)
 1327 
 1328  
 1329  locK_P_RHO.resize(row_nb_dofs, col_nb_dofs, false);
 1330  locK_P_RHO.clear();
 1331 
 1333 
 1334  const double n = commonData.n;
 
 1335  const double rho_ref = commonData.rHo_ref;
 1336 
 1339 
 1340  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
 
 1342 
 1343  const int row_nb_gauss_pts = row_data.getN().size1();
 1344  const int row_nb_base_functions = row_data.getN().size2();
 1345  auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
 1346 
 1347  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
 1348 
 1349    
 1350    double w = getVolume() * getGaussPts()(3, gg);
 
 1351    double a = 
w * (
n / 
rho) * pow(
rho / rho_ref, 
n);
 
 1352 
 1353    int row_bb = 0;
 1354    for (; row_bb != row_nb_dofs / 3; row_bb++) {
 1355      t1(
i) = 
a * row_diff_base_functions(
I) * 
P(
i, 
I);
 
 1356      auto base_functions = col_data.getFTensor0N(gg, 0);
 1357      for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
 1358        
 1359        
 1360        
 1362                                        &locK_P_RHO(3 * row_bb + 1, col_bb),
 1363                                        &locK_P_RHO(3 * row_bb + 2, col_bb));
 1364        k(
i) += t1(
i) * base_functions;
 
 1365        ++base_functions;
 1366      }
 1367      ++row_diff_base_functions;
 1368    }
 1369    for (; row_bb != row_nb_base_functions; row_bb++) {
 1370      ++row_diff_base_functions;
 1371    }
 1372 
 1375  }
 1376 
 1378                      &*row_data.getIndices().begin(), col_nb_dofs,
 1379                      &*col_data.getIndices().begin(), &locK_P_RHO(0, 0),
 1380                      ADD_VALUES);
 1381 
 1383}
 1384 
 1386                                 Remodeling::CommonData &common_data)
 1387    : mField(m_field), commonData(common_data), postProc(m_field),
 1388      postProcElastic(m_field),
 1389      
 1390      iNit(false) {}
 1391 
 1395}
 1396 
 1400}
 1401 
 1404 
 1405  PetscBool mass_postproc = PETSC_FALSE;
 1406  PetscBool equilibrium_flg = PETSC_FALSE;
 1407  PetscBool analysis_mesh_flg = PETSC_FALSE;
 1408  rate = 1;
 1409 
 1410  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling post-process",
 1411                           "none");
 1412 
 1414      "-mass_postproc",
 1415      "is used for testing, calculates mass and energy at each time step", "",
 1416      mass_postproc, &mass_postproc, PETSC_NULL);
 1417 
 1418  CHKERR PetscOptionsBool(
"-analysis_mesh",
 
 1419                          "saves analysis mesh at each time step", "",
 1420                          analysis_mesh_flg, &analysis_mesh_flg, PETSC_NULL);
 1421 
 1422  CHKERR PetscOptionsScalar(
 
 1423      "-equilibrium_stop_rate",
 1424      "is used to stop calculations when equilibium state is achieved", "",
 1425      rate, &rate, &equilibrium_flg);
 1426  PetscOptionsEnd();
 1427 
 1429 
 1430    PetscOptionsBegin(PETSC_COMM_WORLD, "",
 1431                             "Bone remodeling post-process", "none");
 1432 
 1434    CHKERR PetscOptionsInt(
"-my_output_prt",
 
 1435                           "frequncy how often results are dumped on hard disk",
 1437    PetscOptionsEnd();
 1438 
 1442    if (!commonData.less_post_proc) {
 1444    }
 1445    
 1446    
 1447 
 1448    {
 1449      auto &list_of_operators = 
postProc.getOpPtrVector();
 
 1450 
 1451      list_of_operators.push_back(
 1454          "RHO", commonData.data.gradRhoPtr));
 1455      
 1457          "DISPLACEMENTS", commonData.data.gradDispPtr));
 1461    }
 1462 
 1464    CHKERR postProcElastic.addFieldValuesPostProc(
"DISPLACEMENTS");
 
 1465    CHKERR postProcElastic.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
 
 1466    CHKERR postProcElastic.addFieldValuesGradientPostProc(
"DISPLACEMENTS");
 
 1467    std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
 1468        commonData.elasticPtr->setOfBlocks.begin();
 1469    for (; sit != commonData.elasticPtr->setOfBlocks.end(); sit++) {
 1471          postProcElastic.postProcMesh, postProcElastic.mapGaussPts,
 1472          "DISPLACEMENTS", sit->second, postProcElastic.commonData));
 1473    }
 1474 
 1476  }
 1477 
 1478  if (mass_postproc) {
 1481    int mass_vec_ghost[] = {0};
 1483                          1, 1, mass_vec_ghost, &mass_vec);
 1484 
 1485    CHKERR VecZeroEntries(mass_vec);
 
 1486    CHKERR VecDuplicate(mass_vec, &energ_vec);
 
 1487 
 1488    Remodeling::Fe energyCalc(
mField);
 
 1490                       true);
 1491    energyCalc.getOpPtrVector().push_back(
 1494        "RHO", commonData.data.gradRhoPtr));
 1495    energyCalc.getOpPtrVector().push_back(
 1497                                                 commonData.data.gradDispPtr));
 1499    energyCalc.getOpPtrVector().push_back(
 1500        new OpMassAndEnergyCalculation("RHO", commonData, energ_vec, mass_vec));
 1502                                    &energyCalc);
 1503 
 1504    CHKERR VecAssemblyBegin(mass_vec);
 
 1505    CHKERR VecAssemblyEnd(mass_vec);
 
 1506    double mass_sum;
 1507    CHKERR VecSum(mass_vec, &mass_sum);
 
 1508 
 1509    CHKERR VecAssemblyBegin(energ_vec);
 
 1510    CHKERR VecAssemblyEnd(energ_vec);
 
 1511    double energ_sum;
 1512    CHKERR VecSum(energ_vec, &energ_sum);
 
 1513 
 1514    CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
 1515                       "Time: %g Mass: %6.5g Elastic energy: %6.5g \n", 
ts_t,
 
 1516                       mass_sum, energ_sum);
 1517    CHKERR VecDestroy(&mass_vec);
 
 1518    CHKERR VecDestroy(&energ_vec);
 
 1519 
 1520    if (equilibrium_flg) {
 1521      double equilibrium_rate =
 1522          fabs(energ_sum - commonData.cUrrent_psi) / energ_sum;
 1523      double equilibrium_mass_rate =
 1524          fabs(mass_sum - commonData.cUrrent_mass) / mass_sum;
 1525      if (equilibrium_rate < rate) {
 1527            PETSC_COMM_WORLD,
 1528            "Energy equilibrium state is achieved! Difference = %0.6g %%. \n",
 1529            equilibrium_rate * 100);
 1530      }
 1531      if (equilibrium_mass_rate < rate) {
 1533            PETSC_COMM_WORLD,
 1534            "Mass equilibrium state is achieved! Difference = %0.6g %%. \n",
 1535            equilibrium_mass_rate * 100);
 1536      }
 1537      commonData.cUrrent_psi = energ_sum;
 1538      commonData.cUrrent_mass = mass_sum;
 1539    }
 1540  }
 1541 
 1543#if PETSC_VERSION_LE(3, 8, 0)
 1545#else
 1547#endif
 1548 
 1550 
 1551    ostringstream sss;
 1552    sss << 
"out_" << 
step << 
".h5m";
 
 1555                                            "PARALLEL=WRITE_PART");
 1556    if (analysis_mesh_flg) {
 1557      ostringstream ttt;
 1558      ttt << 
"analysis_mesh_" << 
step << 
".h5m";
 
 1560                                          "PARALLEL=WRITE_PART");
 1561    }
 1562 
 1563    if (commonData.nOremodellingBlock) {
 1565                                      &postProcElastic);
 1566      ostringstream ss;
 1567      ss << 
"out_elastic_" << 
step << 
".h5m";
 
 1568      CHKERR postProcElastic.postProcMesh.write_file(ss.str().c_str(), 
"MOAB",
 
 1569                                                     "PARALLEL=WRITE_PART");
 1570    }
 1571  }
 1573}
 1574 
 1576 
 1578  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling", "none");
 1579 
 1580  commonData.oRder = 2;
 1581  CHKERR PetscOptionsInt(
"-my_order", 
"default approximation order", 
"", 2,
 
 1582                         &commonData.oRder, PETSC_NULL);
 1583 
 1584  PetscOptionsEnd();
 1585 
 1586  CHKERR mField.get_moab().get_entities_by_type(0, MBTET, commonData.tEts_all);
 
 1588    string name = it->getName();
 1589    if (name.compare(0, 14, "NO_REMODELLING") == 0) {
 1590      commonData.nOremodellingBlock = true;
 1592      CHKERR this->mField.get_moab().get_entities_by_type(
 
 1593          meshset, MBTET, commonData.tEts_block, true);
 1594      commonData.tEts_all =
 1595          subtract(commonData.tEts_all, commonData.tEts_block);
 1596    }
 1597  }
 1598 
 1600}
 1601 
 1603 
 1605 
 1606  
 1607  
 1608  
 1609  commonData.bitLevel.set(0);
 1611      0, 3, commonData.bitLevel);
 1612  int order = commonData.oRder;
 
 1613 
 1614  
 1615  CHKERR mField.add_field(
"DISPLACEMENTS", 
H1,
 
 1618  
 1621 
 1622  
 1623  bool add_rho_field = false;
 1624  if (!mField.check_field("RHO")) {
 1627    
 1628    CHKERR mField.add_ents_to_field_by_type(commonData.tEts_all, MBTET, 
"RHO");
 
 1629    
 1630    CHKERR mField.synchronise_field_entities(
"RHO");
 
 1631    add_rho_field = true;
 1632 
 1633    CHKERR mField.set_field_order(0, MBVERTEX, 
"RHO", 1);
 
 1634    CHKERR mField.set_field_order(0, MBEDGE, 
"RHO", 
order - 1);
 
 1635    CHKERR mField.set_field_order(0, MBTRI, 
"RHO", 
order - 1);
 
 1636    CHKERR mField.set_field_order(0, MBTET, 
"RHO", 
order - 1);
 
 1637  }
 1638 
 1639  
 1640  CHKERR mField.add_ents_to_field_by_type(0, MBTET, 
"DISPLACEMENTS");
 
 1641  CHKERR mField.add_ents_to_field_by_type(0, MBTET, 
"MESH_NODE_POSITIONS");
 
 1642 
 1643  
 1644  CHKERR mField.set_field_order(0, MBVERTEX, 
"DISPLACEMENTS", 1);
 
 1645  CHKERR mField.set_field_order(0, MBEDGE, 
"DISPLACEMENTS", 
order);
 
 1646  CHKERR mField.set_field_order(0, MBTRI, 
"DISPLACEMENTS", 
order);
 
 1647  CHKERR mField.set_field_order(0, MBTET, 
"DISPLACEMENTS", 
order);
 
 1648 
 1649  
 1650 
 1651  
 1652  {
 1653    
 1654    
 1655    
 1656    
 1657    
 1658    
 1659    
 1660    
 1661    CHKERR mField.set_field_order(0, MBEDGE, 
"MESH_NODE_POSITIONS", 2);
 
 1662  }
 1663 
 1664  CHKERR mField.set_field_order(0, MBVERTEX, 
"MESH_NODE_POSITIONS", 1);
 
 1665 
 1666  
 1667  CHKERR mField.build_fields();
 
 1668 
 1669  
 1670  
 1671  if (add_rho_field) {
 1673                                                      MBVERTEX, "RHO");
 1674    
 1675    
 1676    
 1677    
 1678    
 1679  }
 1680 
 1681  
 1683                                                    "MESH_NODE_POSITIONS");
 1684  CHKERR mField.loop_dofs(
"MESH_NODE_POSITIONS", ent_method_material);
 
 1685 
 1687}
 1688 
 1691 
 1692  CHKERR mField.add_finite_element(
"FE_REMODELLING", 
MF_ZERO);
 
 1693  CHKERR mField.modify_finite_element_add_field_row(
"FE_REMODELLING",
 
 1694                                                    "DISPLACEMENTS");
 1695  CHKERR mField.modify_finite_element_add_field_col(
"FE_REMODELLING",
 
 1696                                                    "DISPLACEMENTS");
 1697  CHKERR mField.modify_finite_element_add_field_data(
"FE_REMODELLING",
 
 1698                                                     "DISPLACEMENTS");
 1699  CHKERR mField.modify_finite_element_add_field_data(
"FE_REMODELLING", 
"RHO");
 
 1700  CHKERR mField.modify_finite_element_add_field_data(
"FE_REMODELLING",
 
 1701                                                     "MESH_NODE_POSITIONS");
 1702  CHKERR mField.modify_finite_element_add_field_row(
"FE_REMODELLING", 
"RHO");
 
 1703  CHKERR mField.modify_finite_element_add_field_col(
"FE_REMODELLING", 
"RHO");
 
 1704  CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_all, MBTET,
 
 1705                                                   "FE_REMODELLING");
 1706 
 1707  CHKERR mField.add_finite_element(
"ELASTIC");
 
 1708  CHKERR mField.modify_finite_element_add_field_row(
"ELASTIC", 
"DISPLACEMENTS");
 
 1709  CHKERR mField.modify_finite_element_add_field_col(
"ELASTIC", 
"DISPLACEMENTS");
 
 1710  CHKERR mField.modify_finite_element_add_field_data(
"ELASTIC",
 
 1711                                                     "DISPLACEMENTS");
 1712  CHKERR mField.modify_finite_element_add_field_data(
"ELASTIC",
 
 1713                                                     "MESH_NODE_POSITIONS");
 1714 
 1715  if (commonData.nOremodellingBlock) {
 1716    CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_block,
 
 1717                                                     MBTET, "ELASTIC");
 1718  }
 1719  
 1720  commonData.elasticMaterialsPtr =
 1722  commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
 1724  CHKERR commonData.elasticMaterialsPtr->setBlocks(
 
 1725      commonData.elasticPtr->setOfBlocks);
 1726  CHKERR commonData.elasticPtr->addElement(
"ELASTIC", 
"DISPLACEMENTS");
 
 1727  CHKERR commonData.elasticPtr->setOperators(
 
 1728      "DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
 1729 
 1730  CHKERR mField.build_finite_elements();
 
 1731  CHKERR mField.build_adjacencies(commonData.bitLevel);
 
 1732 
 1733  
 1734  
 1735  commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(
new VectorDouble());
 
 1736  commonData.data.gradRhoPtr =
 1738  commonData.data.gradDispPtr =
 1740 
 1741  
 1743                     false, false);
 1744  auto &list_of_operators_rhs = commonData.feRhs->getOpPtrVector();
 1745  
 1746  list_of_operators_rhs.push_back(
 1748  list_of_operators_rhs.push_back(
 1750  list_of_operators_rhs.push_back(new OpGetRhoTimeDirevative(commonData));
 1751  
 1753      "DISPLACEMENTS", commonData.data.gradDispPtr));
 1755  list_of_operators_rhs.push_back(new OpAssmbleStressRhs(commonData));
 1756  list_of_operators_rhs.push_back(new OpAssmbleRhoRhs(commonData));
 1757 
 1758  
 1760                     false, false);
 1761  auto &list_of_operators_lhs = commonData.feLhs->getOpPtrVector();
 1762  
 1763  list_of_operators_lhs.push_back(
 1765  list_of_operators_rhs.push_back(
 1767  
 1769      "DISPLACEMENTS", commonData.data.gradDispPtr));
 1770  if (commonData.with_adol_c) {
 1771    list_of_operators_lhs.push_back(
 1772        new OpCalculateStressTangentWithAdolc(commonData));
 1773    list_of_operators_lhs.push_back(
 1774        new OpAssmbleStressLhs_dF<true>(commonData));
 1775  } else {
 1776    list_of_operators_lhs.push_back(new OpCalculateStressTangent(commonData));
 1777    list_of_operators_lhs.push_back(
 1778        new OpAssmbleStressLhs_dF<false>(commonData));
 1779  }
 1780  list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dRho(commonData));
 1781  list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dF(commonData));
 1782  list_of_operators_lhs.push_back(new OpAssmbleStressLhs_dRho(commonData));
 1783 
 1785}
 1786 
 1788 
 1790  
 1791  CHKERR MetaNeummanForces::addNeumannBCElements(mField, 
"DISPLACEMENTS");
 
 1794 
 1795  
 1796  CHKERR MetaNeummanForces::setMomentumFluxOperators(
 
 1797      mField, commonData.neumannForces, PETSC_NULL, "DISPLACEMENTS");
 1798  
 1800                                       PETSC_NULL, "DISPLACEMENTS");
 1801  
 1803                                      "DISPLACEMENTS");
 1804 
 1805  for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
 1806           commonData.neumannForces.begin();
 1807       mit != commonData.neumannForces.end(); mit++) {
 1808    mit->second->methodsOp.push_back(
 1810  }
 1811  for (boost::ptr_map<string, NodalForce>::iterator mit =
 1812           commonData.nodalForces.begin();
 1813       mit != commonData.nodalForces.end(); mit++) {
 1814    mit->second->methodsOp.push_back(
 1816  }
 1817  for (boost::ptr_map<string, EdgeForce>::iterator mit =
 1818           commonData.edgeForces.begin();
 1819       mit != commonData.edgeForces.end(); mit++) {
 1820    mit->second->methodsOp.push_back(
 1822  }
 1823 
 1825}
 1826 
 1828 
 1830  commonData.dm_name = "DM_REMODELING";
 1831  
 1833  CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
 
 1834  CHKERR DMSetType(commonData.dm, commonData.dm_name);
 
 1836  
 1838                            commonData.bitLevel);
 1839  
 1840  
 1841  CHKERR DMSetFromOptions(commonData.dm);
 
 1842  
 1845 
 1846  if (mField.check_finite_element("FORCE_FE")) {
 1848  }
 1849  if (mField.check_finite_element("PRESSURE_FE")) {
 1851  }
 1852  mField.getInterface<
ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
 
 1853  CHKERR DMSetUp(commonData.dm);
 
 1854 
 1856}
 1857 
 1859 
 1861 
 1862  Vec F = commonData.F;
 
 1863  Vec D = commonData.D;
 
 1864  Mat 
A = commonData.A;
 
 1865  DM dm = commonData.dm;
 1866  TS ts = commonData.ts;
 1867 
 1868  CHKERR TSSetIFunction(ts, 
F, PETSC_NULL, PETSC_NULL);
 
 1869  CHKERR TSSetIJacobian(ts, A, A, PETSC_NULL, PETSC_NULL);
 
 1870  double ftime = 1;
 1871#if PETSC_VERSION_GE(3, 8, 0)
 1872  CHKERR TSSetMaxTime(ts, ftime);
 
 1873#endif
 1874  CHKERR TSSetFromOptions(ts);
 
 1876 
 1877  SNES snes;
 1878  CHKERR TSGetSNES(ts, &snes);
 
 1879  KSP ksp;
 1880  CHKERR SNESGetKSP(snes, &ksp);
 
 1881  PC pc;
 1882  CHKERR KSPGetPC(ksp, &pc);
 
 1883  PetscBool is_pcfs = PETSC_FALSE;
 1884  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
 1885 
 1886  
 1887  
 1888  if (is_pcfs == PETSC_TRUE) {
 1889    IS is_disp, is_rho;
 1893        problem_ptr->
getName(), 
ROW, 
"DISPLACEMENTS", 0, 3, &is_disp);
 
 1895        problem_ptr->
getName(), 
ROW, 
"RHO", 0, 1, &is_rho);
 
 1898    CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
 
 1899    CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
 
 1900    CHKERR ISDestroy(&is_disp);
 
 1901    CHKERR ISDestroy(&is_rho);
 
 1902  }
 1903 
 1904  
 1908  {
 1911  }
 1912 
 1913#if PETSC_VERSION_GE(3, 7, 0)
 1914  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
 
 1915#endif
 1917  CHKERR TSGetTime(ts, &ftime);
 
 1918  PetscInt steps, snesfails, rejects, nonlinits, linits;
 1919#if PETSC_VERSION_LE(3, 8, 0)
 1920  CHKERR TSGetTimeStepNumber(ts, &steps);
 
 1921#else
 1922  CHKERR TSGetStepNumber(ts, &steps);
 
 1923#endif
 1924 
 1925  CHKERR TSGetSNESFailures(ts, &snesfails);
 
 1926  CHKERR TSGetStepRejections(ts, &rejects);
 
 1927  CHKERR TSGetSNESIterations(ts, &nonlinits);
 
 1928  CHKERR TSGetKSPIterations(ts, &linits);
 
 1929  PetscPrintf(PETSC_COMM_WORLD,
 1930              "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
 1931              "linits %D\n",
 1932              steps, rejects, snesfails, ftime, nonlinits, linits);
 1933 
 1934  if (commonData.is_atom_testing) {
 1935    if (commonData.cUrrent_psi < 1.67 || commonData.cUrrent_psi > 1.68)
 1937  }
 1938 
 1940}
 1941 
 1943    int row_side, EntityType row_type,
 1944    DataForcesAndSourcesCore::EntData &row_data) {
 1946 
 1947  if (row_type != MBVERTEX)
 1949 
 1950  const int nb_gauss_pts = row_data.getN().size1();
 1951  
 1954 
 1955  double energy = 0;
 1956  double mass = 0;
 1957  for (int gg = 0; gg < nb_gauss_pts; gg++) {
 1958 
 1959    double vol = getVolume() * getGaussPts()(3, gg);
 1960    energy += vol * psi;
 1962 
 1963    ++psi;
 1965  }
 1966 
 1967  CHKERR VecSetValue(energVec, 0, energy, ADD_VALUES);
 
 1968  CHKERR VecSetValue(massVec, 0, mass, ADD_VALUES);
 
 1970}
 1971 
 1972} 
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT)
#define TENSOR2_MAT_PTR(MAT)
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT)
#define TENSOR4_MAT_PTR(MAT)
@ 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 ...
@ 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 ...
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.
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.
#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
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
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
boost::shared_ptr< MatrixDouble > gradDispPtr
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
virtual moab::Interface & get_moab()=0
virtual int get_comm_rank() const =0
Deprecated interface functions.
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.
TS ts
PETSc time stepping solver object.
PetscReal ts_t
Current time value.
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.
Volume finite element base.
MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
MoFEMErrorCode operator()()
Main operator function executed for each loop iteration.
MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
MoFEM::Interface & mField
MonitorPostProc(MoFEM::Interface &m_field, std::map< int, NonlinearElasticElement::BlockData > &set_of_blocks, NonlinearElasticElement::MyVolumeFE &fe_elastic_energy, ConvectiveMassElement::MyVolumeFE &fe_kinetic_energy)
PostProcVolumeOnRefinedMesh postProc
structure grouping operators and data used for calculation of nonlinear elastic element
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateStress(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts)
std::vector< EntityHandle > & mapGaussPts
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, DataAtIntegrationPts &common_data, BlockData &data)
DataAtIntegrationPts & commonData
moab::Interface & postProcMesh
PetscErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Force scale operator for reading two columns.