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