v0.14.0
Loading...
Searching...
No Matches
Header.hpp
Go to the documentation of this file.
1/* This file is part of MoFEM.
2 * MoFEM is free software: you can redistribute it and/or modify it under
3 * the terms of the GNU Lesser General Public License as published by the
4 * Free Software Foundation, either version 3 of the License, or (at your
5 * option) any later version.
6 *
7 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10 * License for more details.
11 *
12 * You should have received a copy of the GNU Lesser General Public
13 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14
15#ifndef __HEADER_HPP__
16#define __HEADER_HPP__
17
18// declaration
19struct CommonData {
20
23 Tag thDens;
24
27
28 boost::shared_ptr<VectorDouble> dataRhotGaussPtr;
29 boost::shared_ptr<MatrixDouble> dispGradAtGaussPts;
30
33 vector<VectorDouble> sTress;
34 double volumeFrac;
39
40 double edgeLength;
41 double lAmbda;
42 double mU;
43 double lAmbda0;
44 double rHo;
45 double yOung;
46 double pOisson;
48 double mOve;
50 double cHange;
51 // for approx
52 double cUrrent_mass;
53 double L_mid; // langange multiplier
55
56 CommonData(MoFEM::Interface &m_field) : mField(m_field) {
57
58 dataRhotGaussPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
59 dispGradAtGaussPts = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
60 volumeFrac = 0.5;
61 yOung = 10;
62 fIlter_radius = 0;
63 penalty = 3;
64 cHange = 0;
65 edgeLength = 1.;
66 lAmbda0 = 0.;
67 cUrrent_mass = 0;
68 L_mid = 0; // langange multiplier
69 uPdate_density = false;
70
71 anisLambdaMat.resize(3,3,false);
72 anisLambdaMat.clear();
73
74 ierr = createTags(m_field.get_moab());
75 }
78 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Topo Config", "none");
79 CHKERR(ierr);
80
81 CHKERR PetscOptionsInt("-penalty", "default penalty ", "", penalty, &penalty,
82 PETSC_NULL);
83
84 CHKERR PetscOptionsScalar("-young", "default yOung modulus", "", yOung, &yOung,
85 PETSC_NULL);
86
87 CHKERR PetscOptionsScalar("-poisson", "poisson ratio", "", 0.1, &pOisson,
88 PETSC_NULL);
89
90 CHKERR PetscOptionsScalar("-move", "max move", "", 0.1, &mOve, PETSC_NULL);
91 // CHKERR PetscOptionsScalar("-lambda", "lambda smooth", "", lAmbda0, &lAmbda0, PETSC_NULL);
92
93 double lambda_coef[3] = {0.,0.,0.};
94 int nmax = 3;
95 PetscBool is_lambda_defined = PETSC_FALSE;
96 PetscOptionsGetRealArray(NULL, "", "-lambda", lambda_coef,
97 &nmax, &is_lambda_defined);
98
99 switch (nmax) {
100
101 case 0:
102 case 1:
103 anisLambdaMat(0, 0) = anisLambdaMat(1, 1) = anisLambdaMat(2, 2) =
104 lambda_coef[0];
105 break;
106 case 3:
107 anisLambdaMat(0, 0) = lambda_coef[0];
108 anisLambdaMat(1, 1) = lambda_coef[1];
109 anisLambdaMat(2, 2) = lambda_coef[2];
110 break;
111 default:
112
113 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
114 "Wrong number of lambda parameters, should be 1 or 3, data "
115 "inconsistency");
116 }
117
118 CHKERR PetscOptionsScalar("-volume_frac", "volume frac", "", volumeFrac,
119 &volumeFrac, PETSC_NULL);
120
121 // CHKERR PetscOptionsScalar("-filter_r", "default filter radius", "", 0,
122 // &fIlter_radius, PETSC_NULL);
123 ierr = PetscOptionsEnd();
124 CHKERRG(ierr);
125
126 lAmbda = (yOung * pOisson) / ((1. + pOisson) * (1. - 2. * pOisson));
127 mU = yOung / (2. * (1. + pOisson));
128
129 const double coefficient = yOung / ((1 + pOisson) * (1 - 2 * pOisson));
130 D.resize(36, 1, false);
131 D.clear();
132
134
139
140 t_D(i, j, k, l) = 0.;
141
142 t_D(0, 0, 0, 0) = 1 - pOisson;
143 t_D(1, 1, 1, 1) = 1 - pOisson;
144 t_D(2, 2, 2, 2) = 1 - pOisson;
145
146 t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * pOisson);
147 t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * pOisson);
148 t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * pOisson);
149
150 t_D(0, 0, 1, 1) = pOisson;
151 t_D(1, 1, 0, 0) = pOisson;
152 t_D(0, 0, 2, 2) = pOisson;
153 t_D(2, 2, 0, 0) = pOisson;
154 t_D(1, 1, 2, 2) = pOisson;
155 t_D(2, 2, 1, 1) = pOisson;
156
157 t_D(i, j, k, l) *= coefficient;
158
160 }
161
162
163 MoFEMErrorCode getSensitivityFromTag(const EntityHandle fe_ent, const int nb_gauss_pts) {
165 double *tag_data;
166 int tag_size;
167 CHKERR mField.get_moab().tag_get_by_ptr(
168 thSensitivity, &fe_ent, 1, (const void **)&tag_data, &tag_size);
169
170 if (tag_size == 1) {
171 sEnsitivityVec.resize(nb_gauss_pts , false);
172 sEnsitivityVec.clear();
173 void const *tag_data[] = {&*sEnsitivityVec.data().begin()};
174 const int tag_size = sEnsitivityVec.data().size();
175 CHKERR mField.get_moab().tag_set_by_ptr(thSensitivity, &fe_ent, 1,
176 tag_data, &tag_size);
177 } else if (tag_size != nb_gauss_pts) {
178 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
179 "Wrong size of the tag, data inconsistency");
180 } else {
182 tag_size, ublas::shallow_array_adaptor<double>(tag_size, tag_data));
183 sEnsitivityVec = tag_vec;
184 }
186 }
187
188
189 MoFEMErrorCode setSensitivityToTag(const EntityHandle fe_ent, const double sens) {
191 void const *tag_data[] = {&*sEnsitivityVec.data().begin()};
192 const int tag_size = sEnsitivityVec.data().size();
193 CHKERR mField.get_moab().tag_set_by_ptr(thSensitivity, &fe_ent, 1,
194 tag_data, &tag_size);
196 }
197
198 MoFEMErrorCode getDensityFromTag(const EntityHandle fe_ent, const int nb_gauss_pts) {
200 double *tag_data;
201 int tag_size;
202 CHKERR mField.get_moab().tag_get_by_ptr(
203 thDens, &fe_ent, 1, (const void **)&tag_data, &tag_size);
204
205 if (tag_size == 1) {
206 rhoVecAtTag.resize(nb_gauss_pts , false);
207 rhoVecAtTag.clear();
208 void const *tag_data[] = {&*rhoVecAtTag.data().begin()};
209 const int tag_size = rhoVecAtTag.data().size();
210 CHKERR mField.get_moab().tag_set_by_ptr(thDens, &fe_ent, 1,
211 tag_data, &tag_size);
212 } else if (tag_size != nb_gauss_pts) {
213 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
214 "Wrong size of the tag, data inconsistency");
215 } else {
217 tag_size, ublas::shallow_array_adaptor<double>(tag_size, tag_data));
218 rhoVecAtTag = tag_vec;
219 }
221 }
222
223 MoFEMErrorCode setDensityToTag(const EntityHandle fe_ent, const double rho) {
225 void const *tag_data[] = {&*rhoVecAtTag.data().begin()};
226 const int tag_size = rhoVecAtTag.data().size();
227 CHKERR mField.get_moab().tag_set_by_ptr(thDens, &fe_ent, 1,
228 tag_data, &tag_size);
230 }
231
232 MoFEMErrorCode createTags(moab::Interface &moab) {
234 std::vector<double> def_val(9, 0);
235 CHKERR mField.get_moab().tag_get_handle(
236 "SENSITIVITY", 1, MB_TYPE_DOUBLE, thSensitivity,
237 MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, &def_val[0]);
238 CHKERR mField.get_moab().tag_get_handle(
239 "RHO", 1, MB_TYPE_DOUBLE, thDens,
240 MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, &def_val[0]);
242 }
243
244 MoFEMErrorCode build_main_dm(DM &dm, DM &sub_dm_elastic, DM &sub_dm_rho,
245 MoFEM::Interface &m_field, BitRefLevel &bit_level0,
246 PetscBool is_partitioned) {
248 DMType dm_name = "DM_TOPO";
249 // Register DM problem
250 CHKERR DMRegister_MoFEM(dm_name);
251 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
252 CHKERR DMSetType(dm, dm_name);
253 CHKERR DMMoFEMCreateMoFEM(dm, &m_field, dm_name, bit_level0);
254 CHKERR DMSetFromOptions(dm);
255 CHKERR DMMoFEMAddElement(dm, "ELASTIC");
256 CHKERR DMMoFEMAddElement(dm, "RHO_APPROX");
257 if (m_field.check_finite_element("PRESSURE_FE"))
258 CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
259 if (m_field.check_finite_element("FORCE_FE"))
260 CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
261 CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
262 m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
263 CHKERR DMSetUp(dm);
264 m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
265
266
268 }
269
270 MoFEMErrorCode build_sub_dms(DM &dm, DM &sub_dm_elastic, DM &sub_dm_rho,
271 MoFEM::Interface &m_field,
272 BitRefLevel &bit_level0,
273 PetscBool is_partitioned) {
274
276 DMType dm_name_sub_elast = "DM_ELASTIC";
277 CHKERR DMRegister_MoFEM(dm_name_sub_elast);
278 CHKERR DMCreate(PETSC_COMM_WORLD, &sub_dm_elastic);
279 CHKERR DMSetType(sub_dm_elastic, dm_name_sub_elast);
280 // m_field.getInterface<ProblemsManager>()->buildProblemFromFields =
281 // PETSC_TRUE;
282 CHKERR DMMoFEMCreateSubDM(sub_dm_elastic, dm, "ELASTIC_PROB");
283 CHKERR DMMoFEMSetDestroyProblem(sub_dm_elastic, PETSC_TRUE);
284 CHKERR DMMoFEMAddSubFieldRow(sub_dm_elastic, "DISPLACEMENTS");
285 CHKERR DMMoFEMAddSubFieldCol(sub_dm_elastic, "DISPLACEMENTS");
286 CHKERR DMMoFEMAddElement(sub_dm_elastic, "ELASTIC");
287 if (m_field.check_finite_element("PRESSURE_FE"))
288 CHKERR DMMoFEMAddElement(sub_dm_elastic, "PRESSURE_FE");
289 if (m_field.check_finite_element("FORCE_FE"))
290 CHKERR DMMoFEMAddElement(sub_dm_elastic, "FORCE_FE");
291
292 CHKERR DMMoFEMSetSquareProblem(sub_dm_elastic, PETSC_TRUE);
293 CHKERR DMMoFEMSetIsPartitioned(sub_dm_elastic, is_partitioned);
294 CHKERR DMSetUp(sub_dm_elastic);
295
296 DMType dm_name_sub_rho = "DM_DENSITY";
297 CHKERR DMRegister_MoFEM(dm_name_sub_rho);
298 CHKERR DMCreate(PETSC_COMM_WORLD, &sub_dm_rho);
299 CHKERR DMSetType(sub_dm_rho, dm_name_sub_rho);
300 // m_field.getInterface<ProblemsManager>()->buildProblemFromFields =
301 // PETSC_TRUE;
302 CHKERR DMMoFEMCreateSubDM(sub_dm_rho, dm, "DENSITY_PROB");
303 CHKERR DMMoFEMSetDestroyProblem(sub_dm_rho, PETSC_TRUE);
304 CHKERR DMMoFEMAddSubFieldRow(sub_dm_rho, "RHO");
305 CHKERR DMMoFEMAddSubFieldCol(sub_dm_rho, "RHO");
306 CHKERR DMMoFEMAddElement(sub_dm_rho, "RHO_APPROX");
307 CHKERR DMMoFEMSetSquareProblem(sub_dm_rho, PETSC_TRUE);
308 CHKERR DMMoFEMSetIsPartitioned(sub_dm_rho, is_partitioned);
309 CHKERR DMSetUp(sub_dm_rho);
311 }
312
313 MoFEMErrorCode getBoundBox(const double material_coords[3],
314 std::vector<double> &coords, double &boxSize) {
316 FTensor::Tensor1<double, 3> t_material_coords(
317 material_coords[0], material_coords[1], material_coords[2]);
318 FTensor::Tensor1<double *, 3> t_coords(&coords[0], &coords[1], &coords[2],
319 3);
321 double max_xyz[] = {0, 0, 0};
322 for (int bb = 0; bb != coords.size() / 3; bb++) {
323 t_coords(i) -= t_material_coords(i);
324 for (int dd = 0; dd != 3; dd++) {
325 max_xyz[dd] = (max_xyz[dd] < fabs(t_coords(dd))) ? fabs(t_coords(dd))
326 : max_xyz[dd];
327 }
328 ++t_coords;
329 }
330 //take the max of the box
331 boxSize = (max_xyz[0] > max_xyz[1]) ? max_xyz[0] : max_xyz[1];
332 boxSize = (boxSize > max_xyz[2]) ? boxSize : max_xyz[2];
334 }
335};
336
337
339 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
340
341 MatrixDouble rowB;
342 MatrixDouble colB;
343 MatrixDouble CB;
344 MatrixDouble K, transK;
345
347
348 OpAssembleK(CommonData &common_data, bool symm = true)
349 : VolumeElementForcesAndSourcesCore::UserDataOperator(
350 "DISPLACEMENTS", "DISPLACEMENTS", OPROWCOL, symm),
351 commonData(common_data) {}
352
353 /**
354 * \brief Do calculations for give operator
355 * @param row_side row side number (local number) of entity on element
356 * @param col_side column side number (local number) of entity on element
357 * @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
358 * @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
359 * @param row_data data for row
360 * @param col_data data for column
361 * @return error code
362 */
363 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
364 EntityType col_type,
365 DataForcesAndSourcesCore::EntData &row_data,
366 DataForcesAndSourcesCore::EntData &col_data) {
368 // get number of dofs on row
369 nbRows = row_data.getIndices().size();
370 // if no dofs on row, exit that work, nothing to do here
371 if (!nbRows)
373 // get number of dofs on column
374 nbCols = col_data.getIndices().size();
375 // if no dofs on Columbia, exit nothing to do here
376 if (!nbCols)
378 // get number of integration points
379 nbIntegrationPts = getGaussPts().size2();
380 // chekk if entity block is on matrix diagonal
381 if (row_side == col_side && row_type == col_type) { // ASK
382 isDiag = true; // yes, it is
383 } else {
384 isDiag = false;
385 }
386
387 // integrate local matrix for entity block
388 CHKERR iNtegrate(row_data, col_data);
389 // asseble local matrix
390 CHKERR aSsemble(row_data, col_data);
391
393 }
394
395protected:
396 ///< error code
397
398 int nbRows; ///< number of dofs on rows
399 int nbCols; ///< number if dof on column
400 int nbIntegrationPts; ///< number of integration points
401 bool isDiag; ///< true if this block is on diagonal
402
403 /**
404 * \brief Integrate grad-grad operator
405 * @param row_data row data (consist base functions on row entity)
406 * @param col_data column data (consist base functions on column entity)
407 * @return error code
408 */
409 virtual MoFEMErrorCode
410 iNtegrate(DataForcesAndSourcesCore::EntData &row_data,
411 DataForcesAndSourcesCore::EntData &col_data) {
413
414 int nb_dofs_row = row_data.getFieldData().size();
415 if (nb_dofs_row == 0)
417 int nb_dofs_col = col_data.getFieldData().size();
418 if (nb_dofs_col == 0)
420
421 K.resize(nb_dofs_row, nb_dofs_col, false);
422 K.clear();
423 auto max = [](double a, double b) { return a > b ? a : b; };
424 auto min = [](double a, double b) { return a < b ? a : b; };
425 auto rho = getFTensor0FromVec(*commonData.dataRhotGaussPtr);
426
427 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
429 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
430 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
431 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
432 };
433
438
439 // get element volume
440 double vol = getVolume();
441
442 // get intergrayion weights
443 auto t_w = getFTensor0IntegrationWeight();
444
445 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
446 const int row_nb_base_fun = row_data.getN().size2();
447
450
451 for (int gg = 0; gg != nbIntegrationPts; gg++) {
452
453 const double coeff = pow(rho, commonData.penalty);
454 double a = t_w * vol;
455 if (getHoGaussPtsDetJac().size()) {
456 a *= getHoGaussPtsDetJac()[gg];
457 }
458 a *= coeff;
459 int rr = 0;
460 for (; rr != nbRows / 3; ++rr) {
461 auto t_m = get_tensor2(K, 3 * rr, 0);
462 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
464 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
465 for (int cc = 0; cc != nbCols / 3; ++cc) {
466 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
467 ++t_col_diff_base;
468 ++t_m;
469 }
470 ++t_row_diff_base;
471 }
472 for (; rr != row_nb_base_fun; ++rr)
473 ++t_row_diff_base;
474 ++t_w;
475 ++rho;
476 }
478 }
479
480 /**
481 * \brief Assemble local entity block matrix
482 * @param row_data row data (consist base functions on row entity)
483 * @param col_data column data (consist base functions on column entity)
484 * @return error code
485 */
486 virtual MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data,
487 DataForcesAndSourcesCore::EntData &col_data) {
489
490 int nb_dofs_row = row_data.getFieldData().size();
491 if (nb_dofs_row == 0)
493 int nb_dofs_col = col_data.getFieldData().size();
494 if (nb_dofs_col == 0)
496
497 // get pointer to first global index on row
498 const int *row_indices = &*row_data.getIndices().data().begin();
499 // get pointer to first global index on column
500 const int *col_indices = &*col_data.getIndices().data().begin();
501 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
502 : getFEMethod()->snes_B;
503
504 CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
505 &*K.data().begin(), ADD_VALUES);
506 if (!isDiag && sYmm) {
507 K = trans(K);
508 CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
509 &*K.data().begin(), ADD_VALUES);
510 }
512 }
513};
514
515
516
520 // double &volumeVec;
522 : VolumeElementForcesAndSourcesCore::UserDataOperator(
523 "RHO", UserDataOperator::OPROW),
524 commonData(common_data) {}
525
526 MoFEMErrorCode doWork(int side, EntityType type,
527 DataForcesAndSourcesCore::EntData &data) {
529 if (type != MBVERTEX)
531 int nb_dofs = data.getFieldData().size();
532 if (nb_dofs == 0) {
534 }
535 const double p = commonData.penalty;
536 const EntityHandle fe_ent = getFEEntityHandle();
537 const int nb_gauss_pts = getGaussPts().size2();
538
541
542 auto rho = getFTensor0FromVec(*commonData.dataRhotGaussPtr);
543 commonData.rhoVecAtTag.resize(nb_gauss_pts);
544 commonData.rhoVecAtTag.clear();
545 auto rho_tag = getFTensor0FromVec(commonData.rhoVecAtTag);
546 auto grad = getFTensor2FromMat<3, 3>(*commonData.dispGradAtGaussPts);
547
549 const double lambda = commonData.lAmbda;
550 const double mu = commonData.mU;
551
552 commonData.getSensitivityFromTag(fe_ent, nb_gauss_pts);
553 commonData.getDensityFromTag(fe_ent, nb_gauss_pts);
554 auto t_sensitivity = getFTensor0FromVec(commonData.sEnsitivityVec);
555
556 double newSensitivity = 0;
557 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
558
559 double vol = getVolume() * getGaussPts()(3, gg);
560 strain(i, j) = 0.5 * (grad(i, j) || grad(j, i));
561 const double trace = strain(i, i);
562
563 double energy =
564 0.5 * trace * trace * lambda + mu * strain(i, j) * strain(i, j);
565 rho = rho > 1. ? 1. : (rho < 0.001 ? 0.001 : rho); //FIXME: use lobatto quadrature to avoid this
566 t_sensitivity = (-p) * pow(rho, p - 1) * energy * vol; //
567 // if(rho < 0.001 || rho > 1.00001)
568 // SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
569 // "rho is not within boundaries, that should not happen");
570 if (t_sensitivity > 0) //FIXME:
571 SETERRQ(
572 PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
573 "t_sensitivity is not within boundaries, that should not happen");
574 rho_tag = rho;
575
576 ++t_sensitivity;
577 ++rho;
578 ++grad;
579 ++rho_tag;
580 }
581
582 commonData.setDensityToTag(fe_ent, nb_gauss_pts);
583 commonData.setSensitivityToTag(fe_ent, nb_gauss_pts);
584
586 }
587};
588
589
590template <bool UPDATE = true>
595 : VolumeElementForcesAndSourcesCore::UserDataOperator(
596 "RHO", UserDataOperator::OPROW),
597 commonData(common_data) {}
598
599 MoFEMErrorCode doWork(int side, EntityType type,
600 DataForcesAndSourcesCore::EntData &data) {
602 if (type != MBVERTEX)
604 int nb_dofs = data.getFieldData().size();
605 if (nb_dofs == 0) {
607 }
608
609 int nb_gauss_pts = data.getN().size1();
610
612 commonData.getDensityFromTag(fe_ent, nb_gauss_pts);
613 const int &p = commonData.penalty;
614 auto rho_old = getFTensor0FromVec(commonData.rhoVecAtTag);
615 auto rho_new = getFTensor0FromVec(*commonData.dataRhotGaussPtr);
616 const double &mv = commonData.mOve;
617 commonData.getSensitivityFromTag(fe_ent, nb_gauss_pts);
618 auto t_sens = getFTensor0FromVec(commonData.sEnsitivityVec);
619
620 auto max = [](double a, double b) { return a > b ? a : b; };
621 auto min = [](double a, double b) { return a < b ? a : b; };
622
623
624 double change = 0;
625 for (int gg = 0; gg != nb_gauss_pts; gg++) {
626
627 VectorAdaptor shape_function = data.getN(gg, nb_dofs);
628
629 const double vol = getVolume() * getGaussPts()(3,gg);
630 const double sensitivity = t_sens;
631 // double rho = rho_old > 1.0 ? 1. : (rho_old < 0.001 ? 0.001 : rho_old);
632 double rho = rho_old ;
633
634 const double &x = rho;
635
636 //FROM ZEGARD PAPER
637 const double den = commonData.L_mid * vol;
638 const double xnew = rho * sqrt(-sensitivity/den);
639
640 rho_new = max(0.001, min(1.0, max(x - mv, min(1.0, min(x + mv, xnew)))));
641 if (rho_new < 0.001 || rho_new > 1.00001)
642 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
643 "rho_new is not within boundaries, that should not happen");
644 if (rho_old < 0.001 || rho_old > 1.00001)
645 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
646 "rho_new is not within boundaries, that should not happen");
647 const double &check = rho_new;
648 if (!std::isnormal(check))
649 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
650 "rho_new is not normal, that should not happen"); //FIXME:
651 if (UPDATE)
652 change += vol * fabs(rho_old - rho_new);
653
654 ++t_sens;
655 ++rho_new;
656 ++rho_old;
657 }
658 if (UPDATE)
659 CHKERR VecSetValue(commonData.changeVec,0,change,ADD_VALUES);
660
662 }
663};
664
665
668
670 OpVolumeCalculation(Vec &volume_vec)
671 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
672 "RHO", UserDataOperator::OPROW),
673 volumeVec(volume_vec) {}
674
675 MoFEMErrorCode doWork(int row_side, EntityType row_type,
676 DataForcesAndSourcesCore::EntData &row_data) {
678 if (row_type != MBVERTEX)
680 double vol = getVolume();
681 CHKERR VecSetValue(volumeVec, 0, vol, ADD_VALUES);
682
684 }
685};
686
689
692
693 OpMassCalculation(const string &field_name, CommonData &common_data,
694 Vec mass_vec)
695 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
697 commonData(common_data), massVec(mass_vec) {}
698
699 MoFEMErrorCode doWork(int row_side, EntityType row_type,
700 DataForcesAndSourcesCore::EntData &row_data) {
702 if (row_type != MBVERTEX)
704 int nb_gauss_pts = row_data.getN().size1();
705 auto rho = getFTensor0FromVec(*commonData.dataRhotGaussPtr);
706 double mass = 0;
707 for (int gg = 0; gg < nb_gauss_pts; gg++) {
708
709 double vol = getVolume() * getGaussPts()(3, gg);
710 if (getHoGaussPtsDetJac().size() > 0) {
711 vol *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
712 }
713 mass += rho * vol;
714 // commonData.cUrrent_mass += rho * vol;
715 ++rho;
716 }
717 CHKERR VecSetValue(commonData.massVec,0,mass,ADD_VALUES);
719 }
720};
721
724
726 MatrixDouble locK;
727 MatrixDouble transK;
728
729 OpAssembleRhoLhs(const std::string &row_field, const std::string &col_field,
730 CommonData &common_data)
731 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
732 row_field, col_field, UserDataOperator::OPROWCOL),
733 commonData(common_data) {
734 sYmm = true;
735 }
736
737 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
738 EntityType col_type,
739 DataForcesAndSourcesCore::EntData &row_data,
740 DataForcesAndSourcesCore::EntData &col_data) {
741
743
744 int nb_rows = row_data.getIndices().size(); // number of rows
745 int nb_cols = col_data.getIndices().size(); // number of columns
746 if (nb_rows == 0)
748 if (nb_cols == 0)
750
751 locK.resize(nb_rows, nb_cols, false);
752 locK.clear();
753
754 // double lsquared = commonData.lAmbda0 * commonData.lAmbda0 *
755 // commonData.edgeLength * commonData.edgeLength;
756 MatrixDouble anisotropic_filter = commonData.anisLambdaMat;
757 anisotropic_filter *= commonData.edgeLength;
758 // anisotropic_filter.resize(3,3);
759 // anisotropic_filter.clear();
760 // anisotropic_filter(0,0) = sqrt(lsquared);
761 // anisotropic_filter(1,1) = sqrt(lsquared);
762 // anisotropic_filter(2,2) = 1e2 * commonData.edgeLength;
763
764 // lsquared = commonData.lAmbda0; //FIXME:
765 int nb_gauss_pts = row_data.getN().size1();
766 for (int gg = 0; gg < nb_gauss_pts; gg++) {
767 double vol = getVolume() * getGaussPts()(3, gg);
768 if (getHoGaussPtsDetJac().size() > 0) {
769 vol *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
770 }
771 VectorAdaptor row_base_function = row_data.getN(gg);
772 VectorAdaptor col_base_function = col_data.getN(gg);
773 locK += vol * outer_prod(row_base_function, col_base_function);
774
775 MatrixAdaptor row_diff_base_function = row_data.getDiffN(gg);
776 MatrixAdaptor col_diff_base_function = col_data.getDiffN(gg);
777
778 row_diff_base_function = prod(row_diff_base_function, anisotropic_filter);
779 col_diff_base_function = prod(col_diff_base_function, anisotropic_filter);
780
781 // locK += vol * lsquared * //FIXME:
782 locK += vol * prod(row_diff_base_function, trans(col_diff_base_function));
783 }
784 CHKERR MatSetValues(getFEMethod()->ksp_B, row_data.getIndices().size(),
785 &row_data.getIndices()[0], col_data.getIndices().size(),
786 &col_data.getIndices()[0], &locK(0, 0), ADD_VALUES);
787
788 // Assemble off symmetric side
789 if (row_side != col_side || row_type != col_type) {
790 transK.resize(locK.size2(), locK.size1(), false);
791 noalias(transK) = trans(locK);
792 CHKERR MatSetValues(getFEMethod()->ksp_B, col_data.getIndices().size(),
793 &col_data.getIndices()[0],
794 row_data.getIndices().size(),
795 &row_data.getIndices()[0], &transK(0, 0), ADD_VALUES);
796 }
797
799 }
800
801 };
802
805
807 VectorDouble nF;
808 OpAssembleDensityRhs(const string &row_field, CommonData &common_data)
809 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
810 row_field, row_field, UserDataOperator::OPROW),
811 commonData(common_data) {}
812
813 MoFEMErrorCode doWork(int row_side, EntityType row_type,
814 DataForcesAndSourcesCore::EntData &row_data) {
816
817 int nb_rows = row_data.getIndices().size();
818 if (nb_rows == 0)
820 nF.resize(nb_rows, false);
821 nF.clear();
822
823 auto rho = getFTensor0FromVec(*commonData.dataRhotGaussPtr);
824 // auto rho_old = getFTensor0FromVec(commonData.rhoVecAtTag);
825 int nb_gauss_pts = row_data.getN().size1();
826
827 for (int gg = 0; gg < nb_gauss_pts; gg++) {
828
829 double vol = getVolume() * getGaussPts()(3, gg);
830 if (getHoGaussPtsDetJac().size() > 0) {
831 vol *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
832 }
833 VectorAdaptor base_function = row_data.getN(gg);
834
835 nF += vol * rho * base_function;
836
837 if (rho > 1.000001 || rho < 0.001) // for debug only
838 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, //FIXME:
839 "wrong rho. data inconsistency");
840
841 ++rho;
842 }
843
844 CHKERR VecSetValues(getFEMethod()->ksp_f, row_data.getIndices().size(),
845 &row_data.getIndices()[0], &nF[0], ADD_VALUES);
846
848 }
849 };
850
853 moab::Interface &postProcMesh;
854 std::vector<EntityHandle> &mapGaussPts;
856
857 OpPostProcCompliance(moab::Interface &post_proc_mesh,
858 std::vector<EntityHandle> &map_gauss_pts,
859 CommonData &common_data)
860 : VolumeElementForcesAndSourcesCore::UserDataOperator(
861 "DISPLACEMENTS", UserDataOperator::OPROW),
862 postProcMesh(post_proc_mesh),mapGaussPts(map_gauss_pts),commonData(common_data) {}
863
864 MoFEMErrorCode doWork(int side, EntityType type,
865 DataForcesAndSourcesCore::EntData &data) {
867 if (type != MBVERTEX)
869 double def_VAL[9];
870 bzero(def_VAL, 9 * sizeof(double));
874 Tag th_psi;
875
876 const double p = commonData.penalty;
877 const double lambda = commonData.lAmbda;
878 const double mu = commonData.mU;
879
880 CHKERR postProcMesh.tag_get_handle("compliance", 1, MB_TYPE_DOUBLE, th_psi,
881 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
882 auto rho = getFTensor0FromVec(*commonData.dataRhotGaussPtr);
883 auto grad = getFTensor2FromMat<3, 3>(*commonData.dispGradAtGaussPts);
884
885 const int nb_gauss_pts = mapGaussPts.size();
886
887 for (int gg = 0; gg != nb_gauss_pts; gg++) {
888
889 strain(i, j) = 0.5 * (grad(i, j) || grad(j, i));
890 const double trace = grad(i, i);
891 const double energy =
892 // const double energy = pow(rho, p) *
893 0.5 * lambda * trace * trace + mu * strain(i, j) * strain(i, j);
894
895 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &energy);
896
897 ++rho;
898 ++grad;
899 }
900
902 }
903};
904
906
907 DirichletDensityBc(MoFEM::Interface &m_field, const std::string &field_name,
908 Mat aij, Vec x, Vec f)
909 : DirichletDisplacementBc(m_field, field_name, aij, x, f) {}
910
911 DirichletDensityBc(MoFEM::Interface &m_field, const std::string &field_name)
913
914 MoFEMErrorCode iNitalize() {
916 if (mapZeroRows.empty() || !methodsOp.empty()) {
917 ParallelComm *pcomm =
918 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
919
921 string name = it->getName();
922 if (name.compare(0, 7, "DENSITY") == 0) {
923 // if (name.compare(0, 7, "fdsfdsfsd") == 0) {
924
925 VectorDouble scaled_values(1);
926 scaled_values[0] = 1;
927
928 for (int dim = 0; dim < 3; dim++) {
929 Range ents;
930 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), dim,
931 ents, true);
932 if (dim > 1) {
933 Range _edges;
934 CHKERR mField.get_moab().get_adjacencies(ents, 1, false, _edges,
935 moab::Interface::UNION);
936 ents.insert(_edges.begin(), _edges.end());
937 }
938 if (dim > 0) {
939 Range _nodes;
940 CHKERR mField.get_moab().get_connectivity(ents, _nodes, true);
941 ents.insert(_nodes.begin(), _nodes.end());
942 }
943 auto &dofs_by_uid = problemPtr->getNumeredRowDofs()->get<Unique_mi_tag>();
944 for (auto eit = ents.pair_begin(); eit != ents.pair_end(); ++eit) {
945 auto bit_number = mField.get_field_bit_number(fieldName);
946 auto lo_dit = dofs_by_uid.lower_bound(
947 DofEntity::getLoFieldEntityUId(bit_number, eit->first));
948 auto hi_dit = dofs_by_uid.upper_bound(
949 DofEntity::getHiFieldEntityUId(bit_number, eit->second));
950 for (; lo_dit != hi_dit; ++lo_dit) {
951 auto &dof = *lo_dit;
952 if (dof->getEntType() == MBVERTEX) {
953 mapZeroRows[dof->getPetscGlobalDofIdx()] = scaled_values[0];
954 } else {
955 mapZeroRows[dof->getPetscGlobalDofIdx()] = 0;
956 }
957 }
958 }
959 // for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
960 // for (_IT_NUMEREDDOF_ROW_BY_NAME_ENT_PART_FOR_LOOP_(
961 // problemPtr, fieldName, *eit, pcomm->rank(), dof_ptr)) {
962 // NumeredDofEntity *dof = dof_ptr->get();
963 // if (dof->getEntType() == MBVERTEX) {
964 // mapZeroRows[dof->getPetscGlobalDofIdx()] = scaled_values[0];
965 // // dof->getFieldData() = scaled_values[0];
966 // } else {
967 // mapZeroRows[dof->getPetscGlobalDofIdx()] = 0;
968 // // dof->getFieldData() = 0;
969 // }
970 // }
971 // }
972 // }
973 }
974 }
975 }
976 dofsIndices.resize(mapZeroRows.size());
977 dofsValues.resize(mapZeroRows.size());
978 int ii = 0;
979 std::map<DofIdx, FieldData>::iterator mit = mapZeroRows.begin();
980 for (; mit != mapZeroRows.end(); mit++, ii++) {
981 dofsIndices[ii] = mit->first;
982 dofsValues[ii] = mit->second;
983 }
984 }
986 }
987
988 MoFEMErrorCode postProcess() {
990
991 switch (ts_ctx) {
992 case CTX_TSSETIFUNCTION: {
994 snes_x = ts_u;
995 snes_f = ts_F;
996 break;
997 }
998 case CTX_TSSETIJACOBIAN: {
1000 snes_B = ts_B;
1001 break;
1002 }
1003 default:
1004 break;
1005 }
1006
1007 if (snes_B) {
1008 CHKERR MatAssemblyBegin(snes_B, MAT_FINAL_ASSEMBLY);
1009 CHKERR MatAssemblyEnd(snes_B, MAT_FINAL_ASSEMBLY);
1010 CHKERR MatZeroRowsColumns(snes_B, dofsIndices.size(),
1011 dofsIndices.empty() ? PETSC_NULL
1012 : &dofsIndices[0],
1013 dIag, PETSC_NULL, PETSC_NULL);
1014 }
1015 if (snes_f) {
1016 CHKERR VecAssemblyBegin(snes_f);
1017 CHKERR VecAssemblyEnd(snes_f);
1018 for (std::vector<int>::iterator vit = dofsIndices.begin();
1019 vit != dofsIndices.end(); vit++) {
1020 CHKERR VecSetValue(snes_f, *vit, 0., INSERT_VALUES);
1021 }
1022 CHKERR VecAssemblyBegin(snes_f);
1023 CHKERR VecAssemblyEnd(snes_f);
1024 }
1025
1027 }
1028};
1029
1030#endif //__HEADER_HPP__
static Index< 'p', 3 > p
#define MAT_TO_DDG(SM)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
static PetscErrorCode ierr
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
double trace(FTensor::Tensor2< T, 2, 2 > &t_stress)
const int dim
FTensor::Index< 'm', SPACE_DIM > m
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
double rho
Definition: plastic.cpp:191
constexpr auto field_name
boost::shared_ptr< VectorDouble > dataRhotGaussPtr
Definition: Header.hpp:28
Tag thDens
Definition: Header.hpp:23
MoFEMErrorCode build_main_dm(DM &dm, DM &sub_dm_elastic, DM &sub_dm_rho, MoFEM::Interface &m_field, BitRefLevel &bit_level0, PetscBool is_partitioned)
Definition: Header.hpp:244
double rHo
Definition: Header.hpp:44
MoFEM::Interface & mField
Definition: Header.hpp:21
VectorDouble sEnsitivityVec
Definition: Header.hpp:36
CommonData(MoFEM::Interface &m_field)
Definition: Header.hpp:56
vector< VectorDouble > sTress
Definition: Header.hpp:33
MoFEMErrorCode getParameters()
Definition: Header.hpp:76
Vec massVec
Definition: Header.hpp:25
MoFEMErrorCode setSensitivityToTag(const EntityHandle fe_ent, const double sens)
Definition: Header.hpp:189
double volumeFrac
Definition: Header.hpp:34
double L_mid
Definition: Header.hpp:53
MatrixDouble anisLambdaMat
Definition: Header.hpp:38
Tag thSensitivity
Definition: Header.hpp:22
double yOung
Definition: Header.hpp:45
double cHange
Definition: Header.hpp:50
MoFEMErrorCode build_sub_dms(DM &dm, DM &sub_dm_elastic, DM &sub_dm_rho, MoFEM::Interface &m_field, BitRefLevel &bit_level0, PetscBool is_partitioned)
Definition: Header.hpp:270
MoFEMErrorCode getSensitivityFromTag(const EntityHandle fe_ent, const int nb_gauss_pts)
Definition: Header.hpp:163
double lAmbda
Definition: Header.hpp:41
MoFEMErrorCode getBoundBox(const double material_coords[3], std::vector< double > &coords, double &boxSize)
Definition: Header.hpp:313
bool uPdate_density
Definition: Header.hpp:54
boost::shared_ptr< MatrixDouble > dispGradAtGaussPts
Definition: Header.hpp:29
VectorDouble sTrain
Definition: Header.hpp:32
double cUrrent_mass
Definition: Header.hpp:52
MatrixDouble cOmplianceMat
Definition: Header.hpp:35
double mOve
Definition: Header.hpp:48
double lAmbda0
Definition: Header.hpp:43
MoFEMErrorCode setDensityToTag(const EntityHandle fe_ent, const double rho)
Definition: Header.hpp:223
double fIlter_radius
Definition: Header.hpp:49
double mU
Definition: Header.hpp:42
VectorDouble rhoVecAtTag
Definition: Header.hpp:37
int penalty
Definition: Header.hpp:47
double edgeLength
Definition: Header.hpp:40
MoFEMErrorCode createTags(moab::Interface &moab)
Definition: Header.hpp:232
Vec changeVec
Definition: Header.hpp:26
double pOisson
Definition: Header.hpp:46
MoFEMErrorCode getDensityFromTag(const EntityHandle fe_ent, const int nb_gauss_pts)
Definition: Header.hpp:198
MatrixDouble D
Definition: Header.hpp:31
DirichletDensityBc(MoFEM::Interface &m_field, const std::string &field_name, Mat aij, Vec x, Vec f)
Definition: Header.hpp:907
MoFEMErrorCode iNitalize()
Definition: Header.hpp:914
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: Header.hpp:988
DirichletDensityBc(MoFEM::Interface &m_field, const std::string &field_name)
Definition: Header.hpp:911
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
std::map< DofIdx, FieldData > mapZeroRows
Definition: DirichletBC.hpp:70
const std::string fieldName
field name to set Dirichlet BC
Definition: DirichletBC.hpp:60
double dIag
diagonal value set on zeroed column and rows
Definition: DirichletBC.hpp:61
std::vector< int > dofsIndices
Definition: DirichletBC.hpp:71
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: DirichletBC.hpp:76
std::vector< double > dofsValues
Definition: DirichletBC.hpp:72
MoFEM::Interface & mField
Definition: DirichletBC.hpp:59
const Problem * problemPtr
raw pointer to problem
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Problem manager is used to build and partition problems.
Vec & snes_f
residual
Vec & snes_x
state vector
Mat & snes_B
preconditioner of jacobian matrix
SNESContext snes_ctx
Vec & ts_F
residual vector
TSContext ts_ctx
Vec & ts_u
state vector
Mat & ts_B
Preconditioner for ts_A.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
CommonData & commonData
Definition: Header.hpp:806
VectorDouble nF
Definition: Header.hpp:807
OpAssembleDensityRhs(const string &row_field, CommonData &common_data)
Definition: Header.hpp:808
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: Header.hpp:813
Assemble K *.
Definition: Header.hpp:339
virtual MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate grad-grad operator.
Definition: Header.hpp:410
int nbIntegrationPts
number of integration points
Definition: Header.hpp:400
CommonData & commonData
Definition: Header.hpp:346
MatrixDouble colB
Definition: Header.hpp:342
OpAssembleK(CommonData &common_data, bool symm=true)
Definition: Header.hpp:348
int nbCols
number if dof on column
Definition: Header.hpp:399
MatrixDouble transK
Definition: Header.hpp:344
MatrixDouble rowB
Definition: Header.hpp:341
MatrixDouble CB
Definition: Header.hpp:343
bool isDiag
true if this block is on diagonal
Definition: Header.hpp:401
MatrixDouble K
Definition: Header.hpp:344
virtual MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble local entity block matrix.
Definition: Header.hpp:486
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Do calculations for give operator.
Definition: Header.hpp:363
int nbRows
< error code
Definition: Header.hpp:398
MatrixDouble transK
Definition: Header.hpp:727
CommonData & commonData
Definition: Header.hpp:725
OpAssembleRhoLhs(const std::string &row_field, const std::string &col_field, CommonData &common_data)
Definition: Header.hpp:729
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Header.hpp:737
MatrixDouble locK
Definition: Header.hpp:726
OpCalculateSensitivities(CommonData &common_data)
Definition: Header.hpp:521
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Header.hpp:526
CommonData & commonData
Definition: Header.hpp:519
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: Header.hpp:699
OpMassCalculation(const string &field_name, CommonData &common_data, Vec mass_vec)
Definition: Header.hpp:693
CommonData & commonData
Definition: Header.hpp:691
Definition: Header.hpp:852
CommonData & commonData
Definition: Header.hpp:855
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Header.hpp:864
OpPostProcCompliance(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, CommonData &common_data)
Definition: Header.hpp:857
std::vector< EntityHandle > & mapGaussPts
Definition: Header.hpp:854
moab::Interface & postProcMesh
Definition: Header.hpp:853
OpUpdateDensity(CommonData &common_data)
Definition: Header.hpp:594
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Header.hpp:599
CommonData & commonData
Definition: Header.hpp:593
OpVolumeCalculation(Vec &volume_vec)
Definition: Header.hpp:670
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: Header.hpp:675