v0.15.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 : mField(m_field), fePtr(fe_ptr), vRhs(rhs), sevLevel(sev),
420 printBlockName(PETSC_FALSE), reactionPtr(reaction_ptr),
421 reactionBlockName("") {
422
423 CHK_THROW_MESSAGE(PetscOptionsGetBool(PETSC_NULLPTR, PETSC_NULLPTR,
424 "-reaction_print_block_name",
425 &printBlockName, PETSC_NULLPTR),
426 "can not get option");
427 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "-reaction_block_name",
428 reactionBlockName, 255, PETSC_NULLPTR);
429}
430
431MoFEMErrorCode EssentialPreProcReaction<DisplacementCubitBcData>::operator()() {
432 MOFEM_LOG_CHANNEL("WORLD");
433 MoFEMFunctionBegin;
434
435 enum { X = 0, Y, Z, MX, MY, MZ, LAST };
436
437 if (auto fe_ptr = fePtr.lock()) {
438
439 SmartPetscObj<Vec> f = vRhs ? vRhs : SmartPetscObj<Vec>(fe_ptr->f, true);
440
441 if (fe_ptr->vecAssembleSwitch) {
442 if ((*fe_ptr->vecAssembleSwitch) && !vRhs) {
443 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
444 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
445 CHKERR VecAssemblyBegin(f);
446 CHKERR VecAssemblyEnd(f);
447 *fe_ptr->vecAssembleSwitch = false;
448 }
449 }
450
451 auto get_low_hi_uid_by_entities = [](auto bit_number, auto f, auto s) {
452 return std::make_pair(DofEntity::getLoFieldEntityUId(bit_number, f),
453 DofEntity::getHiFieldEntityUId(bit_number, s));
454 };
455
456 auto get_low_hi = [fe_ptr](auto lo_uid, auto hi_uid) {
457 auto it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
458 .lower_bound(lo_uid);
459 auto hi_it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
460 .upper_bound(hi_uid);
461 return std::make_pair(it, hi_it);
462 };
463
464 auto mpi_array_reduce = [this](auto &array) {
465 std::array<double, LAST> array_sum{0, 0, 0, 0, 0, 0};
466 MPI_Allreduce(&array[0], &array_sum[0], LAST, MPI_DOUBLE, MPI_SUM,
467 mField.get_comm());
468 return array_sum;
469 };
470
471 const double *a;
472 CHKERR VecGetArrayRead(f, &a);
473
474 auto bc_mng = mField.getInterface<BcManager>();
475 const auto problem_name = fe_ptr->problemPtr->getName();
476 const auto nb_local_dofs = fe_ptr->problemPtr->nbLocDofsRow;
477
478 std::array<double, LAST> total_reactions{0, 0, 0, 0, 0, 0};
479
480 for (auto bc : bc_mng->getBcMapByBlockName()) {
481 if (auto disp_bc = bc.second->dispBcPtr) {
482
483 auto &bc_id = bc.first;
484
485 auto regex_str = (boost::format("%s_(.*)") % problem_name).str();
486 if (std::regex_match(bc_id, std::regex(regex_str))) {
487
488 auto [field_name, block_name] =
489 BcManager::extractStringFromBlockId(bc_id, problem_name);
490
491 MOFEM_TAG_AND_LOG("WORLD", sevLevel, "Essential")
492 << "EssentialPreProc<DisplacementCubitBcData>: " << problem_name
493 << "_" << field_name << "_" << block_name;
494 auto bit_number = mField.get_field_bit_number(field_name);
495
496 FTensor::Tensor1<double, 3> t_off{0., 0., 0.};
497 if (auto ext_disp_bc =
498 dynamic_cast<DisplacementCubitBcDataWithRotation const *>(
499 disp_bc.get())) {
500 for (int a = 0; a != 3; ++a)
501 t_off(a) = ext_disp_bc->rotOffset[a];
502 }
503
504 auto verts = bc.second->bcEnts.subset_by_type(MBVERTEX);
505
506 FTensor::Index<'i', 3> i;
507 FTensor::Index<'j', 3> j;
508 FTensor::Index<'k', 3> k;
509
510 auto get_coords_vec = [&]() {
511 VectorDouble coords_vec(3 * verts.size());
512 if (verts.size()) {
513 CHKERR mField.get_moab().get_coords(verts, &*coords_vec.begin());
514 }
515 return coords_vec;
516 };
517
518 auto coords_vec = get_coords_vec();
519 std::array<double, LAST> reactions{0, 0, 0, 0, 0, 0};
520
521 for (auto pit = verts.const_pair_begin();
522 pit != verts.const_pair_end(); ++pit) {
523 auto [lo_uid, hi_uid] =
524 get_low_hi_uid_by_entities(bit_number, pit->first, pit->second);
525 auto [lo, hi] = get_low_hi(lo_uid, hi_uid);
526
527 for (; lo != hi; ++lo) {
528 const auto loc_dof = (*lo)->getPetscLocalDofIdx();
529 if (loc_dof < nb_local_dofs) {
530
531 const auto coeff = (*lo)->getDofCoeffIdx();
532
533 if (
534
535 ((disp_bc->data.flag1 || disp_bc->data.flag4) &&
536 coeff == 0) ||
537 ((disp_bc->data.flag2 || disp_bc->data.flag5) &&
538 coeff == 1) ||
539 ((disp_bc->data.flag3 || disp_bc->data.flag6) &&
540 coeff == 2)) {
541
542 const auto ent = (*lo)->getEnt();
543 reactions[coeff] += a[loc_dof];
544
545 auto force = [&]() {
546 FTensor::Tensor1<double, 3> t_force{0., 0., 0.};
547 t_force(coeff) = a[loc_dof];
548 return t_force;
549 };
550
551 auto coord = [&]() {
552 const auto idx = verts.index(ent);
553 FTensor::Tensor1<double, 3> t_coords{
554 coords_vec[3 * idx + X], coords_vec[3 * idx + Y],
555 coords_vec[3 * idx + Z]};
556 t_coords(i) -= t_off(i);
557 return t_coords;
558 };
559
560 auto moment = [&](auto &&t_force, auto &&t_coords) {
561 FTensor::Tensor1<double, 3> t_moment;
562 t_moment(i) =
563 (FTensor::levi_civita<double>(i, j, k) * t_coords(k)) *
564 t_force(j);
565 return t_moment;
566 };
567
568 auto t_moment = moment(force(), coord());
569 reactions[MX] += t_moment(X);
570 reactions[MY] += t_moment(Y);
571 reactions[MZ] += t_moment(Z);
572 }
573 }
574 }
575 }
576
577 FTensor::Tensor1<double, 3> t_force{reactions[X], reactions[Y],
578 reactions[Z]};
579 FTensor::Tensor1<double, 3> t_moment{reactions[MX], reactions[MY],
580 reactions[MZ]};
581 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_total_force{
582 &total_reactions[X], &total_reactions[Y], &total_reactions[Z]};
583 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_total_moment{
584 &total_reactions[MX], &total_reactions[MY], &total_reactions[MZ]};
585 t_total_force(i) += t_force(i);
586 t_total_moment(i) +=
587 t_moment(i) +
588 (FTensor::levi_civita<double>(i, j, k) * t_off(k)) * t_force(j);
589
590 auto mpi_reactions = mpi_array_reduce(reactions);
591 if (printBlockName) {
592 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
593 "Block %s Offset: %6.4e %6.4e %6.4e",
594 block_name.c_str(), t_off(X), t_off(Y),
595 t_off(Z));
596 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
597 "Block %s Force: %6.4e %6.4e %6.4e",
598 block_name.c_str(), mpi_reactions[X],
599 mpi_reactions[Y], mpi_reactions[Z]);
600 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
601 "Block %s Moment: %6.4e %6.4e %6.4e",
602 block_name.c_str(), mpi_reactions[MX],
603 mpi_reactions[MY], mpi_reactions[MZ]);
604 if ((reactionPtr != nullptr) &&
605 (block_name.compare(reactionBlockName) == 0)) {
606 (*reactionPtr)[X] = mpi_reactions[X];
607 (*reactionPtr)[Y] = mpi_reactions[Y];
608 (*reactionPtr)[Z] = mpi_reactions[Z];
609 }
610 } else {
611 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
612 "Offset: %6.4e %6.4e %6.4e", t_off(X), t_off(Y),
613 t_off(Z));
614 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
615 "Force: %6.4e %6.4e %6.4e", mpi_reactions[X],
616 mpi_reactions[Y], mpi_reactions[Z]);
617 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
618 "Moment: %6.4e %6.4e %6.4e", mpi_reactions[MX],
619 mpi_reactions[MY], mpi_reactions[MZ]);
620 }
621 }
622 }
623 }
624
625 CHKERR VecRestoreArrayRead(f, &a);
626
627 auto mpi_total_reactions = mpi_array_reduce(total_reactions);
628 MOFEM_TAG_AND_LOG_C(
629 "WORLD", sevLevel, "Essential", "Total force: %6.4e %6.4e %6.4e",
630 mpi_total_reactions[X], mpi_total_reactions[Y], mpi_total_reactions[Z]);
631 MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
632 "Total moment: %6.4e %6.4e %6.4e",
633 mpi_total_reactions[MX], mpi_total_reactions[MY],
634 mpi_total_reactions[MZ]);
635
636 } else {
637 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
638 "Can not lock shared pointer");
639 }
640
641 MoFEMFunctionReturn(0);
642}
643
644} // 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:1932
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition Levi_Civita.hpp:617
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
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
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 form 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::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.12.0 and hosted at University of Glasgow