v0.16.0
Loading...
Searching...
No Matches
EssentialDisplacementCubitBcData.cpp
Go to the documentation of this file.
1/**
2 * @file EssentialDisplacementCubitBcData.cpp
3 * @brief Essential boundary conditions
4 * @version 13.1
5 * @date 2022-09-03
6 *
7 * @copyright Copyright (c) 2022
8 *
9 */
10
11#include <MoFEM.hpp>
12
13namespace MoFEM {
14
15EssentialPreProc<DisplacementCubitBcData>::EssentialPreProc(
16 MoFEM::Interface &m_field, boost::shared_ptr<FEMethod> fe_ptr,
17 std::vector<boost::shared_ptr<ScalingMethod>> smv, bool get_coords)
18 : mField(m_field), fePtr(fe_ptr), vecOfTimeScalingMethods(smv),
19 getCoords(get_coords) {}
20
21MoFEMErrorCode EssentialPreProc<DisplacementCubitBcData>::operator()() {
22 MOFEM_LOG_CHANNEL("WORLD");
23 MoFEMFunctionBegin;
24
25 if (auto fe_method_ptr = fePtr.lock()) {
26
27 auto bc_mng = mField.getInterface<BcManager>();
28 auto fb = mField.getInterface<FieldBlas>();
29 const auto problem_name = fe_method_ptr->problemPtr->getName();
30
31 for (auto bc : bc_mng->getBcMapByBlockName()) {
32 if (auto disp_bc = bc.second->dispBcPtr) {
33
34 auto &bc_id = bc.first;
35
36 auto regex_str = (boost::format("%s_(.*)") % problem_name).str();
37 if (std::regex_match(bc_id, std::regex(regex_str))) {
38
39 auto [field_name, block_name] =
40 BcManager::extractStringFromBlockId(bc_id, problem_name);
41
42 auto get_field_coeffs = [&](auto field_name) {
43 auto field_ptr = mField.get_field_structure(field_name);
44 return field_ptr->getNbOfCoeffs();
45 };
46 const auto nb_field_coeffs = get_field_coeffs(field_name);
47
48 MOFEM_LOG("WORLD", Sev::noisy)
49 << "Apply EssentialPreProc<DisplacementCubitBcData>: "
50 << problem_name << "_" << field_name << "_" << block_name;
51
52 FTensor::Tensor1<double, 3> t_angles{0., 0., 0.};
53 FTensor::Tensor1<double, 3> t_vals{0., 0., 0.};
54 FTensor::Tensor1<double, 3> t_off{0., 0., 0.};
55
56 if (auto ext_disp_bc =
57 dynamic_cast<DisplacementCubitBcDataWithRotation const *>(
58 disp_bc.get())) {
59 for (int a = 0; a != 3; ++a)
60 t_off(a) = ext_disp_bc->rotOffset[a];
61 }
62
63 auto scale_value = [&](const double &c) {
64 double val = c;
65 for (auto s : vecOfTimeScalingMethods) {
66 val *= s->getScale(fe_method_ptr->ts_t);
67 }
68 return val;
69 };
70
71 if (disp_bc->data.flag1 == 1)
72 t_vals(0) = scale_value(-disp_bc->data.value1);
73 if (disp_bc->data.flag2 == 1)
74 t_vals(1) = scale_value(-disp_bc->data.value2);
75 if (disp_bc->data.flag3 == 1)
76 t_vals(2) = scale_value(-disp_bc->data.value3);
77 if (disp_bc->data.flag4 == 1)
78 t_angles(0) = scale_value(-disp_bc->data.value4);
79 if (disp_bc->data.flag5 == 1)
80 t_angles(1) = scale_value(-disp_bc->data.value5);
81 if (disp_bc->data.flag6 == 1)
82 t_angles(2) = scale_value(-disp_bc->data.value6);
83
84 int coeff;
85 std::array<std::vector<double>, 3> coords;
86 int idx;
87
88 const bool is_rotation =
89 disp_bc->data.flag4 || disp_bc->data.flag5 || disp_bc->data.flag6;
90
91 auto lambda = [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
92 MoFEMFunctionBegin;
93
94 auto v = t_vals(coeff);
95 if (is_rotation) {
96 FTensor::Tensor1<double, 3> t_coords(
97 coords[0][idx], coords[1][idx], coords[2][idx]);
98 v += DisplacementCubitBcDataWithRotation::GetRotDisp(
99 t_angles, t_coords, t_off)(coeff);
100 }
101 if (getCoords) {
102 v += coords[coeff][idx];
103 }
104
105 field_entity_ptr->getEntFieldData()[coeff] = v;
106 ++idx;
107
108 MoFEMFunctionReturn(0);
109 };
110
111 auto zero_lambda =
112 [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
113 MoFEMFunctionBegin;
114 auto size = field_entity_ptr->getEntFieldData().size();
115 for (int i = coeff; i < size; i += nb_field_coeffs)
116 field_entity_ptr->getEntFieldData()[i] = 0;
117 MoFEMFunctionReturn(0);
118 };
119
120 auto verts = bc.second->bcEnts.subset_by_type(MBVERTEX);
121 auto not_verts = subtract(bc.second->bcEnts, verts);
122
123 if (getCoords || is_rotation) {
124 for (auto d : {0, 1, 2})
125 coords[d].resize(verts.size());
126 CHKERR mField.get_moab().get_coords(verts, &*coords[0].begin(),
127 &*coords[1].begin(),
128 &*coords[2].begin());
129 }
130
131 if (disp_bc->data.flag1 || disp_bc->data.flag5 ||
132 disp_bc->data.flag6) {
133 idx = 0;
134 coeff = 0;
135 CHKERR fb->fieldLambdaOnEntities(lambda, field_name, &verts);
136 CHKERR fb->fieldLambdaOnEntities(zero_lambda, field_name,
137 &not_verts);
138 }
139 if (disp_bc->data.flag2 || disp_bc->data.flag4 ||
140 disp_bc->data.flag6) {
141 idx = 0;
142 coeff = 1;
143 if (nb_field_coeffs > 1) {
144 CHKERR fb->fieldLambdaOnEntities(lambda, field_name, &verts);
145 CHKERR fb->fieldLambdaOnEntities(zero_lambda, field_name,
146 &not_verts);
147 }
148 }
149 if (disp_bc->data.flag3 || disp_bc->data.flag4 ||
150 disp_bc->data.flag5 || is_rotation) {
151 idx = 0;
152 coeff = 2;
153 if (nb_field_coeffs > 2) {
154 CHKERR fb->fieldLambdaOnEntities(lambda, field_name, &verts);
155 CHKERR fb->fieldLambdaOnEntities(zero_lambda, field_name,
156 &not_verts);
157 }
158 }
159 }
160 }
161 }
162
163 } else {
164 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
165 "Can not lock shared pointer");
166 }
167
168 MoFEMFunctionReturn(0);
169}
170
171EssentialPostProcRhs<DisplacementCubitBcData>::EssentialPostProcRhs(
172 MoFEM::Interface &m_field, boost::shared_ptr<FEMethod> fe_ptr, double diag,
173 SmartPetscObj<Vec> rhs)
174 : mField(m_field), fePtr(fe_ptr), vDiag(diag), vRhs(rhs) {}
175
176MoFEMErrorCode EssentialPostProcRhs<DisplacementCubitBcData>::operator()() {
177 MOFEM_LOG_CHANNEL("WORLD");
178 MoFEMFunctionBegin;
179
180 if (auto fe_method_ptr = fePtr.lock()) {
181
182 auto bc_mng = mField.getInterface<BcManager>();
183 auto vec_mng = mField.getInterface<VecManager>();
184 auto is_mng = mField.getInterface<ISManager>();
185
186 const auto problem_name = fe_method_ptr->problemPtr->getName();
187
188 SmartPetscObj<IS> is_sum;
189
190 for (auto bc : bc_mng->getBcMapByBlockName()) {
191 if (auto disp_bc = bc.second->dispBcPtr) {
192
193 auto &bc_id = bc.first;
194
195 auto regex_str = (boost::format("%s_(.*)") % problem_name).str();
196 if (std::regex_match(bc_id, std::regex(regex_str))) {
197
198 auto [field_name, block_name] =
199 BcManager::extractStringFromBlockId(bc_id, problem_name);
200
201 MOFEM_LOG("WORLD", Sev::noisy)
202 << "Apply EssentialPreProc<DisplacementCubitBcData>: "
203 << problem_name << "_" << field_name << "_" << block_name;
204
205 const bool is_rotation =
206 disp_bc->data.flag4 || disp_bc->data.flag5 || disp_bc->data.flag6;
207
208 auto ents = bc.second->bcEnts;
209
210 std::array<SmartPetscObj<IS>, 3> is_xyz;
211 auto prb_name = fe_method_ptr->problemPtr->getName();
212
213 if (disp_bc->data.flag1 || is_rotation) {
214 CHKERR is_mng->isCreateProblemFieldAndRankLocal(
215 prb_name, ROW, field_name, 0, 0, is_xyz[0], &ents);
216 }
217 if (disp_bc->data.flag2 || is_rotation) {
218 CHKERR is_mng->isCreateProblemFieldAndRankLocal(
219 prb_name, ROW, field_name, 1, 1, is_xyz[1], &ents);
220 }
221 if (disp_bc->data.flag3 || is_rotation) {
222 CHKERR is_mng->isCreateProblemFieldAndRankLocal(
223 prb_name, ROW, field_name, 2, 2, is_xyz[2], &ents);
224 }
225
226 auto get_is_sum = [](auto is1, auto is2) {
227 IS is;
228 CHK_THROW_MESSAGE(ISExpand(is1, is2, &is), "is sum");
229 return SmartPetscObj<IS>(is);
230 };
231
232 for (auto &is : is_xyz) {
233 if (is) {
234 if (!is_sum) {
235 is_sum = is;
236 } else {
237 is_sum = get_is_sum(is_sum, is);
238 }
239 }
240 }
241 }
242 }
243 }
244
245 if (is_sum) {
246 if (auto fe_ptr = fePtr.lock()) {
247
248 auto snes_ctx = fe_ptr->snes_ctx;
249 auto ts_ctx = fe_ptr->ts_ctx;
250
251 SmartPetscObj<Vec> f =
252 vRhs ? vRhs : SmartPetscObj<Vec>(fe_ptr->f, true);
253
254 if (fe_ptr->vecAssembleSwitch) {
255 if ((*fe_ptr->vecAssembleSwitch) && !vRhs) {
256 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
257 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
258 CHKERR VecAssemblyBegin(f);
259 CHKERR VecAssemblyEnd(f);
260 *fe_ptr->vecAssembleSwitch = false;
261 }
262 }
263
264 const int *index_ptr;
265 CHKERR ISGetIndices(is_sum, &index_ptr);
266 int size;
267 CHKERR ISGetLocalSize(is_sum, &size);
268 double *a;
269 CHKERR VecGetArray(f, &a);
270
271 auto tmp_x = vectorDuplicate(f);
272 CHKERR vec_mng->setLocalGhostVector(problem_name, ROW, tmp_x,
273 INSERT_VALUES, SCATTER_FORWARD);
274 const double *u;
275 CHKERR VecGetArrayRead(tmp_x, &u);
276
277 if (snes_ctx != FEMethod::CTX_SNESNONE ||
278 ts_ctx != FEMethod::CTX_TSNONE) {
279
280 auto x = fe_ptr->x;
281 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
282 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
283
284 const double *v;
285 CHKERR VecGetArrayRead(x, &v);
286
287 for (auto i = 0; i != size; ++i) {
288 a[index_ptr[i]] = vDiag * (v[index_ptr[i]] - u[index_ptr[i]]);
289 }
290
291 CHKERR VecRestoreArrayRead(x, &v);
292
293 } else {
294 for (auto i = 0; i != size; ++i) {
295 a[index_ptr[i]] = vDiag * u[index_ptr[i]];
296 }
297 }
298
299 CHKERR VecRestoreArrayRead(tmp_x, &u);
300 CHKERR VecRestoreArray(f, &a);
301 CHKERR ISRestoreIndices(is_sum, &index_ptr);
302 }
303 }
304
305 } else {
306 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
307 "Can not lock shared pointer");
308 }
309
310 MoFEMFunctionReturn(0);
311}
312
313EssentialPostProcLhs<DisplacementCubitBcData>::EssentialPostProcLhs(
314 MoFEM::Interface &m_field, boost::shared_ptr<FEMethod> fe_ptr, double diag,
315 SmartPetscObj<Mat> lhs, SmartPetscObj<AO> ao)
316 : mField(m_field), fePtr(fe_ptr), vDiag(diag), vLhs(lhs), vAO(ao) {}
317
318MoFEMErrorCode EssentialPostProcLhs<DisplacementCubitBcData>::operator()() {
319 MOFEM_LOG_CHANNEL("WORLD");
320 MoFEMFunctionBegin;
321
322 if (auto fe_method_ptr = fePtr.lock()) {
323
324 auto bc_mng = mField.getInterface<BcManager>();
325 auto is_mng = mField.getInterface<ISManager>();
326
327 const auto problem_name = fe_method_ptr->problemPtr->getName();
328
329 SmartPetscObj<IS> is_sum;
330
331 for (auto bc : bc_mng->getBcMapByBlockName()) {
332 if (auto disp_bc = bc.second->dispBcPtr) {
333
334 auto &bc_id = bc.first;
335
336 auto regex_str = (boost::format("%s_(.*)") % problem_name).str();
337 if (std::regex_match(bc_id, std::regex(regex_str))) {
338
339 auto [field_name, block_name] =
340 BcManager::extractStringFromBlockId(bc_id, problem_name);
341
342 MOFEM_LOG("WORLD", Sev::noisy)
343 << "Apply EssentialPreProc<DisplacementCubitBcData>: "
344 << problem_name << "_" << field_name << "_" << block_name;
345
346 const bool is_rotation =
347 disp_bc->data.flag4 || disp_bc->data.flag5 || disp_bc->data.flag6;
348
349 auto ents = bc.second->bcEnts;
350
351 std::array<SmartPetscObj<IS>, 3> is_xyz;
352 auto prb_name = fe_method_ptr->problemPtr->getName();
353
354 if (disp_bc->data.flag1 || is_rotation) {
355 CHKERR is_mng->isCreateProblemFieldAndRank(
356 prb_name, ROW, field_name, 0, 0, is_xyz[0], &ents);
357 }
358 if (disp_bc->data.flag2 || is_rotation) {
359 CHKERR is_mng->isCreateProblemFieldAndRank(
360 prb_name, ROW, field_name, 1, 1, is_xyz[1], &ents);
361 }
362 if (disp_bc->data.flag3 || is_rotation) {
363 CHKERR is_mng->isCreateProblemFieldAndRank(
364 prb_name, ROW, field_name, 2, 2, is_xyz[2], &ents);
365 }
366
367 auto get_is_sum = [](auto is1, auto is2) {
368 IS is;
369 CHK_THROW_MESSAGE(ISExpand(is1, is2, &is), "is sum");
370 return SmartPetscObj<IS>(is);
371 };
372
373 for (auto &is : is_xyz) {
374 if (is) {
375 if (!is_sum) {
376 is_sum = is;
377 } else {
378 is_sum = get_is_sum(is_sum, is);
379 }
380 }
381 }
382 }
383 }
384 }
385
386 if (is_sum) {
387 if (auto fe_ptr = fePtr.lock()) {
388 SmartPetscObj<Mat> B =
389 vLhs ? vLhs : SmartPetscObj<Mat>(fe_ptr->B, true);
390 // User is responsible for assembly if vLhs is provided
391 if ((*fe_ptr->matAssembleSwitch) && !vLhs) {
392 if (*fe_ptr->matAssembleSwitch) {
393 CHKERR MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
394 CHKERR MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
395 *fe_ptr->matAssembleSwitch = false;
396 }
397 }
398 if (vAO) {
399 MOFEM_LOG("WORLD", Sev::noisy) << "Apply AO to IS";
400 CHKERR AOApplicationToPetscIS(vAO, is_sum);
401 }
402 // ISView(is_sum, PETSC_VIEWER_STDOUT_WORLD);
403 CHKERR MatZeroRowsColumnsIS(B, is_sum, vDiag, PETSC_NULLPTR, PETSC_NULLPTR);
404 }
405 }
406
407 } else {
408 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
409 "Can not lock shared pointer");
410 }
411
412 MoFEMFunctionReturn(0);
413}
414
415EssentialPreProcReaction<DisplacementCubitBcData>::EssentialPreProcReaction(
416 MoFEM::Interface &m_field, boost::shared_ptr<FEMethod> fe_ptr,
417 SmartPetscObj<Vec> rhs, LogManager::SeverityLevel sev,
418 boost::shared_ptr<std::vector<double>> reaction_ptr,
419 boost::shared_ptr<BlockReactions> block_reactions)
420 : mField(m_field), fePtr(fe_ptr), vRhs(rhs), sevLevel(sev),
421 printBlockName(PETSC_FALSE), reactionPtr(reaction_ptr),
422 blockReactions(block_reactions),
423 reactionBlockName("") {
424
425 CHK_THROW_MESSAGE(PetscOptionsGetBool(PETSC_NULLPTR, PETSC_NULLPTR,
426 "-reaction_print_block_name",
427 &printBlockName, PETSC_NULLPTR),
428 "can not get option");
429 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "-reaction_block_name",
430 reactionBlockName, 255, PETSC_NULLPTR);
431}
432
433MoFEMErrorCode EssentialPreProcReaction<DisplacementCubitBcData>::operator()() {
434 MOFEM_LOG_CHANNEL("WORLD");
435 MoFEMFunctionBegin;
436
437 enum { X = 0, Y, Z, MX, MY, MZ, LAST };
438
439 if (auto fe_ptr = fePtr.lock()) {
440
441 SmartPetscObj<Vec> f = vRhs ? vRhs : SmartPetscObj<Vec>(fe_ptr->f, true);
442
443 if (fe_ptr->vecAssembleSwitch) {
444 if ((*fe_ptr->vecAssembleSwitch) && !vRhs) {
445 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
446 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
447 CHKERR VecAssemblyBegin(f);
448 CHKERR VecAssemblyEnd(f);
449 *fe_ptr->vecAssembleSwitch = false;
450 }
451 }
452
453 auto get_low_hi_uid_by_entities = [](auto bit_number, auto f, auto s) {
454 return std::make_pair(DofEntity::getLoFieldEntityUId(bit_number, f),
455 DofEntity::getHiFieldEntityUId(bit_number, s));
456 };
457
458 auto get_low_hi = [fe_ptr](auto lo_uid, auto hi_uid) {
459 auto it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
460 .lower_bound(lo_uid);
461 auto hi_it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
462 .upper_bound(hi_uid);
463 return std::make_pair(it, hi_it);
464 };
465
466 auto mpi_array_reduce = [this](auto &array) {
467 std::array<double, LAST> array_sum{0, 0, 0, 0, 0, 0};
468 MPI_Allreduce(&array[0], &array_sum[0], LAST, MPI_DOUBLE, MPI_SUM,
469 mField.get_comm());
470 return array_sum;
471 };
472
473 const double *a;
474 CHKERR VecGetArrayRead(f, &a);
475
476 auto bc_mng = mField.getInterface<BcManager>();
477 const auto problem_name = fe_ptr->problemPtr->getName();
478 const auto nb_local_dofs = fe_ptr->problemPtr->nbLocDofsRow;
479
480 std::array<double, LAST> total_reactions{0, 0, 0, 0, 0, 0};
481 if (blockReactions) {
482 blockReactions->clear();
483 }
484
485 for (auto bc : bc_mng->getBcMapByBlockName()) {
486 if (auto disp_bc = bc.second->dispBcPtr) {
487
488 auto &bc_id = bc.first;
489
490 auto regex_str = (boost::format("%s_(.*)") % problem_name).str();
491 if (std::regex_match(bc_id, std::regex(regex_str))) {
492
493 auto [field_name, block_name] =
494 BcManager::extractStringFromBlockId(bc_id, problem_name);
495
496 MOFEM_TAG_AND_LOG("WORLD", sevLevel, "Essential")
497 << "EssentialPreProc<DisplacementCubitBcData>: " << problem_name
498 << "_" << field_name << "_" << block_name;
499 auto bit_number = mField.get_field_bit_number(field_name);
500
501 FTensor::Tensor1<double, 3> t_off{0., 0., 0.};
502 if (auto ext_disp_bc =
503 dynamic_cast<DisplacementCubitBcDataWithRotation const *>(
504 disp_bc.get())) {
505 for (int a = 0; a != 3; ++a)
506 t_off(a) = ext_disp_bc->rotOffset[a];
507 }
508
509 auto verts = bc.second->bcEnts.subset_by_type(MBVERTEX);
510
511 FTensor::Index<'i', 3> i;
512 FTensor::Index<'j', 3> j;
513 FTensor::Index<'k', 3> k;
514
515 auto get_coords_vec = [&]() {
516 VectorDouble coords_vec(3 * verts.size());
517 if (verts.size()) {
518 CHKERR mField.get_moab().get_coords(verts, &*coords_vec.begin());
519 }
520 return coords_vec;
521 };
522
523 auto coords_vec = get_coords_vec();
524 std::array<double, LAST> reactions{0, 0, 0, 0, 0, 0};
525
526 for (auto pit = verts.const_pair_begin();
527 pit != verts.const_pair_end(); ++pit) {
528 auto [lo_uid, hi_uid] =
529 get_low_hi_uid_by_entities(bit_number, pit->first, pit->second);
530 auto [lo, hi] = get_low_hi(lo_uid, hi_uid);
531
532 for (; lo != hi; ++lo) {
533 const auto loc_dof = (*lo)->getPetscLocalDofIdx();
534 if (loc_dof < nb_local_dofs) {
535
536 const auto coeff = (*lo)->getDofCoeffIdx();
537
538 if (
539
540 ((disp_bc->data.flag1 || disp_bc->data.flag4) &&
541 coeff == 0) ||
542 ((disp_bc->data.flag2 || disp_bc->data.flag5) &&
543 coeff == 1) ||
544 ((disp_bc->data.flag3 || disp_bc->data.flag6) &&
545 coeff == 2)) {
546
547 const auto ent = (*lo)->getEnt();
548 reactions[coeff] += a[loc_dof];
549
550 auto force = [&]() {
551 FTensor::Tensor1<double, 3> t_force{0., 0., 0.};
552 t_force(coeff) = a[loc_dof];
553 return t_force;
554 };
555
556 auto coord = [&]() {
557 const auto idx = verts.index(ent);
558 FTensor::Tensor1<double, 3> t_coords{
559 coords_vec[3 * idx + X], coords_vec[3 * idx + Y],
560 coords_vec[3 * idx + Z]};
561 t_coords(i) -= t_off(i);
562 return t_coords;
563 };
564
565 auto moment = [&](auto &&t_force, auto &&t_coords) {
566 FTensor::Tensor1<double, 3> t_moment;
567 t_moment(i) =
568 (FTensor::levi_civita<double>(i, j, k) * t_coords(k)) *
569 t_force(j);
570 return t_moment;
571 };
572
573 auto t_moment = moment(force(), coord());
574 reactions[MX] += t_moment(X);
575 reactions[MY] += t_moment(Y);
576 reactions[MZ] += t_moment(Z);
577 }
578 }
579 }
580 }
581
582 FTensor::Tensor1<double, 3> t_force{reactions[X], reactions[Y],
583 reactions[Z]};
584 FTensor::Tensor1<double, 3> t_moment{reactions[MX], reactions[MY],
585 reactions[MZ]};
586 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_total_force{
587 &total_reactions[X], &total_reactions[Y], &total_reactions[Z]};
588 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_total_moment{
589 &total_reactions[MX], &total_reactions[MY], &total_reactions[MZ]};
590 t_total_force(i) += t_force(i);
591 t_total_moment(i) +=
592 t_moment(i) +
593 (FTensor::levi_civita<double>(i, j, k) * t_off(k)) * t_force(j);
594
595 auto mpi_reactions = mpi_array_reduce(reactions);
596 if (blockReactions) {
597 (*blockReactions)[block_name] = mpi_reactions;
598 }
599 if (printBlockName) {
600 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
601 "Block %s Offset: %6.4e %6.4e %6.4e",
602 block_name.c_str(), t_off(X), t_off(Y),
603 t_off(Z));
604 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
605 "Block %s Force: %6.4e %6.4e %6.4e",
606 block_name.c_str(), mpi_reactions[X],
607 mpi_reactions[Y], mpi_reactions[Z]);
608 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
609 "Block %s Moment: %6.4e %6.4e %6.4e",
610 block_name.c_str(), mpi_reactions[MX],
611 mpi_reactions[MY], mpi_reactions[MZ]);
612 if ((reactionPtr != nullptr) &&
613 (block_name.compare(reactionBlockName) == 0)) {
614 (*reactionPtr)[X] = mpi_reactions[X];
615 (*reactionPtr)[Y] = mpi_reactions[Y];
616 (*reactionPtr)[Z] = mpi_reactions[Z];
617 }
618 } else {
619 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
620 "Offset: %6.4e %6.4e %6.4e", t_off(X), t_off(Y),
621 t_off(Z));
622 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
623 "Force: %6.4e %6.4e %6.4e", mpi_reactions[X],
624 mpi_reactions[Y], mpi_reactions[Z]);
625 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
626 "Moment: %6.4e %6.4e %6.4e", mpi_reactions[MX],
627 mpi_reactions[MY], mpi_reactions[MZ]);
628 }
629 }
630 }
631 }
632
633 CHKERR VecRestoreArrayRead(f, &a);
634
635 auto mpi_total_reactions = mpi_array_reduce(total_reactions);
636 MOFEM_TAG_AND_LOG_C(
637 "WORLD", sevLevel, "Essential", "Total force: %6.4e %6.4e %6.4e",
638 mpi_total_reactions[X], mpi_total_reactions[Y], mpi_total_reactions[Z]);
639 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
640 "Total moment: %6.4e %6.4e %6.4e",
641 mpi_total_reactions[MX], mpi_total_reactions[MY],
642 mpi_total_reactions[MZ]);
643
644 } else {
645 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
646 "Can not lock shared pointer");
647 }
648
649 MoFEMFunctionReturn(0);
650}
651
652} // namespace MoFEM
MOFEM_TAG_AND_LOG_C
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
Definition LogManager.hpp:384
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition LogManager.hpp:376
a
constexpr double a
Definition approx_sphere.cpp:30
ROW
@ ROW
Definition definitions.h:136
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition definitions.h:612
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
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
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MoFEM::LogManager::SeverityLevel
SeverityLevel
Severity levels.
Definition LogManager.hpp:33
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition LogManager.hpp:292
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
v
const double v
phase velocity of light in medium (cm/ns)
Definition initial_diffusion.cpp:40
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition level_set.cpp:1923
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition DeprecatedPetsc.hpp:182
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition PetscSmartObj.hpp:223
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
MoFEM::BcManager
Boundary condition manager for finite element problem setup.
Definition BcManager.hpp:384
MoFEM::BcManager::extractStringFromBlockId
static std::pair< std::string, std::string > extractStringFromBlockId(const std::string block_id, const std::string prb_name)
Extract block name and block name from block id.
Definition BcManager.cpp:1223
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::DisplacementCubitBcDataWithRotation
A specialized version of DisplacementCubitBcData that includes an additional rotation offset.
Definition EssentialDisplacementCubitBcData.hpp:24
MoFEM::DisplacementCubitBcDataWithRotation::GetRotDisp
static FTensor::Tensor1< double, 3 > GetRotDisp(const FTensor::Tensor1< double, 3 > &angles, FTensor::Tensor1< double, 3 > coordinates, FTensor::Tensor1< double, 3 > offset=FTensor::Tensor1< double, 3 >{ 0., 0., 0.})
Calculates the rotated displacement given the rotation angles, coordinates, and an optional offset.
Definition EssentialDisplacementCubitBcData.hpp:40
MoFEM::DofEntity::getLoFieldEntityUId
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Definition DofsMultiIndices.hpp:69
MoFEM::DofEntity::getHiFieldEntityUId
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Definition DofsMultiIndices.hpp:75
MoFEM::EssentialPostProcLhs
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition Essential.hpp:33
MoFEM::EssentialPostProcRhs
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition Essential.hpp:41
MoFEM::EssentialPreProcReaction
Class (Function) to calculate residual side diagonal.
Definition Essential.hpp:49
MoFEM::EssentialPreProc
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
MoFEM::FieldBlas
Basic algebra on fields.
Definition FieldBlas.hpp:21
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85
MoFEM::SnesMethod::CTX_SNESNONE
@ CTX_SNESNONE
No specific SNES context.
Definition LoopMethods.hpp:215
MoFEM::TSMethod::CTX_TSNONE
@ CTX_TSNONE
No specific TS context.
Definition LoopMethods.hpp:277
MoFEM::Unique_mi_tag
Definition TagMultiIndices.hpp:18
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition VecManager.hpp:23
Generated by Doxygen 1.9.8 and hosted at University of Glasgow