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