v0.15.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 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 PetscOptionsEnd();
213
215}
216
217inline double Remodeling::CommonData::getCFromDensity(const double &rho) {
218 double mid_rho = 0.5 * (rHo_max + rHo_min);
219 double var_h = (rho - mid_rho) / (rHo_max - mid_rho);
220
221 return 1 / (1 + (b != 0) * pow(var_h, 2 * b));
222}
223
225 double mid_rho = 0.5 * (rHo_max + rHo_min);
226 double var_h = (rho - mid_rho) / (rHo_max - mid_rho);
227 return (b != 0) * (-2) * b * pow(var_h, 2 * b - 1) /
228 ((rHo_max - mid_rho) * pow(pow(var_h, 2 * b) + 1, 2));
229}
230
236
238OpGetRhoTimeDirevative::doWork(int side, EntityType type,
239 DataForcesAndSourcesCore::EntData &data) {
240
242
243 const int nb_gauss_pts = data.getN().size1();
244 if (!nb_gauss_pts)
246 const int nb_base_functions = data.getN().size2();
247
248 auto base_function = data.getFTensor0N();
249 commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
250 if (type == MBVERTEX) {
251 commonData.data.vecRhoDt.clear();
252 }
253
254 int nb_dofs = data.getIndices().size();
255 if (nb_dofs == 0)
257 double *array;
258 CHKERR VecGetArray(getFEMethod()->ts_u_t, &array);
260
261 for (int gg = 0; gg < nb_gauss_pts; gg++) {
262 int bb = 0;
263 for (; bb < nb_dofs; bb++) {
264 const double field_data = array[data.getLocalIndices()[bb]];
265 drho_dt += field_data * base_function;
266 ++base_function;
267 }
268 // It is possible to have more base functions than dofs
269 for (; bb != nb_base_functions; bb++)
270 ++base_function;
271 ++drho_dt;
272 }
273
274 CHKERR VecRestoreArray(getFEMethod()->ts_u_t, &array);
275
277}
278
281 "DISPLACEMENTS", UserDataOperator::OPROW),
282 commonData(common_data) {
283 // Set that this operator is executed only for vertices on element
284 doVertices = true;
285 doEdges = false;
286 doQuads = false;
287 doTris = false;
288 doTets = false;
289 doPrisms = false;
290}
291
293OpCalculateStress::doWork(int side, EntityType type,
294 DataForcesAndSourcesCore::EntData &data) {
295
297
298 if (type != MBVERTEX)
300
302 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
303 "Gradient at integration points of displacement is not calculated");
304 }
306
307 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
308 if (!commonData.data.rhoPtr) {
309 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
310 "Density at integration points is not calculated");
311 }
312
314
315 commonData.data.vecPsi.resize(nb_gauss_pts, false);
317 commonData.data.vecR.resize(nb_gauss_pts, false);
319
320 commonData.data.matInvF.resize(9, nb_gauss_pts, false);
322 commonData.data.vecDetF.resize(nb_gauss_pts, false);
324 MatrixDouble &matS = commonData.data.matS;
325 matS.resize(nb_gauss_pts, 6, false);
327 commonData.data.matP.resize(9, nb_gauss_pts, false);
329 MatrixDouble &matF = commonData.data.matGradientOfDeformation;
330 matF.resize(9, nb_gauss_pts, false);
332
333 const double c = commonData.c;
334 const double n = commonData.n;
335 const double m = commonData.m;
336 const double rho_ref = commonData.rHo_ref;
337 const double psi_ref = commonData.pSi_ref;
338 const double lambda = commonData.lambda;
339 const double mu = commonData.mu;
340 FTensor::Index<'i', 3> i;
341 FTensor::Index<'I', 3> I;
342 FTensor::Index<'J', 3> J;
343
344 for (int gg = 0; gg != nb_gauss_pts; gg++) {
345
346 // Get deformation gradient
347 F(i, I) = grad(i, I);
348 F(0, 0) += 1;
349 F(1, 1) += 1;
350 F(2, 2) += 1;
351
353 CHKERR invertTensor3by3(F, det_f, invF);
354
355 const double log_det = log(det_f);
356 psi = F(i, J) * F(i, J) - 3 - 2 * log_det;
357 psi *= 0.5 * mu;
358 psi += 0.5 * lambda * log_det * log_det;
359
360 const double coef = lambda * log_det - mu;
361 P(i, J) = mu * F(i, J) + coef * invF(J, i);
362 S(i, I) = invF(i, J) ^ P(J, I);
363
364 const double kp = commonData.getCFromDensity(rho);
365 R = c * kp * (pow(rho / rho_ref, n - m) * psi - psi_ref);
366
367 ++det_f;
368 ++invF;
369 ++grad;
370 ++rho;
371 ++psi;
372 ++S;
373 ++P;
374 ++R;
375 ++F;
376 }
377
379}
380
382 Remodeling::CommonData &common_data)
384 "DISPLACEMENTS", UserDataOperator::OPROW),
385 commonData(common_data) {
386 // Set that this operator is executed only for vertices on element
387 doVertices = true;
388 doEdges = false;
389 doQuads = false;
390 doTris = false;
391 doTets = false;
392 doPrisms = false;
393}
394
396 int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
397
399
400 if (type != MBVERTEX)
402
404 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
405
406 vecC.resize(6, false);
408 MatrixDouble &dS_dC = commonData.data.matMaterialTangent;
409 dS_dC.resize(6, 6 * nb_gauss_pts, false);
410 dS_dC.clear();
412
413 MatrixDouble &matS = commonData.data.matS;
414 matS.resize(nb_gauss_pts, 6, false);
416 MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
417 dP_dF.resize(81, nb_gauss_pts, false);
418 MatrixDouble &matF = commonData.data.matGradientOfDeformation;
419 matF.resize(9, nb_gauss_pts, false);
422 commonData.data.matP.resize(9, nb_gauss_pts, false);
424 commonData.data.vecPsi.resize(nb_gauss_pts, false);
426
427 FTensor::Index<'i', 3> i;
428 FTensor::Index<'k', 3> k;
429 FTensor::Index<'I', 3> I;
430 FTensor::Index<'J', 3> J;
431 FTensor::Index<'K', 3> K;
432 FTensor::Index<'L', 3> L;
433
434 for (int gg = 0; gg != nb_gauss_pts; gg++) {
435
436 // Get deformation gradient
437 F(i, I) = grad(i, I);
438 F(0, 0) += 1;
439 F(1, 1) += 1;
440 F(2, 2) += 1;
441
442 // Calculate Green defamation Tensor0
443 // ^ that means multiplication of Tensor 2 which result in symmetric Tensor
444 // rank 2.
445 C(I, J) = F(i, I) ^ F(i, J);
446
447 // Calculate free energy
450 commonData, C, psi);
451
452 int r_g = ::gradient(commonData.tAg, 6, &vecC[0], &matS(gg, 0));
453 if (r_g < 0) {
454 // That means that energy function is not smooth and derivative
455 // can not be calculated,
456 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
457 "ADOL-C function evaluation with error");
458 }
459
460 // Scale energy density and 2nd Piola Stress
461 S(0, 0) *= 2;
462 S(1, 1) *= 2;
463 S(2, 2) *= 2;
464
465 // Calculate 1st Piola-Stress tensor
466 P(i, J) = F(i, I) * S(I, J);
467
468 const int is = 6 * gg;
469 double *hessian_ptr[6] = {&dS_dC(0, is), &dS_dC(1, is), &dS_dC(2, is),
470 &dS_dC(3, is), &dS_dC(4, is), &dS_dC(5, is)};
471 int r_h = ::hessian(commonData.tAg, 6, &vecC[0], hessian_ptr);
472 if (r_h < 0) {
473 // That means that energy function is not smooth and derivative
474 // can not be calculated,
475 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
476 "ADOL-C function evaluation with error");
477 }
478
479 // FIXME: TODO: Here tensor 4 with two minor symmetries is used, since
480 // \psi_{,ijkl} = \psi_{,klij} = \psi_{,lkji} = ... there is additional
481 // major symmetry. That is 21 independent components not 36 how is now
482 // implemented using Tenso4_ddg.
483 //
484 // That way we have to fill off-diagonal part for dS_dC. In future that
485 // need to be fixed by implementation of Tenso4 with additional major
486 // symmetries. Pleas consider doing such implementation by yourself. You
487 // have support of MoFEM developing team to guide you through this process.
488 //
489 // This deteriorate efficiency.
490 //
491 // TODO: Implement matrix with more symmetries.
492 //
493 for (int ii = 0; ii < 6; ii++) {
494 for (int jj = 0; jj < ii; jj++) {
495 dS_dC(jj, is + ii) = dS_dC(ii, is + jj);
496 }
497 }
498
499 // As result that symmetry of deformation gradient C is used, terms which
500 // have mixed indices are two time to big, and terms which has two times
501 // midxed index are four times to big. Those terms are not multiplied.
502 D1(0, 0, 0, 0) *= 4;
503 D1(1, 1, 1, 1) *= 4;
504 D1(2, 2, 2, 2) *= 4;
505
506 D1(0, 0, 1, 1) *= 4;
507 D1(1, 1, 0, 0) *= 4;
508 D1(1, 1, 2, 2) *= 4;
509 D1(2, 2, 1, 1) *= 4;
510 D1(0, 0, 2, 2) *= 4;
511 D1(2, 2, 0, 0) *= 4;
512
513 D1(0, 0, 0, 1) *= 2;
514 D1(0, 0, 0, 2) *= 2;
515 D1(0, 0, 1, 2) *= 2;
516 D1(0, 1, 0, 0) *= 2;
517 D1(0, 2, 0, 0) *= 2;
518 D1(1, 2, 0, 0) *= 2;
519
520 D1(1, 1, 0, 1) *= 2;
521 D1(1, 1, 0, 2) *= 2;
522 D1(1, 1, 1, 2) *= 2;
523 D1(0, 1, 1, 1) *= 2;
524 D1(0, 2, 1, 1) *= 2;
525 D1(1, 2, 1, 1) *= 2;
526
527 D1(2, 2, 0, 1) *= 2;
528 D1(2, 2, 0, 2) *= 2;
529 D1(2, 2, 1, 2) *= 2;
530 D1(0, 1, 2, 2) *= 2;
531 D1(0, 2, 2, 2) *= 2;
532 D1(1, 2, 2, 2) *= 2;
533
534 D2(i, J, k, L) = (F(i, I) * (D1(I, J, K, L) * F(k, K)));
535 // D2(i, J, k, L) += I(i, k) * S(J, L);
536
537 ++grad;
538 ++S;
539 ++F;
540 ++D1;
541 ++D2;
542 ++P;
543 ++psi;
544 }
545
547}
548
549OpPostProcStress::OpPostProcStress(moab::Interface &post_proc_mesh,
550 std::vector<EntityHandle> &map_gauss_pts,
551 Remodeling::CommonData &common_data)
553 "DISPLACEMENTS", UserDataOperator::OPROW),
554 commonData(common_data), postProcMesh(post_proc_mesh),
555 mapGaussPts(map_gauss_pts) {
556 // Set that this operator is executed only for vertices on element
557 doVertices = true;
558 doEdges = false;
559 doQuads = false;
560 doTris = false;
561 doTets = false;
562 doPrisms = false;
563}
564
566OpPostProcStress::doWork(int side, EntityType type,
567 DataForcesAndSourcesCore::EntData &data) {
568
570
571 if (type != MBVERTEX)
573 double def_VAL[9];
574 bzero(def_VAL, 9 * sizeof(double));
575
576 Tag th_piola2;
577 CHKERR postProcMesh.tag_get_handle("Piola_Stress2", 9, MB_TYPE_DOUBLE,
578 th_piola2, MB_TAG_CREAT | MB_TAG_SPARSE,
579 def_VAL);
580
581 Tag th_psi;
582 CHKERR postProcMesh.tag_get_handle("psi", 1, MB_TYPE_DOUBLE, th_psi,
583 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
584
585 Tag th_rho;
586 CHKERR postProcMesh.tag_get_handle("RHO", 1, MB_TYPE_DOUBLE, th_rho,
587 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
588
589 Tag th_piola1;
590 Tag principal;
591 Tag th_prin_stress_vect1, th_prin_stress_vect2, th_prin_stress_vect3;
592 Tag th_R;
593 Tag th_grad_rho;
595
596 CHKERR postProcMesh.tag_get_handle("Piola_Stress1", 9, MB_TYPE_DOUBLE,
597 th_piola1, MB_TAG_CREAT | MB_TAG_SPARSE,
598 def_VAL);
599
600 CHKERR postProcMesh.tag_get_handle("Principal_stress", 3, MB_TYPE_DOUBLE,
601 principal, MB_TAG_CREAT | MB_TAG_SPARSE,
602 def_VAL);
603
604 CHKERR postProcMesh.tag_get_handle("Principal_stress_S1", 3, MB_TYPE_DOUBLE,
605 th_prin_stress_vect1,
606 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
607 CHKERR postProcMesh.tag_get_handle("Principal_stress_S2", 3, MB_TYPE_DOUBLE,
608 th_prin_stress_vect2,
609 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
610 CHKERR postProcMesh.tag_get_handle("Principal_stress_S3", 3, MB_TYPE_DOUBLE,
611 th_prin_stress_vect3,
612 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
613
614 CHKERR postProcMesh.tag_get_handle("R", 1, MB_TYPE_DOUBLE, th_R,
615 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
616
617 CHKERR postProcMesh.tag_get_handle("grad_rho", 3, MB_TYPE_DOUBLE,
618 th_grad_rho,
619 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
620 }
621
622 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
623 MatrixDouble &matS = commonData.data.matS;
628
631
632 MatrixDouble3by3 stress(3, 3);
633 VectorDouble3 eigen_values(3);
634 MatrixDouble3by3 eigen_vectors(3, 3);
635
636 for (int gg = 0; gg != nb_gauss_pts; gg++) {
637
638 double val = psi;
639 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &val);
640 val = rho;
641 CHKERR postProcMesh.tag_set_data(th_rho, &mapGaussPts[gg], 1, &val);
642
643 // S is symmetric tensor, need to map it to non-symmetric tensor to save it
644 // on moab tag to be viable in paraview.
645 for (int ii = 0; ii != 3; ii++) {
646 for (int jj = 0; jj != 3; jj++) {
647 stress(ii, jj) = S(ii, jj);
648 }
649 }
650 CHKERR postProcMesh.tag_set_data(th_piola2, &mapGaussPts[gg], 1,
651 &stress(0, 0));
652
654
655 for (int ii = 0; ii != 3; ii++) {
656 for (int jj = 0; jj != 3; jj++) {
657 stress(ii, jj) = P(ii, jj);
658 }
659 }
660 CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
661 &stress(0, 0));
662 val = R;
663 CHKERR postProcMesh.tag_set_data(th_R, &mapGaussPts[gg], 1, &val);
664 const double t2[] = {grad_rho(0), grad_rho(1), grad_rho(2)};
665 CHKERR postProcMesh.tag_set_data(th_grad_rho, &mapGaussPts[gg], 1, t2);
666
667 // Add calculation of eigen values, i.e. prinple stresses.
668 noalias(eigen_vectors) = stress;
669 VectorDouble3 prin_stress_vect1(3);
670 VectorDouble3 prin_stress_vect2(3);
671 VectorDouble3 prin_stress_vect3(3);
672 // LAPACK - eigenvalues and vectors. Applied twice for initial creates
673 // memory space
674 int n = 3, lda = 3, info, lwork = -1;
675 double wkopt;
676 info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
677 &(eigen_values.data()[0]), &wkopt, lwork);
678 if (info != 0)
679 SETERRQ1(PETSC_COMM_SELF, 1,
680 "something is wrong with lapack_dsyev info = %d", info);
681 lwork = (int)wkopt;
682 double work[lwork];
683 info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
684 &(eigen_values.data()[0]), work, lwork);
685 if (info != 0)
686 SETERRQ1(PETSC_COMM_SELF, 1,
687 "something is wrong with lapack_dsyev info = %d", info);
688
689 for (int ii = 0; ii < 3; ii++) {
690 prin_stress_vect1[ii] = eigen_vectors(0, ii);
691 prin_stress_vect2[ii] = eigen_vectors(1, ii);
692 prin_stress_vect3[ii] = eigen_vectors(2, ii);
693 }
694
695 CHKERR postProcMesh.tag_set_data(principal, &mapGaussPts[gg], 1,
696 &eigen_values[0]);
697 CHKERR postProcMesh.tag_set_data(th_prin_stress_vect1, &mapGaussPts[gg],
698 1, &*prin_stress_vect1.begin());
699 CHKERR postProcMesh.tag_set_data(th_prin_stress_vect2, &mapGaussPts[gg],
700 1, &*prin_stress_vect2.begin());
701 CHKERR postProcMesh.tag_set_data(th_prin_stress_vect3, &mapGaussPts[gg],
702 1, &*prin_stress_vect3.begin());
703 }
704
705 ++S;
706 ++P;
707 ++psi;
708 ++R;
709 ++rho;
710 ++grad_rho;
711 }
712
714}
715
717 Remodeling::CommonData &common_data)
719 "DISPLACEMENTS", UserDataOperator::OPROW),
720 commonData(common_data) {
721 // Set that this operator is executed only for vertices on element
722 doVertices = true;
723 doEdges = false;
724 doQuads = false;
725 doTris = false;
726 doTets = false;
727 doPrisms = false;
728}
729
731OpCalculateStressTangent::doWork(int side, EntityType type,
732 DataForcesAndSourcesCore::EntData &data) {
733
735
736 if (type != MBVERTEX)
738
740 const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
741
742 commonData.data.matInvF.resize(9, nb_gauss_pts, false);
744 commonData.data.vecDetF.resize(nb_gauss_pts, false);
746
747 MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
748 dP_dF.resize(81, nb_gauss_pts, false);
749 MatrixDouble &matF = commonData.data.matGradientOfDeformation;
750 matF.resize(9, nb_gauss_pts, false);
753 commonData.data.matP.resize(9, nb_gauss_pts, false);
755 commonData.data.vecPsi.resize(nb_gauss_pts, false);
757 FTensor::Index<'i', 3> i;
758 FTensor::Index<'k', 3> k;
759 FTensor::Index<'I', 3> I;
760 FTensor::Index<'J', 3> J;
761 FTensor::Index<'K', 3> K;
762 FTensor::Index<'L', 3> L;
763 const double lambda = commonData.lambda;
764 const double mu = commonData.mu;
765
766 for (int gg = 0; gg != nb_gauss_pts; gg++) {
767
768 // Get deformation gradient
769 F(i, I) = grad(i, I);
770 F(0, 0) += 1;
771 F(1, 1) += 1;
772 F(2, 2) += 1;
773
775 CHKERR invertTensor3by3(F, det_f, invF);
776
777 const double log_det = log(det_f);
778 psi = F(i, J) * F(i, J) - 3 - 2 * log_det;
779 psi *= 0.5 * mu;
780 psi += 0.5 * lambda * log_det * log_det;
781
782 const double var = lambda * log_det - mu;
783 P(i, J) = mu * F(i, J) + var * invF(J, i);
784
785 const double coef = mu - lambda * log_det;
786 // TODO: This can be improved by utilising the symmetries of the D2 tensor
787
788 D2(i, J, k, L) =
789 invF(J, i) * (invF(L, k) * lambda) + invF(L, i) * (invF(J, k) * coef);
790
791 D2(0, 0, 0, 0) += mu;
792 D2(0, 1, 0, 1) += mu;
793 D2(0, 2, 0, 2) += mu;
794 D2(1, 0, 1, 0) += mu;
795 D2(1, 1, 1, 1) += mu;
796 D2(1, 2, 1, 2) += mu;
797 D2(2, 0, 2, 0) += mu;
798 D2(2, 1, 2, 1) += mu;
799 D2(2, 2, 2, 2) += mu;
800
801 ++det_f;
802 ++invF;
803 ++grad;
804 ++F;
805 ++D2;
806 ++P;
807 ++psi;
808 }
809
811}
812
815 "DISPLACEMENTS", UserDataOperator::OPROW),
816 commonData(common_data) {}
817
819OpAssmbleStressRhs::doWork(int side, EntityType type,
820 DataForcesAndSourcesCore::EntData &data) {
821
823
824 const int nb_dofs = data.getIndices().size();
825 if (!nb_dofs)
827 nF.resize(nb_dofs, false);
828 nF.clear();
829
830 const int nb_gauss_pts = data.getN().size1();
831 const int nb_base_functions = data.getN().size2();
832 auto diff_base_functions = data.getFTensor1DiffN<3>();
835
836 // MatrixDouble &matS = commonData.data.matS;
837 // MatrixDouble &matF = commonData.data.matGradientOfDeformation;
838 // FTensor::Tensor2<double*,3,3> F(TENSOR2_MAT_PTR(matF));
839 // FTensor::Tensor2_symmetric<double*,3> S(SYMMETRIC_TENSOR2_MAT_PTR(matS),6);
840
841 const double n = commonData.n;
842 const double rho_ref = commonData.rHo_ref;
843
844 FTensor::Index<'I', 3> I;
845 FTensor::Index<'J', 3> J;
846 FTensor::Index<'i', 3> i;
847
848 for (int gg = 0; gg != nb_gauss_pts; gg++) {
849
850 // Get volume and integration weight
851 double w = getVolume() * getGaussPts()(3, gg);
852 double a = w * pow(rho / rho_ref, n);
853
854 int bb = 0;
855 for (; bb != nb_dofs / 3; bb++) {
856 double *ptr = &nF[3 * bb];
857 FTensor::Tensor1<double *, 3> t1(ptr, &ptr[1], &ptr[2]);
858 t1(i) += a * P(i, I) * diff_base_functions(I);
859 // t1(i) += 0.5*a*diff_base_functions(J)*(F(i,I)*S(I,J));
860 // t1(i) += 0.5*a*diff_base_functions(I)*(F(i,J)*S(I,J));
861 ++diff_base_functions;
862 }
863 // Could be that we base more base functions than needed to approx. disp.
864 // field.
865 for (; bb != nb_base_functions; bb++) {
866 ++diff_base_functions;
867 }
868 ++P;
869 // ++S;
870 // ++F;
871 ++rho;
872 }
873
874 // Assemble internal force vector for time solver (TS)
876 getFEMethod()->ts_F, /// Get the right hand side vector for time solver
877 nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
878
880}
881
886
888OpAssmbleRhoRhs::doWork(int side, EntityType type,
889 DataForcesAndSourcesCore::EntData &data) {
890
892
893 const int nb_dofs = data.getIndices().size();
894 if (!nb_dofs)
896 nF.resize(nb_dofs, false);
897 nF.clear();
898
899 const int nb_gauss_pts = data.getN().size1();
900 const int nb_base_functions = data.getN().size2();
901 auto base_functions = data.getFTensor0N();
902 auto diff_base_functions = data.getFTensor1DiffN<3>();
903 if (!commonData.data.rhoPtr) {
904 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "rhoPtsr is null");
905 }
908 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "gradRhoPtr is null");
909 }
912
913 if (commonData.data.vecRhoDt.size() != nb_gauss_pts) {
914 commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
915 commonData.data.vecRhoDt.clear();
916 // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Time derivative of
917 // density at integration not calculated");
918 }
920
921 const double R0 = commonData.R0;
922
923 FTensor::Index<'I', 3> I;
924
925 for (int gg = 0; gg != nb_gauss_pts; gg++) {
926
927 // Get volume and integration weight
928 double w = getVolume() * getGaussPts()(3, gg);
929
930 int bb = 0;
931 for (; bb != nb_dofs; bb++) {
932 nF[bb] += w * base_functions * drho_dt; // Density rate term
933 nF[bb] -= w * base_functions * R; // Source term
934 nF[bb] -= w * R0 * diff_base_functions(I) * grad_rho(I); // Gradient term
935 ++base_functions;
936 ++diff_base_functions;
937 }
938 // Could be that we base more base functions than needed to approx. density
939 // field.
940 for (; bb != nb_base_functions; bb++) {
941 ++base_functions;
942 ++diff_base_functions;
943 }
944
945 // Increment quantities for integration pts.
946 ++rho;
947 ++grad_rho;
948 ++drho_dt;
949 ++R;
950 }
951
952 // Assemble internal force vector for time solver (TS)
954 getFEMethod()->ts_F, /// Get the right hand side vector for time solver
955 nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
956
958}
959
962 "RHO", "RHO", UserDataOperator::OPROWCOL),
963 commonData(common_data) {
964 sYmm = true;
965}
966
968OpAssmbleRhoLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
969 EntityType col_type,
970 DataForcesAndSourcesCore::EntData &row_data,
971 DataForcesAndSourcesCore::EntData &col_data) {
972
974
975 const int row_nb_dofs = row_data.getIndices().size();
976 if (!row_nb_dofs)
978 const int col_nb_dofs = col_data.getIndices().size();
979 if (!col_nb_dofs)
981
982 // Set size can clear local tangent matrix
983 locK_rho_rho.resize(row_nb_dofs, col_nb_dofs, false);
984 locK_rho_rho.clear();
985
986 const int row_nb_gauss_pts = row_data.getN().size1();
987 if (!row_nb_gauss_pts)
989 const int row_nb_base_functions = row_data.getN().size2();
990
991 const double ts_a = getFEMethod()->ts_a;
992 const double R0 = commonData.R0;
993 const double c = commonData.c;
994 const double n = commonData.n;
995 const double m = commonData.m;
996 const double rho_ref = commonData.rHo_ref;
997 const double psi_ref = commonData.pSi_ref;
998
999 FTensor::Index<'I', 3> I;
1000
1003 auto row_base_functions = row_data.getFTensor0N();
1004 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1005
1006 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1007
1008 // Get volume and integration weight
1009 double w = getVolume() * getGaussPts()(3, gg);
1010
1011 const double kp = commonData.getCFromDensity(rho);
1012 const double kp_diff = commonData.getCFromDensityDiff(rho);
1013 const double a = c * kp * ((n - m) / rho) * pow(rho / rho_ref, n - m) * psi;
1014 const double a_diff =
1015 c * kp_diff * (pow(rho / rho_ref, n - m) * psi - psi_ref);
1016
1017 int row_bb = 0;
1018 for (; row_bb != row_nb_dofs; row_bb++) {
1019 auto col_base_functions = col_data.getFTensor0N(gg, 0);
1020 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1021 for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1022 locK_rho_rho(row_bb, col_bb) +=
1023 ts_a * w * row_base_functions * col_base_functions;
1024 locK_rho_rho(row_bb, col_bb) -=
1025 w * R0 * row_diff_base_functions(I) * col_diff_base_functions(I);
1026 locK_rho_rho(row_bb, col_bb) -=
1027 a * w * row_base_functions * col_base_functions;
1028 locK_rho_rho(row_bb, col_bb) -=
1029 a_diff * w * row_base_functions * col_base_functions;
1030
1031 ++col_base_functions;
1032 ++col_diff_base_functions;
1033 }
1034 ++row_base_functions;
1035 ++row_diff_base_functions;
1036 }
1037 for (; row_bb != row_nb_base_functions; row_bb++) {
1038 ++row_base_functions;
1039 ++row_diff_base_functions;
1040 }
1041 ++psi;
1042 ++rho;
1043 }
1044
1045 // cerr << locK_rho_rho << endl;
1046
1047 CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1048 &*row_data.getIndices().begin(), col_nb_dofs,
1049 &*col_data.getIndices().begin(), &locK_rho_rho(0, 0),
1050 ADD_VALUES);
1051
1052 // is symmetric
1053 if (row_side != col_side || row_type != col_type) {
1054 transLocK_rho_rho.resize(col_nb_dofs, row_nb_dofs, false);
1055 noalias(transLocK_rho_rho) = trans(locK_rho_rho);
1056 CHKERR MatSetValues(getFEMethod()->ts_B, col_nb_dofs,
1057 &*col_data.getIndices().begin(), row_nb_dofs,
1058 &*row_data.getIndices().begin(),
1059 &transLocK_rho_rho(0, 0), ADD_VALUES);
1060 }
1061
1063}
1064
1067 "RHO", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
1068 commonData(common_data) {
1069 sYmm = false;
1070}
1071
1073OpAssmbleRhoLhs_dF::doWork(int row_side, int col_side, EntityType row_type,
1074 EntityType col_type,
1075 DataForcesAndSourcesCore::EntData &row_data,
1076 DataForcesAndSourcesCore::EntData &col_data) {
1077
1079
1080 const int row_nb_dofs = row_data.getIndices().size();
1081 if (!row_nb_dofs)
1083 const int col_nb_dofs = col_data.getIndices().size();
1084 if (!col_nb_dofs)
1086
1087 // Set size can clear local tangent matrix
1088 locK_rho_F.resize(row_nb_dofs, col_nb_dofs, false);
1089 locK_rho_F.clear();
1090
1091 const int row_nb_gauss_pts = row_data.getN().size1();
1092 const int col_nb_base_functions = col_data.getN().size2();
1093
1094 const double c = commonData.c;
1095 const double n = commonData.n;
1096 const double m = commonData.m;
1097 const double rho_ref = commonData.rHo_ref;
1098
1099 FTensor::Index<'I', 3> I;
1100 FTensor::Index<'i', 3> i;
1101
1104
1106 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>();
1107
1108 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1109
1110 // Get volume and integration weight
1111 double w = getVolume() * getGaussPts()(3, gg);
1112 const double kp = commonData.getCFromDensity(rho);
1113 const double a = w * c * kp * pow(rho / rho_ref, n - m);
1114
1115 int col_bb = 0;
1116 for (; col_bb != col_nb_dofs / 3; col_bb++) {
1117 t1(i) = a * P(i, I) * col_diff_base_functions(I);
1118 auto row_base_functions = row_data.getFTensor0N(gg, 0);
1119 for (int row_bb = 0; row_bb != row_nb_dofs; row_bb++) {
1120 FTensor::Tensor1<double *, 3> k(&locK_rho_F(row_bb, 3 * col_bb + 0),
1121 &locK_rho_F(row_bb, 3 * col_bb + 1),
1122 &locK_rho_F(row_bb, 3 * col_bb + 2));
1123 k(i) -= row_base_functions * t1(i);
1124 ++row_base_functions;
1125 }
1126 ++col_diff_base_functions;
1127 }
1128 for (; col_bb != col_nb_base_functions; col_bb++) {
1129 ++col_diff_base_functions;
1130 }
1131
1132 ++P;
1133 ++rho;
1134 }
1135
1136 CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1137 &*row_data.getIndices().begin(), col_nb_dofs,
1138 &*col_data.getIndices().begin(), &locK_rho_F(0, 0),
1139 ADD_VALUES);
1140
1142}
1143template <bool ADOLC>
1145 Remodeling::CommonData &common_data)
1147 "DISPLACEMENTS", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
1148 commonData(common_data) {
1149 sYmm = true;
1150}
1151template <bool ADOLC>
1153 int row_side, int col_side, EntityType row_type, EntityType col_type,
1154 DataForcesAndSourcesCore::EntData &row_data,
1155 DataForcesAndSourcesCore::EntData &col_data) {
1156
1158
1159 const int row_nb_dofs = row_data.getIndices().size();
1160 if (!row_nb_dofs)
1162 const int col_nb_dofs = col_data.getIndices().size();
1163 if (!col_nb_dofs)
1165
1166 const bool diagonal_block = (row_type == col_type) && (row_side == col_side);
1167
1168 // Set size can clear local tangent matrix
1169 locK_P_F.resize(row_nb_dofs, col_nb_dofs, false);
1170 locK_P_F.clear();
1171
1172 const int row_nb_gauss_pts = row_data.getN().size1();
1173 const int row_nb_base_functions = row_data.getN().size2();
1174
1175 MatrixDouble &matS = commonData.data.matS;
1176 MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
1177 // MatrixDouble &dS_dC = commonData.data.matMaterialTangent;
1178 // MatrixDouble &matF = commonData.data.matGradientOfDeformation;
1179
1180 double t0;
1181
1182 const double n = commonData.n;
1183 const double rho_ref = commonData.rHo_ref;
1184
1185 FTensor::Index<'i', 3> i;
1186 FTensor::Index<'j', 3> j;
1187 FTensor::Index<'k', 3> k;
1188 FTensor::Index<'I', 3> I;
1189 FTensor::Index<'J', 3> J;
1190 FTensor::Index<'K', 3> K;
1191 FTensor::Index<'L', 3> L;
1192
1193 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1196 auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
1197
1198 // FTensor::Tensor4_ddg<double*,3,3> D1(SYMMETRIC_TENSOR4_MAT_PTR(dS_dC),6);
1199 // FTensor::Tensor2<double*,3,3> F(TENSOR2_MAT_PTR(matF));
1200
1201 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1202
1203 // Get volume and integration weight
1204 double w = getVolume() * getGaussPts()(3, gg);
1205 const double a = w * pow(rho / rho_ref, n);
1206
1207 int row_bb = 0;
1208 for (; row_bb != row_nb_dofs / 3; row_bb++) {
1209
1210 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
1211
1212 const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
1213 int col_bb = 0;
1214 for (; col_bb != final_bb; col_bb++) {
1215
1217 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1218 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1219 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1220 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1221 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1222 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1223 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1224 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1225 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1226
1227 diffDiff(J, L) =
1228 a * row_diff_base_functions(J) * col_diff_base_functions(L);
1229
1230 // TODO: FIXME: Tensor D2 has major symmetry, we do not have such tensor
1231 // implemented at the moment. Using tensor with major symmetry will
1232 // improve code efficiency.
1233 t_assemble(i, k) += diffDiff(J, L) * D2(i, J, k, L);
1234
1235 if (ADOLC) {
1236 // Stress stiffening part (only with adolc since it is using dS / dC
1237 // derrivative. Check OpCalculateStressTangentWithAdolc operator)
1238 t0 = diffDiff(J, L) * S(J, L);
1239 t_assemble(0, 0) += t0;
1240 t_assemble(1, 1) += t0;
1241 t_assemble(2, 2) += t0;
1242 }
1243
1244 // Next base function for column
1245 ++col_diff_base_functions;
1246 }
1247
1248 ++row_diff_base_functions;
1249 }
1250 for (; row_bb != row_nb_base_functions; row_bb++) {
1251 ++row_diff_base_functions;
1252 }
1253
1254 // ++D1;
1255 ++D2;
1256 ++S;
1257 // ++F;
1258 ++rho;
1259 }
1260
1261 // Copy of symmetric part for assemble
1262 if (diagonal_block) {
1263 for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
1264 int col_bb = 0;
1265 for (; col_bb != row_bb + 1; col_bb++) {
1267 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
1268 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
1269 &locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
1270 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
1271 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
1272 &locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
1273 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
1274 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
1275 &locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
1277 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 0),
1278 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 1),
1279 &locK_P_F(3 * col_bb + 0, 3 * row_bb + 2),
1280 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 0),
1281 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 1),
1282 &locK_P_F(3 * col_bb + 1, 3 * row_bb + 2),
1283 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 0),
1284 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 1),
1285 &locK_P_F(3 * col_bb + 2, 3 * row_bb + 2));
1286 t_off_side(i, k) = t_assemble(k, i);
1287 }
1288 }
1289 }
1290
1291 const int *row_ind = &*row_data.getIndices().begin();
1292 const int *col_ind = &*col_data.getIndices().begin();
1293 CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs, row_ind, col_nb_dofs,
1294 col_ind, &locK_P_F(0, 0), ADD_VALUES);
1295
1296 if (row_type != col_type || row_side != col_side) {
1297 transLocK_P_F.resize(col_nb_dofs, row_nb_dofs, false);
1298 noalias(transLocK_P_F) = trans(locK_P_F);
1299 CHKERR MatSetValues(getFEMethod()->ts_B, col_nb_dofs, col_ind, row_nb_dofs,
1300 row_ind, &transLocK_P_F(0, 0), ADD_VALUES);
1301 }
1302
1304}
1305
1307 Remodeling::CommonData &common_data)
1309 "DISPLACEMENTS", "RHO", UserDataOperator::OPROWCOL),
1310 commonData(common_data) {
1311 sYmm = false; // This is off-diagonal (upper-left), so no symmetric
1312}
1313
1315OpAssmbleStressLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
1316 EntityType col_type,
1317 DataForcesAndSourcesCore::EntData &row_data,
1318 DataForcesAndSourcesCore::EntData &col_data) {
1319
1321
1322 const int row_nb_dofs = row_data.getIndices().size();
1323 if (!row_nb_dofs)
1325 const int col_nb_dofs = col_data.getIndices().size();
1326 if (!col_nb_dofs)
1328
1329 // Set size can clear local tangent matrix
1330 locK_P_RHO.resize(row_nb_dofs, col_nb_dofs, false);
1331 locK_P_RHO.clear();
1332
1334
1335 const double n = commonData.n;
1336 const double rho_ref = commonData.rHo_ref;
1337
1338 FTensor::Index<'I', 3> I;
1339 FTensor::Index<'i', 3> i;
1340
1343
1344 const int row_nb_gauss_pts = row_data.getN().size1();
1345 const int row_nb_base_functions = row_data.getN().size2();
1346 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
1347
1348 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
1349
1350 // Get volume and integration weight
1351 double w = getVolume() * getGaussPts()(3, gg);
1352 double a = w * (n / rho) * pow(rho / rho_ref, n);
1353
1354 int row_bb = 0;
1355 for (; row_bb != row_nb_dofs / 3; row_bb++) {
1356 t1(i) = a * row_diff_base_functions(I) * P(i, I);
1357 auto base_functions = col_data.getFTensor0N(gg, 0);
1358 for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
1359 // if(base_functions!=col_data.getN(gg)[col_bb]) {
1360 // cerr << "Error" << endl;
1361 // }
1362 FTensor::Tensor1<double *, 3> k(&locK_P_RHO(3 * row_bb + 0, col_bb),
1363 &locK_P_RHO(3 * row_bb + 1, col_bb),
1364 &locK_P_RHO(3 * row_bb + 2, col_bb));
1365 k(i) += t1(i) * base_functions;
1366 ++base_functions;
1367 }
1368 ++row_diff_base_functions;
1369 }
1370 for (; row_bb != row_nb_base_functions; row_bb++) {
1371 ++row_diff_base_functions;
1372 }
1373
1374 ++P;
1375 ++rho;
1376 }
1377
1378 CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
1379 &*row_data.getIndices().begin(), col_nb_dofs,
1380 &*col_data.getIndices().begin(), &locK_P_RHO(0, 0),
1381 ADD_VALUES);
1382
1384}
1385
1387 Remodeling::CommonData &common_data)
1388 : mField(m_field), commonData(common_data), postProc(m_field),
1389 postProcElastic(m_field),
1390 // densityMaps(m_field),
1391 iNit(false) {}
1392
1397
1402
1405
1406 PetscBool mass_postproc = PETSC_FALSE;
1407 PetscBool equilibrium_flg = PETSC_FALSE;
1408 PetscBool analysis_mesh_flg = PETSC_FALSE;
1409 rate = 1;
1410
1411 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling post-process",
1412 "none");
1413
1414 CHKERR PetscOptionsBool(
1415 "-mass_postproc",
1416 "is used for testing, calculates mass and energy at each time step", "",
1417 mass_postproc, &mass_postproc, PETSC_NULL);
1418
1419 CHKERR PetscOptionsBool("-analysis_mesh",
1420 "saves analysis mesh at each time step", "",
1421 analysis_mesh_flg, &analysis_mesh_flg, PETSC_NULL);
1422
1423 CHKERR PetscOptionsScalar(
1424 "-equilibrium_stop_rate",
1425 "is used to stop calculations when equilibium state is achieved", "",
1427 PetscOptionsEnd();
1428
1429 if (!iNit) {
1430
1431 PetscOptionsBegin(PETSC_COMM_WORLD, "",
1432 "Bone remodeling post-process", "none");
1433
1434 pRT = 1;
1435 CHKERR PetscOptionsInt("-my_output_prt",
1436 "frequncy how often results are dumped on hard disk",
1437 "", 1, &pRT, NULL);
1438 PetscOptionsEnd();
1439
1441 CHKERR postProc.addFieldValuesPostProc("DISPLACEMENTS");
1442 CHKERR postProc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
1445 }
1446 // CHKERR postProc.addFieldValuesPostProc("RHO");
1447 // CHKERR postProc.addFieldValuesGradientPostProc("RHO");
1448
1449 {
1450 auto &list_of_operators = postProc.getOpPtrVector();
1451
1452 list_of_operators.push_back(
1454 list_of_operators.push_back(new OpCalculateScalarFieldGradient<3>(
1455 "RHO", commonData.data.gradRhoPtr));
1456 // Get displacement gradient at integration points
1457 list_of_operators.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1458 "DISPLACEMENTS", commonData.data.gradDispPtr));
1459 list_of_operators.push_back(new OpCalculateStress(commonData));
1460 list_of_operators.push_back(new OpPostProcStress(
1462 }
1463
1466 CHKERR postProcElastic.addFieldValuesPostProc("MESH_NODE_POSITIONS");
1468 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
1469 commonData.elasticPtr->setOfBlocks.begin();
1470 for (; sit != commonData.elasticPtr->setOfBlocks.end(); sit++) {
1473 "DISPLACEMENTS", sit->second, postProcElastic.commonData));
1474 }
1475
1476 iNit = true;
1477 }
1478
1479 if (mass_postproc) {
1480 Vec mass_vec;
1481 Vec energ_vec;
1482 int mass_vec_ghost[] = {0};
1483 CHKERR VecCreateGhost(PETSC_COMM_WORLD, (!mField.get_comm_rank()) ? 1 : 0,
1484 1, 1, mass_vec_ghost, &mass_vec);
1485
1486 CHKERR VecZeroEntries(mass_vec);
1487 CHKERR VecDuplicate(mass_vec, &energ_vec);
1488
1489 Remodeling::Fe energyCalc(mField);
1490 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", energyCalc, true, false, false,
1491 true);
1492 energyCalc.getOpPtrVector().push_back(
1494 energyCalc.getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<3>(
1495 "RHO", commonData.data.gradRhoPtr));
1496 energyCalc.getOpPtrVector().push_back(
1497 new OpCalculateVectorFieldGradient<3, 3>("DISPLACEMENTS",
1499 energyCalc.getOpPtrVector().push_back(new OpCalculateStress(commonData));
1500 energyCalc.getOpPtrVector().push_back(
1501 new OpMassAndEnergyCalculation("RHO", commonData, energ_vec, mass_vec));
1502 CHKERR DMoFEMLoopFiniteElements(commonData.dm, "FE_REMODELLING",
1503 &energyCalc);
1504
1505 CHKERR VecAssemblyBegin(mass_vec);
1506 CHKERR VecAssemblyEnd(mass_vec);
1507 double mass_sum;
1508 CHKERR VecSum(mass_vec, &mass_sum);
1509
1510 CHKERR VecAssemblyBegin(energ_vec);
1511 CHKERR VecAssemblyEnd(energ_vec);
1512 double energ_sum;
1513 CHKERR VecSum(energ_vec, &energ_sum);
1514
1515 CHKERR PetscPrintf(PETSC_COMM_WORLD,
1516 "Time: %g Mass: %6.5g Elastic energy: %6.5g \n", ts_t,
1517 mass_sum, energ_sum);
1518 CHKERR VecDestroy(&mass_vec);
1519 CHKERR VecDestroy(&energ_vec);
1520
1521 if (equilibrium_flg) {
1522 double equilibrium_rate =
1523 fabs(energ_sum - commonData.cUrrent_psi) / energ_sum;
1524 double equilibrium_mass_rate =
1525 fabs(mass_sum - commonData.cUrrent_mass) / mass_sum;
1526 if (equilibrium_rate < rate) {
1527 CHKERR PetscPrintf(
1528 PETSC_COMM_WORLD,
1529 "Energy equilibrium state is achieved! Difference = %0.6g %%. \n",
1530 equilibrium_rate * 100);
1531 }
1532 if (equilibrium_mass_rate < rate) {
1533 CHKERR PetscPrintf(
1534 PETSC_COMM_WORLD,
1535 "Mass equilibrium state is achieved! Difference = %0.6g %%. \n",
1536 equilibrium_mass_rate * 100);
1537 }
1538 commonData.cUrrent_psi = energ_sum;
1539 commonData.cUrrent_mass = mass_sum;
1540 }
1541 }
1542
1543 int step;
1544#if PETSC_VERSION_LE(3, 8, 0)
1545 CHKERR TSGetTimeStepNumber(ts, &step);
1546#else
1547 CHKERR TSGetStepNumber(ts, &step);
1548#endif
1549
1550 if ((step) % pRT == 0) {
1551
1552 ostringstream sss;
1553 sss << "out_" << step << ".h5m";
1555 CHKERR postProc.postProcMesh.write_file(sss.str().c_str(), "MOAB",
1556 "PARALLEL=WRITE_PART");
1557 if (analysis_mesh_flg) {
1558 ostringstream ttt;
1559 ttt << "analysis_mesh_" << step << ".h5m";
1560 CHKERR mField.get_moab().write_file(ttt.str().c_str(), "MOAB",
1561 "PARALLEL=WRITE_PART");
1562 }
1563
1567 ostringstream ss;
1568 ss << "out_elastic_" << step << ".h5m";
1569 CHKERR postProcElastic.postProcMesh.write_file(ss.str().c_str(), "MOAB",
1570 "PARALLEL=WRITE_PART");
1571 }
1572 }
1574}
1575
1577
1579 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling", "none");
1580
1581 commonData.oRder = 2;
1582 CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
1583 &commonData.oRder, PETSC_NULL);
1584
1585 PetscOptionsEnd();
1586
1587 CHKERR mField.get_moab().get_entities_by_type(0, MBTET, commonData.tEts_all);
1589 string name = it->getName();
1590 if (name.compare(0, 14, "NO_REMODELLING") == 0) {
1592 EntityHandle meshset = it->getMeshset();
1593 CHKERR this->mField.get_moab().get_entities_by_type(
1594 meshset, MBTET, commonData.tEts_block, true);
1597 }
1598 }
1599
1601}
1602
1604
1606
1607 // Seed all mesh entities to MoFEM database, those entities can be potentially
1608 // used as finite elements or as entities which carry some approximation
1609 // field.
1610 commonData.bitLevel.set(0);
1611 CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
1612 0, 3, commonData.bitLevel);
1613 int order = commonData.oRder;
1614
1615 // Add displacement field
1616 CHKERR mField.add_field("DISPLACEMENTS", H1,
1617 /*AINSWORTH_LOBATTO_BASE*/ AINSWORTH_LEGENDRE_BASE, 3,
1618 MB_TAG_SPARSE, MF_ZERO);
1619 // Add field representing ho-geometry
1620 CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3,
1621 MB_TAG_SPARSE, MF_ZERO);
1622
1623 // Check if density is available, if not add density field.
1624 bool add_rho_field = false;
1625 if (!mField.check_field("RHO")) {
1627 MB_TAG_SPARSE, MF_ZERO);
1628 // FIXME
1630 // CHKERR mField.add_ents_to_field_by_type(0,MBTET,"RHO");
1632 add_rho_field = true;
1633
1634 CHKERR mField.set_field_order(0, MBVERTEX, "RHO", 1);
1635 CHKERR mField.set_field_order(0, MBEDGE, "RHO", order - 1);
1636 CHKERR mField.set_field_order(0, MBTRI, "RHO", order - 1);
1637 CHKERR mField.set_field_order(0, MBTET, "RHO", order - 1);
1638 }
1639
1640 // Add entities to field
1641 CHKERR mField.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENTS");
1642 CHKERR mField.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
1643
1644 // Set approximation order to entities
1645 CHKERR mField.set_field_order(0, MBVERTEX, "DISPLACEMENTS", 1);
1646 CHKERR mField.set_field_order(0, MBEDGE, "DISPLACEMENTS", order);
1647 CHKERR mField.set_field_order(0, MBTRI, "DISPLACEMENTS", order);
1648 CHKERR mField.set_field_order(0, MBTET, "DISPLACEMENTS", order);
1649
1650 // Assumes that geometry is approximated using 2nd order polynomials.
1651
1652 // Apply 2nd order only on skin
1653 {
1654 // Skinner skin(&mField.get_moab());
1655 // Range faces,tets;
1656 // CHKERR mField.get_moab().get_entities_by_type(0,MBTET,tets);
1657 // CHKERR skin.find_skin(0,tets,false,faces);
1658 // Range edges;
1659 // CHKERR mField.get_moab().get_adjacencies(
1660 // faces,1,false,edges,moab::Interface::UNION
1661 // );
1662 CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
1663 }
1664
1665 CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
1666
1667 // Build fields
1669
1670 // If order was not set from CT scans set homogenous order equal to
1671 // reference bone density
1672 if (add_rho_field) {
1674 MBVERTEX, "RHO");
1675 // const DofEntity_multiIndex *dofs_ptr;
1676 // CHKERR mField.get_dofs(&dofs_ptr);
1677 // for(_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(mField,"RHO",dit)) {
1678 // cerr << (*dit)->getFieldData() << endl;
1679 // }
1680 }
1681
1682 // Project geometry given on 10-node tets on ho-geometry
1683 Projection10NodeCoordsOnField ent_method_material(mField,
1684 "MESH_NODE_POSITIONS");
1685 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
1686
1688}
1689
1692
1693 CHKERR mField.add_finite_element("FE_REMODELLING", MF_ZERO);
1695 "DISPLACEMENTS");
1697 "DISPLACEMENTS");
1699 "DISPLACEMENTS");
1700 CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING", "RHO");
1702 "MESH_NODE_POSITIONS");
1703 CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING", "RHO");
1704 CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING", "RHO");
1706 "FE_REMODELLING");
1707
1708 CHKERR mField.add_finite_element("ELASTIC");
1709 CHKERR mField.modify_finite_element_add_field_row("ELASTIC", "DISPLACEMENTS");
1710 CHKERR mField.modify_finite_element_add_field_col("ELASTIC", "DISPLACEMENTS");
1712 "DISPLACEMENTS");
1714 "MESH_NODE_POSITIONS");
1715
1718 MBTET, "ELASTIC");
1719 }
1720
1722 boost::shared_ptr<ElasticMaterials>(new ElasticMaterials(mField));
1723 commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
1726 commonData.elasticPtr->setOfBlocks);
1727 CHKERR commonData.elasticPtr->addElement("ELASTIC", "DISPLACEMENTS");
1728 CHKERR commonData.elasticPtr->setOperators(
1729 "DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
1730
1733
1734 // Allocate memory for density and gradient of displacements at integration
1735 // points
1736 commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
1738 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1740 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1741
1742 // Add operators Rhs
1743 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feRhs), true, false,
1744 false, false);
1745 auto &list_of_operators_rhs = commonData.feRhs->getOpPtrVector();
1746 // Get density at integration points
1747 list_of_operators_rhs.push_back(
1749 list_of_operators_rhs.push_back(
1751 list_of_operators_rhs.push_back(new OpGetRhoTimeDirevative(commonData));
1752 // Get displacement gradient at integration points
1753 list_of_operators_rhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1754 "DISPLACEMENTS", commonData.data.gradDispPtr));
1755 list_of_operators_rhs.push_back(new OpCalculateStress(commonData));
1756 list_of_operators_rhs.push_back(new OpAssmbleStressRhs(commonData));
1757 list_of_operators_rhs.push_back(new OpAssmbleRhoRhs(commonData));
1758
1759 // Add operators Lhs
1760 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feLhs), true, false,
1761 false, false);
1762 auto &list_of_operators_lhs = commonData.feLhs->getOpPtrVector();
1763 // Get density at integration points
1764 list_of_operators_lhs.push_back(
1766 list_of_operators_rhs.push_back(
1768 // Get displacement gradient at integration points
1769 list_of_operators_lhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1770 "DISPLACEMENTS", commonData.data.gradDispPtr));
1771 if (commonData.with_adol_c) {
1772 list_of_operators_lhs.push_back(
1774 list_of_operators_lhs.push_back(
1776 } else {
1777 list_of_operators_lhs.push_back(new OpCalculateStressTangent(commonData));
1778 list_of_operators_lhs.push_back(
1780 }
1781 list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dRho(commonData));
1782 list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dF(commonData));
1783 list_of_operators_lhs.push_back(new OpAssmbleStressLhs_dRho(commonData));
1784
1786}
1787
1789
1791 // Add Neumann forces elements
1792 CHKERR MetaNeummanForces::addNeumannBCElements(mField, "DISPLACEMENTS");
1793 CHKERR MetaNodalForces::addElement(mField, "DISPLACEMENTS");
1794 CHKERR MetaEdgeForces::addElement(mField, "DISPLACEMENTS");
1795
1796 // forces and pressures on surface
1797 CHKERR MetaNeummanForces::setMomentumFluxOperators(
1798 mField, commonData.neumannForces, PETSC_NULL, "DISPLACEMENTS");
1799 // noadl forces
1801 PETSC_NULL, "DISPLACEMENTS");
1802 // edge forces
1804 "DISPLACEMENTS");
1805
1806 for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
1807 commonData.neumannForces.begin();
1808 mit != commonData.neumannForces.end(); mit++) {
1809 mit->second->methodsOp.push_back(
1810 new TimeForceScale("-my_load_history", false));
1811 }
1812 for (boost::ptr_map<string, NodalForce>::iterator mit =
1813 commonData.nodalForces.begin();
1814 mit != commonData.nodalForces.end(); mit++) {
1815 mit->second->methodsOp.push_back(
1816 new TimeForceScale("-my_load_history", false));
1817 }
1818 for (boost::ptr_map<string, EdgeForce>::iterator mit =
1819 commonData.edgeForces.begin();
1820 mit != commonData.edgeForces.end(); mit++) {
1821 mit->second->methodsOp.push_back(
1822 new TimeForceScale("-my_load_history", false));
1823 }
1824
1826}
1827
1829
1831 commonData.dm_name = "DM_REMODELING";
1832 // Register DM problem
1834 CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
1837 // Create DM instance
1840 // Configure DM form line command options (DM itself, solvers,
1841 // pre-conditioners, ... )
1842 CHKERR DMSetFromOptions(commonData.dm);
1843 // Add elements to dm (only one here)
1844 CHKERR DMMoFEMAddElement(commonData.dm, "FE_REMODELLING");
1846
1847 if (mField.check_finite_element("FORCE_FE")) {
1849 }
1850 if (mField.check_finite_element("PRESSURE_FE")) {
1851 CHKERR DMMoFEMAddElement(commonData.dm, "PRESSURE_FE");
1852 }
1853 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
1854 CHKERR DMSetUp(commonData.dm);
1855
1857}
1858
1860
1862
1863 Vec F = commonData.F;
1864 Vec D = commonData.D;
1865 Mat A = commonData.A;
1866 DM dm = commonData.dm;
1867 TS ts = commonData.ts;
1868
1869 CHKERR TSSetIFunction(ts, F, PETSC_NULL, PETSC_NULL);
1870 CHKERR TSSetIJacobian(ts, A, A, PETSC_NULL, PETSC_NULL);
1871 double ftime = 1;
1872#if PETSC_VERSION_GE(3, 8, 0)
1873 CHKERR TSSetMaxTime(ts, ftime);
1874#endif
1875 CHKERR TSSetFromOptions(ts);
1876 CHKERR TSSetDM(ts, dm);
1877
1878 SNES snes;
1879 CHKERR TSGetSNES(ts, &snes);
1880 KSP ksp;
1881 CHKERR SNESGetKSP(snes, &ksp);
1882 PC pc;
1883 CHKERR KSPGetPC(ksp, &pc);
1884 PetscBool is_pcfs = PETSC_FALSE;
1885 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1886
1887 // Set up FIELDSPLIT
1888 // Only is user set -pc_type fieldsplit
1889 if (is_pcfs == PETSC_TRUE) {
1890 IS is_disp, is_rho;
1891 const MoFEM::Problem *problem_ptr;
1892 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
1893 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1894 problem_ptr->getName(), ROW, "DISPLACEMENTS", 0, 3, &is_disp);
1895 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1896 problem_ptr->getName(), ROW, "RHO", 0, 1, &is_rho);
1897 CHKERR ISSort(is_disp);
1898 CHKERR ISSort(is_rho);
1899 CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
1900 CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
1901 CHKERR ISDestroy(&is_disp);
1902 CHKERR ISDestroy(&is_rho);
1903 }
1904
1905 // Monitor
1906 MonitorPostProc post_proc(mField, commonData);
1907 TsCtx *ts_ctx;
1908 DMMoFEMGetTsCtx(dm, &ts_ctx);
1909 {
1910 ts_ctx->getPostProcessMonitor().push_back(&post_proc);
1911 CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
1912 }
1913
1914#if PETSC_VERSION_GE(3, 7, 0)
1915 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1916#endif
1917 CHKERR TSSolve(ts, D);
1918 CHKERR TSGetTime(ts, &ftime);
1919 PetscInt steps, snesfails, rejects, nonlinits, linits;
1920#if PETSC_VERSION_LE(3, 8, 0)
1921 CHKERR TSGetTimeStepNumber(ts, &steps);
1922#else
1923 CHKERR TSGetStepNumber(ts, &steps);
1924#endif
1925
1926 CHKERR TSGetSNESFailures(ts, &snesfails);
1927 CHKERR TSGetStepRejections(ts, &rejects);
1928 CHKERR TSGetSNESIterations(ts, &nonlinits);
1929 CHKERR TSGetKSPIterations(ts, &linits);
1930 PetscPrintf(PETSC_COMM_WORLD,
1931 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
1932 "linits %D\n",
1933 steps, rejects, snesfails, ftime, nonlinits, linits);
1934
1936 if (commonData.cUrrent_psi < 1.67 || commonData.cUrrent_psi > 1.68)
1937 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "atom test diverged");
1938 }
1939
1941}
1942
1944 int row_side, EntityType row_type,
1945 DataForcesAndSourcesCore::EntData &row_data) {
1947
1948 if (row_type != MBVERTEX)
1950
1951 const int nb_gauss_pts = row_data.getN().size1();
1952 // commonData.data.vecPsi.resize(nb_gauss_pts,false);
1955
1956 double energy = 0;
1957 double mass = 0;
1958 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1959
1960 double vol = getVolume() * getGaussPts()(3, gg);
1961 energy += vol * psi;
1962 mass += vol * rho;
1963
1964 ++psi;
1965 ++rho;
1966 }
1967
1968 CHKERR VecSetValue(energVec, 0, energy, ADD_VALUES);
1969 CHKERR VecSetValue(massVec, 0, mass, ADD_VALUES);
1971}
1972
1973} // namespace BoneRemodeling
Elastic materials.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)
@ ROW
@ MF_ZERO
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
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
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT)
@ 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 ...
#define TENSOR4_MAT_PTR(MAT)
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
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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
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 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.
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
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
MoFEMErrorCode recordFreeEnergy_dC(Remodeling::CommonData &common_data, T &dC, B1 &psi)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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: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.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
constexpr IntegrationType I
double young_modulus
Young modulus.
Definition plastic.cpp:125
double rho
Definition plastic.cpp:144
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:126
FTensor::Index< 'm', 3 > m
Remodeling::CommonData & commonData
PostProcVolumeOnRefinedMesh postProcElastic
PostProcVolumeOnRefinedMesh postProc
MonitorPostProc(MoFEM::Interface &m_field, Remodeling::CommonData &common_data)
Off diagonal block of tangent matrix .
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
Diagonal block of tangent matrix .
OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
Remodeling::CommonData & commonData
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble residual for conservation of mass (density)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
VectorDouble nF
Vector of the right hand side (internal forces)
Remodeling::CommonData & commonData
Off diagonal block of tangent matrix .
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}...
Remodeling::CommonData & commonData
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)
Remodeling::CommonData & commonData
VectorDouble nF
Vector of the right hand side (internal forces)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpAssmbleStressRhs(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpCalculateStressTangentWithAdolc(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpCalculateStressTangent(Remodeling::CommonData &common_data)
Evaluate physical equations at integration points.
Remodeling::CommonData & commonData
OpCalculateStress(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Evaluate density derivative with respect to time in case of Backward Euler Method.
OpGetRhoTimeDirevative(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Remodeling::CommonData & commonData
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Used to post proc stresses, energy, source term.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, Remodeling::CommonData &common_data)
Remodeling::CommonData & commonData
std::vector< EntityHandle > & mapGaussPts
boost::shared_ptr< MatrixDouble > gradDispPtr
Ptr to gradient of displacements matrix container.
boost::shared_ptr< VectorDouble > rhoPtr
Ptr to density matrix container.
boost::shared_ptr< MatrixDouble > gradRhoPtr
Gradient of density field.
MatrixDouble matInvF
inverse of deformation gradient
VectorDouble vecRhoDt
Time derivative of density.
MatrixDouble matGradientOfDeformation
Gradient of deformation.
double getCFromDensity(const double &rho)
double pSi_ref
reference free energy
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
double n
porosity exponent [-]
int b
b exponent for bell function
double getCFromDensityDiff(const double &rho)
double m
algorithmic exponent [-]
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
PetscBool is_atom_testing
for atom tests
double c
density evolution (growth) velocity [d/m^2]
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
PetscBool less_post_proc
reduce file size
boost::shared_ptr< ElasticMaterials > elasticMaterialsPtr
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
double R0
mass conduction coefficient
double cUrrent_psi
current free energy for evaluating equilibrium state
boost::shared_ptr< NonlinearElasticElement > elasticPtr
double cUrrent_mass
current free energy for evaluating equilibrium state
Volume finite element.
MoFEM::Interface & mField
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.
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 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 field rank 0, i.e. vector field.
Get value at integration points for scalar field.
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.
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:148
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
structure grouping operators and data used for calculation of nonlinear elastic element
std::vector< EntityHandle > mapGaussPts
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Force scale operator for reading two columns.