17 std::vector<boost::shared_ptr<ScalingMethod>> smv,
bool get_coords)
18 : mField(m_field), fePtr(fe_ptr), vecOfTimeScalingMethods(smv),
19 getCoords(get_coords) {}
25 if (
auto fe_method_ptr = fePtr.lock()) {
27 auto bc_mng = mField.getInterface<
BcManager>();
29 const auto problem_name = fe_method_ptr->problemPtr->getName();
31 for (
auto bc : bc_mng->getBcMapByBlockName()) {
32 if (
auto disp_bc = bc.second->dispBcPtr) {
34 auto &bc_id = bc.first;
36 auto regex_str = (boost::format(
"%s_(.*)") % problem_name).str();
37 if (std::regex_match(bc_id, std::regex(regex_str))) {
43 auto field_ptr = mField.get_field_structure(
field_name);
44 return field_ptr->getNbOfCoeffs();
46 const auto nb_field_coeffs = get_field_coeffs(
field_name);
49 <<
"Apply EssentialPreProc<DisplacementCubitBcData>: "
50 << problem_name <<
"_" <<
field_name <<
"_" << block_name;
56 if (
auto ext_disp_bc =
59 for (
int a = 0;
a != 3; ++
a)
60 t_off(
a) = ext_disp_bc->rotOffset[
a];
63 auto scale_value = [&](
const double &
c) {
65 for (
auto s : vecOfTimeScalingMethods) {
66 val *= s->getScale(fe_method_ptr->ts_t);
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);
85 std::array<std::vector<double>, 3> coords;
88 const bool is_rotation =
89 disp_bc->data.flag4 || disp_bc->data.flag5 || disp_bc->data.flag6;
91 auto lambda = [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
94 auto v = t_vals(coeff);
97 coords[0][idx], coords[1][idx], coords[2][idx]);
99 t_angles, t_coords, t_off)(coeff);
102 v += coords[coeff][idx];
105 field_entity_ptr->getEntFieldData()[coeff] =
v;
112 [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
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;
120 auto verts = bc.second->bcEnts.subset_by_type(MBVERTEX);
121 auto not_verts = subtract(bc.second->bcEnts, verts);
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(),
128 &*coords[2].begin());
131 if (disp_bc->data.flag1 || disp_bc->data.flag5 ||
132 disp_bc->data.flag6) {
139 if (disp_bc->data.flag2 || disp_bc->data.flag4 ||
140 disp_bc->data.flag6) {
143 if (nb_field_coeffs > 1) {
149 if (disp_bc->data.flag3 || disp_bc->data.flag4 ||
150 disp_bc->data.flag5 || is_rotation) {
153 if (nb_field_coeffs > 2) {
165 "Can not lock shared pointer");
174 : mField(m_field), fePtr(fe_ptr), vDiag(diag), vRhs(rhs) {}
180 if (
auto fe_method_ptr = fePtr.lock()) {
182 auto bc_mng = mField.getInterface<
BcManager>();
186 const auto problem_name = fe_method_ptr->problemPtr->getName();
190 for (
auto bc : bc_mng->getBcMapByBlockName()) {
191 if (
auto disp_bc = bc.second->dispBcPtr) {
193 auto &bc_id = bc.first;
195 auto regex_str = (boost::format(
"%s_(.*)") % problem_name).str();
196 if (std::regex_match(bc_id, std::regex(regex_str))) {
202 <<
"Apply EssentialPreProc<DisplacementCubitBcData>: "
203 << problem_name <<
"_" <<
field_name <<
"_" << block_name;
205 const bool is_rotation =
206 disp_bc->data.flag4 || disp_bc->data.flag5 || disp_bc->data.flag6;
208 auto ents = bc.second->bcEnts;
210 std::array<SmartPetscObj<IS>, 3> is_xyz;
211 auto prb_name = fe_method_ptr->problemPtr->getName();
213 if (disp_bc->data.flag1 || is_rotation) {
214 CHKERR is_mng->isCreateProblemFieldAndRankLocal(
217 if (disp_bc->data.flag2 || is_rotation) {
218 CHKERR is_mng->isCreateProblemFieldAndRankLocal(
221 if (disp_bc->data.flag3 || is_rotation) {
222 CHKERR is_mng->isCreateProblemFieldAndRankLocal(
226 auto get_is_sum = [](
auto is1,
auto is2) {
232 for (
auto &is : is_xyz) {
237 is_sum = get_is_sum(is_sum, is);
246 if (
auto fe_ptr = fePtr.lock()) {
248 auto snes_ctx = fe_ptr->snes_ctx;
249 auto ts_ctx = fe_ptr->ts_ctx;
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);
260 *fe_ptr->vecAssembleSwitch =
false;
264 const int *index_ptr;
265 CHKERR ISGetIndices(is_sum, &index_ptr);
267 CHKERR ISGetLocalSize(is_sum, &size);
272 CHKERR vec_mng->setLocalGhostVector(problem_name,
ROW, tmp_x,
273 INSERT_VALUES, SCATTER_FORWARD);
275 CHKERR VecGetArrayRead(tmp_x, &u);
281 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
282 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
287 for (
auto i = 0;
i != size; ++
i) {
288 a[index_ptr[
i]] = vDiag * (
v[index_ptr[
i]] - u[index_ptr[
i]]);
291 CHKERR VecRestoreArrayRead(x, &
v);
294 for (
auto i = 0;
i != size; ++
i) {
295 a[index_ptr[
i]] = vDiag * u[index_ptr[
i]];
299 CHKERR VecRestoreArrayRead(tmp_x, &u);
301 CHKERR ISRestoreIndices(is_sum, &index_ptr);
307 "Can not lock shared pointer");
316 : mField(m_field), fePtr(fe_ptr), vDiag(diag), vLhs(lhs), vAO(ao) {}
322 if (
auto fe_method_ptr = fePtr.lock()) {
324 auto bc_mng = mField.getInterface<
BcManager>();
327 const auto problem_name = fe_method_ptr->problemPtr->getName();
331 for (
auto bc : bc_mng->getBcMapByBlockName()) {
332 if (
auto disp_bc = bc.second->dispBcPtr) {
334 auto &bc_id = bc.first;
336 auto regex_str = (boost::format(
"%s_(.*)") % problem_name).str();
337 if (std::regex_match(bc_id, std::regex(regex_str))) {
343 <<
"Apply EssentialPreProc<DisplacementCubitBcData>: "
344 << problem_name <<
"_" <<
field_name <<
"_" << block_name;
346 const bool is_rotation =
347 disp_bc->data.flag4 || disp_bc->data.flag5 || disp_bc->data.flag6;
349 auto ents = bc.second->bcEnts;
351 std::array<SmartPetscObj<IS>, 3> is_xyz;
352 auto prb_name = fe_method_ptr->problemPtr->getName();
354 if (disp_bc->data.flag1 || is_rotation) {
355 CHKERR is_mng->isCreateProblemFieldAndRank(
358 if (disp_bc->data.flag2 || is_rotation) {
359 CHKERR is_mng->isCreateProblemFieldAndRank(
362 if (disp_bc->data.flag3 || is_rotation) {
363 CHKERR is_mng->isCreateProblemFieldAndRank(
367 auto get_is_sum = [](
auto is1,
auto is2) {
373 for (
auto &is : is_xyz) {
378 is_sum = get_is_sum(is_sum, is);
387 if (
auto fe_ptr = fePtr.lock()) {
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;
399 MOFEM_LOG(
"WORLD", Sev::noisy) <<
"Apply AO to IS";
400 CHKERR AOApplicationToPetscIS(vAO, is_sum);
403 CHKERR MatZeroRowsColumnsIS(
B, is_sum, vDiag, PETSC_NULL, PETSC_NULL);
409 "Can not lock shared pointer");
418 : mField(m_field), fePtr(fe_ptr), vRhs(rhs), sevLevel(sev) {}
424 enum { X = 0, Y, Z, MX, MY, MZ, LAST };
426 if (
auto fe_ptr = fePtr.lock()) {
430 if (fe_ptr->vecAssembleSwitch) {
431 if ((*fe_ptr->vecAssembleSwitch) && !vRhs) {
432 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
433 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
434 CHKERR VecAssemblyBegin(f);
436 *fe_ptr->vecAssembleSwitch =
false;
440 auto get_low_hi_uid_by_entities = [](
auto bit_number,
auto f,
auto s) {
445 auto get_low_hi = [fe_ptr](
auto lo_uid,
auto hi_uid) {
446 auto it = fe_ptr->problemPtr->numeredRowDofsPtr->get<
Unique_mi_tag>()
447 .lower_bound(lo_uid);
448 auto hi_it = fe_ptr->problemPtr->numeredRowDofsPtr->get<
Unique_mi_tag>()
449 .upper_bound(hi_uid);
450 return std::make_pair(it, hi_it);
453 auto mpi_array_reduce = [
this](
auto &array) {
454 std::array<double, LAST> array_sum{0, 0, 0, 0, 0, 0};
455 MPI_Allreduce(&array[0], &array_sum[0], LAST, MPI_DOUBLE, MPI_SUM,
463 auto bc_mng = mField.getInterface<
BcManager>();
464 const auto problem_name = fe_ptr->problemPtr->getName();
465 const auto nb_local_dofs = fe_ptr->problemPtr->nbLocDofsRow;
467 std::array<double, LAST> total_reactions{0, 0, 0, 0, 0, 0};
469 for (
auto bc : bc_mng->getBcMapByBlockName()) {
470 if (
auto disp_bc = bc.second->dispBcPtr) {
472 auto &bc_id = bc.first;
474 auto regex_str = (boost::format(
"%s_(.*)") % problem_name).str();
475 if (std::regex_match(bc_id, std::regex(regex_str))) {
481 <<
"EssentialPreProc<DisplacementCubitBcData>: " << problem_name
483 auto bit_number = mField.get_field_bit_number(
field_name);
486 if (
auto ext_disp_bc =
489 for (
int a = 0;
a != 3; ++
a)
490 t_off(
a) = ext_disp_bc->rotOffset[
a];
493 auto verts = bc.second->bcEnts.subset_by_type(MBVERTEX);
499 auto get_coords_vec = [&]() {
502 CHKERR mField.get_moab().get_coords(verts, &*coords_vec.begin());
507 auto coords_vec = get_coords_vec();
508 std::array<double, LAST> reactions{0, 0, 0, 0, 0, 0};
510 for (
auto pit = verts.const_pair_begin();
511 pit != verts.const_pair_end(); ++pit) {
512 auto [lo_uid, hi_uid] =
513 get_low_hi_uid_by_entities(bit_number, pit->first, pit->second);
514 auto [lo, hi] = get_low_hi(lo_uid, hi_uid);
516 for (; lo != hi; ++lo) {
517 const auto loc_dof = (*lo)->getPetscLocalDofIdx();
518 if (loc_dof < nb_local_dofs) {
520 const auto coeff = (*lo)->getDofCoeffIdx();
524 ((disp_bc->data.flag1 || disp_bc->data.flag4) &&
526 ((disp_bc->data.flag2 || disp_bc->data.flag5) &&
528 ((disp_bc->data.flag3 || disp_bc->data.flag6) &&
531 const auto ent = (*lo)->getEnt();
532 reactions[coeff] +=
a[loc_dof];
536 t_force(coeff) =
a[loc_dof];
541 const auto idx = verts.index(ent);
543 coords_vec[3 * idx + X], coords_vec[3 * idx + Y],
544 coords_vec[3 * idx + Z]};
545 t_coords(
i) -= t_off(
i);
549 auto moment = [&](
auto &&t_force,
auto &&t_coords) {
552 (FTensor::levi_civita<double>(
i,
j,
k) * t_coords(
k)) *
557 auto t_moment = moment(force(), coord());
558 reactions[MX] += t_moment(X);
559 reactions[MY] += t_moment(Y);
560 reactions[MZ] += t_moment(Z);
571 &total_reactions[X], &total_reactions[Y], &total_reactions[Z]};
573 &total_reactions[MX], &total_reactions[MY], &total_reactions[MZ]};
574 t_total_force(
i) += t_force(
i);
577 (FTensor::levi_civita<double>(
i,
j,
k) * t_off(
k)) * t_force(
j);
579 auto mpi_reactions = mpi_array_reduce(reactions);
581 "Offset: %4.3e %4.3e %4.3e", t_off(X), t_off(Y),
584 "Force: %4.3e %4.3e %4.3e", mpi_reactions[X],
585 mpi_reactions[Y], mpi_reactions[Z]);
587 "Moment: %4.3e %4.3e %4.3e", mpi_reactions[MX],
588 mpi_reactions[MY], mpi_reactions[MZ]);
595 CHKERR VecRestoreArrayRead(f, &
a);
597 auto mpi_total_reactions = mpi_array_reduce(total_reactions);
599 "WORLD", sevLevel,
"Essential",
"Total force: %4.3e %4.3e %4.3e",
600 mpi_total_reactions[X], mpi_total_reactions[Y], mpi_total_reactions[Z]);
602 "Total moment: %4.3e %4.3e %4.3e",
603 mpi_total_reactions[MX], mpi_total_reactions[MY],
604 mpi_total_reactions[MZ]);
608 "Can not lock shared pointer");
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
constexpr auto field_name
Simple interface for fast problem set-up.
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.
Deprecated interface functions.
A specialized version of DisplacementCubitBcData that includes an additional rotation offset.
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.
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Class (Function) to enforce essential constrains on the left hand side diagonal.
Class (Function) to enforce essential constrains on the right hand side diagonal.
Class (Function) to enforce essential constrains.
Class (Function) to calculate residual side diagonal.
Section manager is used to create indexes and sections.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.