v0.9.1
Remodeling.cpp
Go to the documentation of this file.
1 /**
2  * \file Remodeling.cpp
3  * \example Remodeling.cpp
4  *
5  * \brief Integration points are taken from OOFEM
6  * <http://www.oofem.org/resources/doc/oofemrefman/microplanematerial_8C_source.html>
7  *
8  */
9 
10 /*
11  * This file is part of MoFEM.
12  * MoFEM is free software: you can redistribute it and/or modify it under
13  * the terms of the GNU Lesser General Public License as published by the
14  * Free Software Foundation, either version 3 of the License, or (at your
15  * option) any later version.
16  *
17  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
18  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
19  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
20  * License for more details.
21  *
22  * You should have received a copy of the GNU Lesser General Public
23  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
24 
25 #include <boost/program_options.hpp>
26 using namespace std;
27 namespace po = boost::program_options;
28 
29 #include <BasicFiniteElements.hpp>
30 using namespace MoFEM;
31 #include <ElasticMaterials.hpp>
32 #include <Remodeling.hpp>
33 
34 #ifndef WITH_ADOL_C
35 #error "MoFEM need to be compiled with ADOL-C"
36 #endif
37 
38 #define TENSOR1_VEC_PTR(VEC) &VEC[0], &VEC[1], &VEC[2]
39 
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)
47 
48 #define TENSOR4_MAT_PTR(MAT) &MAT(0, 0), MAT.size2()
49 
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)
53 
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)
56 
57 #define SYMMETRIC_TENSOR2_VEC_PTR(VEC) \
58  &VEC[0], &VEC[1], &VEC[2], &VEC[3], &VEC[4], &VEC[5]
59 
60 namespace BoneRemodeling {
61 
63 
64  PetscBool flg_config = PETSC_FALSE;
65  is_atom_testing = PETSC_FALSE;
66  with_adol_c = PETSC_FALSE;
67  char my_config_file_name[255];
68  double young_modulus = 5.0e+9; // Set here some typical value for bone 5 GPa
69  double poisson_ratio = 0.3; // Set here some typical value for bone
70  c = 4.0e-5; // Set here some typical value for bone. days/ m2 // this
71  // parameters governs how fast we can achieve equilibrium
72  m = 3.25; // Set some parameter typical for bone
73  n = 2.25; // Set Some parameter typical for bone
74  rHo_ref = 1000; // Set Some parameter typical for bone. kg/m3
75  pSi_ref = 2.75e4; // Set Some parameter typical for bone. energy / volume unit
76  // J/m3 = Pa
77  rHo_max = rHo_ref + 0.5 * rHo_ref;
78  rHo_min = rHo_ref - 0.5 * rHo_ref;
79  b = 0;
80  R0 = 0.0; // Set Some parameter typical for bone.
81  cUrrent_psi = 0.0;
82  cUrrent_mass = 0.0;
83  nOremodellingBlock = false;
84  less_post_proc = PETSC_FALSE;
85 
87  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling parameters",
88  "none");
89 
90  // config file
91  CHKERR PetscOptionsString("-my_config", "configuration file name", "",
92  "my_config.in", my_config_file_name, 255,
93  &flg_config);
94 
95  if (flg_config) {
96  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Config file: %s loaded.\n",
97  my_config_file_name);
98  try {
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 ***");
102  }
103  po::variables_map vm;
104  po::options_description config_file_options;
105  // std::cout << my_config_file_name << std::endl;
106 
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));
123 
124  po::parsed_options parsed =
125  parse_config_file(ini_file, config_file_options, true);
126  store(parsed, vm);
127  po::notify(vm);
128 
129  } catch (const std::exception &ex) {
130  std::ostringstream ss;
131  ss << ex.what() << std::endl;
132  SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
133  }
134  }
135 
136  CHKERR PetscOptionsScalar("-young_modulus", "get young modulus", "",
137  young_modulus, &young_modulus, PETSC_NULL);
138 
139  CHKERR PetscOptionsScalar("-poisson_ratio", "get poisson_ratio", "",
140  poisson_ratio, &poisson_ratio, PETSC_NULL);
141 
142  CHKERR PetscOptionsScalar("-c", "density evolution (growth) velocity [d/m^2]",
143  "", c, &c, PETSC_NULL);
144 
145  CHKERR PetscOptionsScalar("-m", "algorithmic exponent", "", m, &m,
146  PETSC_NULL);
147 
148  CHKERR PetscOptionsScalar("-n", "porosity exponent", "", n, &n, PETSC_NULL);
149 
150  CHKERR PetscOptionsInt("-b", "bell function exponent", "", b, &b, PETSC_NULL);
151 
152  CHKERR PetscOptionsScalar("-rho_ref", "reference bone density", "", rHo_ref,
153  &rHo_ref, PETSC_NULL);
154 
155  CHKERR PetscOptionsScalar("-rho_max", "reference bone density", "", rHo_max,
156  &rHo_max, PETSC_NULL);
157  CHKERR PetscOptionsScalar("-rho_min", "reference bone density", "", rHo_min,
158  &rHo_min, PETSC_NULL);
159 
160  CHKERR PetscOptionsScalar("-psi_ref", "reference energy density", "", pSi_ref,
161  &pSi_ref, PETSC_NULL);
162 
163  CHKERR PetscOptionsScalar("-r0", "mass source", "", R0, &R0, PETSC_NULL);
164 
165  CHKERR PetscOptionsBool("-my_is_atom_test",
166  "is used with testing, exit with error when diverged",
167  "", is_atom_testing, &is_atom_testing, PETSC_NULL);
168 
169  CHKERR PetscOptionsBool("-less_post_proc",
170  "is used to reduce output file size", "",
171  less_post_proc, &less_post_proc, PETSC_NULL);
172 
173  CHKERR PetscOptionsBool(
174  "-with_adolc",
175  "calculate the material tangent with automatic differentiation", "",
176  with_adol_c, &with_adol_c, PETSC_NULL);
177 
180 
181  CHKERR PetscPrintf(PETSC_COMM_WORLD,
182  "Young's modulus E[Pa]: %4.3g\n",
183  young_modulus);
184  CHKERR PetscPrintf(PETSC_COMM_WORLD,
185  "Poisson ratio nu[-] %4.3g\n",
186  poisson_ratio);
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);
205  if (b != 0) {
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);
210  }
211 
212  ierr = PetscOptionsEnd();
213  CHKERRQ(ierr);
214 
216 }
217 
218 inline double Remodeling::CommonData::getCFromDensity(const double &rho) {
219  double mid_rho = 0.5 * (rHo_max + rHo_min);
220  double var_h = (rho - mid_rho) / (rHo_max - mid_rho);
221 
222  return 1 / (1 + (b != 0) * pow(var_h, 2 * b));
223 }
224 
225 inline double Remodeling::CommonData::getCFromDensityDiff(const double &rho) {
226  double mid_rho = 0.5 * (rHo_max + rHo_min);
227  double var_h = (rho - mid_rho) / (rHo_max - mid_rho);
228  return (b != 0) * (-2) * b * pow(var_h, 2 * b - 1) /
229  ((rHo_max - mid_rho) * pow(pow(var_h, 2 * b) + 1, 2));
230 }
231 
232 OpGetRhoTimeDirevative::OpGetRhoTimeDirevative(
233  Remodeling::CommonData &common_data)
235  "RHO", UserDataOperator::OPROW),
236  commonData(common_data) {}
237 
241 
243 
244  const int nb_gauss_pts = data.getN().size1();
245  if (!nb_gauss_pts)
247  const int nb_base_functions = data.getN().size2();
248 
249  auto base_function = data.getFTensor0N();
250  commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
251  if (type == MBVERTEX) {
252  commonData.data.vecRhoDt.clear();
253  }
254 
255  int nb_dofs = data.getIndices().size();
256  if (nb_dofs == 0)
258  double *array;
259  CHKERR VecGetArray(getFEMethod()->ts_u_t, &array);
260  auto drho_dt = getFTensor0FromVec(commonData.data.vecRhoDt);
261 
262  for (int gg = 0; gg < nb_gauss_pts; gg++) {
263  int bb = 0;
264  for (; bb < nb_dofs; bb++) {
265  const double field_data = array[data.getLocalIndices()[bb]];
266  drho_dt += field_data * base_function;
267  ++base_function;
268  }
269  // It is possible to have more base functions than dofs
270  for (; bb != nb_base_functions; bb++)
271  ++base_function;
272  ++drho_dt;
273  }
274 
275  CHKERR VecRestoreArray(getFEMethod()->ts_u_t, &array);
276 
278 }
279 
282  "DISPLACEMENTS", UserDataOperator::OPROW),
283  commonData(common_data) {
284  // Set that this operator is executed only for vertices on element
285  doVertices = true;
286  doEdges = false;
287  doQuads = false;
288  doTris = false;
289  doTets = false;
290  doPrisms = false;
291 }
292 
294 OpCalculateStress::doWork(int side, EntityType type,
296 
298 
299  if (type != MBVERTEX)
301 
302  if (!commonData.data.gradDispPtr) {
303  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
304  "Gradient at integration points of displacement is not calculated");
305  }
306  auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
307 
308  const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
309  if (!commonData.data.rhoPtr) {
310  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
311  "Density at integration points is not calculated");
312  }
313 
315 
316  commonData.data.vecPsi.resize(nb_gauss_pts, false);
318  commonData.data.vecR.resize(nb_gauss_pts, false);
320 
321  commonData.data.matInvF.resize(9, nb_gauss_pts, false);
322  auto invF = getFTensor2FromMat<3, 3>(commonData.data.matInvF);
323  commonData.data.vecDetF.resize(nb_gauss_pts, false);
326  matS.resize(nb_gauss_pts, 6, false);
328  commonData.data.matP.resize(9, nb_gauss_pts, false);
329  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
331  matF.resize(9, nb_gauss_pts, false);
333 
334  const double c = commonData.c;
335  const double n = commonData.n;
336  const double m = commonData.m;
337  const double rho_ref = commonData.rHo_ref;
338  const double psi_ref = commonData.pSi_ref;
339  const double lambda = commonData.lambda;
340  const double mu = commonData.mu;
344 
345  for (int gg = 0; gg != nb_gauss_pts; gg++) {
346 
347  // Get deformation gradient
348  F(i, I) = grad(i, I);
349  F(0, 0) += 1;
350  F(1, 1) += 1;
351  F(2, 2) += 1;
352 
353  CHKERR determinantTensor3by3(F, det_f);
354  CHKERR invertTensor3by3(F, det_f, invF);
355 
356  const double log_det = log(det_f);
357  psi = F(i, J) * F(i, J) - 3 - 2 * log_det;
358  psi *= 0.5 * mu;
359  psi += 0.5 * lambda * log_det * log_det;
360 
361  const double coef = lambda * log_det - mu;
362  P(i, J) = mu * F(i, J) + coef * invF(J, i);
363  S(i, I) = invF(i, J) ^ P(J, I);
364 
365  const double kp = commonData.getCFromDensity(rho);
366  R = c * kp * (pow(rho / rho_ref, n - m) * psi - psi_ref);
367 
368  ++det_f;
369  ++invF;
370  ++grad;
371  ++rho;
372  ++psi;
373  ++S;
374  ++P;
375  ++R;
376  ++F;
377  }
378 
380 }
381 
383  Remodeling::CommonData &common_data)
385  "DISPLACEMENTS", UserDataOperator::OPROW),
386  commonData(common_data) {
387  // Set that this operator is executed only for vertices on element
388  doVertices = true;
389  doEdges = false;
390  doQuads = false;
391  doTris = false;
392  doTets = false;
393  doPrisms = false;
394 }
395 
397  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
398 
400 
401  if (type != MBVERTEX)
403 
404  auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
405  const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
406 
407  vecC.resize(6, false);
410  dS_dC.resize(6, 6 * nb_gauss_pts, false);
411  dS_dC.clear();
413 
415  matS.resize(nb_gauss_pts, 6, false);
418  dP_dF.resize(81, nb_gauss_pts, false);
420  matF.resize(9, nb_gauss_pts, false);
423  commonData.data.matP.resize(9, nb_gauss_pts, false);
424  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
425  commonData.data.vecPsi.resize(nb_gauss_pts, false);
427 
434 
435  for (int gg = 0; gg != nb_gauss_pts; gg++) {
436 
437  // Get deformation gradient
438  F(i, I) = grad(i, I);
439  F(0, 0) += 1;
440  F(1, 1) += 1;
441  F(2, 2) += 1;
442 
443  // Calculate Green defamation Tensor0
444  // ^ that means multiplication of Tensor 2 which result in symmetric Tensor
445  // rank 2.
446  C(I, J) = F(i, I) ^ F(i, J);
447 
448  // Calculate free energy
449  CHKERR recordFreeEnergy_dC<FTensor::Tensor0<FTensor::PackPtr<double *, 1>>,
451  commonData, C, psi);
452 
453  int r_g = ::gradient(commonData.tAg, 6, &vecC[0], &matS(gg, 0));
454  if (r_g < 0) {
455  // That means that energy function is not smooth and derivative
456  // can not be calculated,
457  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
458  "ADOL-C function evaluation with error");
459  }
460 
461  // Scale energy density and 2nd Piola Stress
462  S(0, 0) *= 2;
463  S(1, 1) *= 2;
464  S(2, 2) *= 2;
465 
466  // Calculate 1st Piola-Stress tensor
467  P(i, J) = F(i, I) * S(I, J);
468 
469  const int is = 6 * gg;
470  double *hessian_ptr[6] = {&dS_dC(0, is), &dS_dC(1, is), &dS_dC(2, is),
471  &dS_dC(3, is), &dS_dC(4, is), &dS_dC(5, is)};
472  int r_h = ::hessian(commonData.tAg, 6, &vecC[0], hessian_ptr);
473  if (r_h < 0) {
474  // That means that energy function is not smooth and derivative
475  // can not be calculated,
476  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
477  "ADOL-C function evaluation with error");
478  }
479 
480  // FIXME: TODO: Here tensor 4 with two minor symmetries is used, since
481  // \psi_{,ijkl} = \psi_{,klij} = \psi_{,lkji} = ... there is additional
482  // major symmetry. That is 21 independent components not 36 how is now
483  // implemented using Tenso4_ddg.
484  //
485  // That way we have to fill off-diagonal part for dS_dC. In future that
486  // need to be fixed by implementation of Tenso4 with additional major
487  // symmetries. Pleas consider doing such implementation by yourself. You
488  // have support of MoFEM developing team to guide you through this process.
489  //
490  // This deteriorate efficiency.
491  //
492  // TODO: Implement matrix with more symmetries.
493  //
494  for (int ii = 0; ii < 6; ii++) {
495  for (int jj = 0; jj < ii; jj++) {
496  dS_dC(jj, is + ii) = dS_dC(ii, is + jj);
497  }
498  }
499 
500  // As result that symmetry of deformation gradient C is used, terms which
501  // have mixed indices are two time to big, and terms which has two times
502  // midxed index are four times to big. Those terms are not multiplied.
503  D1(0, 0, 0, 0) *= 4;
504  D1(1, 1, 1, 1) *= 4;
505  D1(2, 2, 2, 2) *= 4;
506 
507  D1(0, 0, 1, 1) *= 4;
508  D1(1, 1, 0, 0) *= 4;
509  D1(1, 1, 2, 2) *= 4;
510  D1(2, 2, 1, 1) *= 4;
511  D1(0, 0, 2, 2) *= 4;
512  D1(2, 2, 0, 0) *= 4;
513 
514  D1(0, 0, 0, 1) *= 2;
515  D1(0, 0, 0, 2) *= 2;
516  D1(0, 0, 1, 2) *= 2;
517  D1(0, 1, 0, 0) *= 2;
518  D1(0, 2, 0, 0) *= 2;
519  D1(1, 2, 0, 0) *= 2;
520 
521  D1(1, 1, 0, 1) *= 2;
522  D1(1, 1, 0, 2) *= 2;
523  D1(1, 1, 1, 2) *= 2;
524  D1(0, 1, 1, 1) *= 2;
525  D1(0, 2, 1, 1) *= 2;
526  D1(1, 2, 1, 1) *= 2;
527 
528  D1(2, 2, 0, 1) *= 2;
529  D1(2, 2, 0, 2) *= 2;
530  D1(2, 2, 1, 2) *= 2;
531  D1(0, 1, 2, 2) *= 2;
532  D1(0, 2, 2, 2) *= 2;
533  D1(1, 2, 2, 2) *= 2;
534 
535  D2(i, J, k, L) = (F(i, I) * (D1(I, J, K, L) * F(k, K)));
536  // D2(i, J, k, L) += I(i, k) * S(J, L);
537 
538  ++grad;
539  ++S;
540  ++F;
541  ++D1;
542  ++D2;
543  ++P;
544  ++psi;
545  }
546 
548 }
549 
551  std::vector<EntityHandle> &map_gauss_pts,
552  Remodeling::CommonData &common_data)
554  "DISPLACEMENTS", UserDataOperator::OPROW),
555  commonData(common_data), postProcMesh(post_proc_mesh),
556  mapGaussPts(map_gauss_pts) {
557  // Set that this operator is executed only for vertices on element
558  doVertices = true;
559  doEdges = false;
560  doQuads = false;
561  doTris = false;
562  doTets = false;
563  doPrisms = false;
564 }
565 
567 OpPostProcStress::doWork(int side, EntityType type,
569 
571 
572  if (type != MBVERTEX)
574  double def_VAL[9];
575  bzero(def_VAL, 9 * sizeof(double));
576 
577  Tag th_piola2;
578  CHKERR postProcMesh.tag_get_handle("Piola_Stress2", 9, MB_TYPE_DOUBLE,
579  th_piola2, MB_TAG_CREAT | MB_TAG_SPARSE,
580  def_VAL);
581 
582  Tag th_psi;
583  CHKERR postProcMesh.tag_get_handle("psi", 1, MB_TYPE_DOUBLE, th_psi,
584  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
585 
586  Tag th_rho;
587  CHKERR postProcMesh.tag_get_handle("RHO", 1, MB_TYPE_DOUBLE, th_rho,
588  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
589 
590  Tag th_piola1;
591  Tag principal;
592  Tag th_prin_stress_vect1, th_prin_stress_vect2, th_prin_stress_vect3;
593  Tag th_R;
594  Tag th_grad_rho;
595  if (!commonData.less_post_proc) {
596 
597  CHKERR postProcMesh.tag_get_handle("Piola_Stress1", 9, MB_TYPE_DOUBLE,
598  th_piola1, MB_TAG_CREAT | MB_TAG_SPARSE,
599  def_VAL);
600 
601  CHKERR postProcMesh.tag_get_handle("Principal_stress", 3, MB_TYPE_DOUBLE,
602  principal, MB_TAG_CREAT | MB_TAG_SPARSE,
603  def_VAL);
604 
605  CHKERR postProcMesh.tag_get_handle("Principal_stress_S1", 3, MB_TYPE_DOUBLE,
606  th_prin_stress_vect1,
607  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
608  CHKERR postProcMesh.tag_get_handle("Principal_stress_S2", 3, MB_TYPE_DOUBLE,
609  th_prin_stress_vect2,
610  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
611  CHKERR postProcMesh.tag_get_handle("Principal_stress_S3", 3, MB_TYPE_DOUBLE,
612  th_prin_stress_vect3,
613  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
614 
615  CHKERR postProcMesh.tag_get_handle("R", 1, MB_TYPE_DOUBLE, th_R,
616  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
617 
618  CHKERR postProcMesh.tag_get_handle("grad_rho", 3, MB_TYPE_DOUBLE,
619  th_grad_rho,
620  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
621  }
622 
623  const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
626  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
629 
631  auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
632 
633  MatrixDouble3by3 stress(3, 3);
634  VectorDouble3 eigen_values(3);
635  MatrixDouble3by3 eigen_vectors(3, 3);
636 
637  for (int gg = 0; gg != nb_gauss_pts; gg++) {
638 
639  double val = psi;
640  CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &val);
641  val = rho;
642  CHKERR postProcMesh.tag_set_data(th_rho, &mapGaussPts[gg], 1, &val);
643 
644  // S is symmetric tensor, need to map it to non-symmetric tensor to save it
645  // on moab tag to be viable in paraview.
646  for (int ii = 0; ii != 3; ii++) {
647  for (int jj = 0; jj != 3; jj++) {
648  stress(ii, jj) = S(ii, jj);
649  }
650  }
651  CHKERR postProcMesh.tag_set_data(th_piola2, &mapGaussPts[gg], 1,
652  &stress(0, 0));
653 
654  if (!commonData.less_post_proc) {
655 
656  for (int ii = 0; ii != 3; ii++) {
657  for (int jj = 0; jj != 3; jj++) {
658  stress(ii, jj) = P(ii, jj);
659  }
660  }
661  CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
662  &stress(0, 0));
663  val = R;
664  CHKERR postProcMesh.tag_set_data(th_R, &mapGaussPts[gg], 1, &val);
665  const double t2[] = {grad_rho(0), grad_rho(1), grad_rho(2)};
666  CHKERR postProcMesh.tag_set_data(th_grad_rho, &mapGaussPts[gg], 1, t2);
667 
668  // Add calculation of eigen values, i.e. prinple stresses.
669  noalias(eigen_vectors) = stress;
670  VectorDouble3 prin_stress_vect1(3);
671  VectorDouble3 prin_stress_vect2(3);
672  VectorDouble3 prin_stress_vect3(3);
673  // LAPACK - eigenvalues and vectors. Applied twice for initial creates
674  // memory space
675  int n = 3, lda = 3, info, lwork = -1;
676  double wkopt;
677  info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
678  &(eigen_values.data()[0]), &wkopt, lwork);
679  if (info != 0)
680  SETERRQ1(PETSC_COMM_SELF, 1,
681  "something is wrong with lapack_dsyev info = %d", info);
682  lwork = (int)wkopt;
683  double work[lwork];
684  info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
685  &(eigen_values.data()[0]), work, lwork);
686  if (info != 0)
687  SETERRQ1(PETSC_COMM_SELF, 1,
688  "something is wrong with lapack_dsyev info = %d", info);
689 
690  for (int ii = 0; ii < 3; ii++) {
691  prin_stress_vect1[ii] = eigen_vectors(0, ii);
692  prin_stress_vect2[ii] = eigen_vectors(1, ii);
693  prin_stress_vect3[ii] = eigen_vectors(2, ii);
694  }
695 
696  CHKERR postProcMesh.tag_set_data(principal, &mapGaussPts[gg], 1,
697  &eigen_values[0]);
698  CHKERR postProcMesh.tag_set_data(th_prin_stress_vect1, &mapGaussPts[gg],
699  1, &*prin_stress_vect1.begin());
700  CHKERR postProcMesh.tag_set_data(th_prin_stress_vect2, &mapGaussPts[gg],
701  1, &*prin_stress_vect2.begin());
702  CHKERR postProcMesh.tag_set_data(th_prin_stress_vect3, &mapGaussPts[gg],
703  1, &*prin_stress_vect3.begin());
704  }
705 
706  ++S;
707  ++P;
708  ++psi;
709  ++R;
710  ++rho;
711  ++grad_rho;
712  }
713 
715 }
716 
718  Remodeling::CommonData &common_data)
720  "DISPLACEMENTS", UserDataOperator::OPROW),
721  commonData(common_data) {
722  // Set that this operator is executed only for vertices on element
723  doVertices = true;
724  doEdges = false;
725  doQuads = false;
726  doTris = false;
727  doTets = false;
728  doPrisms = false;
729 }
730 
734 
736 
737  if (type != MBVERTEX)
739 
740  auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
741  const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
742 
743  commonData.data.matInvF.resize(9, nb_gauss_pts, false);
744  auto invF = getFTensor2FromMat<3, 3>(commonData.data.matInvF);
745  commonData.data.vecDetF.resize(nb_gauss_pts, false);
747 
749  dP_dF.resize(81, nb_gauss_pts, false);
751  matF.resize(9, nb_gauss_pts, false);
754  commonData.data.matP.resize(9, nb_gauss_pts, false);
755  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
756  commonData.data.vecPsi.resize(nb_gauss_pts, false);
764  const double lambda = commonData.lambda;
765  const double mu = commonData.mu;
766 
767  for (int gg = 0; gg != nb_gauss_pts; gg++) {
768 
769  // Get deformation gradient
770  F(i, I) = grad(i, I);
771  F(0, 0) += 1;
772  F(1, 1) += 1;
773  F(2, 2) += 1;
774 
775  CHKERR determinantTensor3by3(F, det_f);
776  CHKERR invertTensor3by3(F, det_f, invF);
777 
778  const double log_det = log(det_f);
779  psi = F(i, J) * F(i, J) - 3 - 2 * log_det;
780  psi *= 0.5 * mu;
781  psi += 0.5 * lambda * log_det * log_det;
782 
783  const double var = lambda * log_det - mu;
784  P(i, J) = mu * F(i, J) + var * invF(J, i);
785 
786  const double coef = mu - lambda * log_det;
787  // TODO: This can be improved by utilising the symmetries of the D2 tensor
788 
789  D2(i, J, k, L) =
790  invF(J, i) * (invF(L, k) * lambda) + invF(L, i) * (invF(J, k) * coef);
791 
792  D2(0, 0, 0, 0) += mu;
793  D2(0, 1, 0, 1) += mu;
794  D2(0, 2, 0, 2) += mu;
795  D2(1, 0, 1, 0) += mu;
796  D2(1, 1, 1, 1) += mu;
797  D2(1, 2, 1, 2) += mu;
798  D2(2, 0, 2, 0) += mu;
799  D2(2, 1, 2, 1) += mu;
800  D2(2, 2, 2, 2) += mu;
801 
802  ++det_f;
803  ++invF;
804  ++grad;
805  ++F;
806  ++D2;
807  ++P;
808  ++psi;
809  }
810 
812 }
813 
816  "DISPLACEMENTS", UserDataOperator::OPROW),
817  commonData(common_data) {}
818 
820 OpAssmbleStressRhs::doWork(int side, EntityType type,
822 
824 
825  const int nb_dofs = data.getIndices().size();
826  if (!nb_dofs)
828  nF.resize(nb_dofs, false);
829  nF.clear();
830 
831  const int nb_gauss_pts = data.getN().size1();
832  const int nb_base_functions = data.getN().size2();
833  auto diff_base_functions = data.getFTensor1DiffN<3>();
834  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
836 
837  // MatrixDouble &matS = commonData.data.matS;
838  // MatrixDouble &matF = commonData.data.matGradientOfDeformation;
839  // FTensor::Tensor2<double*,3,3> F(TENSOR2_MAT_PTR(matF));
840  // FTensor::Tensor2_symmetric<double*,3> S(SYMMETRIC_TENSOR2_MAT_PTR(matS),6);
841 
842  const double n = commonData.n;
843  const double rho_ref = commonData.rHo_ref;
844 
848 
849  for (int gg = 0; gg != nb_gauss_pts; gg++) {
850 
851  // Get volume and integration weight
852  double w = getVolume() * getGaussPts()(3, gg);
853  if (getHoGaussPtsDetJac().size() > 0) {
854  w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
855  }
856  double a = w * pow(rho / rho_ref, n);
857 
858  int bb = 0;
859  for (; bb != nb_dofs / 3; bb++) {
860  double *ptr = &nF[3 * bb];
861  FTensor::Tensor1<double *, 3> t1(ptr, &ptr[1], &ptr[2]);
862  t1(i) += a * P(i, I) * diff_base_functions(I);
863  // t1(i) += 0.5*a*diff_base_functions(J)*(F(i,I)*S(I,J));
864  // t1(i) += 0.5*a*diff_base_functions(I)*(F(i,J)*S(I,J));
865  ++diff_base_functions;
866  }
867  // Could be that we base more base functions than needed to approx. disp.
868  // field.
869  for (; bb != nb_base_functions; bb++) {
870  ++diff_base_functions;
871  }
872  ++P;
873  // ++S;
874  // ++F;
875  ++rho;
876  }
877 
878  // Assemble internal force vector for time solver (TS)
880  getFEMethod()->ts_F, /// Get the right hand side vector for time solver
881  nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
882 
884 }
885 
888  "RHO", UserDataOperator::OPROW),
889  commonData(common_data) {}
890 
892 OpAssmbleRhoRhs::doWork(int side, EntityType type,
894 
896 
897  const int nb_dofs = data.getIndices().size();
898  if (!nb_dofs)
900  nF.resize(nb_dofs, false);
901  nF.clear();
902 
903  const int nb_gauss_pts = data.getN().size1();
904  const int nb_base_functions = data.getN().size2();
905  auto base_functions = data.getFTensor0N();
906  auto diff_base_functions = data.getFTensor1DiffN<3>();
907  if (!commonData.data.rhoPtr) {
908  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "rhoPtsr is null");
909  }
911  if (!commonData.data.gradRhoPtr) {
912  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "gradRhoPtr is null");
913  }
914  auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
916 
917  if (commonData.data.vecRhoDt.size() != nb_gauss_pts) {
918  commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
919  commonData.data.vecRhoDt.clear();
920  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Time derivative of
921  // density at integration not calculated");
922  }
923  auto drho_dt = getFTensor0FromVec(commonData.data.vecRhoDt);
924 
925  const double R0 = commonData.R0;
926 
928 
929  for (int gg = 0; gg != nb_gauss_pts; gg++) {
930 
931  // Get volume and integration weight
932  double w = getVolume() * getGaussPts()(3, gg);
933  if (getHoGaussPtsDetJac().size() > 0) {
934  w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
935  }
936 
937  // cerr << R << endl;
938 
939  int bb = 0;
940  for (; bb != nb_dofs; bb++) {
941  nF[bb] += w * base_functions * drho_dt; // Density rate term
942  nF[bb] -= w * base_functions * R; // Source term
943  nF[bb] -= w * R0 * diff_base_functions(I) * grad_rho(I); // Gradient term
944  ++base_functions;
945  ++diff_base_functions;
946  }
947  // Could be that we base more base functions than needed to approx. density
948  // field.
949  for (; bb != nb_base_functions; bb++) {
950  ++base_functions;
951  ++diff_base_functions;
952  }
953 
954  // Increment quantities for integration pts.
955  ++rho;
956  ++grad_rho;
957  ++drho_dt;
958  ++R;
959  }
960 
961  // Assemble internal force vector for time solver (TS)
963  getFEMethod()->ts_F, /// Get the right hand side vector for time solver
964  nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
965 
967 }
968 
971  "RHO", "RHO", UserDataOperator::OPROWCOL),
972  commonData(common_data) {
973  sYmm = true;
974 }
975 
977 OpAssmbleRhoLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
978  EntityType col_type,
981 
983 
984  const int row_nb_dofs = row_data.getIndices().size();
985  if (!row_nb_dofs)
987  const int col_nb_dofs = col_data.getIndices().size();
988  if (!col_nb_dofs)
990 
991  // Set size can clear local tangent matrix
992  locK_rho_rho.resize(row_nb_dofs, col_nb_dofs, false);
993  locK_rho_rho.clear();
994 
995  const int row_nb_gauss_pts = row_data.getN().size1();
996  if (!row_nb_gauss_pts)
998  const int row_nb_base_functions = row_data.getN().size2();
999 
1000  const double ts_a = getFEMethod()->ts_a;
1001  const double R0 = commonData.R0;
1002  const double c = commonData.c;
1003  const double n = commonData.n;
1004  const double m = commonData.m;
1005  const double rho_ref = commonData.rHo_ref;
1006  const double psi_ref = commonData.pSi_ref;
1007 
1009 
1012  auto row_base_functions = row_data.getFTensor0N();
1013  auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1014 
1015  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1016 
1017  // Get volume and integration weight
1018  double w = getVolume() * getGaussPts()(3, gg);
1019  if (getHoGaussPtsDetJac().size() > 0) {
1020  w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
1021  }
1022 
1023  const double kp = commonData.getCFromDensity(rho);
1024  const double kp_diff = commonData.getCFromDensityDiff(rho);
1025  const double a = c * kp * ((n - m) / rho) * pow(rho / rho_ref, n - m) * psi;
1026  const double a_diff =
1027  c * kp_diff * (pow(rho / rho_ref, n - m) * psi - psi_ref);
1028 
1029  int row_bb = 0;
1030  for (; row_bb != row_nb_dofs; row_bb++) {
1031  auto col_base_functions = col_data.getFTensor0N(gg, 0);
1032  auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1033  for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1034  locK_rho_rho(row_bb, col_bb) +=
1035  ts_a * w * row_base_functions * col_base_functions;
1036  locK_rho_rho(row_bb, col_bb) -=
1037  w * R0 * row_diff_base_functions(I) * col_diff_base_functions(I);
1038  locK_rho_rho(row_bb, col_bb) -=
1039  a * w * row_base_functions * col_base_functions;
1040  locK_rho_rho(row_bb, col_bb) -=
1041  a_diff * w * row_base_functions * col_base_functions;
1042 
1043  ++col_base_functions;
1044  ++col_diff_base_functions;
1045  }
1046  ++row_base_functions;
1047  ++row_diff_base_functions;
1048  }
1049  for (; row_bb != row_nb_base_functions; row_bb++) {
1050  ++row_base_functions;
1051  ++row_diff_base_functions;
1052  }
1053  ++psi;
1054  ++rho;
1055  }
1056 
1057  // cerr << locK_rho_rho << endl;
1058 
1059  CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1060  &*row_data.getIndices().begin(), col_nb_dofs,
1061  &*col_data.getIndices().begin(), &locK_rho_rho(0, 0),
1062  ADD_VALUES);
1063 
1064  // is symmetric
1065  if (row_side != col_side || row_type != col_type) {
1066  transLocK_rho_rho.resize(col_nb_dofs, row_nb_dofs, false);
1067  noalias(transLocK_rho_rho) = trans(locK_rho_rho);
1068  CHKERR MatSetValues(getFEMethod()->ts_B, col_nb_dofs,
1069  &*col_data.getIndices().begin(), row_nb_dofs,
1070  &*row_data.getIndices().begin(),
1071  &transLocK_rho_rho(0, 0), ADD_VALUES);
1072  }
1073 
1075 }
1076 
1079  "RHO", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
1080  commonData(common_data) {
1081  sYmm = false;
1082 }
1083 
1085 OpAssmbleRhoLhs_dF::doWork(int row_side, int col_side, EntityType row_type,
1086  EntityType col_type,
1089 
1091 
1092  const int row_nb_dofs = row_data.getIndices().size();
1093  if (!row_nb_dofs)
1095  const int col_nb_dofs = col_data.getIndices().size();
1096  if (!col_nb_dofs)
1098 
1099  // Set size can clear local tangent matrix
1100  locK_rho_F.resize(row_nb_dofs, col_nb_dofs, false);
1101  locK_rho_F.clear();
1102 
1103  const int row_nb_gauss_pts = row_data.getN().size1();
1104  const int col_nb_base_functions = col_data.getN().size2();
1105 
1106  const double c = commonData.c;
1107  const double n = commonData.n;
1108  const double m = commonData.m;
1109  const double rho_ref = commonData.rHo_ref;
1110 
1113 
1115  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
1116 
1118  auto col_diff_base_functions = col_data.getFTensor1DiffN<3>();
1119 
1120  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1121 
1122  // Get volume and integration weight
1123  double w = getVolume() * getGaussPts()(3, gg);
1124  if (getHoGaussPtsDetJac().size() > 0) {
1125  w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
1126  }
1127 
1128  const double kp = commonData.getCFromDensity(rho);
1129  const double a = w * c * kp * pow(rho / rho_ref, n - m);
1130 
1131  int col_bb = 0;
1132  for (; col_bb != col_nb_dofs / 3; col_bb++) {
1133  t1(i) = a * P(i, I) * col_diff_base_functions(I);
1134  auto row_base_functions = row_data.getFTensor0N(gg, 0);
1135  for (int row_bb = 0; row_bb != row_nb_dofs; row_bb++) {
1136  FTensor::Tensor1<double *, 3> k(&locK_rho_F(row_bb, 3 * col_bb + 0),
1137  &locK_rho_F(row_bb, 3 * col_bb + 1),
1138  &locK_rho_F(row_bb, 3 * col_bb + 2));
1139  k(i) -= row_base_functions * t1(i);
1140  ++row_base_functions;
1141  }
1142  ++col_diff_base_functions;
1143  }
1144  for (; col_bb != col_nb_base_functions; col_bb++) {
1145  ++col_diff_base_functions;
1146  }
1147 
1148  ++P;
1149  ++rho;
1150  }
1151 
1152  CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1153  &*row_data.getIndices().begin(), col_nb_dofs,
1154  &*col_data.getIndices().begin(), &locK_rho_F(0, 0),
1155  ADD_VALUES);
1156 
1158 }
1159 template <bool ADOLC>
1161  Remodeling::CommonData &common_data)
1163  "DISPLACEMENTS", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
1164  commonData(common_data) {
1165  sYmm = true;
1166 }
1167 template <bool ADOLC>
1169  int row_side, int col_side, EntityType row_type, EntityType col_type,
1172 
1174 
1175  const int row_nb_dofs = row_data.getIndices().size();
1176  if (!row_nb_dofs)
1178  const int col_nb_dofs = col_data.getIndices().size();
1179  if (!col_nb_dofs)
1181 
1182  const bool diagonal_block = (row_type == col_type) && (row_side == col_side);
1183 
1184  // Set size can clear local tangent matrix
1185  locK_P_F.resize(row_nb_dofs, col_nb_dofs, false);
1186  locK_P_F.clear();
1187 
1188  const int row_nb_gauss_pts = row_data.getN().size1();
1189  const int row_nb_base_functions = row_data.getN().size2();
1190 
1191  MatrixDouble &matS = commonData.data.matS;
1192  MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
1193  // MatrixDouble &dS_dC = commonData.data.matMaterialTangent;
1194  // MatrixDouble &matF = commonData.data.matGradientOfDeformation;
1195 
1196  double t0;
1197 
1198  const double n = commonData.n;
1199  const double rho_ref = commonData.rHo_ref;
1200 
1208 
1209  auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1212  auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1213 
1214  // FTensor::Tensor4_ddg<double*,3,3> D1(SYMMETRIC_TENSOR4_MAT_PTR(dS_dC),6);
1215  // FTensor::Tensor2<double*,3,3> F(TENSOR2_MAT_PTR(matF));
1216 
1217  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1218 
1219  // Get volume and integration weight
1220  double w = getVolume() * getGaussPts()(3, gg);
1221  if (getHoGaussPtsDetJac().size() > 0) {
1222  w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
1223  }
1224  const double a = w * pow(rho / rho_ref, n);
1225 
1226  int row_bb = 0;
1227  for (; row_bb != row_nb_dofs / 3; row_bb++) {
1228 
1229  auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1230 
1231  const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
1232  int col_bb = 0;
1233  for (; col_bb != final_bb; col_bb++) {
1234 
1236  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1237  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1238  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1239  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1240  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1241  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1242  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1243  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1244  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1245 
1246  diffDiff(J, L) =
1247  a * row_diff_base_functions(J) * col_diff_base_functions(L);
1248 
1249  // TODO: FIXME: Tensor D2 has major symmetry, we do not have such tensor
1250  // implemented at the moment. Using tensor with major symmetry will
1251  // improve code efficiency.
1252  t_assemble(i, k) += diffDiff(J, L) * D2(i, J, k, L);
1253 
1254  if (ADOLC) {
1255  // Stress stiffening part (only with adolc since it is using dS / dC
1256  // derrivative. Check OpCalculateStressTangentWithAdolc operator)
1257  t0 = diffDiff(J, L) * S(J, L);
1258  t_assemble(0, 0) += t0;
1259  t_assemble(1, 1) += t0;
1260  t_assemble(2, 2) += t0;
1261  }
1262 
1263  // Next base function for column
1264  ++col_diff_base_functions;
1265  }
1266 
1267  ++row_diff_base_functions;
1268  }
1269  for (; row_bb != row_nb_base_functions; row_bb++) {
1270  ++row_diff_base_functions;
1271  }
1272 
1273  // ++D1;
1274  ++D2;
1275  ++S;
1276  // ++F;
1277  ++rho;
1278  }
1279 
1280  // Copy of symmetric part for assemble
1281  if (diagonal_block) {
1282  for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
1283  int col_bb = 0;
1284  for (; col_bb != row_bb + 1; col_bb++) {
1286  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1287  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1288  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1289  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1290  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1291  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1292  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1293  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1294  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1296  &locK_P_F(3 * col_bb + 0, 3 * row_bb + 0),
1297  &locK_P_F(3 * col_bb + 0, 3 * row_bb + 1),
1298  &locK_P_F(3 * col_bb + 0, 3 * row_bb + 2),
1299  &locK_P_F(3 * col_bb + 1, 3 * row_bb + 0),
1300  &locK_P_F(3 * col_bb + 1, 3 * row_bb + 1),
1301  &locK_P_F(3 * col_bb + 1, 3 * row_bb + 2),
1302  &locK_P_F(3 * col_bb + 2, 3 * row_bb + 0),
1303  &locK_P_F(3 * col_bb + 2, 3 * row_bb + 1),
1304  &locK_P_F(3 * col_bb + 2, 3 * row_bb + 2));
1305  t_off_side(i, k) = t_assemble(k, i);
1306  }
1307  }
1308  }
1309 
1310  const int *row_ind = &*row_data.getIndices().begin();
1311  const int *col_ind = &*col_data.getIndices().begin();
1312  CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs, row_ind, col_nb_dofs,
1313  col_ind, &locK_P_F(0, 0), ADD_VALUES);
1314 
1315  if (row_type != col_type || row_side != col_side) {
1316  transLocK_P_F.resize(col_nb_dofs, row_nb_dofs, false);
1317  noalias(transLocK_P_F) = trans(locK_P_F);
1318  CHKERR MatSetValues(getFEMethod()->ts_B, col_nb_dofs, col_ind, row_nb_dofs,
1319  row_ind, &transLocK_P_F(0, 0), ADD_VALUES);
1320  }
1321 
1323 }
1324 
1326  Remodeling::CommonData &common_data)
1328  "DISPLACEMENTS", "RHO", UserDataOperator::OPROWCOL),
1329  commonData(common_data) {
1330  sYmm = false; // This is off-diagonal (upper-left), so no symmetric
1331 }
1332 
1334 OpAssmbleStressLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
1335  EntityType col_type,
1338 
1340 
1341  const int row_nb_dofs = row_data.getIndices().size();
1342  if (!row_nb_dofs)
1344  const int col_nb_dofs = col_data.getIndices().size();
1345  if (!col_nb_dofs)
1347 
1348  // Set size can clear local tangent matrix
1349  locK_P_RHO.resize(row_nb_dofs, col_nb_dofs, false);
1350  locK_P_RHO.clear();
1351 
1353 
1354  const double n = commonData.n;
1355  const double rho_ref = commonData.rHo_ref;
1356 
1359 
1360  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
1362 
1363  const int row_nb_gauss_pts = row_data.getN().size1();
1364  const int row_nb_base_functions = row_data.getN().size2();
1365  auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1366 
1367  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1368 
1369  // Get volume and integration weight
1370  double w = getVolume() * getGaussPts()(3, gg);
1371  if (getHoGaussPtsDetJac().size() > 0) {
1372  w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
1373  }
1374  double a = w * (n / rho) * pow(rho / rho_ref, n);
1375 
1376  int row_bb = 0;
1377  for (; row_bb != row_nb_dofs / 3; row_bb++) {
1378  t1(i) = a * row_diff_base_functions(I) * P(i, I);
1379  auto base_functions = col_data.getFTensor0N(gg, 0);
1380  for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1381  // if(base_functions!=col_data.getN(gg)[col_bb]) {
1382  // cerr << "Error" << endl;
1383  // }
1384  FTensor::Tensor1<double *, 3> k(&locK_P_RHO(3 * row_bb + 0, col_bb),
1385  &locK_P_RHO(3 * row_bb + 1, col_bb),
1386  &locK_P_RHO(3 * row_bb + 2, col_bb));
1387  k(i) += t1(i) * base_functions;
1388  ++base_functions;
1389  }
1390  ++row_diff_base_functions;
1391  }
1392  for (; row_bb != row_nb_base_functions; row_bb++) {
1393  ++row_diff_base_functions;
1394  }
1395 
1396  ++P;
1397  ++rho;
1398  }
1399 
1400  CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1401  &*row_data.getIndices().begin(), col_nb_dofs,
1402  &*col_data.getIndices().begin(), &locK_P_RHO(0, 0),
1403  ADD_VALUES);
1404 
1406 }
1407 
1409  Remodeling::CommonData &common_data)
1410  : mField(m_field), commonData(common_data), postProc(m_field),
1411  postProcElastic(m_field),
1412  // densityMaps(m_field),
1413  iNit(false) {}
1414 
1418 }
1419 
1423 }
1424 
1427 
1428  PetscBool mass_postproc = PETSC_FALSE;
1429  PetscBool equilibrium_flg = PETSC_FALSE;
1430  PetscBool analysis_mesh_flg = PETSC_FALSE;
1431  rate = 1;
1432 
1433  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling post-process",
1434  "none");
1435 
1436  CHKERR PetscOptionsBool(
1437  "-mass_postproc",
1438  "is used for testing, calculates mass and energy at each time step", "",
1439  mass_postproc, &mass_postproc, PETSC_NULL);
1440 
1441  CHKERR PetscOptionsBool("-analysis_mesh",
1442  "saves analysis mesh at each time step", "",
1443  analysis_mesh_flg, &analysis_mesh_flg, PETSC_NULL);
1444 
1445  CHKERR PetscOptionsScalar(
1446  "-equilibrium_stop_rate",
1447  "is used to stop calculations when equilibium state is achieved", "",
1448  rate, &rate, &equilibrium_flg);
1449  ierr = PetscOptionsEnd();
1450  CHKERRQ(ierr);
1451 
1452  if (!iNit) {
1453 
1454  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "",
1455  "Bone remodeling post-process", "none");
1456 
1457  pRT = 1;
1458  CHKERR PetscOptionsInt("-my_output_prt",
1459  "frequncy how often results are dumped on hard disk",
1460  "", 1, &pRT, NULL);
1461  ierr = PetscOptionsEnd();
1462  CHKERRQ(ierr);
1463 
1465  CHKERR postProc.addFieldValuesPostProc("DISPLACEMENTS");
1466  CHKERR postProc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
1467  if (!commonData.less_post_proc) {
1469  }
1470  // CHKERR postProc.addFieldValuesPostProc("RHO");
1471  // CHKERR postProc.addFieldValuesGradientPostProc("RHO");
1472 
1473  {
1474  boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
1475  &list_of_operators = postProc.getOpPtrVector();
1476 
1477  list_of_operators.push_back(
1479  list_of_operators.push_back(new OpCalculateScalarFieldGradient<3>(
1480  "RHO", commonData.data.gradRhoPtr));
1481  // Get displacement gradient at integration points
1482  list_of_operators.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1483  "DISPLACEMENTS", commonData.data.gradDispPtr));
1484  list_of_operators.push_back(new OpCalculateStress(commonData));
1485  list_of_operators.push_back(new OpPostProcStress(
1487  }
1488 
1491  CHKERR postProcElastic.addFieldValuesPostProc("MESH_NODE_POSITIONS");
1493  std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
1494  commonData.elasticPtr->setOfBlocks.begin();
1495  for (; sit != commonData.elasticPtr->setOfBlocks.end(); sit++) {
1498  "DISPLACEMENTS", sit->second, postProcElastic.commonData));
1499  }
1500 
1501  iNit = true;
1502  }
1503 
1504  if (mass_postproc) {
1505  Vec mass_vec;
1506  Vec energ_vec;
1507  int mass_vec_ghost[] = {0};
1508  CHKERR VecCreateGhost(PETSC_COMM_WORLD, (!mField.get_comm_rank()) ? 1 : 0,
1509  1, 1, mass_vec_ghost, &mass_vec);
1510 
1511  CHKERR VecZeroEntries(mass_vec);
1512  CHKERR VecDuplicate(mass_vec, &energ_vec);
1513 
1514  Remodeling::Fe energyCalc(mField);
1515 
1516  energyCalc.getOpPtrVector().push_back(
1518  energyCalc.getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<3>(
1519  "RHO", commonData.data.gradRhoPtr));
1520  energyCalc.getOpPtrVector().push_back(
1521  new OpCalculateVectorFieldGradient<3, 3>("DISPLACEMENTS",
1523  energyCalc.getOpPtrVector().push_back(new OpCalculateStress(commonData));
1524  energyCalc.getOpPtrVector().push_back(
1525  new OpMassAndEnergyCalculation("RHO", commonData, energ_vec, mass_vec));
1526  CHKERR DMoFEMLoopFiniteElements(commonData.dm, "FE_REMODELLING",
1527  &energyCalc);
1528 
1529  CHKERR VecAssemblyBegin(mass_vec);
1530  CHKERR VecAssemblyEnd(mass_vec);
1531  double mass_sum;
1532  CHKERR VecSum(mass_vec, &mass_sum);
1533 
1534  CHKERR VecAssemblyBegin(energ_vec);
1535  CHKERR VecAssemblyEnd(energ_vec);
1536  double energ_sum;
1537  CHKERR VecSum(energ_vec, &energ_sum);
1538 
1539  CHKERR PetscPrintf(PETSC_COMM_WORLD,
1540  "Time: %g Mass: %6.5g Elastic energy: %6.5g \n", ts_t,
1541  mass_sum, energ_sum);
1542  CHKERR VecDestroy(&mass_vec);
1543  CHKERR VecDestroy(&energ_vec);
1544 
1545  if (equilibrium_flg) {
1546  double equilibrium_rate =
1547  fabs(energ_sum - commonData.cUrrent_psi) / energ_sum;
1548  double equilibrium_mass_rate =
1549  fabs(mass_sum - commonData.cUrrent_mass) / mass_sum;
1550  if (equilibrium_rate < rate) {
1551  CHKERR PetscPrintf(
1552  PETSC_COMM_WORLD,
1553  "Energy equilibrium state is achieved! Difference = %0.6g %%. \n",
1554  equilibrium_rate * 100);
1555  }
1556  if (equilibrium_mass_rate < rate) {
1557  CHKERR PetscPrintf(
1558  PETSC_COMM_WORLD,
1559  "Mass equilibrium state is achieved! Difference = %0.6g %%. \n",
1560  equilibrium_mass_rate * 100);
1561  }
1562  commonData.cUrrent_psi = energ_sum;
1563  commonData.cUrrent_mass = mass_sum;
1564  }
1565  }
1566 
1567  int step;
1568 #if PETSC_VERSION_LE(3, 8, 0)
1569  CHKERR TSGetTimeStepNumber(ts, &step);
1570 #else
1571  CHKERR TSGetStepNumber(ts, &step);
1572 #endif
1573 
1574  if ((step) % pRT == 0) {
1575 
1576  ostringstream sss;
1577  sss << "out_" << step << ".h5m";
1578  CHKERR DMoFEMLoopFiniteElements(commonData.dm, "FE_REMODELLING", &postProc);
1579  CHKERR postProc.postProcMesh.write_file(sss.str().c_str(), "MOAB",
1580  "PARALLEL=WRITE_PART");
1581  if (analysis_mesh_flg) {
1582  ostringstream ttt;
1583  ttt << "analysis_mesh_" << step << ".h5m";
1584  CHKERR mField.get_moab().write_file(ttt.str().c_str(), "MOAB",
1585  "PARALLEL=WRITE_PART");
1586  }
1587 
1590  &postProcElastic);
1591  ostringstream ss;
1592  ss << "out_elastic_" << step << ".h5m";
1593  CHKERR postProcElastic.postProcMesh.write_file(ss.str().c_str(), "MOAB",
1594  "PARALLEL=WRITE_PART");
1595  }
1596  }
1598 }
1599 
1601 
1603  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling", "none");
1604 
1605  commonData.oRder = 2;
1606  CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
1607  &commonData.oRder, PETSC_NULL);
1608 
1609  ierr = PetscOptionsEnd();
1610  CHKERRQ(ierr);
1611 
1612  CHKERR mField.get_moab().get_entities_by_type(0, MBTET, commonData.tEts_all);
1614  string name = it->getName();
1615  if (name.compare(0, 14, "NO_REMODELLING") == 0) {
1617  EntityHandle meshset = it->getMeshset();
1618  CHKERR this->mField.get_moab().get_entities_by_type(
1619  meshset, MBTET, commonData.tEts_block, true);
1622  }
1623  }
1624 
1626 }
1627 
1629 
1631 
1632  // Seed all mesh entities to MoFEM database, those entities can be potentially
1633  // used as finite elements or as entities which carry some approximation
1634  // field.
1635  commonData.bitLevel.set(0);
1636  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
1637  0, 3, commonData.bitLevel);
1638  int order = commonData.oRder;
1639 
1640  // Add displacement field
1641  CHKERR mField.add_field("DISPLACEMENTS", H1,
1642  /*AINSWORTH_LOBATTO_BASE*/ AINSWORTH_LEGENDRE_BASE, 3,
1643  MB_TAG_SPARSE, MF_ZERO);
1644  // Add field representing ho-geometry
1645  CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3,
1646  MB_TAG_SPARSE, MF_ZERO);
1647 
1648  // Check if density is available, if not add density field.
1649  bool add_rho_field = false;
1650  if (!mField.check_field("RHO")) {
1652  MB_TAG_SPARSE, MF_ZERO);
1653  // FIXME
1655  // CHKERR mField.add_ents_to_field_by_type(0,MBTET,"RHO");
1657  add_rho_field = true;
1658 
1659  CHKERR mField.set_field_order(0, MBVERTEX, "RHO", 1);
1660  CHKERR mField.set_field_order(0, MBEDGE, "RHO", order - 1);
1661  CHKERR mField.set_field_order(0, MBTRI, "RHO", order - 1);
1662  CHKERR mField.set_field_order(0, MBTET, "RHO", order - 1);
1663  }
1664 
1665  // Add entities to field
1666  CHKERR mField.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENTS");
1667  CHKERR mField.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
1668 
1669  // Set approximation order to entities
1670  CHKERR mField.set_field_order(0, MBVERTEX, "DISPLACEMENTS", 1);
1671  CHKERR mField.set_field_order(0, MBEDGE, "DISPLACEMENTS", order);
1672  CHKERR mField.set_field_order(0, MBTRI, "DISPLACEMENTS", order);
1673  CHKERR mField.set_field_order(0, MBTET, "DISPLACEMENTS", order);
1674 
1675  // Assumes that geometry is approximated using 2nd order polynomials.
1676 
1677  // Apply 2nd order only on skin
1678  {
1679  // Skinner skin(&mField.get_moab());
1680  // Range faces,tets;
1681  // CHKERR mField.get_moab().get_entities_by_type(0,MBTET,tets);
1682  // CHKERR skin.find_skin(0,tets,false,faces);
1683  // Range edges;
1684  // CHKERR mField.get_moab().get_adjacencies(
1685  // faces,1,false,edges,moab::Interface::UNION
1686  // );
1687  CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
1688  }
1689 
1690  CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
1691 
1692  // Build fields
1694 
1695  // If order was not set from CT scans set homogenous order equal to
1696  // reference bone density
1697  if (add_rho_field) {
1699  MBVERTEX, "RHO");
1700  // const DofEntity_multiIndex *dofs_ptr;
1701  // CHKERR mField.get_dofs(&dofs_ptr);
1702  // for(_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(mField,"RHO",dit)) {
1703  // cerr << (*dit)->getFieldData() << endl;
1704  // }
1705  }
1706 
1707  // Project geometry given on 10-node tets on ho-geometry
1708  Projection10NodeCoordsOnField ent_method_material(mField,
1709  "MESH_NODE_POSITIONS");
1710  CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
1711 
1713 }
1714 
1717 
1718  CHKERR mField.add_finite_element("FE_REMODELLING", MF_ZERO);
1720  "DISPLACEMENTS");
1722  "DISPLACEMENTS");
1724  "DISPLACEMENTS");
1725  CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING", "RHO");
1727  "MESH_NODE_POSITIONS");
1728  CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING", "RHO");
1729  CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING", "RHO");
1731  "FE_REMODELLING");
1732 
1733  CHKERR mField.add_finite_element("ELASTIC");
1734  CHKERR mField.modify_finite_element_add_field_row("ELASTIC", "DISPLACEMENTS");
1735  CHKERR mField.modify_finite_element_add_field_col("ELASTIC", "DISPLACEMENTS");
1737  "DISPLACEMENTS");
1739  "MESH_NODE_POSITIONS");
1740 
1743  MBTET, "ELASTIC");
1744  }
1745 
1747  boost::shared_ptr<ElasticMaterials>(new ElasticMaterials(mField));
1748  commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
1749  new NonlinearElasticElement(mField, 10));
1751  commonData.elasticPtr->setOfBlocks);
1752  CHKERR commonData.elasticPtr->addElement("ELASTIC", "DISPLACEMENTS");
1753  CHKERR commonData.elasticPtr->setOperators(
1754  "DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
1755 
1758 
1759  // Allocate memory for density and gradient of displacements at integration
1760  // points
1761  commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
1763  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1765  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1766 
1767  // Add operators Rhs
1768  boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
1769  &list_of_operators_rhs = commonData.feRhs->getOpPtrVector();
1770 
1771  // Get density at integration points
1772  list_of_operators_rhs.push_back(
1774  list_of_operators_rhs.push_back(
1776  list_of_operators_rhs.push_back(new OpGetRhoTimeDirevative(commonData));
1777  // Get displacement gradient at integration points
1778  list_of_operators_rhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1779  "DISPLACEMENTS", commonData.data.gradDispPtr));
1780  list_of_operators_rhs.push_back(new OpCalculateStress(commonData));
1781  list_of_operators_rhs.push_back(new OpAssmbleStressRhs(commonData));
1782  list_of_operators_rhs.push_back(new OpAssmbleRhoRhs(commonData));
1783 
1784  // Add operators Lhs
1785  boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
1786  &list_of_operators_lhs = commonData.feLhs->getOpPtrVector();
1787  // Get density at integration points
1788  list_of_operators_lhs.push_back(
1790  list_of_operators_rhs.push_back(
1792  // Get displacement gradient at integration points
1793  list_of_operators_lhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1794  "DISPLACEMENTS", commonData.data.gradDispPtr));
1795  if (commonData.with_adol_c) {
1796  list_of_operators_lhs.push_back(
1798  list_of_operators_lhs.push_back(
1800  } else {
1801  list_of_operators_lhs.push_back(new OpCalculateStressTangent(commonData));
1802  list_of_operators_lhs.push_back(
1804  }
1805  list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dRho(commonData));
1806  list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dF(commonData));
1807  list_of_operators_lhs.push_back(new OpAssmbleStressLhs_dRho(commonData));
1808 
1810 }
1811 
1813 
1815  // Add Neumann forces elements
1816  CHKERR MetaNeummanForces::addNeumannBCElements(mField, "DISPLACEMENTS");
1817  CHKERR MetaNodalForces::addElement(mField, "DISPLACEMENTS");
1818  CHKERR MetaEdgeForces::addElement(mField, "DISPLACEMENTS");
1819 
1820  // forces and pressures on surface
1821  CHKERR MetaNeummanForces::setMomentumFluxOperators(
1822  mField, commonData.neumannForces, PETSC_NULL, "DISPLACEMENTS");
1823  // noadl forces
1825  PETSC_NULL, "DISPLACEMENTS");
1826  // edge forces
1828  "DISPLACEMENTS");
1829 
1830  for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
1831  commonData.neumannForces.begin();
1832  mit != commonData.neumannForces.end(); mit++) {
1833  mit->second->methodsOp.push_back(
1834  new TimeForceScale("-my_load_history", false));
1835  }
1836  for (boost::ptr_map<string, NodalForce>::iterator mit =
1837  commonData.nodalForces.begin();
1838  mit != commonData.nodalForces.end(); mit++) {
1839  mit->second->methodsOp.push_back(
1840  new TimeForceScale("-my_load_history", false));
1841  }
1842  for (boost::ptr_map<string, EdgeForce>::iterator mit =
1843  commonData.edgeForces.begin();
1844  mit != commonData.edgeForces.end(); mit++) {
1845  mit->second->methodsOp.push_back(
1846  new TimeForceScale("-my_load_history", false));
1847  }
1848 
1850 }
1851 
1853 
1855  commonData.dm_name = "DM_REMODELING";
1856  // Register DM problem
1858  CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
1859  CHKERR DMSetType(commonData.dm, commonData.dm_name);
1861  // Create DM instance
1864  // Configure DM form line command options (DM itself, solvers,
1865  // pre-conditioners, ... )
1866  CHKERR DMSetFromOptions(commonData.dm);
1867  // Add elements to dm (only one here)
1868  CHKERR DMMoFEMAddElement(commonData.dm, "FE_REMODELLING");
1869  CHKERR DMMoFEMAddElement(commonData.dm, "ELASTIC");
1870 
1871  if (mField.check_finite_element("FORCE_FE")) {
1872  CHKERR DMMoFEMAddElement(commonData.dm, "FORCE_FE");
1873  }
1874  if (mField.check_finite_element("PRESSURE_FE")) {
1875  CHKERR DMMoFEMAddElement(commonData.dm, "PRESSURE_FE");
1876  }
1877  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
1878  CHKERR DMSetUp(commonData.dm);
1879 
1881 }
1882 
1884 
1886 
1887  Vec F = commonData.F;
1888  Vec D = commonData.D;
1889  Mat A = commonData.A;
1890  DM dm = commonData.dm;
1891  TS ts = commonData.ts;
1892 
1893  CHKERR TSSetIFunction(ts, F, PETSC_NULL, PETSC_NULL);
1894  CHKERR TSSetIJacobian(ts, A, A, PETSC_NULL, PETSC_NULL);
1895  double ftime = 1;
1896 #if PETSC_VERSION_GE(3, 8, 0)
1897  CHKERR TSSetMaxTime(ts, ftime);
1898 #endif
1899  CHKERR TSSetFromOptions(ts);
1900  CHKERR TSSetDM(ts, dm);
1901 
1902  SNES snes;
1903  CHKERR TSGetSNES(ts, &snes);
1904  KSP ksp;
1905  CHKERR SNESGetKSP(snes, &ksp);
1906  PC pc;
1907  CHKERR KSPGetPC(ksp, &pc);
1908  PetscBool is_pcfs = PETSC_FALSE;
1909  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1910 
1911  // Set up FIELDSPLIT
1912  // Only is user set -pc_type fieldsplit
1913  if (is_pcfs == PETSC_TRUE) {
1914  IS is_disp, is_rho;
1915  const MoFEM::Problem *problem_ptr;
1916  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
1917  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1918  problem_ptr->getName(), ROW, "DISPLACEMENTS", 0, 3, &is_disp);
1919  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1920  problem_ptr->getName(), ROW, "RHO", 0, 1, &is_rho);
1921  CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
1922  CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
1923  CHKERR ISDestroy(&is_disp);
1924  CHKERR ISDestroy(&is_rho);
1925  }
1926 
1927  // Monitor
1928  MonitorPostProc post_proc(mField, commonData);
1929  TsCtx *ts_ctx;
1930  DMMoFEMGetTsCtx(dm, &ts_ctx);
1931  {
1932  ts_ctx->get_postProcess_to_do_Monitor().push_back(&post_proc);
1933  CHKERR TSMonitorSet(ts, f_TSMonitorSet, ts_ctx, PETSC_NULL);
1934  }
1935 
1936 #if PETSC_VERSION_GE(3, 7, 0)
1937  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1938 #endif
1939  CHKERR TSSolve(ts, D);
1940  CHKERR TSGetTime(ts, &ftime);
1941  PetscInt steps, snesfails, rejects, nonlinits, linits;
1942 #if PETSC_VERSION_LE(3, 8, 0)
1943  CHKERR TSGetTimeStepNumber(ts, &steps);
1944 #else
1945  CHKERR TSGetStepNumber(ts, &steps);
1946 #endif
1947 
1948  CHKERR TSGetSNESFailures(ts, &snesfails);
1949  CHKERR TSGetStepRejections(ts, &rejects);
1950  CHKERR TSGetSNESIterations(ts, &nonlinits);
1951  CHKERR TSGetKSPIterations(ts, &linits);
1952  PetscPrintf(PETSC_COMM_WORLD,
1953  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
1954  "linits %D\n",
1955  steps, rejects, snesfails, ftime, nonlinits, linits);
1956 
1958  if (commonData.cUrrent_psi < 1.67 || commonData.cUrrent_psi > 1.68)
1959  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "atom test diverged");
1960  }
1961 
1963 }
1964 
1966  int row_side, EntityType row_type,
1969 
1970  if (row_type != MBVERTEX)
1972 
1973  const int nb_gauss_pts = row_data.getN().size1();
1974  // commonData.data.vecPsi.resize(nb_gauss_pts,false);
1977 
1978  double energy = 0;
1979  double mass = 0;
1980  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1981 
1982  double vol = getVolume() * getGaussPts()(3, gg);
1983  if (getHoGaussPtsDetJac().size() > 0) {
1984  vol *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
1985  }
1986 
1987  energy += vol * psi;
1988  mass += vol * rho;
1989 
1990  ++psi;
1991  ++rho;
1992  }
1993 
1994  CHKERR VecSetValue(energVec, 0, energy, ADD_VALUES);
1995  CHKERR VecSetValue(massVec, 0, mass, ADD_VALUES);
1997 }
1998 
1999 } // namespace BoneRemodeling
PetscBool equilibrium_flg
Definition: Remodeling.hpp:561
Evaluate density derivative with respect to time in case of Backward Euler Method The density within...
Definition: Remodeling.hpp:332
moab::Interface & postProcMesh
Definition: Remodeling.hpp:378
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:377
bool & doTris
\deprectaed
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:820
MoFEMErrorCode getParameters()
Get parameters form line command or config file.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
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
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:141
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
boost::shared_ptr< VectorDouble > rhoPtr
Ptr to density matrix container.
Definition: Remodeling.hpp:187
#define LAMBDA(E, NU)
Definition: fem_tools.h:32
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
Deprecated interface functions.
Problem manager is used to build and partition problems.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
CommonData commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:456
OpCalculateStress(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:280
Force scale operator for reading two columns.
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT)
Definition: Remodeling.cpp:40
PetscBool less_post_proc
reduce file size
Definition: Remodeling.hpp:163
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:482
virtual moab::Interface & get_moab()=0
OpCalculateStressTangent(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:717
double cUrrent_psi
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:158
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:892
Off diagonal block of tangent matrix .
Definition: Remodeling.hpp:493
Assemble residual for conservation of mass (density) .
Definition: Remodeling.hpp:453
bool & doVertices
\deprectaed If false skip vertices
MoFEMErrorCode getParameters()
Get the material parameters from the command line.
boost::shared_ptr< MatrixDouble > gradRhoPtr
Gradient of density field.
Definition: Remodeling.hpp:188
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:475
Definition: Remodeling.hpp:553
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
Diagonal block of tangent matrix .
Definition: Remodeling.hpp:472
VectorDouble vecRhoDt
Time derivative of density.
Definition: Remodeling.hpp:201
Off diagonal block of tangent matrix .
Definition: Remodeling.hpp:513
Section manager is used to create indexes and sectionsFIXME: ISManager is not properly testsed by ato...
Definition: ISManager.hpp:36
moab::Interface & postProcMesh
double pSi_ref
reference free energy
Definition: Remodeling.hpp:154
bool & doPrisms
\deprectaed
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
BasicMethodsSequence & get_postProcess_to_do_Monitor()
Get the postProcess to do Monitor object.
Definition: TsCtx.hpp:187
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:27
MoFEMErrorCode addMomentumFluxes()
Finite elements to calculate tractions.
#define TENSOR4_MAT_PTR(MAT)
Definition: Remodeling.cpp:48
bool iNit
Definition: Remodeling.hpp:563
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:425
double cUrrent_mass
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:160
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:732
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
double rate
Definition: Remodeling.hpp:562
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Remodeling.cpp:977
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
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.
PetscBool mass_postproc
Definition: Remodeling.hpp:560
FTensor::Index< 'I', 3 > I
Definition: PlasticOps.hpp:70
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:550
OpCalculateStressTangentWithAdolc(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:382
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.
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:74
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
MoFEMErrorCode solveDM()
Solve problem set up in DM.
MoFEMErrorCode addElements()
Set and add finite elements.
keeps basic data about problemThis is low level structure with information about problem,...
VectorDouble nF
Vector of the right hand side (internal forces)
Definition: Remodeling.hpp:457
double getCFromDensity(const double &rho)
Definition: Remodeling.cpp:218
MonitorPostProc(MoFEM::Interface &m_field, Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
Off diagonal block of tangent matrix /f[ K_{u \rho}=\intop_{V} \left[\frac{n}{\rho_{0}}\right] \left...
Definition: Remodeling.hpp:536
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:407
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:71
Projection of edge entities with one mid-node on hierarchical basis.
boost::shared_ptr< ElasticMaterials > elasticMaterialsPtr
Definition: Remodeling.hpp:141
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:335
MoFEMErrorCode preProcess()
constexpr double young_modulus
MatrixDouble matGradientOfDeformation
Gradient of deformation.
Definition: Remodeling.hpp:198
OpAssmbleStressLhs_dRho(Remodeling::CommonData &common_data)
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT)
Definition: Remodeling.cpp:54
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:539
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.
Definition: Templates.hpp:463
constexpr double poisson_ratio
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)
Definition: lapack_wrap.h:273
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
OpAssmbleStressLhs_dF(Remodeling::CommonData &common_data)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Assemble residual for conservation of momentum (stresses)
Definition: Remodeling.hpp:435
boost::shared_ptr< MatrixDouble > gradDispPtr
Ptr to gradient of displacements matrix container.
Definition: Remodeling.hpp:185
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode operator()()
bool sYmm
If true assume that matrix is symmetric structure.
virtual int get_comm_rank() const =0
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEM::Interface & mField
Definition: Remodeling.hpp:555
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:239
PetscReal ts_a
shift for U_tt (see PETSc Time Solver)
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:140
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:359
MoFEM::Interface & mField
Definition: Remodeling.hpp:206
Evaluate physical equations at integration points.
Definition: Remodeling.hpp:356
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Basic algebra on fields.
Definition: FieldBlas.hpp:34
CHKERRQ(ierr)
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
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.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::vector< EntityHandle > & mapGaussPts
Definition: Remodeling.hpp:379
double m
algorithmic exponent [-]
Definition: Remodeling.hpp:147
Managing BitRefLevels.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:104
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
VectorDouble vecPsi
Elastic energy density.
Definition: Remodeling.hpp:189
bool & doEdges
\deprectaed If false skip edges
double getCFromDensityDiff(const double &rho)
Definition: Remodeling.cpp:225
MoFEMErrorCode addFields()
Set and add entities to approximation fields.
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.
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:496
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:101
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
Definition: Remodeling.hpp:136
int pRT
Definition: Remodeling.hpp:564
bool & doQuads
\deprectaed
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.
Definition: DMMMoFEM.cpp:105
Elastic materials.
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
Definition: Remodeling.hpp:137
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
PetscBool is_atom_testing
for atom tests
Definition: Remodeling.hpp:162
#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.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:602
#define MU(E, NU)
Definition: fem_tools.h:33
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
PostProcVolumeOnRefinedMesh postProc
Definition: Remodeling.hpp:557
OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:969
std::vector< EntityHandle > mapGaussPts
constexpr int order
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
MoFEMErrorCode buildDM()
Set problem and DM.
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data)
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
double n
porosity exponent [-]
Definition: Remodeling.hpp:148
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:438
DM dm
Discretization manager.
Definition: Remodeling.hpp:127
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:1001
OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:886
#define TENSOR2_MAT_PTR(MAT)
Definition: Remodeling.cpp:50
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: EdgeForce.hpp:109
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
double c
density evolution (growth) velocity [d/m^2]
Definition: Remodeling.hpp:146
virtual bool check_field(const std::string &name) const =0
check if field is in database
Used to post proc stresses, energy, source term.
Definition: Remodeling.hpp:374
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
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
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
Definition: Remodeling.hpp:130
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:348
MoFEMErrorCode postProcess()
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
continuous field
Definition: definitions.h:177
Volume finite element.
Definition: Remodeling.hpp:41
DEPRECATED PetscErrorCode f_TSMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Definition: TsCtx.hpp:327
std::string getName() const
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
PostProcVolumeOnRefinedMesh postProcElastic
Definition: Remodeling.hpp:558
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)
Definition: Remodeling.cpp:57
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:556
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Definition: Remodeling.hpp:140
VectorDouble nF
Vector of the right hand side (internal forces)
Definition: Remodeling.hpp:439
Get value at integration points for scalar field.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:982
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
Definition: Remodeling.hpp:138
MatrixDouble matInvF
inverse of deformation gradient
Definition: Remodeling.hpp:199
bool & doTets
\deprectaed
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
double R0
mass conduction coefficient
Definition: Remodeling.hpp:155
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:396
double w(const double g, const double t)
Definition: ContactOps.hpp:160
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...
OpAssmbleStressRhs(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:814
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
Definition: Remodeling.hpp:131