v0.15.0
Loading...
Searching...
No Matches
Namespaces | Macros
Remodeling.cpp File Reference
#include <boost/program_options.hpp>
#include <BasicFiniteElements.hpp>
#include <ElasticMaterials.hpp>
#include <Remodeling.hpp>

Go to the source code of this file.

Namespaces

namespace  BoneRemodeling
 

Macros

#define TENSOR1_VEC_PTR(VEC)   &VEC[0], &VEC[1], &VEC[2]
 
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT)
 
#define TENSOR4_MAT_PTR(MAT)   &MAT(0, 0), MAT.size2()
 
#define TENSOR2_MAT_PTR(MAT)
 
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT)    &MAT(0, 0), &MAT(0, 1), &MAT(0, 2), &MAT(0, 3), &MAT(0, 4), &MAT(0, 5)
 
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)    &VEC[0], &VEC[1], &VEC[2], &VEC[3], &VEC[4], &VEC[5]
 

Macro Definition Documentation

◆ SYMMETRIC_TENSOR2_MAT_PTR

#define SYMMETRIC_TENSOR2_MAT_PTR (   MAT)     &MAT(0, 0), &MAT(0, 1), &MAT(0, 2), &MAT(0, 3), &MAT(0, 4), &MAT(0, 5)

Definition at line 54 of file Remodeling.cpp.

59 {
60
61MoFEMErrorCode Remodeling::CommonData::getParameters() {
62
63 PetscBool flg_config = PETSC_FALSE;
64 is_atom_testing = PETSC_FALSE;
65 with_adol_c = PETSC_FALSE;
66 char my_config_file_name[255];
67 double young_modulus = 5.0e+9; // Set here some typical value for bone 5 GPa
68 double poisson_ratio = 0.3; // Set here some typical value for bone
69 c = 4.0e-5; // Set here some typical value for bone. days/ m2 // this
70 // parameters governs how fast we can achieve equilibrium
71 m = 3.25; // Set some parameter typical for bone
72 n = 2.25; // Set Some parameter typical for bone
73 rHo_ref = 1000; // Set Some parameter typical for bone. kg/m3
74 pSi_ref = 2.75e4; // Set Some parameter typical for bone. energy / volume unit
75 // J/m3 = Pa
76 rHo_max = rHo_ref + 0.5 * rHo_ref;
77 rHo_min = rHo_ref - 0.5 * rHo_ref;
78 b = 0;
79 R0 = 0.0; // Set Some parameter typical for bone.
80 cUrrent_psi = 0.0;
81 cUrrent_mass = 0.0;
82 nOremodellingBlock = false;
83 less_post_proc = PETSC_FALSE;
84
86 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling parameters",
87 "none");
88
89 // config file
90 CHKERR PetscOptionsString("-my_config", "configuration file name", "",
91 "my_config.in", my_config_file_name, 255,
92 &flg_config);
93
94 if (flg_config) {
95 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Config file: %s loaded.\n",
96 my_config_file_name);
97 try {
98 ifstream ini_file(my_config_file_name);
99 if (!ini_file.is_open()) {
100 SETERRQ(PETSC_COMM_SELF, 1, "*** -my_config does not exist ***");
101 }
102 po::variables_map vm;
103 po::options_description config_file_options;
104 // std::cout << my_config_file_name << std::endl;
105
106 config_file_options.add_options()(
107 "young_modulus", po::value<double>(&young_modulus)->default_value(1));
108 config_file_options.add_options()(
109 "poisson_ratio", po::value<double>(&poisson_ratio)->default_value(1));
110 config_file_options.add_options()(
111 "c", po::value<double>(&c)->default_value(1));
112 config_file_options.add_options()(
113 "m", po::value<double>(&m)->default_value(1));
114 config_file_options.add_options()(
115 "n", po::value<double>(&n)->default_value(1));
116 config_file_options.add_options()(
117 "rHo_ref", po::value<double>(&rHo_ref)->default_value(1));
118 config_file_options.add_options()(
119 "pSi_ref", po::value<double>(&pSi_ref)->default_value(1));
120 config_file_options.add_options()(
121 "R0", po::value<double>(&R0)->default_value(1));
122
123 po::parsed_options parsed =
124 parse_config_file(ini_file, config_file_options, true);
125 store(parsed, vm);
126 po::notify(vm);
127
128 } catch (const std::exception &ex) {
129 std::ostringstream ss;
130 ss << ex.what() << std::endl;
131 SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
132 }
133 }
134
135 CHKERR PetscOptionsScalar("-young_modulus", "get young modulus", "",
136 young_modulus, &young_modulus, PETSC_NULL);
137
138 CHKERR PetscOptionsScalar("-poisson_ratio", "get poisson_ratio", "",
139 poisson_ratio, &poisson_ratio, PETSC_NULL);
140
141 CHKERR PetscOptionsScalar("-c", "density evolution (growth) velocity [d/m^2]",
142 "", c, &c, PETSC_NULL);
143
144 CHKERR PetscOptionsScalar("-m", "algorithmic exponent", "", m, &m,
145 PETSC_NULL);
146
147 CHKERR PetscOptionsScalar("-n", "porosity exponent", "", n, &n, PETSC_NULL);
148
149 CHKERR PetscOptionsInt("-b", "bell function exponent", "", b, &b, PETSC_NULL);
150
151 CHKERR PetscOptionsScalar("-rho_ref", "reference bone density", "", rHo_ref,
152 &rHo_ref, PETSC_NULL);
153
154 CHKERR PetscOptionsScalar("-rho_max", "reference bone density", "", rHo_max,
155 &rHo_max, PETSC_NULL);
156 CHKERR PetscOptionsScalar("-rho_min", "reference bone density", "", rHo_min,
157 &rHo_min, PETSC_NULL);
158
159 CHKERR PetscOptionsScalar("-psi_ref", "reference energy density", "", pSi_ref,
160 &pSi_ref, PETSC_NULL);
161
162 CHKERR PetscOptionsScalar("-r0", "mass source", "", R0, &R0, PETSC_NULL);
163
164 CHKERR PetscOptionsBool("-my_is_atom_test",
165 "is used with testing, exit with error when diverged",
166 "", is_atom_testing, &is_atom_testing, PETSC_NULL);
167
168 CHKERR PetscOptionsBool("-less_post_proc",
169 "is used to reduce output file size", "",
170 less_post_proc, &less_post_proc, PETSC_NULL);
171
172 CHKERR PetscOptionsBool(
173 "-with_adolc",
174 "calculate the material tangent with automatic differentiation", "",
175 with_adol_c, &with_adol_c, PETSC_NULL);
176
179
180 CHKERR PetscPrintf(PETSC_COMM_WORLD,
181 "Young's modulus E[Pa]: %4.3g\n",
183 CHKERR PetscPrintf(PETSC_COMM_WORLD,
184 "Poisson ratio nu[-] %4.3g\n",
186 CHKERR PetscPrintf(PETSC_COMM_WORLD,
187 "Lame coefficient lambda[Pa]: %4.3g\n", lambda);
188 CHKERR PetscPrintf(PETSC_COMM_WORLD,
189 "Lame coefficient mu[Pa]: %4.3g\n", mu);
190 CHKERR PetscPrintf(PETSC_COMM_WORLD,
191 "Density growth velocity [d/m2] %4.3g\n", c);
192 CHKERR PetscPrintf(PETSC_COMM_WORLD,
193 "Algorithmic exponent m[-]: %4.3g\n", m);
194 CHKERR PetscPrintf(PETSC_COMM_WORLD,
195 "Porosity exponent n[-]: %4.3g\n", n);
196 CHKERR PetscPrintf(PETSC_COMM_WORLD,
197 "Reference density rHo_ref[kg/m3]: %4.3g\n", rHo_ref);
198 CHKERR PetscPrintf(PETSC_COMM_WORLD,
199 "Reference energy pSi_ref[Pa]: %4.3g\n", pSi_ref);
200 CHKERR PetscPrintf(PETSC_COMM_WORLD,
201 "Mass conduction R0[kg/m3s]: %4.3g\n", R0);
202 CHKERR PetscPrintf(PETSC_COMM_WORLD,
203 "Bell function coefficent b[-]: %4.3g\n", b);
204 if (b != 0) {
205 CHKERR PetscPrintf(PETSC_COMM_WORLD,
206 "Max density rHo_max[kg/m3]: %4.3g\n", rHo_max);
207 CHKERR PetscPrintf(PETSC_COMM_WORLD,
208 "Min density rHo_min[kg/m3]: %4.3g\n", rHo_min);
209 }
210
211 PetscOptionsEnd();
212
214}
215
216inline double Remodeling::CommonData::getCFromDensity(const double &rho) {
217 double mid_rho = 0.5 * (rHo_max + rHo_min);
218 double var_h = (rho - mid_rho) / (rHo_max - mid_rho);
219
220 return 1 / (1 + (b != 0) * pow(var_h, 2 * b));
221}
222
223inline double Remodeling::CommonData::getCFromDensityDiff(const double &rho) {
224 double mid_rho = 0.5 * (rHo_max + rHo_min);
225 double var_h = (rho - mid_rho) / (rHo_max - mid_rho);
226 return (b != 0) * (-2) * b * pow(var_h, 2 * b - 1) /
227 ((rHo_max - mid_rho) * pow(pow(var_h, 2 * b) + 1, 2));
228}
229
230OpGetRhoTimeDirevative::OpGetRhoTimeDirevative(
231 Remodeling::CommonData &common_data)
233 "RHO", UserDataOperator::OPROW),
234 commonData(common_data) {}
235
237OpGetRhoTimeDirevative::doWork(int side, EntityType type,
238 DataForcesAndSourcesCore::EntData &data) {
239
241
242 const int nb_gauss_pts = data.getN().size1();
243 if (!nb_gauss_pts)
245 const int nb_base_functions = data.getN().size2();
246
247 auto base_function = data.getFTensor0N();
248 commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
249 if (type == MBVERTEX) {
250 commonData.data.vecRhoDt.clear();
251 }
252
253 int nb_dofs = data.getIndices().size();
254 if (nb_dofs == 0)
256 double *array;
257 CHKERR VecGetArray(getFEMethod()->ts_u_t, &array);
258 auto drho_dt = getFTensor0FromVec(commonData.data.vecRhoDt);
259
260 for (int gg = 0; gg < nb_gauss_pts; gg++) {
261 int bb = 0;
262 for (; bb < nb_dofs; bb++) {
263 const double field_data = array[data.getLocalIndices()[bb]];
264 drho_dt += field_data * base_function;
265 ++base_function;
266 }
267 // It is possible to have more base functions than dofs
268 for (; bb != nb_base_functions; bb++)
269 ++base_function;
270 ++drho_dt;
271 }
272
273 CHKERR VecRestoreArray(getFEMethod()->ts_u_t, &array);
274
276}
277
278OpCalculateStress::OpCalculateStress(Remodeling::CommonData &common_data)
280 "DISPLACEMENTS", UserDataOperator::OPROW),
281 commonData(common_data) {
282 // Set that this operator is executed only for vertices on element
283 doVertices = true;
284 doEdges = false;
285 doQuads = false;
286 doTris = false;
287 doTets = false;
288 doPrisms = false;
289}
290
292OpCalculateStress::doWork(int side, EntityType type,
293 DataForcesAndSourcesCore::EntData &data) {
294
296
297 if (type != MBVERTEX)
299
300 if (!commonData.data.gradDispPtr) {
301 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
302 "Gradient at integration points of displacement is not calculated");
303 }
304 auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
305
306 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
307 if (!commonData.data.rhoPtr) {
308 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
309 "Density at integration points is not calculated");
310 }
311
312 auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
313
314 commonData.data.vecPsi.resize(nb_gauss_pts, false);
315 auto psi = getFTensor0FromVec(commonData.data.vecPsi);
316 commonData.data.vecR.resize(nb_gauss_pts, false);
317 auto R = getFTensor0FromVec(commonData.data.vecR);
318
319 commonData.data.matInvF.resize(9, nb_gauss_pts, false);
320 auto invF = getFTensor2FromMat<3, 3>(commonData.data.matInvF);
321 commonData.data.vecDetF.resize(nb_gauss_pts, false);
322 auto det_f = getFTensor0FromVec(commonData.data.vecDetF);
323 MatrixDouble &matS = commonData.data.matS;
324 matS.resize(nb_gauss_pts, 6, false);
326 commonData.data.matP.resize(9, nb_gauss_pts, false);
327 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
328 MatrixDouble &matF = commonData.data.matGradientOfDeformation;
329 matF.resize(9, nb_gauss_pts, false);
331
332 const double c = commonData.c;
333 const double n = commonData.n;
334 const double m = commonData.m;
335 const double rho_ref = commonData.rHo_ref;
336 const double psi_ref = commonData.pSi_ref;
337 const double lambda = commonData.lambda;
338 const double mu = commonData.mu;
339 FTensor::Index<'i', 3> i;
340 FTensor::Index<'I', 3> I;
341 FTensor::Index<'J', 3> J;
342
343 for (int gg = 0; gg != nb_gauss_pts; gg++) {
344
345 // Get deformation gradient
346 F(i, I) = grad(i, I);
347 F(0, 0) += 1;
348 F(1, 1) += 1;
349 F(2, 2) += 1;
350
352 CHKERR invertTensor3by3(F, det_f, invF);
353
354 const double log_det = log(det_f);
355 psi = F(i, J) * F(i, J) - 3 - 2 * log_det;
356 psi *= 0.5 * mu;
357 psi += 0.5 * lambda * log_det * log_det;
358
359 const double coef = lambda * log_det - mu;
360 P(i, J) = mu * F(i, J) + coef * invF(J, i);
361 S(i, I) = invF(i, J) ^ P(J, I);
362
363 const double kp = commonData.getCFromDensity(rho);
364 R = c * kp * (pow(rho / rho_ref, n - m) * psi - psi_ref);
365
366 ++det_f;
367 ++invF;
368 ++grad;
369 ++rho;
370 ++psi;
371 ++S;
372 ++P;
373 ++R;
374 ++F;
375 }
376
378}
379
380OpCalculateStressTangentWithAdolc::OpCalculateStressTangentWithAdolc(
381 Remodeling::CommonData &common_data)
383 "DISPLACEMENTS", UserDataOperator::OPROW),
384 commonData(common_data) {
385 // Set that this operator is executed only for vertices on element
386 doVertices = true;
387 doEdges = false;
388 doQuads = false;
389 doTris = false;
390 doTets = false;
391 doPrisms = false;
392}
393
394MoFEMErrorCode OpCalculateStressTangentWithAdolc::doWork(
395 int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
396
398
399 if (type != MBVERTEX)
401
402 auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
403 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
404
405 vecC.resize(6, false);
407 MatrixDouble &dS_dC = commonData.data.matMaterialTangent;
408 dS_dC.resize(6, 6 * nb_gauss_pts, false);
409 dS_dC.clear();
411
412 MatrixDouble &matS = commonData.data.matS;
413 matS.resize(nb_gauss_pts, 6, false);
415 MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
416 dP_dF.resize(81, nb_gauss_pts, false);
417 MatrixDouble &matF = commonData.data.matGradientOfDeformation;
418 matF.resize(9, nb_gauss_pts, false);
421 commonData.data.matP.resize(9, nb_gauss_pts, false);
422 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
423 commonData.data.vecPsi.resize(nb_gauss_pts, false);
424 auto psi = getFTensor0FromVec(commonData.data.vecPsi);
425
426 FTensor::Index<'i', 3> i;
427 FTensor::Index<'k', 3> k;
428 FTensor::Index<'I', 3> I;
429 FTensor::Index<'J', 3> J;
430 FTensor::Index<'K', 3> K;
431 FTensor::Index<'L', 3> L;
432
433 for (int gg = 0; gg != nb_gauss_pts; gg++) {
434
435 // Get deformation gradient
436 F(i, I) = grad(i, I);
437 F(0, 0) += 1;
438 F(1, 1) += 1;
439 F(2, 2) += 1;
440
441 // Calculate Green defamation Tensor0
442 // ^ that means multiplication of Tensor 2 which result in symmetric Tensor
443 // rank 2.
444 C(I, J) = F(i, I) ^ F(i, J);
445
446 // Calculate free energy
447 CHKERR recordFreeEnergy_dC<FTensor::Tensor0<FTensor::PackPtr<double *, 1>>,
449 commonData, C, psi);
450
451 int r_g = ::gradient(commonData.tAg, 6, &vecC[0], &matS(gg, 0));
452 if (r_g < 0) {
453 // That means that energy function is not smooth and derivative
454 // can not be calculated,
455 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
456 "ADOL-C function evaluation with error");
457 }
458
459 // Scale energy density and 2nd Piola Stress
460 S(0, 0) *= 2;
461 S(1, 1) *= 2;
462 S(2, 2) *= 2;
463
464 // Calculate 1st Piola-Stress tensor
465 P(i, J) = F(i, I) * S(I, J);
466
467 const int is = 6 * gg;
468 double *hessian_ptr[6] = {&dS_dC(0, is), &dS_dC(1, is), &dS_dC(2, is),
469 &dS_dC(3, is), &dS_dC(4, is), &dS_dC(5, is)};
470 int r_h = ::hessian(commonData.tAg, 6, &vecC[0], hessian_ptr);
471 if (r_h < 0) {
472 // That means that energy function is not smooth and derivative
473 // can not be calculated,
474 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
475 "ADOL-C function evaluation with error");
476 }
477
478 // FIXME: TODO: Here tensor 4 with two minor symmetries is used, since
479 // \psi_{,ijkl} = \psi_{,klij} = \psi_{,lkji} = ... there is additional
480 // major symmetry. That is 21 independent components not 36 how is now
481 // implemented using Tenso4_ddg.
482 //
483 // That way we have to fill off-diagonal part for dS_dC. In future that
484 // need to be fixed by implementation of Tenso4 with additional major
485 // symmetries. Pleas consider doing such implementation by yourself. You
486 // have support of MoFEM developing team to guide you through this process.
487 //
488 // This deteriorate efficiency.
489 //
490 // TODO: Implement matrix with more symmetries.
491 //
492 for (int ii = 0; ii < 6; ii++) {
493 for (int jj = 0; jj < ii; jj++) {
494 dS_dC(jj, is + ii) = dS_dC(ii, is + jj);
495 }
496 }
497
498 // As result that symmetry of deformation gradient C is used, terms which
499 // have mixed indices are two time to big, and terms which has two times
500 // midxed index are four times to big. Those terms are not multiplied.
501 D1(0, 0, 0, 0) *= 4;
502 D1(1, 1, 1, 1) *= 4;
503 D1(2, 2, 2, 2) *= 4;
504
505 D1(0, 0, 1, 1) *= 4;
506 D1(1, 1, 0, 0) *= 4;
507 D1(1, 1, 2, 2) *= 4;
508 D1(2, 2, 1, 1) *= 4;
509 D1(0, 0, 2, 2) *= 4;
510 D1(2, 2, 0, 0) *= 4;
511
512 D1(0, 0, 0, 1) *= 2;
513 D1(0, 0, 0, 2) *= 2;
514 D1(0, 0, 1, 2) *= 2;
515 D1(0, 1, 0, 0) *= 2;
516 D1(0, 2, 0, 0) *= 2;
517 D1(1, 2, 0, 0) *= 2;
518
519 D1(1, 1, 0, 1) *= 2;
520 D1(1, 1, 0, 2) *= 2;
521 D1(1, 1, 1, 2) *= 2;
522 D1(0, 1, 1, 1) *= 2;
523 D1(0, 2, 1, 1) *= 2;
524 D1(1, 2, 1, 1) *= 2;
525
526 D1(2, 2, 0, 1) *= 2;
527 D1(2, 2, 0, 2) *= 2;
528 D1(2, 2, 1, 2) *= 2;
529 D1(0, 1, 2, 2) *= 2;
530 D1(0, 2, 2, 2) *= 2;
531 D1(1, 2, 2, 2) *= 2;
532
533 D2(i, J, k, L) = (F(i, I) * (D1(I, J, K, L) * F(k, K)));
534 // D2(i, J, k, L) += I(i, k) * S(J, L);
535
536 ++grad;
537 ++S;
538 ++F;
539 ++D1;
540 ++D2;
541 ++P;
542 ++psi;
543 }
544
546}
547
548OpPostProcStress::OpPostProcStress(moab::Interface &post_proc_mesh,
549 std::vector<EntityHandle> &map_gauss_pts,
550 Remodeling::CommonData &common_data)
552 "DISPLACEMENTS", UserDataOperator::OPROW),
553 commonData(common_data), postProcMesh(post_proc_mesh),
554 mapGaussPts(map_gauss_pts) {
555 // Set that this operator is executed only for vertices on element
556 doVertices = true;
557 doEdges = false;
558 doQuads = false;
559 doTris = false;
560 doTets = false;
561 doPrisms = false;
562}
563
565OpPostProcStress::doWork(int side, EntityType type,
566 DataForcesAndSourcesCore::EntData &data) {
567
569
570 if (type != MBVERTEX)
572 double def_VAL[9];
573 bzero(def_VAL, 9 * sizeof(double));
574
575 Tag th_piola2;
576 CHKERR postProcMesh.tag_get_handle("Piola_Stress2", 9, MB_TYPE_DOUBLE,
577 th_piola2, MB_TAG_CREAT | MB_TAG_SPARSE,
578 def_VAL);
579
580 Tag th_psi;
581 CHKERR postProcMesh.tag_get_handle("psi", 1, MB_TYPE_DOUBLE, th_psi,
582 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
583
584 Tag th_rho;
585 CHKERR postProcMesh.tag_get_handle("RHO", 1, MB_TYPE_DOUBLE, th_rho,
586 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
587
588 Tag th_piola1;
589 Tag principal;
590 Tag th_prin_stress_vect1, th_prin_stress_vect2, th_prin_stress_vect3;
591 Tag th_R;
592 Tag th_grad_rho;
593 if (!commonData.less_post_proc) {
594
595 CHKERR postProcMesh.tag_get_handle("Piola_Stress1", 9, MB_TYPE_DOUBLE,
596 th_piola1, MB_TAG_CREAT | MB_TAG_SPARSE,
597 def_VAL);
598
599 CHKERR postProcMesh.tag_get_handle("Principal_stress", 3, MB_TYPE_DOUBLE,
600 principal, MB_TAG_CREAT | MB_TAG_SPARSE,
601 def_VAL);
602
603 CHKERR postProcMesh.tag_get_handle("Principal_stress_S1", 3, MB_TYPE_DOUBLE,
604 th_prin_stress_vect1,
605 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
606 CHKERR postProcMesh.tag_get_handle("Principal_stress_S2", 3, MB_TYPE_DOUBLE,
607 th_prin_stress_vect2,
608 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
609 CHKERR postProcMesh.tag_get_handle("Principal_stress_S3", 3, MB_TYPE_DOUBLE,
610 th_prin_stress_vect3,
611 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
612
613 CHKERR postProcMesh.tag_get_handle("R", 1, MB_TYPE_DOUBLE, th_R,
614 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
615
616 CHKERR postProcMesh.tag_get_handle("grad_rho", 3, MB_TYPE_DOUBLE,
617 th_grad_rho,
618 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
619 }
620
621 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
622 MatrixDouble &matS = commonData.data.matS;
624 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
625 auto psi = getFTensor0FromVec(commonData.data.vecPsi);
626 auto R = getFTensor0FromVec(commonData.data.vecR);
627
628 auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
629 auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
630
631 MatrixDouble3by3 stress(3, 3);
632 VectorDouble3 eigen_values(3);
633 MatrixDouble3by3 eigen_vectors(3, 3);
634
635 for (int gg = 0; gg != nb_gauss_pts; gg++) {
636
637 double val = psi;
638 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &val);
639 val = rho;
640 CHKERR postProcMesh.tag_set_data(th_rho, &mapGaussPts[gg], 1, &val);
641
642 // S is symmetric tensor, need to map it to non-symmetric tensor to save it
643 // on moab tag to be viable in paraview.
644 for (int ii = 0; ii != 3; ii++) {
645 for (int jj = 0; jj != 3; jj++) {
646 stress(ii, jj) = S(ii, jj);
647 }
648 }
649 CHKERR postProcMesh.tag_set_data(th_piola2, &mapGaussPts[gg], 1,
650 &stress(0, 0));
651
652 if (!commonData.less_post_proc) {
653
654 for (int ii = 0; ii != 3; ii++) {
655 for (int jj = 0; jj != 3; jj++) {
656 stress(ii, jj) = P(ii, jj);
657 }
658 }
659 CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
660 &stress(0, 0));
661 val = R;
662 CHKERR postProcMesh.tag_set_data(th_R, &mapGaussPts[gg], 1, &val);
663 const double t2[] = {grad_rho(0), grad_rho(1), grad_rho(2)};
664 CHKERR postProcMesh.tag_set_data(th_grad_rho, &mapGaussPts[gg], 1, t2);
665
666 // Add calculation of eigen values, i.e. prinple stresses.
667 noalias(eigen_vectors) = stress;
668 VectorDouble3 prin_stress_vect1(3);
669 VectorDouble3 prin_stress_vect2(3);
670 VectorDouble3 prin_stress_vect3(3);
671 // LAPACK - eigenvalues and vectors. Applied twice for initial creates
672 // memory space
673 int n = 3, lda = 3, info, lwork = -1;
674 double wkopt;
675 info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
676 &(eigen_values.data()[0]), &wkopt, lwork);
677 if (info != 0)
678 SETERRQ1(PETSC_COMM_SELF, 1,
679 "something is wrong with lapack_dsyev info = %d", info);
680 lwork = (int)wkopt;
681 double work[lwork];
682 info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
683 &(eigen_values.data()[0]), work, lwork);
684 if (info != 0)
685 SETERRQ1(PETSC_COMM_SELF, 1,
686 "something is wrong with lapack_dsyev info = %d", info);
687
688 for (int ii = 0; ii < 3; ii++) {
689 prin_stress_vect1[ii] = eigen_vectors(0, ii);
690 prin_stress_vect2[ii] = eigen_vectors(1, ii);
691 prin_stress_vect3[ii] = eigen_vectors(2, ii);
692 }
693
694 CHKERR postProcMesh.tag_set_data(principal, &mapGaussPts[gg], 1,
695 &eigen_values[0]);
696 CHKERR postProcMesh.tag_set_data(th_prin_stress_vect1, &mapGaussPts[gg],
697 1, &*prin_stress_vect1.begin());
698 CHKERR postProcMesh.tag_set_data(th_prin_stress_vect2, &mapGaussPts[gg],
699 1, &*prin_stress_vect2.begin());
700 CHKERR postProcMesh.tag_set_data(th_prin_stress_vect3, &mapGaussPts[gg],
701 1, &*prin_stress_vect3.begin());
702 }
703
704 ++S;
705 ++P;
706 ++psi;
707 ++R;
708 ++rho;
709 ++grad_rho;
710 }
711
713}
714
715OpCalculateStressTangent::OpCalculateStressTangent(
716 Remodeling::CommonData &common_data)
718 "DISPLACEMENTS", UserDataOperator::OPROW),
719 commonData(common_data) {
720 // Set that this operator is executed only for vertices on element
721 doVertices = true;
722 doEdges = false;
723 doQuads = false;
724 doTris = false;
725 doTets = false;
726 doPrisms = false;
727}
728
730OpCalculateStressTangent::doWork(int side, EntityType type,
731 DataForcesAndSourcesCore::EntData &data) {
732
734
735 if (type != MBVERTEX)
737
738 auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
739 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
740
741 commonData.data.matInvF.resize(9, nb_gauss_pts, false);
742 auto invF = getFTensor2FromMat<3, 3>(commonData.data.matInvF);
743 commonData.data.vecDetF.resize(nb_gauss_pts, false);
744 auto det_f = getFTensor0FromVec(commonData.data.vecDetF);
745
746 MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
747 dP_dF.resize(81, nb_gauss_pts, false);
748 MatrixDouble &matF = commonData.data.matGradientOfDeformation;
749 matF.resize(9, nb_gauss_pts, false);
752 commonData.data.matP.resize(9, nb_gauss_pts, false);
753 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
754 commonData.data.vecPsi.resize(nb_gauss_pts, false);
755 auto psi = getFTensor0FromVec(commonData.data.vecPsi);
756 FTensor::Index<'i', 3> i;
757 FTensor::Index<'k', 3> k;
758 FTensor::Index<'I', 3> I;
759 FTensor::Index<'J', 3> J;
760 FTensor::Index<'K', 3> K;
761 FTensor::Index<'L', 3> L;
762 const double lambda = commonData.lambda;
763 const double mu = commonData.mu;
764
765 for (int gg = 0; gg != nb_gauss_pts; gg++) {
766
767 // Get deformation gradient
768 F(i, I) = grad(i, I);
769 F(0, 0) += 1;
770 F(1, 1) += 1;
771 F(2, 2) += 1;
772
774 CHKERR invertTensor3by3(F, det_f, invF);
775
776 const double log_det = log(det_f);
777 psi = F(i, J) * F(i, J) - 3 - 2 * log_det;
778 psi *= 0.5 * mu;
779 psi += 0.5 * lambda * log_det * log_det;
780
781 const double var = lambda * log_det - mu;
782 P(i, J) = mu * F(i, J) + var * invF(J, i);
783
784 const double coef = mu - lambda * log_det;
785 // TODO: This can be improved by utilising the symmetries of the D2 tensor
786
787 D2(i, J, k, L) =
788 invF(J, i) * (invF(L, k) * lambda) + invF(L, i) * (invF(J, k) * coef);
789
790 D2(0, 0, 0, 0) += mu;
791 D2(0, 1, 0, 1) += mu;
792 D2(0, 2, 0, 2) += mu;
793 D2(1, 0, 1, 0) += mu;
794 D2(1, 1, 1, 1) += mu;
795 D2(1, 2, 1, 2) += mu;
796 D2(2, 0, 2, 0) += mu;
797 D2(2, 1, 2, 1) += mu;
798 D2(2, 2, 2, 2) += mu;
799
800 ++det_f;
801 ++invF;
802 ++grad;
803 ++F;
804 ++D2;
805 ++P;
806 ++psi;
807 }
808
810}
811
812OpAssmbleStressRhs::OpAssmbleStressRhs(Remodeling::CommonData &common_data)
814 "DISPLACEMENTS", UserDataOperator::OPROW),
815 commonData(common_data) {}
816
818OpAssmbleStressRhs::doWork(int side, EntityType type,
819 DataForcesAndSourcesCore::EntData &data) {
820
822
823 const int nb_dofs = data.getIndices().size();
824 if (!nb_dofs)
826 nF.resize(nb_dofs, false);
827 nF.clear();
828
829 const int nb_gauss_pts = data.getN().size1();
830 const int nb_base_functions = data.getN().size2();
831 auto diff_base_functions = data.getFTensor1DiffN<3>();
832 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
833 auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
834
835 // MatrixDouble &matS = commonData.data.matS;
836 // MatrixDouble &matF = commonData.data.matGradientOfDeformation;
837 // FTensor::Tensor2<double*,3,3> F(TENSOR2_MAT_PTR(matF));
838 // FTensor::Tensor2_symmetric<double*,3> S(SYMMETRIC_TENSOR2_MAT_PTR(matS),6);
839
840 const double n = commonData.n;
841 const double rho_ref = commonData.rHo_ref;
842
843 FTensor::Index<'I', 3> I;
844 FTensor::Index<'J', 3> J;
845 FTensor::Index<'i', 3> i;
846
847 for (int gg = 0; gg != nb_gauss_pts; gg++) {
848
849 // Get volume and integration weight
850 double w = getVolume() * getGaussPts()(3, gg);
851 double a = w * pow(rho / rho_ref, n);
852
853 int bb = 0;
854 for (; bb != nb_dofs / 3; bb++) {
855 double *ptr = &nF[3 * bb];
856 FTensor::Tensor1<double *, 3> t1(ptr, &ptr[1], &ptr[2]);
857 t1(i) += a * P(i, I) * diff_base_functions(I);
858 // t1(i) += 0.5*a*diff_base_functions(J)*(F(i,I)*S(I,J));
859 // t1(i) += 0.5*a*diff_base_functions(I)*(F(i,J)*S(I,J));
860 ++diff_base_functions;
861 }
862 // Could be that we base more base functions than needed to approx. disp.
863 // field.
864 for (; bb != nb_base_functions; bb++) {
865 ++diff_base_functions;
866 }
867 ++P;
868 // ++S;
869 // ++F;
870 ++rho;
871 }
872
873 // Assemble internal force vector for time solver (TS)
875 getFEMethod()->ts_F, /// Get the right hand side vector for time solver
876 nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
877
879}
880
881OpAssmbleRhoRhs::OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
883 "RHO", UserDataOperator::OPROW),
884 commonData(common_data) {}
885
887OpAssmbleRhoRhs::doWork(int side, EntityType type,
888 DataForcesAndSourcesCore::EntData &data) {
889
891
892 const int nb_dofs = data.getIndices().size();
893 if (!nb_dofs)
895 nF.resize(nb_dofs, false);
896 nF.clear();
897
898 const int nb_gauss_pts = data.getN().size1();
899 const int nb_base_functions = data.getN().size2();
900 auto base_functions = data.getFTensor0N();
901 auto diff_base_functions = data.getFTensor1DiffN<3>();
902 if (!commonData.data.rhoPtr) {
903 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "rhoPtsr is null");
904 }
905 auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
906 if (!commonData.data.gradRhoPtr) {
907 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "gradRhoPtr is null");
908 }
909 auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
910 auto R = getFTensor0FromVec(commonData.data.vecR);
911
912 if (commonData.data.vecRhoDt.size() != nb_gauss_pts) {
913 commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
914 commonData.data.vecRhoDt.clear();
915 // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Time derivative of
916 // density at integration not calculated");
917 }
918 auto drho_dt = getFTensor0FromVec(commonData.data.vecRhoDt);
919
920 const double R0 = commonData.R0;
921
922 FTensor::Index<'I', 3> I;
923
924 for (int gg = 0; gg != nb_gauss_pts; gg++) {
925
926 // Get volume and integration weight
927 double w = getVolume() * getGaussPts()(3, gg);
928
929 int bb = 0;
930 for (; bb != nb_dofs; bb++) {
931 nF[bb] += w * base_functions * drho_dt; // Density rate term
932 nF[bb] -= w * base_functions * R; // Source term
933 nF[bb] -= w * R0 * diff_base_functions(I) * grad_rho(I); // Gradient term
934 ++base_functions;
935 ++diff_base_functions;
936 }
937 // Could be that we base more base functions than needed to approx. density
938 // field.
939 for (; bb != nb_base_functions; bb++) {
940 ++base_functions;
941 ++diff_base_functions;
942 }
943
944 // Increment quantities for integration pts.
945 ++rho;
946 ++grad_rho;
947 ++drho_dt;
948 ++R;
949 }
950
951 // Assemble internal force vector for time solver (TS)
953 getFEMethod()->ts_F, /// Get the right hand side vector for time solver
954 nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
955
957}
958
959OpAssmbleRhoLhs_dRho::OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
961 "RHO", "RHO", UserDataOperator::OPROWCOL),
962 commonData(common_data) {
963 sYmm = true;
964}
965
967OpAssmbleRhoLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
968 EntityType col_type,
969 DataForcesAndSourcesCore::EntData &row_data,
970 DataForcesAndSourcesCore::EntData &col_data) {
971
973
974 const int row_nb_dofs = row_data.getIndices().size();
975 if (!row_nb_dofs)
977 const int col_nb_dofs = col_data.getIndices().size();
978 if (!col_nb_dofs)
980
981 // Set size can clear local tangent matrix
982 locK_rho_rho.resize(row_nb_dofs, col_nb_dofs, false);
983 locK_rho_rho.clear();
984
985 const int row_nb_gauss_pts = row_data.getN().size1();
986 if (!row_nb_gauss_pts)
988 const int row_nb_base_functions = row_data.getN().size2();
989
990 const double ts_a = getFEMethod()->ts_a;
991 const double R0 = commonData.R0;
992 const double c = commonData.c;
993 const double n = commonData.n;
994 const double m = commonData.m;
995 const double rho_ref = commonData.rHo_ref;
996 const double psi_ref = commonData.pSi_ref;
997
998 FTensor::Index<'I', 3> I;
999
1000 auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1001 auto psi = getFTensor0FromVec(commonData.data.vecPsi);
1002 auto row_base_functions = row_data.getFTensor0N();
1003 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1004
1005 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1006
1007 // Get volume and integration weight
1008 double w = getVolume() * getGaussPts()(3, gg);
1009
1010 const double kp = commonData.getCFromDensity(rho);
1011 const double kp_diff = commonData.getCFromDensityDiff(rho);
1012 const double a = c * kp * ((n - m) / rho) * pow(rho / rho_ref, n - m) * psi;
1013 const double a_diff =
1014 c * kp_diff * (pow(rho / rho_ref, n - m) * psi - psi_ref);
1015
1016 int row_bb = 0;
1017 for (; row_bb != row_nb_dofs; row_bb++) {
1018 auto col_base_functions = col_data.getFTensor0N(gg, 0);
1019 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1020 for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1021 locK_rho_rho(row_bb, col_bb) +=
1022 ts_a * w * row_base_functions * col_base_functions;
1023 locK_rho_rho(row_bb, col_bb) -=
1024 w * R0 * row_diff_base_functions(I) * col_diff_base_functions(I);
1025 locK_rho_rho(row_bb, col_bb) -=
1026 a * w * row_base_functions * col_base_functions;
1027 locK_rho_rho(row_bb, col_bb) -=
1028 a_diff * w * row_base_functions * col_base_functions;
1029
1030 ++col_base_functions;
1031 ++col_diff_base_functions;
1032 }
1033 ++row_base_functions;
1034 ++row_diff_base_functions;
1035 }
1036 for (; row_bb != row_nb_base_functions; row_bb++) {
1037 ++row_base_functions;
1038 ++row_diff_base_functions;
1039 }
1040 ++psi;
1041 ++rho;
1042 }
1043
1044 // cerr << locK_rho_rho << endl;
1045
1046 CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1047 &*row_data.getIndices().begin(), col_nb_dofs,
1048 &*col_data.getIndices().begin(), &locK_rho_rho(0, 0),
1049 ADD_VALUES);
1050
1051 // is symmetric
1052 if (row_side != col_side || row_type != col_type) {
1053 transLocK_rho_rho.resize(col_nb_dofs, row_nb_dofs, false);
1054 noalias(transLocK_rho_rho) = trans(locK_rho_rho);
1055 CHKERR MatSetValues(getFEMethod()->ts_B, col_nb_dofs,
1056 &*col_data.getIndices().begin(), row_nb_dofs,
1057 &*row_data.getIndices().begin(),
1058 &transLocK_rho_rho(0, 0), ADD_VALUES);
1059 }
1060
1062}
1063
1064OpAssmbleRhoLhs_dF::OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data)
1066 "RHO", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
1067 commonData(common_data) {
1068 sYmm = false;
1069}
1070
1072OpAssmbleRhoLhs_dF::doWork(int row_side, int col_side, EntityType row_type,
1073 EntityType col_type,
1074 DataForcesAndSourcesCore::EntData &row_data,
1075 DataForcesAndSourcesCore::EntData &col_data) {
1076
1078
1079 const int row_nb_dofs = row_data.getIndices().size();
1080 if (!row_nb_dofs)
1082 const int col_nb_dofs = col_data.getIndices().size();
1083 if (!col_nb_dofs)
1085
1086 // Set size can clear local tangent matrix
1087 locK_rho_F.resize(row_nb_dofs, col_nb_dofs, false);
1088 locK_rho_F.clear();
1089
1090 const int row_nb_gauss_pts = row_data.getN().size1();
1091 const int col_nb_base_functions = col_data.getN().size2();
1092
1093 const double c = commonData.c;
1094 const double n = commonData.n;
1095 const double m = commonData.m;
1096 const double rho_ref = commonData.rHo_ref;
1097
1098 FTensor::Index<'I', 3> I;
1099 FTensor::Index<'i', 3> i;
1100
1101 auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1102 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
1103
1105 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>();
1106
1107 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1108
1109 // Get volume and integration weight
1110 double w = getVolume() * getGaussPts()(3, gg);
1111 const double kp = commonData.getCFromDensity(rho);
1112 const double a = w * c * kp * pow(rho / rho_ref, n - m);
1113
1114 int col_bb = 0;
1115 for (; col_bb != col_nb_dofs / 3; col_bb++) {
1116 t1(i) = a * P(i, I) * col_diff_base_functions(I);
1117 auto row_base_functions = row_data.getFTensor0N(gg, 0);
1118 for (int row_bb = 0; row_bb != row_nb_dofs; row_bb++) {
1119 FTensor::Tensor1<double *, 3> k(&locK_rho_F(row_bb, 3 * col_bb + 0),
1120 &locK_rho_F(row_bb, 3 * col_bb + 1),
1121 &locK_rho_F(row_bb, 3 * col_bb + 2));
1122 k(i) -= row_base_functions * t1(i);
1123 ++row_base_functions;
1124 }
1125 ++col_diff_base_functions;
1126 }
1127 for (; col_bb != col_nb_base_functions; col_bb++) {
1128 ++col_diff_base_functions;
1129 }
1130
1131 ++P;
1132 ++rho;
1133 }
1134
1135 CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1136 &*row_data.getIndices().begin(), col_nb_dofs,
1137 &*col_data.getIndices().begin(), &locK_rho_F(0, 0),
1138 ADD_VALUES);
1139
1141}
1142template <bool ADOLC>
1143OpAssmbleStressLhs_dF<ADOLC>::OpAssmbleStressLhs_dF(
1144 Remodeling::CommonData &common_data)
1146 "DISPLACEMENTS", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
1147 commonData(common_data) {
1148 sYmm = true;
1149}
1150template <bool ADOLC>
1151MoFEMErrorCode OpAssmbleStressLhs_dF<ADOLC>::doWork(
1152 int row_side, int col_side, EntityType row_type, EntityType col_type,
1153 DataForcesAndSourcesCore::EntData &row_data,
1154 DataForcesAndSourcesCore::EntData &col_data) {
1155
1157
1158 const int row_nb_dofs = row_data.getIndices().size();
1159 if (!row_nb_dofs)
1161 const int col_nb_dofs = col_data.getIndices().size();
1162 if (!col_nb_dofs)
1164
1165 const bool diagonal_block = (row_type == col_type) && (row_side == col_side);
1166
1167 // Set size can clear local tangent matrix
1168 locK_P_F.resize(row_nb_dofs, col_nb_dofs, false);
1169 locK_P_F.clear();
1170
1171 const int row_nb_gauss_pts = row_data.getN().size1();
1172 const int row_nb_base_functions = row_data.getN().size2();
1173
1174 MatrixDouble &matS = commonData.data.matS;
1175 MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
1176 // MatrixDouble &dS_dC = commonData.data.matMaterialTangent;
1177 // MatrixDouble &matF = commonData.data.matGradientOfDeformation;
1178
1179 double t0;
1180
1181 const double n = commonData.n;
1182 const double rho_ref = commonData.rHo_ref;
1183
1184 FTensor::Index<'i', 3> i;
1185 FTensor::Index<'j', 3> j;
1186 FTensor::Index<'k', 3> k;
1187 FTensor::Index<'I', 3> I;
1188 FTensor::Index<'J', 3> J;
1189 FTensor::Index<'K', 3> K;
1190 FTensor::Index<'L', 3> L;
1191
1192 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1195 auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1196
1197 // FTensor::Tensor4_ddg<double*,3,3> D1(SYMMETRIC_TENSOR4_MAT_PTR(dS_dC),6);
1198 // FTensor::Tensor2<double*,3,3> F(TENSOR2_MAT_PTR(matF));
1199
1200 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1201
1202 // Get volume and integration weight
1203 double w = getVolume() * getGaussPts()(3, gg);
1204 const double a = w * pow(rho / rho_ref, n);
1205
1206 int row_bb = 0;
1207 for (; row_bb != row_nb_dofs / 3; row_bb++) {
1208
1209 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1210
1211 const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
1212 int col_bb = 0;
1213 for (; col_bb != final_bb; col_bb++) {
1214
1216 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1217 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1218 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1219 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1220 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1221 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1222 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1223 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1224 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1225
1226 diffDiff(J, L) =
1227 a * row_diff_base_functions(J) * col_diff_base_functions(L);
1228
1229 // TODO: FIXME: Tensor D2 has major symmetry, we do not have such tensor
1230 // implemented at the moment. Using tensor with major symmetry will
1231 // improve code efficiency.
1232 t_assemble(i, k) += diffDiff(J, L) * D2(i, J, k, L);
1233
1234 if (ADOLC) {
1235 // Stress stiffening part (only with adolc since it is using dS / dC
1236 // derrivative. Check OpCalculateStressTangentWithAdolc operator)
1237 t0 = diffDiff(J, L) * S(J, L);
1238 t_assemble(0, 0) += t0;
1239 t_assemble(1, 1) += t0;
1240 t_assemble(2, 2) += t0;
1241 }
1242
1243 // Next base function for column
1244 ++col_diff_base_functions;
1245 }
1246
1247 ++row_diff_base_functions;
1248 }
1249 for (; row_bb != row_nb_base_functions; row_bb++) {
1250 ++row_diff_base_functions;
1251 }
1252
1253 // ++D1;
1254 ++D2;
1255 ++S;
1256 // ++F;
1257 ++rho;
1258 }
1259
1260 // Copy of symmetric part for assemble
1261 if (diagonal_block) {
1262 for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
1263 int col_bb = 0;
1264 for (; col_bb != row_bb + 1; col_bb++) {
1266 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1267 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1268 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1269 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1270 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1271 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1272 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1273 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1274 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1276 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 0),
1277 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 1),
1278 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 2),
1279 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 0),
1280 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 1),
1281 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 2),
1282 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 0),
1283 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 1),
1284 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 2));
1285 t_off_side(i, k) = t_assemble(k, i);
1286 }
1287 }
1288 }
1289
1290 const int *row_ind = &*row_data.getIndices().begin();
1291 const int *col_ind = &*col_data.getIndices().begin();
1292 CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs, row_ind, col_nb_dofs,
1293 col_ind, &locK_P_F(0, 0), ADD_VALUES);
1294
1295 if (row_type != col_type || row_side != col_side) {
1296 transLocK_P_F.resize(col_nb_dofs, row_nb_dofs, false);
1297 noalias(transLocK_P_F) = trans(locK_P_F);
1298 CHKERR MatSetValues(getFEMethod()->ts_B, col_nb_dofs, col_ind, row_nb_dofs,
1299 row_ind, &transLocK_P_F(0, 0), ADD_VALUES);
1300 }
1301
1303}
1304
1305OpAssmbleStressLhs_dRho::OpAssmbleStressLhs_dRho(
1306 Remodeling::CommonData &common_data)
1308 "DISPLACEMENTS", "RHO", UserDataOperator::OPROWCOL),
1309 commonData(common_data) {
1310 sYmm = false; // This is off-diagonal (upper-left), so no symmetric
1311}
1312
1314OpAssmbleStressLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
1315 EntityType col_type,
1316 DataForcesAndSourcesCore::EntData &row_data,
1317 DataForcesAndSourcesCore::EntData &col_data) {
1318
1320
1321 const int row_nb_dofs = row_data.getIndices().size();
1322 if (!row_nb_dofs)
1324 const int col_nb_dofs = col_data.getIndices().size();
1325 if (!col_nb_dofs)
1327
1328 // Set size can clear local tangent matrix
1329 locK_P_RHO.resize(row_nb_dofs, col_nb_dofs, false);
1330 locK_P_RHO.clear();
1331
1333
1334 const double n = commonData.n;
1335 const double rho_ref = commonData.rHo_ref;
1336
1337 FTensor::Index<'I', 3> I;
1338 FTensor::Index<'i', 3> i;
1339
1340 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
1341 auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1342
1343 const int row_nb_gauss_pts = row_data.getN().size1();
1344 const int row_nb_base_functions = row_data.getN().size2();
1345 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1346
1347 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1348
1349 // Get volume and integration weight
1350 double w = getVolume() * getGaussPts()(3, gg);
1351 double a = w * (n / rho) * pow(rho / rho_ref, n);
1352
1353 int row_bb = 0;
1354 for (; row_bb != row_nb_dofs / 3; row_bb++) {
1355 t1(i) = a * row_diff_base_functions(I) * P(i, I);
1356 auto base_functions = col_data.getFTensor0N(gg, 0);
1357 for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1358 // if(base_functions!=col_data.getN(gg)[col_bb]) {
1359 // cerr << "Error" << endl;
1360 // }
1361 FTensor::Tensor1<double *, 3> k(&locK_P_RHO(3 * row_bb + 0, col_bb),
1362 &locK_P_RHO(3 * row_bb + 1, col_bb),
1363 &locK_P_RHO(3 * row_bb + 2, col_bb));
1364 k(i) += t1(i) * base_functions;
1365 ++base_functions;
1366 }
1367 ++row_diff_base_functions;
1368 }
1369 for (; row_bb != row_nb_base_functions; row_bb++) {
1370 ++row_diff_base_functions;
1371 }
1372
1373 ++P;
1374 ++rho;
1375 }
1376
1377 CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1378 &*row_data.getIndices().begin(), col_nb_dofs,
1379 &*col_data.getIndices().begin(), &locK_P_RHO(0, 0),
1380 ADD_VALUES);
1381
1383}
1384
1386 Remodeling::CommonData &common_data)
1387 : mField(m_field), commonData(common_data), postProc(m_field),
1388 postProcElastic(m_field),
1389 // densityMaps(m_field),
1390 iNit(false) {}
1391
1395}
1396
1400}
1401
1404
1405 PetscBool mass_postproc = PETSC_FALSE;
1406 PetscBool equilibrium_flg = PETSC_FALSE;
1407 PetscBool analysis_mesh_flg = PETSC_FALSE;
1408 rate = 1;
1409
1410 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling post-process",
1411 "none");
1412
1413 CHKERR PetscOptionsBool(
1414 "-mass_postproc",
1415 "is used for testing, calculates mass and energy at each time step", "",
1416 mass_postproc, &mass_postproc, PETSC_NULL);
1417
1418 CHKERR PetscOptionsBool("-analysis_mesh",
1419 "saves analysis mesh at each time step", "",
1420 analysis_mesh_flg, &analysis_mesh_flg, PETSC_NULL);
1421
1422 CHKERR PetscOptionsScalar(
1423 "-equilibrium_stop_rate",
1424 "is used to stop calculations when equilibium state is achieved", "",
1425 rate, &rate, &equilibrium_flg);
1426 PetscOptionsEnd();
1427
1428 if (!iNit) {
1429
1430 PetscOptionsBegin(PETSC_COMM_WORLD, "",
1431 "Bone remodeling post-process", "none");
1432
1433 pRT = 1;
1434 CHKERR PetscOptionsInt("-my_output_prt",
1435 "frequncy how often results are dumped on hard disk",
1436 "", 1, &pRT, NULL);
1437 PetscOptionsEnd();
1438
1440 CHKERR postProc.addFieldValuesPostProc("DISPLACEMENTS");
1441 CHKERR postProc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
1442 if (!commonData.less_post_proc) {
1444 }
1445 // CHKERR postProc.addFieldValuesPostProc("RHO");
1446 // CHKERR postProc.addFieldValuesGradientPostProc("RHO");
1447
1448 {
1449 auto &list_of_operators = postProc.getOpPtrVector();
1450
1451 list_of_operators.push_back(
1452 new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
1453 list_of_operators.push_back(new OpCalculateScalarFieldGradient<3>(
1454 "RHO", commonData.data.gradRhoPtr));
1455 // Get displacement gradient at integration points
1456 list_of_operators.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1457 "DISPLACEMENTS", commonData.data.gradDispPtr));
1458 list_of_operators.push_back(new OpCalculateStress(commonData));
1459 list_of_operators.push_back(new OpPostProcStress(
1461 }
1462
1463 CHKERR postProcElastic.generateReferenceElementMesh();
1464 CHKERR postProcElastic.addFieldValuesPostProc("DISPLACEMENTS");
1465 CHKERR postProcElastic.addFieldValuesPostProc("MESH_NODE_POSITIONS");
1466 CHKERR postProcElastic.addFieldValuesGradientPostProc("DISPLACEMENTS");
1467 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
1468 commonData.elasticPtr->setOfBlocks.begin();
1469 for (; sit != commonData.elasticPtr->setOfBlocks.end(); sit++) {
1470 postProcElastic.getOpPtrVector().push_back(new PostProcStress(
1471 postProcElastic.postProcMesh, postProcElastic.mapGaussPts,
1472 "DISPLACEMENTS", sit->second, postProcElastic.commonData));
1473 }
1474
1475 iNit = true;
1476 }
1477
1478 if (mass_postproc) {
1479 Vec mass_vec;
1480 Vec energ_vec;
1481 int mass_vec_ghost[] = {0};
1482 CHKERR VecCreateGhost(PETSC_COMM_WORLD, (!mField.get_comm_rank()) ? 1 : 0,
1483 1, 1, mass_vec_ghost, &mass_vec);
1484
1485 CHKERR VecZeroEntries(mass_vec);
1486 CHKERR VecDuplicate(mass_vec, &energ_vec);
1487
1488 Remodeling::Fe energyCalc(mField);
1489 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", energyCalc, true, false, false,
1490 true);
1491 energyCalc.getOpPtrVector().push_back(
1492 new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
1493 energyCalc.getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<3>(
1494 "RHO", commonData.data.gradRhoPtr));
1495 energyCalc.getOpPtrVector().push_back(
1496 new OpCalculateVectorFieldGradient<3, 3>("DISPLACEMENTS",
1497 commonData.data.gradDispPtr));
1498 energyCalc.getOpPtrVector().push_back(new OpCalculateStress(commonData));
1499 energyCalc.getOpPtrVector().push_back(
1500 new OpMassAndEnergyCalculation("RHO", commonData, energ_vec, mass_vec));
1501 CHKERR DMoFEMLoopFiniteElements(commonData.dm, "FE_REMODELLING",
1502 &energyCalc);
1503
1504 CHKERR VecAssemblyBegin(mass_vec);
1505 CHKERR VecAssemblyEnd(mass_vec);
1506 double mass_sum;
1507 CHKERR VecSum(mass_vec, &mass_sum);
1508
1509 CHKERR VecAssemblyBegin(energ_vec);
1510 CHKERR VecAssemblyEnd(energ_vec);
1511 double energ_sum;
1512 CHKERR VecSum(energ_vec, &energ_sum);
1513
1514 CHKERR PetscPrintf(PETSC_COMM_WORLD,
1515 "Time: %g Mass: %6.5g Elastic energy: %6.5g \n", ts_t,
1516 mass_sum, energ_sum);
1517 CHKERR VecDestroy(&mass_vec);
1518 CHKERR VecDestroy(&energ_vec);
1519
1520 if (equilibrium_flg) {
1521 double equilibrium_rate =
1522 fabs(energ_sum - commonData.cUrrent_psi) / energ_sum;
1523 double equilibrium_mass_rate =
1524 fabs(mass_sum - commonData.cUrrent_mass) / mass_sum;
1525 if (equilibrium_rate < rate) {
1526 CHKERR PetscPrintf(
1527 PETSC_COMM_WORLD,
1528 "Energy equilibrium state is achieved! Difference = %0.6g %%. \n",
1529 equilibrium_rate * 100);
1530 }
1531 if (equilibrium_mass_rate < rate) {
1532 CHKERR PetscPrintf(
1533 PETSC_COMM_WORLD,
1534 "Mass equilibrium state is achieved! Difference = %0.6g %%. \n",
1535 equilibrium_mass_rate * 100);
1536 }
1537 commonData.cUrrent_psi = energ_sum;
1538 commonData.cUrrent_mass = mass_sum;
1539 }
1540 }
1541
1542 int step;
1543#if PETSC_VERSION_LE(3, 8, 0)
1544 CHKERR TSGetTimeStepNumber(ts, &step);
1545#else
1546 CHKERR TSGetStepNumber(ts, &step);
1547#endif
1548
1549 if ((step) % pRT == 0) {
1550
1551 ostringstream sss;
1552 sss << "out_" << step << ".h5m";
1553 CHKERR DMoFEMLoopFiniteElements(commonData.dm, "FE_REMODELLING", &postProc);
1554 CHKERR postProc.postProcMesh.write_file(sss.str().c_str(), "MOAB",
1555 "PARALLEL=WRITE_PART");
1556 if (analysis_mesh_flg) {
1557 ostringstream ttt;
1558 ttt << "analysis_mesh_" << step << ".h5m";
1559 CHKERR mField.get_moab().write_file(ttt.str().c_str(), "MOAB",
1560 "PARALLEL=WRITE_PART");
1561 }
1562
1563 if (commonData.nOremodellingBlock) {
1564 CHKERR DMoFEMLoopFiniteElements(commonData.dm, "ELASTIC",
1565 &postProcElastic);
1566 ostringstream ss;
1567 ss << "out_elastic_" << step << ".h5m";
1568 CHKERR postProcElastic.postProcMesh.write_file(ss.str().c_str(), "MOAB",
1569 "PARALLEL=WRITE_PART");
1570 }
1571 }
1573}
1574
1575MoFEMErrorCode Remodeling::getParameters() {
1576
1578 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling", "none");
1579
1580 commonData.oRder = 2;
1581 CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
1582 &commonData.oRder, PETSC_NULL);
1583
1584 PetscOptionsEnd();
1585
1586 CHKERR mField.get_moab().get_entities_by_type(0, MBTET, commonData.tEts_all);
1588 string name = it->getName();
1589 if (name.compare(0, 14, "NO_REMODELLING") == 0) {
1590 commonData.nOremodellingBlock = true;
1591 EntityHandle meshset = it->getMeshset();
1592 CHKERR this->mField.get_moab().get_entities_by_type(
1593 meshset, MBTET, commonData.tEts_block, true);
1594 commonData.tEts_all =
1595 subtract(commonData.tEts_all, commonData.tEts_block);
1596 }
1597 }
1598
1600}
1601
1602MoFEMErrorCode Remodeling::addFields() {
1603
1605
1606 // Seed all mesh entities to MoFEM database, those entities can be potentially
1607 // used as finite elements or as entities which carry some approximation
1608 // field.
1609 commonData.bitLevel.set(0);
1610 CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
1611 0, 3, commonData.bitLevel);
1612 int order = commonData.oRder;
1613
1614 // Add displacement field
1615 CHKERR mField.add_field("DISPLACEMENTS", H1,
1616 /*AINSWORTH_LOBATTO_BASE*/ AINSWORTH_LEGENDRE_BASE, 3,
1617 MB_TAG_SPARSE, MF_ZERO);
1618 // Add field representing ho-geometry
1619 CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3,
1620 MB_TAG_SPARSE, MF_ZERO);
1621
1622 // Check if density is available, if not add density field.
1623 bool add_rho_field = false;
1624 if (!mField.check_field("RHO")) {
1625 CHKERR mField.add_field("RHO", H1, AINSWORTH_LEGENDRE_BASE, 1,
1626 MB_TAG_SPARSE, MF_ZERO);
1627 // FIXME
1628 CHKERR mField.add_ents_to_field_by_type(commonData.tEts_all, MBTET, "RHO");
1629 // CHKERR mField.add_ents_to_field_by_type(0,MBTET,"RHO");
1630 CHKERR mField.synchronise_field_entities("RHO");
1631 add_rho_field = true;
1632
1633 CHKERR mField.set_field_order(0, MBVERTEX, "RHO", 1);
1634 CHKERR mField.set_field_order(0, MBEDGE, "RHO", order - 1);
1635 CHKERR mField.set_field_order(0, MBTRI, "RHO", order - 1);
1636 CHKERR mField.set_field_order(0, MBTET, "RHO", order - 1);
1637 }
1638
1639 // Add entities to field
1640 CHKERR mField.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENTS");
1641 CHKERR mField.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
1642
1643 // Set approximation order to entities
1644 CHKERR mField.set_field_order(0, MBVERTEX, "DISPLACEMENTS", 1);
1645 CHKERR mField.set_field_order(0, MBEDGE, "DISPLACEMENTS", order);
1646 CHKERR mField.set_field_order(0, MBTRI, "DISPLACEMENTS", order);
1647 CHKERR mField.set_field_order(0, MBTET, "DISPLACEMENTS", order);
1648
1649 // Assumes that geometry is approximated using 2nd order polynomials.
1650
1651 // Apply 2nd order only on skin
1652 {
1653 // Skinner skin(&mField.get_moab());
1654 // Range faces,tets;
1655 // CHKERR mField.get_moab().get_entities_by_type(0,MBTET,tets);
1656 // CHKERR skin.find_skin(0,tets,false,faces);
1657 // Range edges;
1658 // CHKERR mField.get_moab().get_adjacencies(
1659 // faces,1,false,edges,moab::Interface::UNION
1660 // );
1661 CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
1662 }
1663
1664 CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
1665
1666 // Build fields
1667 CHKERR mField.build_fields();
1668
1669 // If order was not set from CT scans set homogenous order equal to
1670 // reference bone density
1671 if (add_rho_field) {
1672 CHKERR mField.getInterface<FieldBlas>()->setField(commonData.rHo_ref,
1673 MBVERTEX, "RHO");
1674 // const DofEntity_multiIndex *dofs_ptr;
1675 // CHKERR mField.get_dofs(&dofs_ptr);
1676 // for(_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(mField,"RHO",dit)) {
1677 // cerr << (*dit)->getFieldData() << endl;
1678 // }
1679 }
1680
1681 // Project geometry given on 10-node tets on ho-geometry
1682 Projection10NodeCoordsOnField ent_method_material(mField,
1683 "MESH_NODE_POSITIONS");
1684 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
1685
1687}
1688
1689MoFEMErrorCode Remodeling::addElements() {
1691
1692 CHKERR mField.add_finite_element("FE_REMODELLING", MF_ZERO);
1693 CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING",
1694 "DISPLACEMENTS");
1695 CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING",
1696 "DISPLACEMENTS");
1697 CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING",
1698 "DISPLACEMENTS");
1699 CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING", "RHO");
1700 CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING",
1701 "MESH_NODE_POSITIONS");
1702 CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING", "RHO");
1703 CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING", "RHO");
1704 CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_all, MBTET,
1705 "FE_REMODELLING");
1706
1707 CHKERR mField.add_finite_element("ELASTIC");
1708 CHKERR mField.modify_finite_element_add_field_row("ELASTIC", "DISPLACEMENTS");
1709 CHKERR mField.modify_finite_element_add_field_col("ELASTIC", "DISPLACEMENTS");
1710 CHKERR mField.modify_finite_element_add_field_data("ELASTIC",
1711 "DISPLACEMENTS");
1712 CHKERR mField.modify_finite_element_add_field_data("ELASTIC",
1713 "MESH_NODE_POSITIONS");
1714
1715 if (commonData.nOremodellingBlock) {
1716 CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_block,
1717 MBTET, "ELASTIC");
1718 }
1719
1720 commonData.elasticMaterialsPtr =
1721 boost::shared_ptr<ElasticMaterials>(new ElasticMaterials(mField));
1722 commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
1723 new NonlinearElasticElement(mField, 10));
1724 CHKERR commonData.elasticMaterialsPtr->setBlocks(
1725 commonData.elasticPtr->setOfBlocks);
1726 CHKERR commonData.elasticPtr->addElement("ELASTIC", "DISPLACEMENTS");
1727 CHKERR commonData.elasticPtr->setOperators(
1728 "DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
1729
1730 CHKERR mField.build_finite_elements();
1731 CHKERR mField.build_adjacencies(commonData.bitLevel);
1732
1733 // Allocate memory for density and gradient of displacements at integration
1734 // points
1735 commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
1736 commonData.data.gradRhoPtr =
1737 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1738 commonData.data.gradDispPtr =
1739 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1740
1741 // Add operators Rhs
1742 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feRhs), true, false,
1743 false, false);
1744 auto &list_of_operators_rhs = commonData.feRhs->getOpPtrVector();
1745 // Get density at integration points
1746 list_of_operators_rhs.push_back(
1747 new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
1748 list_of_operators_rhs.push_back(
1749 new OpCalculateScalarFieldGradient<3>("RHO", commonData.data.gradRhoPtr));
1750 list_of_operators_rhs.push_back(new OpGetRhoTimeDirevative(commonData));
1751 // Get displacement gradient at integration points
1752 list_of_operators_rhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1753 "DISPLACEMENTS", commonData.data.gradDispPtr));
1754 list_of_operators_rhs.push_back(new OpCalculateStress(commonData));
1755 list_of_operators_rhs.push_back(new OpAssmbleStressRhs(commonData));
1756 list_of_operators_rhs.push_back(new OpAssmbleRhoRhs(commonData));
1757
1758 // Add operators Lhs
1759 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feLhs), true, false,
1760 false, false);
1761 auto &list_of_operators_lhs = commonData.feLhs->getOpPtrVector();
1762 // Get density at integration points
1763 list_of_operators_lhs.push_back(
1764 new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
1765 list_of_operators_rhs.push_back(
1766 new OpCalculateScalarFieldGradient<3>("RHO", commonData.data.gradRhoPtr));
1767 // Get displacement gradient at integration points
1768 list_of_operators_lhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1769 "DISPLACEMENTS", commonData.data.gradDispPtr));
1770 if (commonData.with_adol_c) {
1771 list_of_operators_lhs.push_back(
1772 new OpCalculateStressTangentWithAdolc(commonData));
1773 list_of_operators_lhs.push_back(
1774 new OpAssmbleStressLhs_dF<true>(commonData));
1775 } else {
1776 list_of_operators_lhs.push_back(new OpCalculateStressTangent(commonData));
1777 list_of_operators_lhs.push_back(
1778 new OpAssmbleStressLhs_dF<false>(commonData));
1779 }
1780 list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dRho(commonData));
1781 list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dF(commonData));
1782 list_of_operators_lhs.push_back(new OpAssmbleStressLhs_dRho(commonData));
1783
1785}
1786
1787MoFEMErrorCode Remodeling::addMomentumFluxes() {
1788
1790 // Add Neumann forces elements
1791 CHKERR MetaNeummanForces::addNeumannBCElements(mField, "DISPLACEMENTS");
1792 CHKERR MetaNodalForces::addElement(mField, "DISPLACEMENTS");
1793 CHKERR MetaEdgeForces::addElement(mField, "DISPLACEMENTS");
1794
1795 // forces and pressures on surface
1796 CHKERR MetaNeummanForces::setMomentumFluxOperators(
1797 mField, commonData.neumannForces, PETSC_NULL, "DISPLACEMENTS");
1798 // noadl forces
1799 CHKERR MetaNodalForces::setOperators(mField, commonData.nodalForces,
1800 PETSC_NULL, "DISPLACEMENTS");
1801 // edge forces
1802 CHKERR MetaEdgeForces::setOperators(mField, commonData.edgeForces, PETSC_NULL,
1803 "DISPLACEMENTS");
1804
1805 for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
1806 commonData.neumannForces.begin();
1807 mit != commonData.neumannForces.end(); mit++) {
1808 mit->second->methodsOp.push_back(
1809 new TimeForceScale("-my_load_history", false));
1810 }
1811 for (boost::ptr_map<string, NodalForce>::iterator mit =
1812 commonData.nodalForces.begin();
1813 mit != commonData.nodalForces.end(); mit++) {
1814 mit->second->methodsOp.push_back(
1815 new TimeForceScale("-my_load_history", false));
1816 }
1817 for (boost::ptr_map<string, EdgeForce>::iterator mit =
1818 commonData.edgeForces.begin();
1819 mit != commonData.edgeForces.end(); mit++) {
1820 mit->second->methodsOp.push_back(
1821 new TimeForceScale("-my_load_history", false));
1822 }
1823
1825}
1826
1827MoFEMErrorCode Remodeling::buildDM() {
1828
1830 commonData.dm_name = "DM_REMODELING";
1831 // Register DM problem
1832 CHKERR DMRegister_MoFEM(commonData.dm_name);
1833 CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
1834 CHKERR DMSetType(commonData.dm, commonData.dm_name);
1835 CHKERR DMMoFEMSetIsPartitioned(commonData.dm, PETSC_TRUE);
1836 // Create DM instance
1837 CHKERR DMMoFEMCreateMoFEM(commonData.dm, &mField, commonData.dm_name,
1838 commonData.bitLevel);
1839 // Configure DM form line command options (DM itself, solvers,
1840 // pre-conditioners, ... )
1841 CHKERR DMSetFromOptions(commonData.dm);
1842 // Add elements to dm (only one here)
1843 CHKERR DMMoFEMAddElement(commonData.dm, "FE_REMODELLING");
1844 CHKERR DMMoFEMAddElement(commonData.dm, "ELASTIC");
1845
1846 if (mField.check_finite_element("FORCE_FE")) {
1847 CHKERR DMMoFEMAddElement(commonData.dm, "FORCE_FE");
1848 }
1849 if (mField.check_finite_element("PRESSURE_FE")) {
1850 CHKERR DMMoFEMAddElement(commonData.dm, "PRESSURE_FE");
1851 }
1852 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
1853 CHKERR DMSetUp(commonData.dm);
1854
1856}
1857
1858MoFEMErrorCode Remodeling::solveDM() {
1859
1861
1862 Vec F = commonData.F;
1863 Vec D = commonData.D;
1864 Mat A = commonData.A;
1865 DM dm = commonData.dm;
1866 TS ts = commonData.ts;
1867
1868 CHKERR TSSetIFunction(ts, F, PETSC_NULL, PETSC_NULL);
1869 CHKERR TSSetIJacobian(ts, A, A, PETSC_NULL, PETSC_NULL);
1870 double ftime = 1;
1871#if PETSC_VERSION_GE(3, 8, 0)
1872 CHKERR TSSetMaxTime(ts, ftime);
1873#endif
1874 CHKERR TSSetFromOptions(ts);
1875 CHKERR TSSetDM(ts, dm);
1876
1877 SNES snes;
1878 CHKERR TSGetSNES(ts, &snes);
1879 KSP ksp;
1880 CHKERR SNESGetKSP(snes, &ksp);
1881 PC pc;
1882 CHKERR KSPGetPC(ksp, &pc);
1883 PetscBool is_pcfs = PETSC_FALSE;
1884 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1885
1886 // Set up FIELDSPLIT
1887 // Only is user set -pc_type fieldsplit
1888 if (is_pcfs == PETSC_TRUE) {
1889 IS is_disp, is_rho;
1890 const MoFEM::Problem *problem_ptr;
1891 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
1892 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1893 problem_ptr->getName(), ROW, "DISPLACEMENTS", 0, 3, &is_disp);
1894 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1895 problem_ptr->getName(), ROW, "RHO", 0, 1, &is_rho);
1896 CHKERR ISSort(is_disp);
1897 CHKERR ISSort(is_rho);
1898 CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
1899 CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
1900 CHKERR ISDestroy(&is_disp);
1901 CHKERR ISDestroy(&is_rho);
1902 }
1903
1904 // Monitor
1905 MonitorPostProc post_proc(mField, commonData);
1906 TsCtx *ts_ctx;
1907 DMMoFEMGetTsCtx(dm, &ts_ctx);
1908 {
1909 ts_ctx->getPostProcessMonitor().push_back(&post_proc);
1910 CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
1911 }
1912
1913#if PETSC_VERSION_GE(3, 7, 0)
1914 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1915#endif
1916 CHKERR TSSolve(ts, D);
1917 CHKERR TSGetTime(ts, &ftime);
1918 PetscInt steps, snesfails, rejects, nonlinits, linits;
1919#if PETSC_VERSION_LE(3, 8, 0)
1920 CHKERR TSGetTimeStepNumber(ts, &steps);
1921#else
1922 CHKERR TSGetStepNumber(ts, &steps);
1923#endif
1924
1925 CHKERR TSGetSNESFailures(ts, &snesfails);
1926 CHKERR TSGetStepRejections(ts, &rejects);
1927 CHKERR TSGetSNESIterations(ts, &nonlinits);
1928 CHKERR TSGetKSPIterations(ts, &linits);
1929 PetscPrintf(PETSC_COMM_WORLD,
1930 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
1931 "linits %D\n",
1932 steps, rejects, snesfails, ftime, nonlinits, linits);
1933
1934 if (commonData.is_atom_testing) {
1935 if (commonData.cUrrent_psi < 1.67 || commonData.cUrrent_psi > 1.68)
1936 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "atom test diverged");
1937 }
1938
1940}
1941
1942MoFEMErrorCode OpMassAndEnergyCalculation::doWork(
1943 int row_side, EntityType row_type,
1944 DataForcesAndSourcesCore::EntData &row_data) {
1946
1947 if (row_type != MBVERTEX)
1949
1950 const int nb_gauss_pts = row_data.getN().size1();
1951 // commonData.data.vecPsi.resize(nb_gauss_pts,false);
1952 auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1953 auto psi = getFTensor0FromVec(commonData.data.vecPsi);
1954
1955 double energy = 0;
1956 double mass = 0;
1957 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1958
1959 double vol = getVolume() * getGaussPts()(3, gg);
1960 energy += vol * psi;
1961 mass += vol * rho;
1962
1963 ++psi;
1964 ++rho;
1965 }
1966
1967 CHKERR VecSetValue(energVec, 0, energy, ADD_VALUES);
1968 CHKERR VecSetValue(massVec, 0, mass, ADD_VALUES);
1970}
1971
1972} // namespace BoneRemodeling
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT)
#define TENSOR2_MAT_PTR(MAT)
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT)
#define TENSOR4_MAT_PTR(MAT)
constexpr double a
@ ROW
@ MF_ZERO
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT)
#define TENSOR2_MAT_PTR(MAT)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ BLOCKSET
@ MOFEM_STD_EXCEPTION_THROW
Definition definitions.h:39
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
#define MU(E, NU)
Definition fem_tools.h:23
#define LAMBDA(E, NU)
Definition fem_tools.h:22
@ R
@ F
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1113
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
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
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition DMMoFEM.cpp:422
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition DMMoFEM.cpp:1132
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
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:576
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2 or H1 field gradient.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
const double c
speed of light (cm/ns)
double D
const double n
refractive index of diffusive medium
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
MoFEM::TsCtx * ts_ctx
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
Definition Types.hpp:92
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
UBlasVector< double > VectorDouble
Definition Types.hpp:68
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition TsCtx.cpp:263
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
constexpr IntegrationType I
constexpr AssemblyType A
double rho
Definition plastic.cpp:144
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:126
double young_modulus
Young modulus.
Definition plastic.cpp:125
FTensor::Index< 'm', 3 > m
boost::shared_ptr< MatrixDouble > gradDispPtr
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
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
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
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.
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual int get_comm_rank() const =0
Deprecated interface functions.
Basic algebra on fields.
Definition FieldBlas.hpp:21
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for double precision scalar field values calculation.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
keeps basic data about problem
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
TS ts
PETSc time stepping solver object.
PetscReal ts_t
Current time value.
Interface for Time Stepping (TS) solver.
Definition TsCtx.hpp:17
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
Definition TsCtx.hpp:148
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
MoFEMErrorCode operator()()
Main operator function executed for each loop iteration.
MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
MoFEM::Interface & mField
MonitorPostProc(MoFEM::Interface &m_field, std::map< int, NonlinearElasticElement::BlockData > &set_of_blocks, NonlinearElasticElement::MyVolumeFE &fe_elastic_energy, ConvectiveMassElement::MyVolumeFE &fe_kinetic_energy)
PostProcVolumeOnRefinedMesh postProc
structure grouping operators and data used for calculation of nonlinear elastic element
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateStress(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts)
std::vector< EntityHandle > & mapGaussPts
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, DataAtIntegrationPts &common_data, BlockData &data)
DataAtIntegrationPts & commonData
PetscErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::vector< EntityHandle > mapGaussPts
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Force scale operator for reading two columns.

◆ SYMMETRIC_TENSOR2_VEC_PTR

#define SYMMETRIC_TENSOR2_VEC_PTR (   VEC)     &VEC[0], &VEC[1], &VEC[2], &VEC[3], &VEC[4], &VEC[5]

Definition at line 57 of file Remodeling.cpp.

◆ SYMMETRIC_TENSOR4_MAT_PTR

#define SYMMETRIC_TENSOR4_MAT_PTR (   MAT)
Value:
&MAT(0, 0), &MAT(0, 1), &MAT(0, 2), &MAT(0, 3), &MAT(0, 4), &MAT(0, 5), \
&MAT(1, 0), &MAT(1, 1), &MAT(1, 2), &MAT(1, 3), &MAT(1, 4), &MAT(1, 5), \
&MAT(2, 0), &MAT(2, 1), &MAT(2, 2), &MAT(2, 3), &MAT(2, 4), &MAT(2, 5), \
&MAT(3, 0), &MAT(3, 1), &MAT(3, 2), &MAT(3, 3), &MAT(3, 4), &MAT(3, 5), \
&MAT(4, 0), &MAT(4, 1), &MAT(4, 2), &MAT(4, 3), &MAT(4, 4), &MAT(4, 5), \
&MAT(5, 0), &MAT(5, 1), &MAT(5, 2), &MAT(5, 3), &MAT(5, 4), &MAT(5, 5)

Definition at line 40 of file Remodeling.cpp.

◆ TENSOR1_VEC_PTR

#define TENSOR1_VEC_PTR (   VEC)    &VEC[0], &VEC[1], &VEC[2]

Definition at line 38 of file Remodeling.cpp.

◆ TENSOR2_MAT_PTR

#define TENSOR2_MAT_PTR (   MAT)
Value:
&MAT(0, 0), &MAT(1, 0), &MAT(2, 0), &MAT(3, 0), &MAT(4, 0), &MAT(5, 0), \
&MAT(6, 0), &MAT(7, 0), &MAT(8, 0)

Definition at line 50 of file Remodeling.cpp.

◆ TENSOR4_MAT_PTR

#define TENSOR4_MAT_PTR (   MAT)    &MAT(0, 0), MAT.size2()

Definition at line 48 of file Remodeling.cpp.