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

Go to the source code of this file.

Namespaces

namespace  BoneRemodeling
 

Macros

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

Macro Definition Documentation

◆ SYMMETRIC_TENSOR2_MAT_PTR

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

Definition at line 54 of file Remodeling.cpp.

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

◆ SYMMETRIC_TENSOR2_VEC_PTR

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

Definition at line 57 of file Remodeling.cpp.

◆ SYMMETRIC_TENSOR4_MAT_PTR

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

Definition at line 40 of file Remodeling.cpp.

◆ TENSOR1_VEC_PTR

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

Definition at line 38 of file Remodeling.cpp.

◆ TENSOR2_MAT_PTR

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

Definition at line 50 of file Remodeling.cpp.

◆ TENSOR4_MAT_PTR

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

Definition at line 48 of file Remodeling.cpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow