v0.14.0
Loading...
Searching...
No Matches
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>
26using namespace std;
27namespace po = boost::program_options;
28
30using 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
60namespace BoneRemodeling {
61
63
64 PetscBool flg_config = PETSC_FALSE;
65 is_atom_testing = PETSC_FALSE;
66 with_adol_c = PETSC_FALSE;
67 char my_config_file_name[255];
68 double young_modulus = 5.0e+9; // Set here some typical value for bone 5 GPa
69 double poisson_ratio = 0.3; // Set here some typical value for bone
70 c = 4.0e-5; // Set here some typical value for bone. days/ m2 // this
71 // parameters governs how fast we can achieve equilibrium
72 m = 3.25; // Set some parameter typical for bone
73 n = 2.25; // Set Some parameter typical for bone
74 rHo_ref = 1000; // Set Some parameter typical for bone. kg/m3
75 pSi_ref = 2.75e4; // Set Some parameter typical for bone. energy / volume unit
76 // J/m3 = Pa
77 rHo_max = rHo_ref + 0.5 * rHo_ref;
78 rHo_min = rHo_ref - 0.5 * rHo_ref;
79 b = 0;
80 R0 = 0.0; // Set Some parameter typical for bone.
81 cUrrent_psi = 0.0;
82 cUrrent_mass = 0.0;
83 nOremodellingBlock = false;
84 less_post_proc = PETSC_FALSE;
85
87 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling parameters",
88 "none");
89
90 // config file
91 CHKERR PetscOptionsString("-my_config", "configuration file name", "",
92 "my_config.in", my_config_file_name, 255,
93 &flg_config);
94
95 if (flg_config) {
96 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Config file: %s loaded.\n",
97 my_config_file_name);
98 try {
99 ifstream ini_file(my_config_file_name);
100 if (!ini_file.is_open()) {
101 SETERRQ(PETSC_COMM_SELF, 1, "*** -my_config does not exist ***");
102 }
103 po::variables_map vm;
104 po::options_description config_file_options;
105 // std::cout << my_config_file_name << std::endl;
106
107 config_file_options.add_options()(
108 "young_modulus", po::value<double>(&young_modulus)->default_value(1));
109 config_file_options.add_options()(
110 "poisson_ratio", po::value<double>(&poisson_ratio)->default_value(1));
111 config_file_options.add_options()(
112 "c", po::value<double>(&c)->default_value(1));
113 config_file_options.add_options()(
114 "m", po::value<double>(&m)->default_value(1));
115 config_file_options.add_options()(
116 "n", po::value<double>(&n)->default_value(1));
117 config_file_options.add_options()(
118 "rHo_ref", po::value<double>(&rHo_ref)->default_value(1));
119 config_file_options.add_options()(
120 "pSi_ref", po::value<double>(&pSi_ref)->default_value(1));
121 config_file_options.add_options()(
122 "R0", po::value<double>(&R0)->default_value(1));
123
124 po::parsed_options parsed =
125 parse_config_file(ini_file, config_file_options, true);
126 store(parsed, vm);
127 po::notify(vm);
128
129 } catch (const std::exception &ex) {
130 std::ostringstream ss;
131 ss << ex.what() << std::endl;
132 SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
133 }
134 }
135
136 CHKERR PetscOptionsScalar("-young_modulus", "get young modulus", "",
137 young_modulus, &young_modulus, PETSC_NULL);
138
139 CHKERR PetscOptionsScalar("-poisson_ratio", "get poisson_ratio", "",
140 poisson_ratio, &poisson_ratio, PETSC_NULL);
141
142 CHKERR PetscOptionsScalar("-c", "density evolution (growth) velocity [d/m^2]",
143 "", c, &c, PETSC_NULL);
144
145 CHKERR PetscOptionsScalar("-m", "algorithmic exponent", "", m, &m,
146 PETSC_NULL);
147
148 CHKERR PetscOptionsScalar("-n", "porosity exponent", "", n, &n, PETSC_NULL);
149
150 CHKERR PetscOptionsInt("-b", "bell function exponent", "", b, &b, PETSC_NULL);
151
152 CHKERR PetscOptionsScalar("-rho_ref", "reference bone density", "", rHo_ref,
153 &rHo_ref, PETSC_NULL);
154
155 CHKERR PetscOptionsScalar("-rho_max", "reference bone density", "", rHo_max,
156 &rHo_max, PETSC_NULL);
157 CHKERR PetscOptionsScalar("-rho_min", "reference bone density", "", rHo_min,
158 &rHo_min, PETSC_NULL);
159
160 CHKERR PetscOptionsScalar("-psi_ref", "reference energy density", "", pSi_ref,
161 &pSi_ref, PETSC_NULL);
162
163 CHKERR PetscOptionsScalar("-r0", "mass source", "", R0, &R0, PETSC_NULL);
164
165 CHKERR PetscOptionsBool("-my_is_atom_test",
166 "is used with testing, exit with error when diverged",
167 "", is_atom_testing, &is_atom_testing, PETSC_NULL);
168
169 CHKERR PetscOptionsBool("-less_post_proc",
170 "is used to reduce output file size", "",
171 less_post_proc, &less_post_proc, PETSC_NULL);
172
173 CHKERR PetscOptionsBool(
174 "-with_adolc",
175 "calculate the material tangent with automatic differentiation", "",
176 with_adol_c, &with_adol_c, PETSC_NULL);
177
180
181 CHKERR PetscPrintf(PETSC_COMM_WORLD,
182 "Young's modulus E[Pa]: %4.3g\n",
184 CHKERR PetscPrintf(PETSC_COMM_WORLD,
185 "Poisson ratio nu[-] %4.3g\n",
187 CHKERR PetscPrintf(PETSC_COMM_WORLD,
188 "Lame coefficient lambda[Pa]: %4.3g\n", lambda);
189 CHKERR PetscPrintf(PETSC_COMM_WORLD,
190 "Lame coefficient mu[Pa]: %4.3g\n", mu);
191 CHKERR PetscPrintf(PETSC_COMM_WORLD,
192 "Density growth velocity [d/m2] %4.3g\n", c);
193 CHKERR PetscPrintf(PETSC_COMM_WORLD,
194 "Algorithmic exponent m[-]: %4.3g\n", m);
195 CHKERR PetscPrintf(PETSC_COMM_WORLD,
196 "Porosity exponent n[-]: %4.3g\n", n);
197 CHKERR PetscPrintf(PETSC_COMM_WORLD,
198 "Reference density rHo_ref[kg/m3]: %4.3g\n", rHo_ref);
199 CHKERR PetscPrintf(PETSC_COMM_WORLD,
200 "Reference energy pSi_ref[Pa]: %4.3g\n", pSi_ref);
201 CHKERR PetscPrintf(PETSC_COMM_WORLD,
202 "Mass conduction R0[kg/m3s]: %4.3g\n", R0);
203 CHKERR PetscPrintf(PETSC_COMM_WORLD,
204 "Bell function coefficent b[-]: %4.3g\n", b);
205 if (b != 0) {
206 CHKERR PetscPrintf(PETSC_COMM_WORLD,
207 "Max density rHo_max[kg/m3]: %4.3g\n", rHo_max);
208 CHKERR PetscPrintf(PETSC_COMM_WORLD,
209 "Min density rHo_min[kg/m3]: %4.3g\n", rHo_min);
210 }
211
212 ierr = PetscOptionsEnd();
213 CHKERRQ(ierr);
214
216}
217
218inline 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
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
233 Remodeling::CommonData &common_data)
235 "RHO", UserDataOperator::OPROW),
236 commonData(common_data) {}
237
240 DataForcesAndSourcesCore::EntData &data) {
241
243
244 const int nb_gauss_pts = data.getN().size1();
245 if (!nb_gauss_pts)
247 const int nb_base_functions = data.getN().size2();
248
249 auto base_function = data.getFTensor0N();
250 commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
251 if (type == MBVERTEX) {
252 commonData.data.vecRhoDt.clear();
253 }
254
255 int nb_dofs = data.getIndices().size();
256 if (nb_dofs == 0)
258 double *array;
259 CHKERR VecGetArray(getFEMethod()->ts_u_t, &array);
261
262 for (int gg = 0; gg < nb_gauss_pts; gg++) {
263 int bb = 0;
264 for (; bb < nb_dofs; bb++) {
265 const double field_data = array[data.getLocalIndices()[bb]];
266 drho_dt += field_data * base_function;
267 ++base_function;
268 }
269 // It is possible to have more base functions than dofs
270 for (; bb != nb_base_functions; bb++)
271 ++base_function;
272 ++drho_dt;
273 }
274
275 CHKERR VecRestoreArray(getFEMethod()->ts_u_t, &array);
276
278}
279
282 "DISPLACEMENTS", UserDataOperator::OPROW),
283 commonData(common_data) {
284 // Set that this operator is executed only for vertices on element
285 doVertices = true;
286 doEdges = false;
287 doQuads = false;
288 doTris = false;
289 doTets = false;
290 doPrisms = false;
291}
292
295 DataForcesAndSourcesCore::EntData &data) {
296
298
299 if (type != MBVERTEX)
301
303 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
304 "Gradient at integration points of displacement is not calculated");
305 }
306 auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
307
308 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
309 if (!commonData.data.rhoPtr) {
310 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
311 "Density at integration points is not calculated");
312 }
313
315
316 commonData.data.vecPsi.resize(nb_gauss_pts, false);
318 commonData.data.vecR.resize(nb_gauss_pts, false);
320
321 commonData.data.matInvF.resize(9, nb_gauss_pts, false);
322 auto invF = getFTensor2FromMat<3, 3>(commonData.data.matInvF);
323 commonData.data.vecDetF.resize(nb_gauss_pts, false);
326 matS.resize(nb_gauss_pts, 6, false);
328 commonData.data.matP.resize(9, nb_gauss_pts, false);
329 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
331 matF.resize(9, nb_gauss_pts, false);
333
334 const double c = commonData.c;
335 const double n = commonData.n;
336 const double m = commonData.m;
337 const double rho_ref = commonData.rHo_ref;
338 const double psi_ref = commonData.pSi_ref;
339 const double lambda = commonData.lambda;
340 const double mu = commonData.mu;
344
345 for (int gg = 0; gg != nb_gauss_pts; gg++) {
346
347 // Get deformation gradient
348 F(i, I) = grad(i, I);
349 F(0, 0) += 1;
350 F(1, 1) += 1;
351 F(2, 2) += 1;
352
354 CHKERR invertTensor3by3(F, det_f, invF);
355
356 const double log_det = log(det_f);
357 psi = F(i, J) * F(i, J) - 3 - 2 * log_det;
358 psi *= 0.5 * mu;
359 psi += 0.5 * lambda * log_det * log_det;
360
361 const double coef = lambda * log_det - mu;
362 P(i, J) = mu * F(i, J) + coef * invF(J, i);
363 S(i, I) = invF(i, J) ^ P(J, I);
364
365 const double kp = commonData.getCFromDensity(rho);
366 R = c * kp * (pow(rho / rho_ref, n - m) * psi - psi_ref);
367
368 ++det_f;
369 ++invF;
370 ++grad;
371 ++rho;
372 ++psi;
373 ++S;
374 ++P;
375 ++R;
376 ++F;
377 }
378
380}
381
383 Remodeling::CommonData &common_data)
385 "DISPLACEMENTS", UserDataOperator::OPROW),
386 commonData(common_data) {
387 // Set that this operator is executed only for vertices on element
388 doVertices = true;
389 doEdges = false;
390 doQuads = false;
391 doTris = false;
392 doTets = false;
393 doPrisms = false;
394}
395
397 int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
398
400
401 if (type != MBVERTEX)
403
404 auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
405 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
406
407 vecC.resize(6, false);
410 dS_dC.resize(6, 6 * nb_gauss_pts, false);
411 dS_dC.clear();
413
415 matS.resize(nb_gauss_pts, 6, false);
418 dP_dF.resize(81, nb_gauss_pts, false);
420 matF.resize(9, nb_gauss_pts, false);
423 commonData.data.matP.resize(9, nb_gauss_pts, false);
424 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
425 commonData.data.vecPsi.resize(nb_gauss_pts, false);
427
434
435 for (int gg = 0; gg != nb_gauss_pts; gg++) {
436
437 // Get deformation gradient
438 F(i, I) = grad(i, I);
439 F(0, 0) += 1;
440 F(1, 1) += 1;
441 F(2, 2) += 1;
442
443 // Calculate Green defamation Tensor0
444 // ^ that means multiplication of Tensor 2 which result in symmetric Tensor
445 // rank 2.
446 C(I, J) = F(i, I) ^ F(i, J);
447
448 // Calculate free energy
449 CHKERR recordFreeEnergy_dC<FTensor::Tensor0<FTensor::PackPtr<double *, 1>>,
451 commonData, C, psi);
452
453 int r_g = ::gradient(commonData.tAg, 6, &vecC[0], &matS(gg, 0));
454 if (r_g < 0) {
455 // That means that energy function is not smooth and derivative
456 // can not be calculated,
457 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
458 "ADOL-C function evaluation with error");
459 }
460
461 // Scale energy density and 2nd Piola Stress
462 S(0, 0) *= 2;
463 S(1, 1) *= 2;
464 S(2, 2) *= 2;
465
466 // Calculate 1st Piola-Stress tensor
467 P(i, J) = F(i, I) * S(I, J);
468
469 const int is = 6 * gg;
470 double *hessian_ptr[6] = {&dS_dC(0, is), &dS_dC(1, is), &dS_dC(2, is),
471 &dS_dC(3, is), &dS_dC(4, is), &dS_dC(5, is)};
472 int r_h = ::hessian(commonData.tAg, 6, &vecC[0], hessian_ptr);
473 if (r_h < 0) {
474 // That means that energy function is not smooth and derivative
475 // can not be calculated,
476 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
477 "ADOL-C function evaluation with error");
478 }
479
480 // FIXME: TODO: Here tensor 4 with two minor symmetries is used, since
481 // \psi_{,ijkl} = \psi_{,klij} = \psi_{,lkji} = ... there is additional
482 // major symmetry. That is 21 independent components not 36 how is now
483 // implemented using Tenso4_ddg.
484 //
485 // That way we have to fill off-diagonal part for dS_dC. In future that
486 // need to be fixed by implementation of Tenso4 with additional major
487 // symmetries. Pleas consider doing such implementation by yourself. You
488 // have support of MoFEM developing team to guide you through this process.
489 //
490 // This deteriorate efficiency.
491 //
492 // TODO: Implement matrix with more symmetries.
493 //
494 for (int ii = 0; ii < 6; ii++) {
495 for (int jj = 0; jj < ii; jj++) {
496 dS_dC(jj, is + ii) = dS_dC(ii, is + jj);
497 }
498 }
499
500 // As result that symmetry of deformation gradient C is used, terms which
501 // have mixed indices are two time to big, and terms which has two times
502 // midxed index are four times to big. Those terms are not multiplied.
503 D1(0, 0, 0, 0) *= 4;
504 D1(1, 1, 1, 1) *= 4;
505 D1(2, 2, 2, 2) *= 4;
506
507 D1(0, 0, 1, 1) *= 4;
508 D1(1, 1, 0, 0) *= 4;
509 D1(1, 1, 2, 2) *= 4;
510 D1(2, 2, 1, 1) *= 4;
511 D1(0, 0, 2, 2) *= 4;
512 D1(2, 2, 0, 0) *= 4;
513
514 D1(0, 0, 0, 1) *= 2;
515 D1(0, 0, 0, 2) *= 2;
516 D1(0, 0, 1, 2) *= 2;
517 D1(0, 1, 0, 0) *= 2;
518 D1(0, 2, 0, 0) *= 2;
519 D1(1, 2, 0, 0) *= 2;
520
521 D1(1, 1, 0, 1) *= 2;
522 D1(1, 1, 0, 2) *= 2;
523 D1(1, 1, 1, 2) *= 2;
524 D1(0, 1, 1, 1) *= 2;
525 D1(0, 2, 1, 1) *= 2;
526 D1(1, 2, 1, 1) *= 2;
527
528 D1(2, 2, 0, 1) *= 2;
529 D1(2, 2, 0, 2) *= 2;
530 D1(2, 2, 1, 2) *= 2;
531 D1(0, 1, 2, 2) *= 2;
532 D1(0, 2, 2, 2) *= 2;
533 D1(1, 2, 2, 2) *= 2;
534
535 D2(i, J, k, L) = (F(i, I) * (D1(I, J, K, L) * F(k, K)));
536 // D2(i, J, k, L) += I(i, k) * S(J, L);
537
538 ++grad;
539 ++S;
540 ++F;
541 ++D1;
542 ++D2;
543 ++P;
544 ++psi;
545 }
546
548}
549
550OpPostProcStress::OpPostProcStress(moab::Interface &post_proc_mesh,
551 std::vector<EntityHandle> &map_gauss_pts,
552 Remodeling::CommonData &common_data)
554 "DISPLACEMENTS", UserDataOperator::OPROW),
555 commonData(common_data), postProcMesh(post_proc_mesh),
556 mapGaussPts(map_gauss_pts) {
557 // Set that this operator is executed only for vertices on element
558 doVertices = true;
559 doEdges = false;
560 doQuads = false;
561 doTris = false;
562 doTets = false;
563 doPrisms = false;
564}
565
568 DataForcesAndSourcesCore::EntData &data) {
569
571
572 if (type != MBVERTEX)
574 double def_VAL[9];
575 bzero(def_VAL, 9 * sizeof(double));
576
577 Tag th_piola2;
578 CHKERR postProcMesh.tag_get_handle("Piola_Stress2", 9, MB_TYPE_DOUBLE,
579 th_piola2, MB_TAG_CREAT | MB_TAG_SPARSE,
580 def_VAL);
581
582 Tag th_psi;
583 CHKERR postProcMesh.tag_get_handle("psi", 1, MB_TYPE_DOUBLE, th_psi,
584 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
585
586 Tag th_rho;
587 CHKERR postProcMesh.tag_get_handle("RHO", 1, MB_TYPE_DOUBLE, th_rho,
588 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
589
590 Tag th_piola1;
591 Tag principal;
592 Tag th_prin_stress_vect1, th_prin_stress_vect2, th_prin_stress_vect3;
593 Tag th_R;
594 Tag th_grad_rho;
596
597 CHKERR postProcMesh.tag_get_handle("Piola_Stress1", 9, MB_TYPE_DOUBLE,
598 th_piola1, MB_TAG_CREAT | MB_TAG_SPARSE,
599 def_VAL);
600
601 CHKERR postProcMesh.tag_get_handle("Principal_stress", 3, MB_TYPE_DOUBLE,
602 principal, MB_TAG_CREAT | MB_TAG_SPARSE,
603 def_VAL);
604
605 CHKERR postProcMesh.tag_get_handle("Principal_stress_S1", 3, MB_TYPE_DOUBLE,
606 th_prin_stress_vect1,
607 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
608 CHKERR postProcMesh.tag_get_handle("Principal_stress_S2", 3, MB_TYPE_DOUBLE,
609 th_prin_stress_vect2,
610 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
611 CHKERR postProcMesh.tag_get_handle("Principal_stress_S3", 3, MB_TYPE_DOUBLE,
612 th_prin_stress_vect3,
613 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
614
615 CHKERR postProcMesh.tag_get_handle("R", 1, MB_TYPE_DOUBLE, th_R,
616 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
617
618 CHKERR postProcMesh.tag_get_handle("grad_rho", 3, MB_TYPE_DOUBLE,
619 th_grad_rho,
620 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
621 }
622
623 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
626 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
629
631 auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
632
633 MatrixDouble3by3 stress(3, 3);
634 VectorDouble3 eigen_values(3);
635 MatrixDouble3by3 eigen_vectors(3, 3);
636
637 for (int gg = 0; gg != nb_gauss_pts; gg++) {
638
639 double val = psi;
640 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &val);
641 val = rho;
642 CHKERR postProcMesh.tag_set_data(th_rho, &mapGaussPts[gg], 1, &val);
643
644 // S is symmetric tensor, need to map it to non-symmetric tensor to save it
645 // on moab tag to be viable in paraview.
646 for (int ii = 0; ii != 3; ii++) {
647 for (int jj = 0; jj != 3; jj++) {
648 stress(ii, jj) = S(ii, jj);
649 }
650 }
651 CHKERR postProcMesh.tag_set_data(th_piola2, &mapGaussPts[gg], 1,
652 &stress(0, 0));
653
655
656 for (int ii = 0; ii != 3; ii++) {
657 for (int jj = 0; jj != 3; jj++) {
658 stress(ii, jj) = P(ii, jj);
659 }
660 }
661 CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
662 &stress(0, 0));
663 val = R;
664 CHKERR postProcMesh.tag_set_data(th_R, &mapGaussPts[gg], 1, &val);
665 const double t2[] = {grad_rho(0), grad_rho(1), grad_rho(2)};
666 CHKERR postProcMesh.tag_set_data(th_grad_rho, &mapGaussPts[gg], 1, t2);
667
668 // Add calculation of eigen values, i.e. prinple stresses.
669 noalias(eigen_vectors) = stress;
670 VectorDouble3 prin_stress_vect1(3);
671 VectorDouble3 prin_stress_vect2(3);
672 VectorDouble3 prin_stress_vect3(3);
673 // LAPACK - eigenvalues and vectors. Applied twice for initial creates
674 // memory space
675 int n = 3, lda = 3, info, lwork = -1;
676 double wkopt;
677 info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
678 &(eigen_values.data()[0]), &wkopt, lwork);
679 if (info != 0)
680 SETERRQ1(PETSC_COMM_SELF, 1,
681 "something is wrong with lapack_dsyev info = %d", info);
682 lwork = (int)wkopt;
683 double work[lwork];
684 info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
685 &(eigen_values.data()[0]), work, lwork);
686 if (info != 0)
687 SETERRQ1(PETSC_COMM_SELF, 1,
688 "something is wrong with lapack_dsyev info = %d", info);
689
690 for (int ii = 0; ii < 3; ii++) {
691 prin_stress_vect1[ii] = eigen_vectors(0, ii);
692 prin_stress_vect2[ii] = eigen_vectors(1, ii);
693 prin_stress_vect3[ii] = eigen_vectors(2, ii);
694 }
695
696 CHKERR postProcMesh.tag_set_data(principal, &mapGaussPts[gg], 1,
697 &eigen_values[0]);
698 CHKERR postProcMesh.tag_set_data(th_prin_stress_vect1, &mapGaussPts[gg],
699 1, &*prin_stress_vect1.begin());
700 CHKERR postProcMesh.tag_set_data(th_prin_stress_vect2, &mapGaussPts[gg],
701 1, &*prin_stress_vect2.begin());
702 CHKERR postProcMesh.tag_set_data(th_prin_stress_vect3, &mapGaussPts[gg],
703 1, &*prin_stress_vect3.begin());
704 }
705
706 ++S;
707 ++P;
708 ++psi;
709 ++R;
710 ++rho;
711 ++grad_rho;
712 }
713
715}
716
718 Remodeling::CommonData &common_data)
720 "DISPLACEMENTS", UserDataOperator::OPROW),
721 commonData(common_data) {
722 // Set that this operator is executed only for vertices on element
723 doVertices = true;
724 doEdges = false;
725 doQuads = false;
726 doTris = false;
727 doTets = false;
728 doPrisms = false;
729}
730
733 DataForcesAndSourcesCore::EntData &data) {
734
736
737 if (type != MBVERTEX)
739
740 auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
741 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
742
743 commonData.data.matInvF.resize(9, nb_gauss_pts, false);
744 auto invF = getFTensor2FromMat<3, 3>(commonData.data.matInvF);
745 commonData.data.vecDetF.resize(nb_gauss_pts, false);
747
749 dP_dF.resize(81, nb_gauss_pts, false);
751 matF.resize(9, nb_gauss_pts, false);
754 commonData.data.matP.resize(9, nb_gauss_pts, false);
755 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
756 commonData.data.vecPsi.resize(nb_gauss_pts, false);
764 const double lambda = commonData.lambda;
765 const double mu = commonData.mu;
766
767 for (int gg = 0; gg != nb_gauss_pts; gg++) {
768
769 // Get deformation gradient
770 F(i, I) = grad(i, I);
771 F(0, 0) += 1;
772 F(1, 1) += 1;
773 F(2, 2) += 1;
774
776 CHKERR invertTensor3by3(F, det_f, invF);
777
778 const double log_det = log(det_f);
779 psi = F(i, J) * F(i, J) - 3 - 2 * log_det;
780 psi *= 0.5 * mu;
781 psi += 0.5 * lambda * log_det * log_det;
782
783 const double var = lambda * log_det - mu;
784 P(i, J) = mu * F(i, J) + var * invF(J, i);
785
786 const double coef = mu - lambda * log_det;
787 // TODO: This can be improved by utilising the symmetries of the D2 tensor
788
789 D2(i, J, k, L) =
790 invF(J, i) * (invF(L, k) * lambda) + invF(L, i) * (invF(J, k) * coef);
791
792 D2(0, 0, 0, 0) += mu;
793 D2(0, 1, 0, 1) += mu;
794 D2(0, 2, 0, 2) += mu;
795 D2(1, 0, 1, 0) += mu;
796 D2(1, 1, 1, 1) += mu;
797 D2(1, 2, 1, 2) += mu;
798 D2(2, 0, 2, 0) += mu;
799 D2(2, 1, 2, 1) += mu;
800 D2(2, 2, 2, 2) += mu;
801
802 ++det_f;
803 ++invF;
804 ++grad;
805 ++F;
806 ++D2;
807 ++P;
808 ++psi;
809 }
810
812}
813
816 "DISPLACEMENTS", UserDataOperator::OPROW),
817 commonData(common_data) {}
818
821 DataForcesAndSourcesCore::EntData &data) {
822
824
825 const int nb_dofs = data.getIndices().size();
826 if (!nb_dofs)
828 nF.resize(nb_dofs, false);
829 nF.clear();
830
831 const int nb_gauss_pts = data.getN().size1();
832 const int nb_base_functions = data.getN().size2();
833 auto diff_base_functions = data.getFTensor1DiffN<3>();
834 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
836
837 // MatrixDouble &matS = commonData.data.matS;
838 // MatrixDouble &matF = commonData.data.matGradientOfDeformation;
839 // FTensor::Tensor2<double*,3,3> F(TENSOR2_MAT_PTR(matF));
840 // FTensor::Tensor2_symmetric<double*,3> S(SYMMETRIC_TENSOR2_MAT_PTR(matS),6);
841
842 const double n = commonData.n;
843 const double rho_ref = commonData.rHo_ref;
844
848
849 for (int gg = 0; gg != nb_gauss_pts; gg++) {
850
851 // Get volume and integration weight
852 double w = getVolume() * getGaussPts()(3, gg);
853 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)
877 getFEMethod()->ts_F, /// Get the right hand side vector for time solver
878 nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
879
881}
882
885 "RHO", UserDataOperator::OPROW),
886 commonData(common_data) {}
887
890 DataForcesAndSourcesCore::EntData &data) {
891
893
894 const int nb_dofs = data.getIndices().size();
895 if (!nb_dofs)
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 }
909 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "gradRhoPtr is null");
910 }
911 auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
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 }
921
922 const double R0 = commonData.R0;
923
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)
955 getFEMethod()->ts_F, /// Get the right hand side vector for time solver
956 nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
957
959}
960
963 "RHO", "RHO", UserDataOperator::OPROWCOL),
964 commonData(common_data) {
965 sYmm = true;
966}
967
969OpAssmbleRhoLhs_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
975
976 const int row_nb_dofs = row_data.getIndices().size();
977 if (!row_nb_dofs)
979 const int col_nb_dofs = col_data.getIndices().size();
980 if (!col_nb_dofs)
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)
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
1001
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
1064}
1065
1068 "RHO", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
1069 commonData(common_data) {
1070 sYmm = false;
1071}
1072
1074OpAssmbleRhoLhs_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
1080
1081 const int row_nb_dofs = row_data.getIndices().size();
1082 if (!row_nb_dofs)
1084 const int col_nb_dofs = col_data.getIndices().size();
1085 if (!col_nb_dofs)
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
1102
1104 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
1105
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
1143}
1144template <bool ADOLC>
1146 Remodeling::CommonData &common_data)
1148 "DISPLACEMENTS", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
1149 commonData(common_data) {
1150 sYmm = true;
1151}
1152template <bool ADOLC>
1154 int row_side, int col_side, EntityType row_type, EntityType col_type,
1155 DataForcesAndSourcesCore::EntData &row_data,
1156 DataForcesAndSourcesCore::EntData &col_data) {
1157
1159
1160 const int row_nb_dofs = row_data.getIndices().size();
1161 if (!row_nb_dofs)
1163 const int col_nb_dofs = col_data.getIndices().size();
1164 if (!col_nb_dofs)
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
1193
1194 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
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
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++) {
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));
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
1305}
1306
1308 Remodeling::CommonData &common_data)
1310 "DISPLACEMENTS", "RHO", UserDataOperator::OPROWCOL),
1311 commonData(common_data) {
1312 sYmm = false; // This is off-diagonal (upper-left), so no symmetric
1313}
1314
1316OpAssmbleStressLhs_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
1322
1323 const int row_nb_dofs = row_data.getIndices().size();
1324 if (!row_nb_dofs)
1326 const int col_nb_dofs = col_data.getIndices().size();
1327 if (!col_nb_dofs)
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
1335
1336 const double n = commonData.n;
1337 const double rho_ref = commonData.rHo_ref;
1338
1341
1342 auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
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
1385}
1386
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
1397}
1398
1402}
1403
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", "",
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
1444 CHKERR postProc.addFieldValuesPostProc("DISPLACEMENTS");
1445 CHKERR postProc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
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(
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(
1465 }
1466
1469 CHKERR postProcElastic.addFieldValuesPostProc("MESH_NODE_POSITIONS");
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(
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(
1497 energyCalc.getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<3>(
1498 "RHO", commonData.data.gradRhoPtr));
1499 energyCalc.getOpPtrVector().push_back(
1500 new OpCalculateVectorFieldGradient<3, 3>("DISPLACEMENTS",
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";
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
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 }
1577}
1578
1580
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);
1593 string name = it->getName();
1594 if (name.compare(0, 14, "NO_REMODELLING") == 0) {
1596 EntityHandle meshset = it->getMeshset();
1597 CHKERR this->mField.get_moab().get_entities_by_type(
1598 meshset, MBTET, commonData.tEts_block, true);
1601 }
1602 }
1603
1605}
1606
1608
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")) {
1631 MB_TAG_SPARSE, MF_ZERO);
1632 // FIXME
1634 // CHKERR mField.add_ents_to_field_by_type(0,MBTET,"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
1673
1674 // If order was not set from CT scans set homogenous order equal to
1675 // reference bone density
1676 if (add_rho_field) {
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
1692}
1693
1696
1697 CHKERR mField.add_finite_element("FE_REMODELLING", MF_ZERO);
1699 "DISPLACEMENTS");
1701 "DISPLACEMENTS");
1703 "DISPLACEMENTS");
1704 CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING", "RHO");
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");
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");
1716 "DISPLACEMENTS");
1718 "MESH_NODE_POSITIONS");
1719
1722 MBTET, "ELASTIC");
1723 }
1724
1726 boost::shared_ptr<ElasticMaterials>(new ElasticMaterials(mField));
1727 commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
1730 commonData.elasticPtr->setOfBlocks);
1731 CHKERR commonData.elasticPtr->addElement("ELASTIC", "DISPLACEMENTS");
1732 CHKERR commonData.elasticPtr->setOperators(
1733 "DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
1734
1737
1738 // Allocate memory for density and gradient of displacements at integration
1739 // points
1740 commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
1742 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
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(
1753 list_of_operators_rhs.push_back(
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(
1770 list_of_operators_rhs.push_back(
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(
1778 list_of_operators_lhs.push_back(
1780 } else {
1781 list_of_operators_lhs.push_back(new OpCalculateStressTangent(commonData));
1782 list_of_operators_lhs.push_back(
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
1790}
1791
1793
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
1805 PETSC_NULL, "DISPLACEMENTS");
1806 // edge forces
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
1830}
1831
1833
1835 commonData.dm_name = "DM_REMODELING";
1836 // Register DM problem
1838 CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
1841 // Create DM instance
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");
1850
1851 if (mField.check_finite_element("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
1861}
1862
1864
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
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
1945}
1946
1948 int row_side, EntityType row_type,
1949 DataForcesAndSourcesCore::EntData &row_data) {
1951
1952 if (row_type != MBVERTEX)
1954
1955 const int nb_gauss_pts = row_data.getN().size1();
1956 // commonData.data.vecPsi.resize(nb_gauss_pts,false);
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);
1975}
1976
1977} // namespace BoneRemodeling
static Index< 'J', 3 > J
Elastic materials.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)
Definition: definitions.h:626
@ ROW
Definition: definitions.h:123
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT)
Definition: definitions.h:623
#define TENSOR2_MAT_PTR(MAT)
Definition: definitions.h:619
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ 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 ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT)
Definition: definitions.h:609
@ 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()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
#define TENSOR4_MAT_PTR(MAT)
Definition: definitions.h:617
constexpr int order
#define MU(E, NU)
Definition: fem_tools.h:23
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
@ R
@ F
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
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:118
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:412
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1128
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
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:572
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
const double c
speed of light (cm/ns)
double D
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
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:259
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:1559
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.
Definition: Templates.hpp:135
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1511
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)
Definition: enable_if.hpp:6
constexpr IntegrationType I
constexpr AssemblyType A
double young_modulus
Young modulus.
Definition: plastic.cpp:172
double rho
Definition: plastic.cpp:191
Definition: Remodeling.hpp:553
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:556
int pRT
Definition: Remodeling.hpp:564
MoFEMErrorCode postProcess()
MoFEMErrorCode operator()()
MoFEMErrorCode preProcess()
MoFEM::Interface & mField
Definition: Remodeling.hpp:555
PetscBool mass_postproc
Definition: Remodeling.hpp:560
PetscBool equilibrium_flg
Definition: Remodeling.hpp:561
double rate
Definition: Remodeling.hpp:562
bool iNit
Definition: Remodeling.hpp:563
PostProcVolumeOnRefinedMesh postProcElastic
Definition: Remodeling.hpp:558
PostProcVolumeOnRefinedMesh postProc
Definition: Remodeling.hpp:557
MonitorPostProc(MoFEM::Interface &m_field, Remodeling::CommonData &common_data)
Off diagonal block of tangent matrix .
Definition: Remodeling.hpp:494
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data)
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:496
Diagonal block of tangent matrix .
Definition: Remodeling.hpp:473
OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:961
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:475
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
Assemble residual for conservation of mass (density)
Definition: Remodeling.hpp:454
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:889
OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:883
VectorDouble nF
Vector of the right hand side (internal forces)
Definition: Remodeling.hpp:457
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:456
Off diagonal block of tangent matrix .
Definition: Remodeling.hpp:514
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssmbleStressLhs_dF(Remodeling::CommonData &common_data)
Off diagonal block of tangent matrix /f[ K_{u \rho}=\intop_{V} \left[\frac{n}{\rho_{0}}\right] \left...
Definition: Remodeling.hpp:537
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:539
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssmbleStressLhs_dRho(Remodeling::CommonData &common_data)
Assemble residual for conservation of momentum (stresses)
Definition: Remodeling.hpp:436
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:438
VectorDouble nF
Vector of the right hand side (internal forces)
Definition: Remodeling.hpp:439
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:820
OpAssmbleStressRhs(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:814
Evaluate physical equations at integration points.
Definition: Remodeling.hpp:357
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:359
OpCalculateStress(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:280
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:294
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:732
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:407
OpCalculateStressTangent(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:717
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:396
OpCalculateStressTangentWithAdolc(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:382
Evaluate density derivative with respect to time in case of Backward Euler Method.
Definition: Remodeling.hpp:333
OpGetRhoTimeDirevative(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:232
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:239
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:335
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Used to post proc stresses, energy, source term.
Definition: Remodeling.hpp:375
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:567
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:550
Remodeling::CommonData & commonData
Definition: Remodeling.hpp:377
moab::Interface & postProcMesh
Definition: Remodeling.hpp:378
std::vector< EntityHandle > & mapGaussPts
Definition: Remodeling.hpp:379
VectorDouble vecPsi
Elastic energy density.
Definition: Remodeling.hpp:189
boost::shared_ptr< MatrixDouble > gradDispPtr
Ptr to gradient of displacements matrix container.
Definition: Remodeling.hpp:185
boost::shared_ptr< VectorDouble > rhoPtr
Ptr to density matrix container.
Definition: Remodeling.hpp:187
boost::shared_ptr< MatrixDouble > gradRhoPtr
Gradient of density field.
Definition: Remodeling.hpp:188
MatrixDouble matInvF
inverse of deformation gradient
Definition: Remodeling.hpp:199
VectorDouble vecRhoDt
Time derivative of density.
Definition: Remodeling.hpp:201
MatrixDouble matGradientOfDeformation
Gradient of deformation.
Definition: Remodeling.hpp:198
double getCFromDensity(const double &rho)
Definition: Remodeling.cpp:218
double pSi_ref
reference free energy
Definition: Remodeling.hpp:154
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
Definition: Remodeling.hpp:130
double n
porosity exponent [-]
Definition: Remodeling.hpp:148
int b
b exponent for bell function
Definition: Remodeling.hpp:153
double getCFromDensityDiff(const double &rho)
Definition: Remodeling.cpp:225
double m
algorithmic exponent [-]
Definition: Remodeling.hpp:147
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
Definition: Remodeling.hpp:131
PetscBool is_atom_testing
for atom tests
Definition: Remodeling.hpp:162
double c
density evolution (growth) velocity [d/m^2]
Definition: Remodeling.hpp:146
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
Definition: Remodeling.hpp:136
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
Definition: Remodeling.hpp:138
PetscBool less_post_proc
reduce file size
Definition: Remodeling.hpp:163
DM dm
Discretization manager.
Definition: Remodeling.hpp:127
boost::shared_ptr< ElasticMaterials > elasticMaterialsPtr
Definition: Remodeling.hpp:141
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
Definition: Remodeling.hpp:137
double R0
mass conduction coefficient
Definition: Remodeling.hpp:155
double cUrrent_psi
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:158
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Definition: Remodeling.hpp:140
double cUrrent_mass
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:160
Volume finite element.
Definition: Remodeling.hpp:41
MoFEM::Interface & mField
Definition: Remodeling.hpp:206
MoFEMErrorCode addElements()
Set and add finite elements.
MoFEMErrorCode solveDM()
Solve problem set up in DM.
MoFEMErrorCode buildDM()
Set problem and DM.
MoFEMErrorCode getParameters()
Get parameters form line command or config file.
MoFEMErrorCode addFields()
Set and add entities to approximation fields.
MoFEMErrorCode addMomentumFluxes()
Finite elements to calculate tractions.
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
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.
Definition: NodalForce.hpp:128
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
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual int get_comm_rank() const =0
bool & doTris
\deprectaed
bool & doPrisms
\deprectaed
bool sYmm
If true assume that matrix is symmetric structure.
bool & doVertices
\deprectaed If false skip vertices
bool & doTets
\deprectaed
bool & doEdges
\deprectaed If false skip edges
bool & doQuads
\deprectaed
Deprecated interface functions.
Basic algebra on fields.
Definition: FieldBlas.hpp:21
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
keeps basic data about problem
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:17
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
Definition: TsCtx.hpp:144
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
structure grouping operators and data used for calculation of nonlinear elastic element
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
CommonData commonData
Force scale operator for reading two columns.