v0.14.0
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 
13 namespace MoFEM {
14 
15 EssentialPreProc<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 
21 MoFEMErrorCode 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 
171 EssentialPostProcRhs<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 
176 MoFEMErrorCode 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 
313 EssentialPostProcLhs<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 
318 MoFEMErrorCode 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_NULL, PETSC_NULL);
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 
415 EssentialPreProcReaction<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_NULL, PETSC_NULL,
424  "-reaction_print_block_name",
425  &printBlockName, PETSC_NULL),
426  "can not get option");
427  CHKERR PetscOptionsGetString(PETSC_NULL, "-reaction_block_name",
428  reactionBlockName, 255, PETSC_NULL);
429 }
430 
431 MoFEMErrorCode 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::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
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:1381
FTensor::Tensor1< double, 3 >
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM.hpp
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:136
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
MoFEM::SnesMethod::CTX_SNESNONE
@ CTX_SNESNONE
Definition: LoopMethods.hpp:107
MoFEM::EssentialPostProcLhs
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition: Essential.hpp:33
MoFEM::TSMethod::CTX_TSNONE
@ CTX_TSNONE
Definition: LoopMethods.hpp:145
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::LogManager::SeverityLevel
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::EssentialPreProcReaction
Class (Function) to calculate residual side diagonal.
Definition: Essential.hpp:49
MoFEM::DofEntity::getLoFieldEntityUId
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Definition: DofsMultiIndices.hpp:69
MoFEM::EssentialPostProcRhs
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition: Essential.hpp:41
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
MoFEM::DofEntity::getHiFieldEntityUId
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Definition: DofsMultiIndices.hpp:75
MOFEM_TAG_AND_LOG_C
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
Definition: LogManager.hpp:370
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
UBlasVector< double >
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MoFEM::MPC::LAST
@ LAST
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::Unique_mi_tag
Definition: TagMultiIndices.hpp:18
MoFEM::DisplacementCubitBcDataWithRotation
A specialized version of DisplacementCubitBcData that includes an additional rotation offset.
Definition: EssentialDisplacementCubitBcData.hpp:24
MoFEM::EssentialPreProc
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:25
MoFEM::SmartPetscObj< Vec >
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow