v0.14.0
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 
62 MoFEMErrorCode Remodeling::CommonData::getParameters() {
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 
86  MoFEMFunctionBegin;
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 
178  lambda = LAMBDA(young_modulus, poisson_ratio);
179  mu = MU(young_modulus, poisson_ratio);
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 
215  MoFEMFunctionReturn(0);
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)
234  : VolumeElementForcesAndSourcesCore::UserDataOperator(
235  "RHO", UserDataOperator::OPROW),
236  commonData(common_data) {}
237 
238 MoFEMErrorCode
239 OpGetRhoTimeDirevative::doWork(int side, EntityType type,
240  DataForcesAndSourcesCore::EntData &data) {
241 
242  MoFEMFunctionBegin;
243 
244  const int nb_gauss_pts = data.getN().size1();
245  if (!nb_gauss_pts)
246  MoFEMFunctionReturnHot(0);
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)
257  MoFEMFunctionReturnHot(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 
277  MoFEMFunctionReturn(0);
278 }
279 
280 OpCalculateStress::OpCalculateStress(Remodeling::CommonData &common_data)
281  : VolumeElementForcesAndSourcesCore::UserDataOperator(
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 
293 MoFEMErrorCode
294 OpCalculateStress::doWork(int side, EntityType type,
295  DataForcesAndSourcesCore::EntData &data) {
296 
297  MoFEMFunctionBegin;
298 
299  if (type != MBVERTEX)
300  MoFEMFunctionReturnHot(9);
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 
314  auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
315 
316  commonData.data.vecPsi.resize(nb_gauss_pts, false);
317  auto psi = getFTensor0FromVec(commonData.data.vecPsi);
318  commonData.data.vecR.resize(nb_gauss_pts, false);
319  auto R = getFTensor0FromVec(commonData.data.vecR);
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);
324  auto det_f = getFTensor0FromVec(commonData.data.vecDetF);
325  MatrixDouble &matS = commonData.data.matS;
326  matS.resize(nb_gauss_pts, 6, false);
327  FTensor::Tensor2_symmetric<double *, 3> S(SYMMETRIC_TENSOR2_MAT_PTR(matS), 6);
328  commonData.data.matP.resize(9, nb_gauss_pts, false);
329  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
330  MatrixDouble &matF = commonData.data.matGradientOfDeformation;
331  matF.resize(9, nb_gauss_pts, false);
332  FTensor::Tensor2<double *, 3, 3> F(TENSOR2_MAT_PTR(matF));
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;
341  FTensor::Index<'i', 3> i;
342  FTensor::Index<'I', 3> I;
343  FTensor::Index<'J', 3> J;
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 
379  MoFEMFunctionReturn(0);
380 }
381 
382 OpCalculateStressTangentWithAdolc::OpCalculateStressTangentWithAdolc(
383  Remodeling::CommonData &common_data)
384  : VolumeElementForcesAndSourcesCore::UserDataOperator(
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 
396 MoFEMErrorCode OpCalculateStressTangentWithAdolc::doWork(
397  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
398 
399  MoFEMFunctionBegin;
400 
401  if (type != MBVERTEX)
402  MoFEMFunctionReturnHot(0);
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);
408  FTensor::Tensor2_symmetric<double *, 3> C(SYMMETRIC_TENSOR2_VEC_PTR(vecC));
409  MatrixDouble &dS_dC = commonData.data.matMaterialTangent;
410  dS_dC.resize(6, 6 * nb_gauss_pts, false);
411  dS_dC.clear();
412  FTensor::Ddg<double *, 3, 3> D1(SYMMETRIC_TENSOR4_MAT_PTR(dS_dC), 6);
413 
414  MatrixDouble &matS = commonData.data.matS;
415  matS.resize(nb_gauss_pts, 6, false);
416  FTensor::Tensor2_symmetric<double *, 3> S(SYMMETRIC_TENSOR2_MAT_PTR(matS), 6);
417  MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
418  dP_dF.resize(81, nb_gauss_pts, false);
419  MatrixDouble &matF = commonData.data.matGradientOfDeformation;
420  matF.resize(9, nb_gauss_pts, false);
421  FTensor::Tensor2<double *, 3, 3> F(TENSOR2_MAT_PTR(matF));
422  FTensor::Tensor4<double *, 3, 3, 3, 3> D2(TENSOR4_MAT_PTR(dP_dF));
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);
426  auto psi = getFTensor0FromVec(commonData.data.vecPsi);
427 
428  FTensor::Index<'i', 3> i;
429  FTensor::Index<'k', 3> k;
430  FTensor::Index<'I', 3> I;
431  FTensor::Index<'J', 3> J;
432  FTensor::Index<'K', 3> K;
433  FTensor::Index<'L', 3> L;
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>>,
450  FTensor::Tensor2_symmetric<double *, 3>>(
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 
547  MoFEMFunctionReturn(0);
548 }
549 
550 OpPostProcStress::OpPostProcStress(moab::Interface &post_proc_mesh,
551  std::vector<EntityHandle> &map_gauss_pts,
552  Remodeling::CommonData &common_data)
553  : VolumeElementForcesAndSourcesCore::UserDataOperator(
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 
566 MoFEMErrorCode
567 OpPostProcStress::doWork(int side, EntityType type,
568  DataForcesAndSourcesCore::EntData &data) {
569 
570  MoFEMFunctionBegin;
571 
572  if (type != MBVERTEX)
573  MoFEMFunctionReturnHot(9);
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();
624  MatrixDouble &matS = commonData.data.matS;
625  FTensor::Tensor2_symmetric<double *, 3> S(SYMMETRIC_TENSOR2_MAT_PTR(matS), 6);
626  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
627  auto psi = getFTensor0FromVec(commonData.data.vecPsi);
628  auto R = getFTensor0FromVec(commonData.data.vecR);
629 
630  auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
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 
714  MoFEMFunctionReturn(0);
715 }
716 
717 OpCalculateStressTangent::OpCalculateStressTangent(
718  Remodeling::CommonData &common_data)
719  : VolumeElementForcesAndSourcesCore::UserDataOperator(
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 
731 MoFEMErrorCode
732 OpCalculateStressTangent::doWork(int side, EntityType type,
733  DataForcesAndSourcesCore::EntData &data) {
734 
735  MoFEMFunctionBegin;
736 
737  if (type != MBVERTEX)
738  MoFEMFunctionReturnHot(0);
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);
746  auto det_f = getFTensor0FromVec(commonData.data.vecDetF);
747 
748  MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
749  dP_dF.resize(81, nb_gauss_pts, false);
750  MatrixDouble &matF = commonData.data.matGradientOfDeformation;
751  matF.resize(9, nb_gauss_pts, false);
752  FTensor::Tensor2<double *, 3, 3> F(TENSOR2_MAT_PTR(matF));
753  FTensor::Tensor4<double *, 3, 3, 3, 3> D2(TENSOR4_MAT_PTR(dP_dF));
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);
757  auto psi = getFTensor0FromVec(commonData.data.vecPsi);
758  FTensor::Index<'i', 3> i;
759  FTensor::Index<'k', 3> k;
760  FTensor::Index<'I', 3> I;
761  FTensor::Index<'J', 3> J;
762  FTensor::Index<'K', 3> K;
763  FTensor::Index<'L', 3> L;
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 
811  MoFEMFunctionReturn(0);
812 }
813 
814 OpAssmbleStressRhs::OpAssmbleStressRhs(Remodeling::CommonData &common_data)
815  : VolumeElementForcesAndSourcesCore::UserDataOperator(
816  "DISPLACEMENTS", UserDataOperator::OPROW),
817  commonData(common_data) {}
818 
819 MoFEMErrorCode
820 OpAssmbleStressRhs::doWork(int side, EntityType type,
821  DataForcesAndSourcesCore::EntData &data) {
822 
823  MoFEMFunctionBegin;
824 
825  const int nb_dofs = data.getIndices().size();
826  if (!nb_dofs)
827  MoFEMFunctionReturnHot(0);
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);
835  auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
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 
845  FTensor::Index<'I', 3> I;
846  FTensor::Index<'J', 3> J;
847  FTensor::Index<'i', 3> i;
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  double a = w * pow(rho / rho_ref, n);
854 
855  int bb = 0;
856  for (; bb != nb_dofs / 3; bb++) {
857  double *ptr = &nF[3 * bb];
858  FTensor::Tensor1<double *, 3> t1(ptr, &ptr[1], &ptr[2]);
859  t1(i) += a * P(i, I) * diff_base_functions(I);
860  // t1(i) += 0.5*a*diff_base_functions(J)*(F(i,I)*S(I,J));
861  // t1(i) += 0.5*a*diff_base_functions(I)*(F(i,J)*S(I,J));
862  ++diff_base_functions;
863  }
864  // Could be that we base more base functions than needed to approx. disp.
865  // field.
866  for (; bb != nb_base_functions; bb++) {
867  ++diff_base_functions;
868  }
869  ++P;
870  // ++S;
871  // ++F;
872  ++rho;
873  }
874 
875  // Assemble internal force vector for time solver (TS)
876  CHKERR VecSetValues(
877  getFEMethod()->ts_F, /// Get the right hand side vector for time solver
878  nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
879 
880  MoFEMFunctionReturn(0);
881 }
882 
883 OpAssmbleRhoRhs::OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
884  : VolumeElementForcesAndSourcesCore::UserDataOperator(
885  "RHO", UserDataOperator::OPROW),
886  commonData(common_data) {}
887 
888 MoFEMErrorCode
889 OpAssmbleRhoRhs::doWork(int side, EntityType type,
890  DataForcesAndSourcesCore::EntData &data) {
891 
892  MoFEMFunctionBegin;
893 
894  const int nb_dofs = data.getIndices().size();
895  if (!nb_dofs)
896  MoFEMFunctionReturnHot(0);
897  nF.resize(nb_dofs, false);
898  nF.clear();
899 
900  const int nb_gauss_pts = data.getN().size1();
901  const int nb_base_functions = data.getN().size2();
902  auto base_functions = data.getFTensor0N();
903  auto diff_base_functions = data.getFTensor1DiffN<3>();
904  if (!commonData.data.rhoPtr) {
905  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "rhoPtsr is null");
906  }
907  auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
908  if (!commonData.data.gradRhoPtr) {
909  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "gradRhoPtr is null");
910  }
911  auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
912  auto R = getFTensor0FromVec(commonData.data.vecR);
913 
914  if (commonData.data.vecRhoDt.size() != nb_gauss_pts) {
915  commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
916  commonData.data.vecRhoDt.clear();
917  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Time derivative of
918  // density at integration not calculated");
919  }
920  auto drho_dt = getFTensor0FromVec(commonData.data.vecRhoDt);
921 
922  const double R0 = commonData.R0;
923 
924  FTensor::Index<'I', 3> I;
925 
926  for (int gg = 0; gg != nb_gauss_pts; gg++) {
927 
928  // Get volume and integration weight
929  double w = getVolume() * getGaussPts()(3, gg);
930 
931  int bb = 0;
932  for (; bb != nb_dofs; bb++) {
933  nF[bb] += w * base_functions * drho_dt; // Density rate term
934  nF[bb] -= w * base_functions * R; // Source term
935  nF[bb] -= w * R0 * diff_base_functions(I) * grad_rho(I); // Gradient term
936  ++base_functions;
937  ++diff_base_functions;
938  }
939  // Could be that we base more base functions than needed to approx. density
940  // field.
941  for (; bb != nb_base_functions; bb++) {
942  ++base_functions;
943  ++diff_base_functions;
944  }
945 
946  // Increment quantities for integration pts.
947  ++rho;
948  ++grad_rho;
949  ++drho_dt;
950  ++R;
951  }
952 
953  // Assemble internal force vector for time solver (TS)
954  CHKERR VecSetValues(
955  getFEMethod()->ts_F, /// Get the right hand side vector for time solver
956  nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
957 
958  MoFEMFunctionReturn(0);
959 }
960 
961 OpAssmbleRhoLhs_dRho::OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
962  : VolumeElementForcesAndSourcesCore::UserDataOperator(
963  "RHO", "RHO", UserDataOperator::OPROWCOL),
964  commonData(common_data) {
965  sYmm = true;
966 }
967 
968 MoFEMErrorCode
969 OpAssmbleRhoLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
970  EntityType col_type,
971  DataForcesAndSourcesCore::EntData &row_data,
972  DataForcesAndSourcesCore::EntData &col_data) {
973 
974  MoFEMFunctionBegin;
975 
976  const int row_nb_dofs = row_data.getIndices().size();
977  if (!row_nb_dofs)
978  MoFEMFunctionReturnHot(0);
979  const int col_nb_dofs = col_data.getIndices().size();
980  if (!col_nb_dofs)
981  MoFEMFunctionReturnHot(0);
982 
983  // Set size can clear local tangent matrix
984  locK_rho_rho.resize(row_nb_dofs, col_nb_dofs, false);
985  locK_rho_rho.clear();
986 
987  const int row_nb_gauss_pts = row_data.getN().size1();
988  if (!row_nb_gauss_pts)
989  MoFEMFunctionReturnHot(0);
990  const int row_nb_base_functions = row_data.getN().size2();
991 
992  const double ts_a = getFEMethod()->ts_a;
993  const double R0 = commonData.R0;
994  const double c = commonData.c;
995  const double n = commonData.n;
996  const double m = commonData.m;
997  const double rho_ref = commonData.rHo_ref;
998  const double psi_ref = commonData.pSi_ref;
999 
1000  FTensor::Index<'I', 3> I;
1001 
1002  auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1003  auto psi = getFTensor0FromVec(commonData.data.vecPsi);
1004  auto row_base_functions = row_data.getFTensor0N();
1005  auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1006 
1007  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1008 
1009  // Get volume and integration weight
1010  double w = getVolume() * getGaussPts()(3, gg);
1011 
1012  const double kp = commonData.getCFromDensity(rho);
1013  const double kp_diff = commonData.getCFromDensityDiff(rho);
1014  const double a = c * kp * ((n - m) / rho) * pow(rho / rho_ref, n - m) * psi;
1015  const double a_diff =
1016  c * kp_diff * (pow(rho / rho_ref, n - m) * psi - psi_ref);
1017 
1018  int row_bb = 0;
1019  for (; row_bb != row_nb_dofs; row_bb++) {
1020  auto col_base_functions = col_data.getFTensor0N(gg, 0);
1021  auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1022  for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1023  locK_rho_rho(row_bb, col_bb) +=
1024  ts_a * w * row_base_functions * col_base_functions;
1025  locK_rho_rho(row_bb, col_bb) -=
1026  w * R0 * row_diff_base_functions(I) * col_diff_base_functions(I);
1027  locK_rho_rho(row_bb, col_bb) -=
1028  a * w * row_base_functions * col_base_functions;
1029  locK_rho_rho(row_bb, col_bb) -=
1030  a_diff * w * row_base_functions * col_base_functions;
1031 
1032  ++col_base_functions;
1033  ++col_diff_base_functions;
1034  }
1035  ++row_base_functions;
1036  ++row_diff_base_functions;
1037  }
1038  for (; row_bb != row_nb_base_functions; row_bb++) {
1039  ++row_base_functions;
1040  ++row_diff_base_functions;
1041  }
1042  ++psi;
1043  ++rho;
1044  }
1045 
1046  // cerr << locK_rho_rho << endl;
1047 
1048  CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1049  &*row_data.getIndices().begin(), col_nb_dofs,
1050  &*col_data.getIndices().begin(), &locK_rho_rho(0, 0),
1051  ADD_VALUES);
1052 
1053  // is symmetric
1054  if (row_side != col_side || row_type != col_type) {
1055  transLocK_rho_rho.resize(col_nb_dofs, row_nb_dofs, false);
1056  noalias(transLocK_rho_rho) = trans(locK_rho_rho);
1057  CHKERR MatSetValues(getFEMethod()->ts_B, col_nb_dofs,
1058  &*col_data.getIndices().begin(), row_nb_dofs,
1059  &*row_data.getIndices().begin(),
1060  &transLocK_rho_rho(0, 0), ADD_VALUES);
1061  }
1062 
1063  MoFEMFunctionReturn(0);
1064 }
1065 
1066 OpAssmbleRhoLhs_dF::OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data)
1067  : VolumeElementForcesAndSourcesCore::UserDataOperator(
1068  "RHO", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
1069  commonData(common_data) {
1070  sYmm = false;
1071 }
1072 
1073 MoFEMErrorCode
1074 OpAssmbleRhoLhs_dF::doWork(int row_side, int col_side, EntityType row_type,
1075  EntityType col_type,
1076  DataForcesAndSourcesCore::EntData &row_data,
1077  DataForcesAndSourcesCore::EntData &col_data) {
1078 
1079  MoFEMFunctionBegin;
1080 
1081  const int row_nb_dofs = row_data.getIndices().size();
1082  if (!row_nb_dofs)
1083  MoFEMFunctionReturnHot(0);
1084  const int col_nb_dofs = col_data.getIndices().size();
1085  if (!col_nb_dofs)
1086  MoFEMFunctionReturnHot(0);
1087 
1088  // Set size can clear local tangent matrix
1089  locK_rho_F.resize(row_nb_dofs, col_nb_dofs, false);
1090  locK_rho_F.clear();
1091 
1092  const int row_nb_gauss_pts = row_data.getN().size1();
1093  const int col_nb_base_functions = col_data.getN().size2();
1094 
1095  const double c = commonData.c;
1096  const double n = commonData.n;
1097  const double m = commonData.m;
1098  const double rho_ref = commonData.rHo_ref;
1099 
1100  FTensor::Index<'I', 3> I;
1101  FTensor::Index<'i', 3> i;
1102 
1103  auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1104  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
1105 
1106  FTensor::Tensor1<double, 3> t1;
1107  auto col_diff_base_functions = col_data.getFTensor1DiffN<3>();
1108 
1109  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1110 
1111  // Get volume and integration weight
1112  double w = getVolume() * getGaussPts()(3, gg);
1113  const double kp = commonData.getCFromDensity(rho);
1114  const double a = w * c * kp * pow(rho / rho_ref, n - m);
1115 
1116  int col_bb = 0;
1117  for (; col_bb != col_nb_dofs / 3; col_bb++) {
1118  t1(i) = a * P(i, I) * col_diff_base_functions(I);
1119  auto row_base_functions = row_data.getFTensor0N(gg, 0);
1120  for (int row_bb = 0; row_bb != row_nb_dofs; row_bb++) {
1121  FTensor::Tensor1<double *, 3> k(&locK_rho_F(row_bb, 3 * col_bb + 0),
1122  &locK_rho_F(row_bb, 3 * col_bb + 1),
1123  &locK_rho_F(row_bb, 3 * col_bb + 2));
1124  k(i) -= row_base_functions * t1(i);
1125  ++row_base_functions;
1126  }
1127  ++col_diff_base_functions;
1128  }
1129  for (; col_bb != col_nb_base_functions; col_bb++) {
1130  ++col_diff_base_functions;
1131  }
1132 
1133  ++P;
1134  ++rho;
1135  }
1136 
1137  CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1138  &*row_data.getIndices().begin(), col_nb_dofs,
1139  &*col_data.getIndices().begin(), &locK_rho_F(0, 0),
1140  ADD_VALUES);
1141 
1142  MoFEMFunctionReturn(0);
1143 }
1144 template <bool ADOLC>
1145 OpAssmbleStressLhs_dF<ADOLC>::OpAssmbleStressLhs_dF(
1146  Remodeling::CommonData &common_data)
1147  : VolumeElementForcesAndSourcesCore::UserDataOperator(
1148  "DISPLACEMENTS", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
1149  commonData(common_data) {
1150  sYmm = true;
1151 }
1152 template <bool ADOLC>
1153 MoFEMErrorCode OpAssmbleStressLhs_dF<ADOLC>::doWork(
1154  int row_side, int col_side, EntityType row_type, EntityType col_type,
1155  DataForcesAndSourcesCore::EntData &row_data,
1156  DataForcesAndSourcesCore::EntData &col_data) {
1157 
1158  MoFEMFunctionBegin;
1159 
1160  const int row_nb_dofs = row_data.getIndices().size();
1161  if (!row_nb_dofs)
1162  MoFEMFunctionReturnHot(0);
1163  const int col_nb_dofs = col_data.getIndices().size();
1164  if (!col_nb_dofs)
1165  MoFEMFunctionReturnHot(0);
1166 
1167  const bool diagonal_block = (row_type == col_type) && (row_side == col_side);
1168 
1169  // Set size can clear local tangent matrix
1170  locK_P_F.resize(row_nb_dofs, col_nb_dofs, false);
1171  locK_P_F.clear();
1172 
1173  const int row_nb_gauss_pts = row_data.getN().size1();
1174  const int row_nb_base_functions = row_data.getN().size2();
1175 
1176  MatrixDouble &matS = commonData.data.matS;
1177  MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
1178  // MatrixDouble &dS_dC = commonData.data.matMaterialTangent;
1179  // MatrixDouble &matF = commonData.data.matGradientOfDeformation;
1180 
1181  double t0;
1182 
1183  const double n = commonData.n;
1184  const double rho_ref = commonData.rHo_ref;
1185 
1186  FTensor::Index<'i', 3> i;
1187  FTensor::Index<'j', 3> j;
1188  FTensor::Index<'k', 3> k;
1189  FTensor::Index<'I', 3> I;
1190  FTensor::Index<'J', 3> J;
1191  FTensor::Index<'K', 3> K;
1192  FTensor::Index<'L', 3> L;
1193 
1194  auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1195  FTensor::Tensor4<double *, 3, 3, 3, 3> D2(TENSOR4_MAT_PTR(dP_dF));
1196  FTensor::Tensor2_symmetric<double *, 3> S(SYMMETRIC_TENSOR2_MAT_PTR(matS), 6);
1197  auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1198 
1199  // FTensor::Tensor4_ddg<double*,3,3> D1(SYMMETRIC_TENSOR4_MAT_PTR(dS_dC),6);
1200  // FTensor::Tensor2<double*,3,3> F(TENSOR2_MAT_PTR(matF));
1201 
1202  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1203 
1204  // Get volume and integration weight
1205  double w = getVolume() * getGaussPts()(3, gg);
1206  const double a = w * pow(rho / rho_ref, n);
1207 
1208  int row_bb = 0;
1209  for (; row_bb != row_nb_dofs / 3; row_bb++) {
1210 
1211  auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1212 
1213  const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
1214  int col_bb = 0;
1215  for (; col_bb != final_bb; col_bb++) {
1216 
1217  FTensor::Tensor2<double *, 3, 3> t_assemble(
1218  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1219  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1220  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1221  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1222  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1223  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1224  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1225  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1226  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1227 
1228  diffDiff(J, L) =
1229  a * row_diff_base_functions(J) * col_diff_base_functions(L);
1230 
1231  // TODO: FIXME: Tensor D2 has major symmetry, we do not have such tensor
1232  // implemented at the moment. Using tensor with major symmetry will
1233  // improve code efficiency.
1234  t_assemble(i, k) += diffDiff(J, L) * D2(i, J, k, L);
1235 
1236  if (ADOLC) {
1237  // Stress stiffening part (only with adolc since it is using dS / dC
1238  // derrivative. Check OpCalculateStressTangentWithAdolc operator)
1239  t0 = diffDiff(J, L) * S(J, L);
1240  t_assemble(0, 0) += t0;
1241  t_assemble(1, 1) += t0;
1242  t_assemble(2, 2) += t0;
1243  }
1244 
1245  // Next base function for column
1246  ++col_diff_base_functions;
1247  }
1248 
1249  ++row_diff_base_functions;
1250  }
1251  for (; row_bb != row_nb_base_functions; row_bb++) {
1252  ++row_diff_base_functions;
1253  }
1254 
1255  // ++D1;
1256  ++D2;
1257  ++S;
1258  // ++F;
1259  ++rho;
1260  }
1261 
1262  // Copy of symmetric part for assemble
1263  if (diagonal_block) {
1264  for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
1265  int col_bb = 0;
1266  for (; col_bb != row_bb + 1; col_bb++) {
1267  FTensor::Tensor2<double *, 3, 3> t_assemble(
1268  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1269  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1270  &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1271  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1272  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1273  &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1274  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1275  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1276  &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1277  FTensor::Tensor2<double *, 3, 3> t_off_side(
1278  &locK_P_F(3 * col_bb + 0, 3 * row_bb + 0),
1279  &locK_P_F(3 * col_bb + 0, 3 * row_bb + 1),
1280  &locK_P_F(3 * col_bb + 0, 3 * row_bb + 2),
1281  &locK_P_F(3 * col_bb + 1, 3 * row_bb + 0),
1282  &locK_P_F(3 * col_bb + 1, 3 * row_bb + 1),
1283  &locK_P_F(3 * col_bb + 1, 3 * row_bb + 2),
1284  &locK_P_F(3 * col_bb + 2, 3 * row_bb + 0),
1285  &locK_P_F(3 * col_bb + 2, 3 * row_bb + 1),
1286  &locK_P_F(3 * col_bb + 2, 3 * row_bb + 2));
1287  t_off_side(i, k) = t_assemble(k, i);
1288  }
1289  }
1290  }
1291 
1292  const int *row_ind = &*row_data.getIndices().begin();
1293  const int *col_ind = &*col_data.getIndices().begin();
1294  CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs, row_ind, col_nb_dofs,
1295  col_ind, &locK_P_F(0, 0), ADD_VALUES);
1296 
1297  if (row_type != col_type || row_side != col_side) {
1298  transLocK_P_F.resize(col_nb_dofs, row_nb_dofs, false);
1299  noalias(transLocK_P_F) = trans(locK_P_F);
1300  CHKERR MatSetValues(getFEMethod()->ts_B, col_nb_dofs, col_ind, row_nb_dofs,
1301  row_ind, &transLocK_P_F(0, 0), ADD_VALUES);
1302  }
1303 
1304  MoFEMFunctionReturn(0);
1305 }
1306 
1307 OpAssmbleStressLhs_dRho::OpAssmbleStressLhs_dRho(
1308  Remodeling::CommonData &common_data)
1309  : VolumeElementForcesAndSourcesCore::UserDataOperator(
1310  "DISPLACEMENTS", "RHO", UserDataOperator::OPROWCOL),
1311  commonData(common_data) {
1312  sYmm = false; // This is off-diagonal (upper-left), so no symmetric
1313 }
1314 
1315 MoFEMErrorCode
1316 OpAssmbleStressLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
1317  EntityType col_type,
1318  DataForcesAndSourcesCore::EntData &row_data,
1319  DataForcesAndSourcesCore::EntData &col_data) {
1320 
1321  MoFEMFunctionBegin;
1322 
1323  const int row_nb_dofs = row_data.getIndices().size();
1324  if (!row_nb_dofs)
1325  MoFEMFunctionReturnHot(0);
1326  const int col_nb_dofs = col_data.getIndices().size();
1327  if (!col_nb_dofs)
1328  MoFEMFunctionReturnHot(0);
1329 
1330  // Set size can clear local tangent matrix
1331  locK_P_RHO.resize(row_nb_dofs, col_nb_dofs, false);
1332  locK_P_RHO.clear();
1333 
1334  FTensor::Tensor1<double, 3> t1;
1335 
1336  const double n = commonData.n;
1337  const double rho_ref = commonData.rHo_ref;
1338 
1339  FTensor::Index<'I', 3> I;
1340  FTensor::Index<'i', 3> i;
1341 
1342  auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
1343  auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1344 
1345  const int row_nb_gauss_pts = row_data.getN().size1();
1346  const int row_nb_base_functions = row_data.getN().size2();
1347  auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1348 
1349  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1350 
1351  // Get volume and integration weight
1352  double w = getVolume() * getGaussPts()(3, gg);
1353  double a = w * (n / rho) * pow(rho / rho_ref, n);
1354 
1355  int row_bb = 0;
1356  for (; row_bb != row_nb_dofs / 3; row_bb++) {
1357  t1(i) = a * row_diff_base_functions(I) * P(i, I);
1358  auto base_functions = col_data.getFTensor0N(gg, 0);
1359  for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1360  // if(base_functions!=col_data.getN(gg)[col_bb]) {
1361  // cerr << "Error" << endl;
1362  // }
1363  FTensor::Tensor1<double *, 3> k(&locK_P_RHO(3 * row_bb + 0, col_bb),
1364  &locK_P_RHO(3 * row_bb + 1, col_bb),
1365  &locK_P_RHO(3 * row_bb + 2, col_bb));
1366  k(i) += t1(i) * base_functions;
1367  ++base_functions;
1368  }
1369  ++row_diff_base_functions;
1370  }
1371  for (; row_bb != row_nb_base_functions; row_bb++) {
1372  ++row_diff_base_functions;
1373  }
1374 
1375  ++P;
1376  ++rho;
1377  }
1378 
1379  CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1380  &*row_data.getIndices().begin(), col_nb_dofs,
1381  &*col_data.getIndices().begin(), &locK_P_RHO(0, 0),
1382  ADD_VALUES);
1383 
1384  MoFEMFunctionReturn(0);
1385 }
1386 
1387 MonitorPostProc::MonitorPostProc(MoFEM::Interface &m_field,
1388  Remodeling::CommonData &common_data)
1389  : mField(m_field), commonData(common_data), postProc(m_field),
1390  postProcElastic(m_field),
1391  // densityMaps(m_field),
1392  iNit(false) {}
1393 
1394 MoFEMErrorCode MonitorPostProc::preProcess() {
1395  MoFEMFunctionBeginHot;
1396  MoFEMFunctionReturnHot(0);
1397 }
1398 
1399 MoFEMErrorCode MonitorPostProc::operator()() {
1400  MoFEMFunctionBeginHot;
1401  MoFEMFunctionReturnHot(0);
1402 }
1403 
1404 MoFEMErrorCode MonitorPostProc::postProcess() {
1405  MoFEMFunctionBegin;
1406 
1407  PetscBool mass_postproc = PETSC_FALSE;
1408  PetscBool equilibrium_flg = PETSC_FALSE;
1409  PetscBool analysis_mesh_flg = PETSC_FALSE;
1410  rate = 1;
1411 
1412  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling post-process",
1413  "none");
1414 
1415  CHKERR PetscOptionsBool(
1416  "-mass_postproc",
1417  "is used for testing, calculates mass and energy at each time step", "",
1418  mass_postproc, &mass_postproc, PETSC_NULL);
1419 
1420  CHKERR PetscOptionsBool("-analysis_mesh",
1421  "saves analysis mesh at each time step", "",
1422  analysis_mesh_flg, &analysis_mesh_flg, PETSC_NULL);
1423 
1424  CHKERR PetscOptionsScalar(
1425  "-equilibrium_stop_rate",
1426  "is used to stop calculations when equilibium state is achieved", "",
1427  rate, &rate, &equilibrium_flg);
1428  ierr = PetscOptionsEnd();
1429  CHKERRQ(ierr);
1430 
1431  if (!iNit) {
1432 
1433  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "",
1434  "Bone remodeling post-process", "none");
1435 
1436  pRT = 1;
1437  CHKERR PetscOptionsInt("-my_output_prt",
1438  "frequncy how often results are dumped on hard disk",
1439  "", 1, &pRT, NULL);
1440  ierr = PetscOptionsEnd();
1441  CHKERRQ(ierr);
1442 
1443  CHKERR postProc.generateReferenceElementMesh();
1444  CHKERR postProc.addFieldValuesPostProc("DISPLACEMENTS");
1445  CHKERR postProc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
1446  if (!commonData.less_post_proc) {
1447  CHKERR postProc.addFieldValuesGradientPostProc("DISPLACEMENTS");
1448  }
1449  // CHKERR postProc.addFieldValuesPostProc("RHO");
1450  // CHKERR postProc.addFieldValuesGradientPostProc("RHO");
1451 
1452  {
1453  auto &list_of_operators = postProc.getOpPtrVector();
1454 
1455  list_of_operators.push_back(
1456  new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
1457  list_of_operators.push_back(new OpCalculateScalarFieldGradient<3>(
1458  "RHO", commonData.data.gradRhoPtr));
1459  // Get displacement gradient at integration points
1460  list_of_operators.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1461  "DISPLACEMENTS", commonData.data.gradDispPtr));
1462  list_of_operators.push_back(new OpCalculateStress(commonData));
1463  list_of_operators.push_back(new OpPostProcStress(
1464  postProc.postProcMesh, postProc.mapGaussPts, commonData));
1465  }
1466 
1467  CHKERR postProcElastic.generateReferenceElementMesh();
1468  CHKERR postProcElastic.addFieldValuesPostProc("DISPLACEMENTS");
1469  CHKERR postProcElastic.addFieldValuesPostProc("MESH_NODE_POSITIONS");
1470  CHKERR postProcElastic.addFieldValuesGradientPostProc("DISPLACEMENTS");
1471  std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
1472  commonData.elasticPtr->setOfBlocks.begin();
1473  for (; sit != commonData.elasticPtr->setOfBlocks.end(); sit++) {
1474  postProcElastic.getOpPtrVector().push_back(new PostProcStress(
1475  postProcElastic.postProcMesh, postProcElastic.mapGaussPts,
1476  "DISPLACEMENTS", sit->second, postProcElastic.commonData));
1477  }
1478 
1479  iNit = true;
1480  }
1481 
1482  if (mass_postproc) {
1483  Vec mass_vec;
1484  Vec energ_vec;
1485  int mass_vec_ghost[] = {0};
1486  CHKERR VecCreateGhost(PETSC_COMM_WORLD, (!mField.get_comm_rank()) ? 1 : 0,
1487  1, 1, mass_vec_ghost, &mass_vec);
1488 
1489  CHKERR VecZeroEntries(mass_vec);
1490  CHKERR VecDuplicate(mass_vec, &energ_vec);
1491 
1492  Remodeling::Fe energyCalc(mField);
1493  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", energyCalc, true, false, false,
1494  true);
1495  energyCalc.getOpPtrVector().push_back(
1496  new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
1497  energyCalc.getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<3>(
1498  "RHO", commonData.data.gradRhoPtr));
1499  energyCalc.getOpPtrVector().push_back(
1500  new OpCalculateVectorFieldGradient<3, 3>("DISPLACEMENTS",
1501  commonData.data.gradDispPtr));
1502  energyCalc.getOpPtrVector().push_back(new OpCalculateStress(commonData));
1503  energyCalc.getOpPtrVector().push_back(
1504  new OpMassAndEnergyCalculation("RHO", commonData, energ_vec, mass_vec));
1505  CHKERR DMoFEMLoopFiniteElements(commonData.dm, "FE_REMODELLING",
1506  &energyCalc);
1507 
1508  CHKERR VecAssemblyBegin(mass_vec);
1509  CHKERR VecAssemblyEnd(mass_vec);
1510  double mass_sum;
1511  CHKERR VecSum(mass_vec, &mass_sum);
1512 
1513  CHKERR VecAssemblyBegin(energ_vec);
1514  CHKERR VecAssemblyEnd(energ_vec);
1515  double energ_sum;
1516  CHKERR VecSum(energ_vec, &energ_sum);
1517 
1518  CHKERR PetscPrintf(PETSC_COMM_WORLD,
1519  "Time: %g Mass: %6.5g Elastic energy: %6.5g \n", ts_t,
1520  mass_sum, energ_sum);
1521  CHKERR VecDestroy(&mass_vec);
1522  CHKERR VecDestroy(&energ_vec);
1523 
1524  if (equilibrium_flg) {
1525  double equilibrium_rate =
1526  fabs(energ_sum - commonData.cUrrent_psi) / energ_sum;
1527  double equilibrium_mass_rate =
1528  fabs(mass_sum - commonData.cUrrent_mass) / mass_sum;
1529  if (equilibrium_rate < rate) {
1530  CHKERR PetscPrintf(
1531  PETSC_COMM_WORLD,
1532  "Energy equilibrium state is achieved! Difference = %0.6g %%. \n",
1533  equilibrium_rate * 100);
1534  }
1535  if (equilibrium_mass_rate < rate) {
1536  CHKERR PetscPrintf(
1537  PETSC_COMM_WORLD,
1538  "Mass equilibrium state is achieved! Difference = %0.6g %%. \n",
1539  equilibrium_mass_rate * 100);
1540  }
1541  commonData.cUrrent_psi = energ_sum;
1542  commonData.cUrrent_mass = mass_sum;
1543  }
1544  }
1545 
1546  int step;
1547 #if PETSC_VERSION_LE(3, 8, 0)
1548  CHKERR TSGetTimeStepNumber(ts, &step);
1549 #else
1550  CHKERR TSGetStepNumber(ts, &step);
1551 #endif
1552 
1553  if ((step) % pRT == 0) {
1554 
1555  ostringstream sss;
1556  sss << "out_" << step << ".h5m";
1557  CHKERR DMoFEMLoopFiniteElements(commonData.dm, "FE_REMODELLING", &postProc);
1558  CHKERR postProc.postProcMesh.write_file(sss.str().c_str(), "MOAB",
1559  "PARALLEL=WRITE_PART");
1560  if (analysis_mesh_flg) {
1561  ostringstream ttt;
1562  ttt << "analysis_mesh_" << step << ".h5m";
1563  CHKERR mField.get_moab().write_file(ttt.str().c_str(), "MOAB",
1564  "PARALLEL=WRITE_PART");
1565  }
1566 
1567  if (commonData.nOremodellingBlock) {
1568  CHKERR DMoFEMLoopFiniteElements(commonData.dm, "ELASTIC",
1569  &postProcElastic);
1570  ostringstream ss;
1571  ss << "out_elastic_" << step << ".h5m";
1572  CHKERR postProcElastic.postProcMesh.write_file(ss.str().c_str(), "MOAB",
1573  "PARALLEL=WRITE_PART");
1574  }
1575  }
1576  MoFEMFunctionReturn(0);
1577 }
1578 
1579 MoFEMErrorCode Remodeling::getParameters() {
1580 
1581  MoFEMFunctionBegin;
1582  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling", "none");
1583 
1584  commonData.oRder = 2;
1585  CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
1586  &commonData.oRder, PETSC_NULL);
1587 
1588  ierr = PetscOptionsEnd();
1589  CHKERRQ(ierr);
1590 
1591  CHKERR mField.get_moab().get_entities_by_type(0, MBTET, commonData.tEts_all);
1592  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
1593  string name = it->getName();
1594  if (name.compare(0, 14, "NO_REMODELLING") == 0) {
1595  commonData.nOremodellingBlock = true;
1596  EntityHandle meshset = it->getMeshset();
1597  CHKERR this->mField.get_moab().get_entities_by_type(
1598  meshset, MBTET, commonData.tEts_block, true);
1599  commonData.tEts_all =
1600  subtract(commonData.tEts_all, commonData.tEts_block);
1601  }
1602  }
1603 
1604  MoFEMFunctionReturn(0);
1605 }
1606 
1607 MoFEMErrorCode Remodeling::addFields() {
1608 
1609  MoFEMFunctionBegin;
1610 
1611  // Seed all mesh entities to MoFEM database, those entities can be potentially
1612  // used as finite elements or as entities which carry some approximation
1613  // field.
1614  commonData.bitLevel.set(0);
1615  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
1616  0, 3, commonData.bitLevel);
1617  int order = commonData.oRder;
1618 
1619  // Add displacement field
1620  CHKERR mField.add_field("DISPLACEMENTS", H1,
1621  /*AINSWORTH_LOBATTO_BASE*/ AINSWORTH_LEGENDRE_BASE, 3,
1622  MB_TAG_SPARSE, MF_ZERO);
1623  // Add field representing ho-geometry
1624  CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3,
1625  MB_TAG_SPARSE, MF_ZERO);
1626 
1627  // Check if density is available, if not add density field.
1628  bool add_rho_field = false;
1629  if (!mField.check_field("RHO")) {
1630  CHKERR mField.add_field("RHO", H1, AINSWORTH_LEGENDRE_BASE, 1,
1631  MB_TAG_SPARSE, MF_ZERO);
1632  // FIXME
1633  CHKERR mField.add_ents_to_field_by_type(commonData.tEts_all, MBTET, "RHO");
1634  // CHKERR mField.add_ents_to_field_by_type(0,MBTET,"RHO");
1635  CHKERR mField.synchronise_field_entities("RHO");
1636  add_rho_field = true;
1637 
1638  CHKERR mField.set_field_order(0, MBVERTEX, "RHO", 1);
1639  CHKERR mField.set_field_order(0, MBEDGE, "RHO", order - 1);
1640  CHKERR mField.set_field_order(0, MBTRI, "RHO", order - 1);
1641  CHKERR mField.set_field_order(0, MBTET, "RHO", order - 1);
1642  }
1643 
1644  // Add entities to field
1645  CHKERR mField.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENTS");
1646  CHKERR mField.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
1647 
1648  // Set approximation order to entities
1649  CHKERR mField.set_field_order(0, MBVERTEX, "DISPLACEMENTS", 1);
1650  CHKERR mField.set_field_order(0, MBEDGE, "DISPLACEMENTS", order);
1651  CHKERR mField.set_field_order(0, MBTRI, "DISPLACEMENTS", order);
1652  CHKERR mField.set_field_order(0, MBTET, "DISPLACEMENTS", order);
1653 
1654  // Assumes that geometry is approximated using 2nd order polynomials.
1655 
1656  // Apply 2nd order only on skin
1657  {
1658  // Skinner skin(&mField.get_moab());
1659  // Range faces,tets;
1660  // CHKERR mField.get_moab().get_entities_by_type(0,MBTET,tets);
1661  // CHKERR skin.find_skin(0,tets,false,faces);
1662  // Range edges;
1663  // CHKERR mField.get_moab().get_adjacencies(
1664  // faces,1,false,edges,moab::Interface::UNION
1665  // );
1666  CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
1667  }
1668 
1669  CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
1670 
1671  // Build fields
1672  CHKERR mField.build_fields();
1673 
1674  // If order was not set from CT scans set homogenous order equal to
1675  // reference bone density
1676  if (add_rho_field) {
1677  CHKERR mField.getInterface<FieldBlas>()->setField(commonData.rHo_ref,
1678  MBVERTEX, "RHO");
1679  // const DofEntity_multiIndex *dofs_ptr;
1680  // CHKERR mField.get_dofs(&dofs_ptr);
1681  // for(_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(mField,"RHO",dit)) {
1682  // cerr << (*dit)->getFieldData() << endl;
1683  // }
1684  }
1685 
1686  // Project geometry given on 10-node tets on ho-geometry
1687  Projection10NodeCoordsOnField ent_method_material(mField,
1688  "MESH_NODE_POSITIONS");
1689  CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
1690 
1691  MoFEMFunctionReturn(0);
1692 }
1693 
1694 MoFEMErrorCode Remodeling::addElements() {
1695  MoFEMFunctionBegin;
1696 
1697  CHKERR mField.add_finite_element("FE_REMODELLING", MF_ZERO);
1698  CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING",
1699  "DISPLACEMENTS");
1700  CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING",
1701  "DISPLACEMENTS");
1702  CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING",
1703  "DISPLACEMENTS");
1704  CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING", "RHO");
1705  CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING",
1706  "MESH_NODE_POSITIONS");
1707  CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING", "RHO");
1708  CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING", "RHO");
1709  CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_all, MBTET,
1710  "FE_REMODELLING");
1711 
1712  CHKERR mField.add_finite_element("ELASTIC");
1713  CHKERR mField.modify_finite_element_add_field_row("ELASTIC", "DISPLACEMENTS");
1714  CHKERR mField.modify_finite_element_add_field_col("ELASTIC", "DISPLACEMENTS");
1715  CHKERR mField.modify_finite_element_add_field_data("ELASTIC",
1716  "DISPLACEMENTS");
1717  CHKERR mField.modify_finite_element_add_field_data("ELASTIC",
1718  "MESH_NODE_POSITIONS");
1719 
1720  if (commonData.nOremodellingBlock) {
1721  CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_block,
1722  MBTET, "ELASTIC");
1723  }
1724 
1725  commonData.elasticMaterialsPtr =
1726  boost::shared_ptr<ElasticMaterials>(new ElasticMaterials(mField));
1727  commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
1728  new NonlinearElasticElement(mField, 10));
1729  CHKERR commonData.elasticMaterialsPtr->setBlocks(
1730  commonData.elasticPtr->setOfBlocks);
1731  CHKERR commonData.elasticPtr->addElement("ELASTIC", "DISPLACEMENTS");
1732  CHKERR commonData.elasticPtr->setOperators(
1733  "DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
1734 
1735  CHKERR mField.build_finite_elements();
1736  CHKERR mField.build_adjacencies(commonData.bitLevel);
1737 
1738  // Allocate memory for density and gradient of displacements at integration
1739  // points
1740  commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
1741  commonData.data.gradRhoPtr =
1742  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1743  commonData.data.gradDispPtr =
1744  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1745 
1746  // Add operators Rhs
1747  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feRhs), true, false,
1748  false, false);
1749  auto &list_of_operators_rhs = commonData.feRhs->getOpPtrVector();
1750  // Get density at integration points
1751  list_of_operators_rhs.push_back(
1752  new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
1753  list_of_operators_rhs.push_back(
1754  new OpCalculateScalarFieldGradient<3>("RHO", commonData.data.gradRhoPtr));
1755  list_of_operators_rhs.push_back(new OpGetRhoTimeDirevative(commonData));
1756  // Get displacement gradient at integration points
1757  list_of_operators_rhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1758  "DISPLACEMENTS", commonData.data.gradDispPtr));
1759  list_of_operators_rhs.push_back(new OpCalculateStress(commonData));
1760  list_of_operators_rhs.push_back(new OpAssmbleStressRhs(commonData));
1761  list_of_operators_rhs.push_back(new OpAssmbleRhoRhs(commonData));
1762 
1763  // Add operators Lhs
1764  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feLhs), true, false,
1765  false, false);
1766  auto &list_of_operators_lhs = commonData.feLhs->getOpPtrVector();
1767  // Get density at integration points
1768  list_of_operators_lhs.push_back(
1769  new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
1770  list_of_operators_rhs.push_back(
1771  new OpCalculateScalarFieldGradient<3>("RHO", commonData.data.gradRhoPtr));
1772  // Get displacement gradient at integration points
1773  list_of_operators_lhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1774  "DISPLACEMENTS", commonData.data.gradDispPtr));
1775  if (commonData.with_adol_c) {
1776  list_of_operators_lhs.push_back(
1777  new OpCalculateStressTangentWithAdolc(commonData));
1778  list_of_operators_lhs.push_back(
1779  new OpAssmbleStressLhs_dF<true>(commonData));
1780  } else {
1781  list_of_operators_lhs.push_back(new OpCalculateStressTangent(commonData));
1782  list_of_operators_lhs.push_back(
1783  new OpAssmbleStressLhs_dF<false>(commonData));
1784  }
1785  list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dRho(commonData));
1786  list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dF(commonData));
1787  list_of_operators_lhs.push_back(new OpAssmbleStressLhs_dRho(commonData));
1788 
1789  MoFEMFunctionReturn(0);
1790 }
1791 
1792 MoFEMErrorCode Remodeling::addMomentumFluxes() {
1793 
1794  MoFEMFunctionBegin;
1795  // Add Neumann forces elements
1796  CHKERR MetaNeummanForces::addNeumannBCElements(mField, "DISPLACEMENTS");
1797  CHKERR MetaNodalForces::addElement(mField, "DISPLACEMENTS");
1798  CHKERR MetaEdgeForces::addElement(mField, "DISPLACEMENTS");
1799 
1800  // forces and pressures on surface
1801  CHKERR MetaNeummanForces::setMomentumFluxOperators(
1802  mField, commonData.neumannForces, PETSC_NULL, "DISPLACEMENTS");
1803  // noadl forces
1804  CHKERR MetaNodalForces::setOperators(mField, commonData.nodalForces,
1805  PETSC_NULL, "DISPLACEMENTS");
1806  // edge forces
1807  CHKERR MetaEdgeForces::setOperators(mField, commonData.edgeForces, PETSC_NULL,
1808  "DISPLACEMENTS");
1809 
1810  for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
1811  commonData.neumannForces.begin();
1812  mit != commonData.neumannForces.end(); mit++) {
1813  mit->second->methodsOp.push_back(
1814  new TimeForceScale("-my_load_history", false));
1815  }
1816  for (boost::ptr_map<string, NodalForce>::iterator mit =
1817  commonData.nodalForces.begin();
1818  mit != commonData.nodalForces.end(); mit++) {
1819  mit->second->methodsOp.push_back(
1820  new TimeForceScale("-my_load_history", false));
1821  }
1822  for (boost::ptr_map<string, EdgeForce>::iterator mit =
1823  commonData.edgeForces.begin();
1824  mit != commonData.edgeForces.end(); mit++) {
1825  mit->second->methodsOp.push_back(
1826  new TimeForceScale("-my_load_history", false));
1827  }
1828 
1829  MoFEMFunctionReturn(0);
1830 }
1831 
1832 MoFEMErrorCode Remodeling::buildDM() {
1833 
1834  MoFEMFunctionBegin;
1835  commonData.dm_name = "DM_REMODELING";
1836  // Register DM problem
1837  CHKERR DMRegister_MoFEM(commonData.dm_name);
1838  CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
1839  CHKERR DMSetType(commonData.dm, commonData.dm_name);
1840  CHKERR DMMoFEMSetIsPartitioned(commonData.dm, PETSC_TRUE);
1841  // Create DM instance
1842  CHKERR DMMoFEMCreateMoFEM(commonData.dm, &mField, commonData.dm_name,
1843  commonData.bitLevel);
1844  // Configure DM form line command options (DM itself, solvers,
1845  // pre-conditioners, ... )
1846  CHKERR DMSetFromOptions(commonData.dm);
1847  // Add elements to dm (only one here)
1848  CHKERR DMMoFEMAddElement(commonData.dm, "FE_REMODELLING");
1849  CHKERR DMMoFEMAddElement(commonData.dm, "ELASTIC");
1850 
1851  if (mField.check_finite_element("FORCE_FE")) {
1852  CHKERR DMMoFEMAddElement(commonData.dm, "FORCE_FE");
1853  }
1854  if (mField.check_finite_element("PRESSURE_FE")) {
1855  CHKERR DMMoFEMAddElement(commonData.dm, "PRESSURE_FE");
1856  }
1857  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
1858  CHKERR DMSetUp(commonData.dm);
1859 
1860  MoFEMFunctionReturn(0);
1861 }
1862 
1863 MoFEMErrorCode Remodeling::solveDM() {
1864 
1865  MoFEMFunctionBegin;
1866 
1867  Vec F = commonData.F;
1868  Vec D = commonData.D;
1869  Mat A = commonData.A;
1870  DM dm = commonData.dm;
1871  TS ts = commonData.ts;
1872 
1873  CHKERR TSSetIFunction(ts, F, PETSC_NULL, PETSC_NULL);
1874  CHKERR TSSetIJacobian(ts, A, A, PETSC_NULL, PETSC_NULL);
1875  double ftime = 1;
1876 #if PETSC_VERSION_GE(3, 8, 0)
1877  CHKERR TSSetMaxTime(ts, ftime);
1878 #endif
1879  CHKERR TSSetFromOptions(ts);
1880  CHKERR TSSetDM(ts, dm);
1881 
1882  SNES snes;
1883  CHKERR TSGetSNES(ts, &snes);
1884  KSP ksp;
1885  CHKERR SNESGetKSP(snes, &ksp);
1886  PC pc;
1887  CHKERR KSPGetPC(ksp, &pc);
1888  PetscBool is_pcfs = PETSC_FALSE;
1889  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1890 
1891  // Set up FIELDSPLIT
1892  // Only is user set -pc_type fieldsplit
1893  if (is_pcfs == PETSC_TRUE) {
1894  IS is_disp, is_rho;
1895  const MoFEM::Problem *problem_ptr;
1896  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
1897  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1898  problem_ptr->getName(), ROW, "DISPLACEMENTS", 0, 3, &is_disp);
1899  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1900  problem_ptr->getName(), ROW, "RHO", 0, 1, &is_rho);
1901  CHKERR ISSort(is_disp);
1902  CHKERR ISSort(is_rho);
1903  CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
1904  CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
1905  CHKERR ISDestroy(&is_disp);
1906  CHKERR ISDestroy(&is_rho);
1907  }
1908 
1909  // Monitor
1910  MonitorPostProc post_proc(mField, commonData);
1911  TsCtx *ts_ctx;
1912  DMMoFEMGetTsCtx(dm, &ts_ctx);
1913  {
1914  ts_ctx->getPostProcessMonitor().push_back(&post_proc);
1915  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
1916  }
1917 
1918 #if PETSC_VERSION_GE(3, 7, 0)
1919  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1920 #endif
1921  CHKERR TSSolve(ts, D);
1922  CHKERR TSGetTime(ts, &ftime);
1923  PetscInt steps, snesfails, rejects, nonlinits, linits;
1924 #if PETSC_VERSION_LE(3, 8, 0)
1925  CHKERR TSGetTimeStepNumber(ts, &steps);
1926 #else
1927  CHKERR TSGetStepNumber(ts, &steps);
1928 #endif
1929 
1930  CHKERR TSGetSNESFailures(ts, &snesfails);
1931  CHKERR TSGetStepRejections(ts, &rejects);
1932  CHKERR TSGetSNESIterations(ts, &nonlinits);
1933  CHKERR TSGetKSPIterations(ts, &linits);
1934  PetscPrintf(PETSC_COMM_WORLD,
1935  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
1936  "linits %D\n",
1937  steps, rejects, snesfails, ftime, nonlinits, linits);
1938 
1939  if (commonData.is_atom_testing) {
1940  if (commonData.cUrrent_psi < 1.67 || commonData.cUrrent_psi > 1.68)
1941  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "atom test diverged");
1942  }
1943 
1944  MoFEMFunctionReturn(0);
1945 }
1946 
1947 MoFEMErrorCode OpMassAndEnergyCalculation::doWork(
1948  int row_side, EntityType row_type,
1949  DataForcesAndSourcesCore::EntData &row_data) {
1950  MoFEMFunctionBegin;
1951 
1952  if (row_type != MBVERTEX)
1953  MoFEMFunctionReturnHot(0);
1954 
1955  const int nb_gauss_pts = row_data.getN().size1();
1956  // commonData.data.vecPsi.resize(nb_gauss_pts,false);
1957  auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1958  auto psi = getFTensor0FromVec(commonData.data.vecPsi);
1959 
1960  double energy = 0;
1961  double mass = 0;
1962  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1963 
1964  double vol = getVolume() * getGaussPts()(3, gg);
1965  energy += vol * psi;
1966  mass += vol * rho;
1967 
1968  ++psi;
1969  ++rho;
1970  }
1971 
1972  CHKERR VecSetValue(energVec, 0, energy, ADD_VALUES);
1973  CHKERR VecSetValue(massVec, 0, mass, ADD_VALUES);
1974  MoFEMFunctionReturn(0);
1975 }
1976 
1977 } // namespace BoneRemodeling
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
BoneRemodeling::OpAssmbleRhoLhs_dRho::doWork
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:969
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
BoneRemodeling::OpAssmbleStressRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:820
BoneRemodeling::OpAssmbleStressLhs_dRho::OpAssmbleStressLhs_dRho
OpAssmbleStressLhs_dRho(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:1307
BoneRemodeling::Remodeling::CommonData::nOremodellingBlock
bool nOremodellingBlock
Definition: Remodeling.hpp:164
BoneRemodeling::Remodeling::buildDM
MoFEMErrorCode buildDM()
Set problem and DM.
Definition: Remodeling.cpp:1832
MoFEM::DataOperator::doTris
bool & doTris
\deprectaed
Definition: DataOperators.hpp:93
BoneRemodeling::Remodeling::addFields
MoFEMErrorCode addFields()
Set and add entities to approximation fields.
Definition: Remodeling.cpp:1607
BoneRemodeling::Remodeling::mField
MoFEM::Interface & mField
Definition: Remodeling.hpp:206
BoneRemodeling::Remodeling::CommonData::feRhs
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
Definition: Remodeling.hpp:131
BoneRemodeling::OpGetRhoTimeDirevative::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:239
MoFEM::CoreInterface::loop_dofs
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.
BoneRemodeling::Remodeling::CommonData::ts
TS ts
Time solver.
Definition: Remodeling.hpp:128
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
MoFEM::TSMethod::ts_a
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
Definition: LoopMethods.hpp:160
BoneRemodeling::Remodeling::CommonData::DataContainers::vecDetF
VectorDouble vecDetF
determinant of F
Definition: Remodeling.hpp:200
BoneRemodeling::Remodeling::CommonData::neumannForces
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
Definition: Remodeling.hpp:136
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
MOFEM_STD_EXCEPTION_THROW
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:39
BoneRemodeling::Remodeling::CommonData::n
double n
porosity exponent [-]
Definition: Remodeling.hpp:148
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
SYMMETRIC_TENSOR4_MAT_PTR
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT)
Definition: Remodeling.cpp:40
BoneRemodeling::MonitorPostProc::mass_postproc
PetscBool mass_postproc
Definition: Remodeling.hpp:560
EntityHandle
BoneRemodeling::OpAssmbleRhoRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:889
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
BoneRemodeling::OpAssmbleStressRhs::nF
VectorDouble nF
Vector of the right hand side (internal forces)
Definition: Remodeling.hpp:439
BoneRemodeling::OpCalculateStressTangent::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:732
BoneRemodeling::Remodeling::CommonData::dm_name
DMType dm_name
dm (problem) name
Definition: Remodeling.hpp:126
MetaNodalForces::addElement
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:92
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:122
BoneRemodeling::OpAssmbleRhoLhs_dF::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:496
rho
double rho
Definition: plastic.cpp:140
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:674
BoneRemodeling::MonitorPostProc::operator()
MoFEMErrorCode operator()()
Definition: Remodeling.cpp:1399
BoneRemodeling::Remodeling::CommonData::nodalForces
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
Definition: Remodeling.hpp:137
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
BoneRemodeling::OpAssmbleRhoLhs_dRho
Diagonal block of tangent matrix .
Definition: Remodeling.hpp:472
BoneRemodeling::OpAssmbleRhoLhs_dF::locK_rho_F
MatrixDouble locK_rho_F
Definition: Remodeling.hpp:497
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
BoneRemodeling::OpMassAndEnergyCalculation::energVec
Vec energVec
Definition: Remodeling.hpp:581
MoFEM::TsCtx::getPostProcessMonitor
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
Definition: TsCtx.hpp:144
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
BoneRemodeling::Remodeling::CommonData::rHo_ref
double rHo_ref
reference density
Definition: Remodeling.hpp:150
BoneRemodeling::OpAssmbleStressLhs_dF
Off diagonal block of tangent matrix .
Definition: Remodeling.hpp:513
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
TENSOR4_MAT_PTR
#define TENSOR4_MAT_PTR(MAT)
Definition: Remodeling.cpp:48
BoneRemodeling::OpPostProcStress::postProcMesh
moab::Interface & postProcMesh
Definition: Remodeling.hpp:378
BoneRemodeling::OpAssmbleRhoLhs_dRho::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:475
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
BoneRemodeling::Remodeling::CommonData::cUrrent_psi
double cUrrent_psi
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:158
BoneRemodeling::Remodeling::CommonData::cUrrent_mass
double cUrrent_mass
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:160
BoneRemodeling::Remodeling::CommonData::DataContainers::vecR
VectorDouble vecR
Mass sorce.
Definition: Remodeling.hpp:192
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::DataOperator::doPrisms
bool & doPrisms
\deprectaed
Definition: DataOperators.hpp:95
BasicFiniteElements.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
BoneRemodeling::OpMassAndEnergyCalculation::massVec
Vec massVec
Definition: Remodeling.hpp:582
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getVolume
double getVolume() const
element volume (linear geometry)
Definition: VolumeElementForcesAndSourcesCore.hpp:161
ElasticMaterials
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
Definition: ElasticMaterials.hpp:24
MoFEM::DataOperator::doTets
bool & doTets
\deprectaed
Definition: DataOperators.hpp:94
BoneRemodeling::OpAssmbleRhoLhs_dRho::OpAssmbleRhoLhs_dRho
OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:961
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
BoneRemodeling::Remodeling::CommonData::DataContainers::rhoPtr
boost::shared_ptr< VectorDouble > rhoPtr
Ptr to density matrix container.
Definition: Remodeling.hpp:187
BoneRemodeling::OpCalculateStressTangent::OpCalculateStressTangent
OpCalculateStressTangent(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:717
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
Remodeling.hpp
MoFEM::DataOperator::doVertices
bool & doVertices
\deprectaed If false skip vertices
Definition: DataOperators.hpp:90
BoneRemodeling::MonitorPostProc::preProcess
MoFEMErrorCode preProcess()
Definition: Remodeling.cpp:1394
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
BoneRemodeling::Remodeling::CommonData::lambda
double lambda
Lame parameter.
Definition: Remodeling.hpp:144
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
FTensor::Tensor2< double *, 3, 3 >
BoneRemodeling::OpCalculateStressTangentWithAdolc::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:396
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
BoneRemodeling::Remodeling::CommonData::DataContainers::gradDispPtr
boost::shared_ptr< MatrixDouble > gradDispPtr
Ptr to gradient of displacements matrix container.
Definition: Remodeling.hpp:185
BoneRemodeling::OpMassAndEnergyCalculation::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:583
BoneRemodeling::OpAssmbleRhoRhs
Assemble residual for conservation of mass (density)
Definition: Remodeling.hpp:453
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
ROW
@ ROW
Definition: definitions.h:136
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::invertTensor3by3
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:1588
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::TsMonitorSet
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:259
MetaNodalForces::setOperators
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:128
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1236
BoneRemodeling::MonitorPostProc
Definition: Remodeling.hpp:553
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
BoneRemodeling::Remodeling::CommonData::less_post_proc
PetscBool less_post_proc
reduce file size
Definition: Remodeling.hpp:163
BoneRemodeling::Remodeling::CommonData::with_adol_c
PetscBool with_adol_c
Definition: Remodeling.hpp:161
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
BoneRemodeling::Remodeling::CommonData::tEts_block
Range tEts_block
Definition: Remodeling.hpp:166
PostProcStress
Definition: PostProcStresses.hpp:17
BoneRemodeling::OpCalculateStressTangent::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:407
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
BoneRemodeling::MonitorPostProc::MonitorPostProc
MonitorPostProc(MoFEM::Interface &m_field, Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:1387
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
BoneRemodeling::Remodeling::CommonData::bitLevel
BitRefLevel bitLevel
Definition: Remodeling.hpp:123
BoneRemodeling::Remodeling::getParameters
MoFEMErrorCode getParameters()
Get parameters form line command or config file.
Definition: Remodeling.cpp:1579
SYMMETRIC_TENSOR2_MAT_PTR
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT)
Definition: Remodeling.cpp:54
BoneRemodeling::OpAssmbleRhoLhs_dRho::transLocK_rho_rho
MatrixDouble transLocK_rho_rho
Definition: Remodeling.hpp:477
BoneRemodeling::OpCalculateStressTangentWithAdolc
Definition: Remodeling.hpp:415
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
BoneRemodeling::OpCalculateStress::OpCalculateStress
OpCalculateStress(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:280
BoneRemodeling::Remodeling::CommonData::mu
double mu
Lame parameter.
Definition: Remodeling.hpp:145
BoneRemodeling::Remodeling::CommonData::D
Vec D
Definition: Remodeling.hpp:120
BoneRemodeling::OpAssmbleRhoLhs_dRho::locK_rho_rho
MatrixDouble locK_rho_rho
Definition: Remodeling.hpp:476
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
a
constexpr double a
Definition: approx_sphere.cpp:30
BoneRemodeling::Remodeling::CommonData::is_atom_testing
PetscBool is_atom_testing
for atom tests
Definition: Remodeling.hpp:162
BoneRemodeling::MonitorPostProc::postProcElastic
PostProcVolumeOnRefinedMesh postProcElastic
Definition: Remodeling.hpp:558
TimeForceScale
Force scale operator for reading two columns.
Definition: TimeForceScale.hpp:18
R
@ R
Definition: free_surface.cpp:394
BoneRemodeling::Remodeling::CommonData::elasticPtr
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Definition: Remodeling.hpp:140
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
BoneRemodeling::OpAssmbleRhoRhs::OpAssmbleRhoRhs
OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:883
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
BoneRemodeling::OpCalculateStressTangentWithAdolc::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:418
BoneRemodeling::MonitorPostProc::equilibrium_flg
PetscBool equilibrium_flg
Definition: Remodeling.hpp:561
BoneRemodeling::Remodeling::CommonData::DataContainers::gradRhoPtr
boost::shared_ptr< MatrixDouble > gradRhoPtr
Gradient of density field.
Definition: Remodeling.hpp:188
convert.type
type
Definition: convert.py:64
PostProcTemplateVolumeOnRefinedMesh::commonData
CommonData commonData
Definition: PostProcOnRefMesh.hpp:287
BoneRemodeling::OpAssmbleRhoLhs_dF
Off diagonal block of tangent matrix .
Definition: Remodeling.hpp:493
BoneRemodeling::MonitorPostProc::postProcess
MoFEMErrorCode postProcess()
Definition: Remodeling.cpp:1404
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:197
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
BoneRemodeling::Remodeling::CommonData::edgeForces
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
Definition: Remodeling.hpp:138
BoneRemodeling
Definition: DensityMaps.hpp:27
MoFEM::TsCtx
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:17
BoneRemodeling::Remodeling::CommonData::DataContainers::vecPsi
VectorDouble vecPsi
Elastic energy density.
Definition: Remodeling.hpp:189
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
BoneRemodeling::MonitorPostProc::rate
double rate
Definition: Remodeling.hpp:562
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
BoneRemodeling::MonitorPostProc::pRT
int pRT
Definition: Remodeling.hpp:564
BoneRemodeling::Remodeling::CommonData::F
Vec F
Definition: Remodeling.hpp:120
BoneRemodeling::Remodeling::Fe
Volume finite element.
Definition: Remodeling.hpp:41
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
PostProcTemplateVolumeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Definition: PostProcOnRefMesh.hpp:301
PostProcTemplateOnRefineMesh::addFieldValuesGradientPostProc
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
Definition: PostProcOnRefMesh.hpp:195
BoneRemodeling::OpAssmbleStressRhs::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:438
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
BoneRemodeling::Remodeling::addMomentumFluxes
MoFEMErrorCode addMomentumFluxes()
Finite elements to calculate tractions.
Definition: Remodeling.cpp:1792
MoFEM::DataOperator::doQuads
bool & doQuads
\deprectaed
Definition: DataOperators.hpp:92
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
BoneRemodeling::Remodeling::CommonData::getCFromDensity
double getCFromDensity(const double &rho)
Definition: Remodeling.cpp:218
BoneRemodeling::OpAssmbleStressRhs
Assemble residual for conservation of momentum (stresses)
Definition: Remodeling.hpp:435
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
ElasticMaterials.hpp
Elastic materials.
MoFEM::DMMoFEMCreateMoFEM
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: DMMoFEM.cpp:114
BoneRemodeling::Remodeling::commonData
CommonData & commonData
Definition: Remodeling.hpp:207
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
BoneRemodeling::Remodeling::CommonData::DataContainers::matGradientOfDeformation
MatrixDouble matGradientOfDeformation
Gradient of deformation.
Definition: Remodeling.hpp:198
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
BoneRemodeling::OpCalculateStress::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:294
BoneRemodeling::OpAssmbleStressLhs_dRho::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:539
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEMethod
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Definition: ForcesAndSourcesCore.hpp:1042
BoneRemodeling::OpAssmbleStressRhs::OpAssmbleStressRhs
OpAssmbleStressRhs(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:814
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
convert.n
n
Definition: convert.py:82
MetaEdgeForces::addElement
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:62
BoneRemodeling::OpCalculateStress
Evaluate physical equations at integration points.
Definition: Remodeling.hpp:356
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
BoneRemodeling::OpCalculateStressTangentWithAdolc::OpCalculateStressTangentWithAdolc
OpCalculateStressTangentWithAdolc(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:382
BoneRemodeling::Remodeling::solveDM
MoFEMErrorCode solveDM()
Solve problem set up in DM.
Definition: Remodeling.cpp:1863
MoFEM::DataOperator::doEdges
bool & doEdges
\deprectaed If false skip edges
Definition: DataOperators.hpp:91
BoneRemodeling::OpGetRhoTimeDirevative
Evaluate density derivative with respect to time in case of Backward Euler Method.
Definition: Remodeling.hpp:332
BoneRemodeling::Remodeling::CommonData::DataContainers::matPushedMaterialTangent
MatrixDouble matPushedMaterialTangent
Definition: Remodeling.hpp:196
BoneRemodeling::OpAssmbleStressLhs_dRho::doWork
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:1316
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
BoneRemodeling::OpCalculateStress::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:359
BoneRemodeling::Remodeling::CommonData::pSi_ref
double pSi_ref
reference free energy
Definition: Remodeling.hpp:154
BoneRemodeling::Remodeling::CommonData::feLhs
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
Definition: Remodeling.hpp:130
BoneRemodeling::Remodeling::CommonData::R0
double R0
mass conduction coefficient
Definition: Remodeling.hpp:155
BoneRemodeling::OpPostProcStress::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: Remodeling.hpp:379
BoneRemodeling::Remodeling::CommonData
Definition: Remodeling.hpp:117
BoneRemodeling::OpAssmbleStressLhs_dRho
Off diagonal block of tangent matrix /f[ K_{u \rho}=\intop_{V} \left[\frac{n}{\rho_{0}}\right] \left...
Definition: Remodeling.hpp:536
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
BoneRemodeling::Remodeling::addElements
MoFEMErrorCode addElements()
Set and add finite elements.
Definition: Remodeling.cpp:1694
BoneRemodeling::MonitorPostProc::postProc
PostProcVolumeOnRefinedMesh postProc
Definition: Remodeling.hpp:557
std
Definition: enable_if.hpp:5
BoneRemodeling::MonitorPostProc::iNit
bool iNit
Definition: Remodeling.hpp:563
PostProcTemplateOnRefineMesh::mapGaussPts
std::vector< EntityHandle > mapGaussPts
Definition: PostProcOnRefMesh.hpp:125
BoneRemodeling::MonitorPostProc::mField
MoFEM::Interface & mField
Definition: Remodeling.hpp:555
BoneRemodeling::Remodeling::CommonData::DataContainers::matP
MatrixDouble matP
1st Piola stress
Definition: Remodeling.hpp:191
BoneRemodeling::Remodeling::CommonData::elasticMaterialsPtr
boost::shared_ptr< ElasticMaterials > elasticMaterialsPtr
Definition: Remodeling.hpp:141
BoneRemodeling::Remodeling::CommonData::getCFromDensityDiff
double getCFromDensityDiff(const double &rho)
Definition: Remodeling.cpp:225
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FTensor::Ddg
Definition: Ddg_value.hpp:7
lapack_dsyev
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:261
TENSOR2_MAT_PTR
#define TENSOR2_MAT_PTR(MAT)
Definition: Remodeling.cpp:50
BoneRemodeling::Remodeling::CommonData::oRder
int oRder
Definition: Remodeling.hpp:122
BoneRemodeling::Remodeling::CommonData::dm
DM dm
Discretization manager.
Definition: Remodeling.hpp:127
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
BoneRemodeling::Remodeling::CommonData::tEts_all
Range tEts_all
Definition: Remodeling.hpp:165
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
BoneRemodeling::OpMassAndEnergyCalculation::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: Remodeling.cpp:1947
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:426
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
BoneRemodeling::OpAssmbleRhoRhs::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:456
MoFEM::Types::MatrixDouble3by3
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
MoFEM::DMMoFEMGetTsCtx
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1142
BoneRemodeling::Remodeling::CommonData::DataContainers::matInvF
MatrixDouble matInvF
inverse of deformation gradient
Definition: Remodeling.hpp:199
MoFEM::DeprecatedCoreInterface::synchronise_field_entities
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
Definition: DeprecatedCoreInterface.cpp:47
BoneRemodeling::Remodeling::CommonData::m
double m
algorithmic exponent [-]
Definition: Remodeling.hpp:147
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
BoneRemodeling::OpPostProcStress::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:377
BoneRemodeling::OpCalculateStressTangentWithAdolc::vecC
VectorDouble vecC
Definition: Remodeling.hpp:419
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
BoneRemodeling::Remodeling::CommonData::A
Mat A
Definition: Remodeling.hpp:119
BoneRemodeling::OpAssmbleRhoLhs_dF::OpAssmbleRhoLhs_dF
OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:1066
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
BoneRemodeling::Remodeling::CommonData::DataContainers::matS
MatrixDouble matS
2nd Piola stress
Definition: Remodeling.hpp:190
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
BoneRemodeling::Remodeling::CommonData::tAg
const int tAg
Definition: Remodeling.hpp:178
BoneRemodeling::Remodeling::CommonData::DataContainers::vecRhoDt
VectorDouble vecRhoDt
Time derivative of density.
Definition: Remodeling.hpp:201
SYMMETRIC_TENSOR2_VEC_PTR
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)
Definition: Remodeling.cpp:57
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
BoneRemodeling::Remodeling::CommonData::DataContainers::matMaterialTangent
MatrixDouble matMaterialTangent
Definition: Remodeling.hpp:193
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
BoneRemodeling::OpAssmbleStressLhs_dF::doWork
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:1153
MoFEM::CoreInterface::check_finite_element
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
BoneRemodeling::OpPostProcStress::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:567
MoFEM::CoreInterface::set_field_order
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.
BoneRemodeling::OpAssmbleStressLhs_dF::OpAssmbleStressLhs_dF
OpAssmbleStressLhs_dF(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:1145
BoneRemodeling::Remodeling::CommonData::c
double c
density evolution (growth) velocity [d/m^2]
Definition: Remodeling.hpp:146
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
BoneRemodeling::Remodeling::CommonData::data
DataContainers data
Definition: Remodeling.hpp:203
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEM::DataOperator::sYmm
bool sYmm
If true assume that matrix is symmetric structure.
Definition: DataOperators.hpp:82
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
BoneRemodeling::OpGetRhoTimeDirevative::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:335
BoneRemodeling::OpMassAndEnergyCalculation
Definition: Remodeling.hpp:578
BoneRemodeling::OpAssmbleStressLhs_dRho::locK_P_RHO
MatrixDouble locK_P_RHO
Definition: Remodeling.hpp:540
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
BoneRemodeling::OpAssmbleRhoLhs_dF::doWork
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:1074
BoneRemodeling::OpAssmbleRhoRhs::nF
VectorDouble nF
Vector of the right hand side (internal forces)
Definition: Remodeling.hpp:457
BoneRemodeling::OpPostProcStress::OpPostProcStress
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:550
BoneRemodeling::MonitorPostProc::commonData
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:556
MoFEM::CoreInterface::add_field
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.
BoneRemodeling::OpPostProcStress
Used to post proc stresses, energy, source term.
Definition: Remodeling.hpp:374
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MU
#define MU(E, NU)
Definition: fem_tools.h:23
F
@ F
Definition: free_surface.cpp:394
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
MetaEdgeForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name, std::string mesh_node_positions="MESH_NODE_POSITIONS")
Set integration point operators.
Definition: EdgeForce.hpp:97
BoneRemodeling::OpCalculateStressTangent
Definition: Remodeling.hpp:404
PostProcTemplateOnRefineMesh::addFieldValuesPostProc
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
Definition: PostProcOnRefMesh.hpp:153
Generated by Doxygen 1.8.17 and hosted at University of Glasgow