Specialization for DisplacementCubitBcData. More...

#include <src/boundary_conditions/EssentialDisplacementCubitBcData.hpp>

Collaboration diagram for MoFEM::EssentialPreProcReaction< DisplacementCubitBcData >:

Public Member Functions
	EssentialPreProcReaction (MoFEM::Interface &m_field, boost::shared_ptr< FEMethod > fe_ptr, SmartPetscObj< Vec > rhs=nullptr, LogManager::SeverityLevel sev=Sev::inform, boost::shared_ptr< std::vector< double >> reaction_ptr=nullptr)

MoFEMErrorCode	operator() ()

Protected Attributes
MoFEM::Interface &	mField

boost::weak_ptr< FEMethod >	fePtr

SmartPetscObj< Vec >	vRhs

LogManager::SeverityLevel	sevLevel

PetscBool	printBlockName

boost::shared_ptr< std::vector< double > >	reactionPtr

char	reactionBlockName [255]

Detailed Description

Specialization for DisplacementCubitBcData.

Template Parameters

Examples: adolc_plasticity.cpp, plastic.cpp, and seepage.cpp.

Definition at line 157 of file EssentialDisplacementCubitBcData.hpp.

Constructor & Destructor Documentation

◆ EssentialPreProcReaction()

MoFEM::EssentialPreProcReaction< DisplacementCubitBcData >::EssentialPreProcReaction	(	MoFEM::Interface &	m_field,
		boost::shared_ptr< FEMethod >	fe_ptr,
		SmartPetscObj< Vec >	rhs = `nullptr`,
		LogManager::SeverityLevel	sev = `Sev::inform`,
		boost::shared_ptr< std::vector< double >>	reaction_ptr = `nullptr`
	)

Definition at line 415 of file EssentialDisplacementCubitBcData.cpp.

     : mField(m_field), fePtr(fe_ptr), vRhs(rhs), sevLevel(sev),
       printBlockName(PETSC_FALSE), reactionPtr(reaction_ptr),
       reactionBlockName("") {
  
   CHK_THROW_MESSAGE(PetscOptionsGetBool(PETSC_NULL, PETSC_NULL,
                                         "-reaction_print_block_name",
                                         &printBlockName, PETSC_NULL),
                     "can not get option");
   CHKERR PetscOptionsGetString(PETSC_NULL, "-reaction_block_name",
                                reactionBlockName, 255, PETSC_NULL);
 }

Member Function Documentation

◆ operator()()

MoFEMErrorCode MoFEM::EssentialPreProcReaction< DisplacementCubitBcData >::operator() ( )

Definition at line 431 of file EssentialDisplacementCubitBcData.cpp.

                                                                              {
   MOFEM_LOG_CHANNEL("WORLD");
   MoFEMFunctionBegin;
  
   enum { X = 0, Y, Z, MX, MY, MZ, LAST };
  
   if (auto fe_ptr = fePtr.lock()) {
  
     SmartPetscObj<Vec> f = vRhs ? vRhs : SmartPetscObj<Vec>(fe_ptr->f, true);
  
     if (fe_ptr->vecAssembleSwitch) {
       if ((*fe_ptr->vecAssembleSwitch) && !vRhs) {
         CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
         CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
         CHKERR VecAssemblyBegin(f);
         CHKERR VecAssemblyEnd(f);
         *fe_ptr->vecAssembleSwitch = false;
       }
     }
  
     auto get_low_hi_uid_by_entities = [](auto bit_number, auto f, auto s) {
       return std::make_pair(DofEntity::getLoFieldEntityUId(bit_number, f),
                             DofEntity::getHiFieldEntityUId(bit_number, s));
     };
  
     auto get_low_hi = [fe_ptr](auto lo_uid, auto hi_uid) {
       auto it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
                     .lower_bound(lo_uid);
       auto hi_it = fe_ptr->problemPtr->numeredRowDofsPtr->get<Unique_mi_tag>()
                        .upper_bound(hi_uid);
       return std::make_pair(it, hi_it);
     };
  
     auto mpi_array_reduce = [this](auto &array) {
       std::array<double, LAST> array_sum{0, 0, 0, 0, 0, 0};
       MPI_Allreduce(&array[0], &array_sum[0], LAST, MPI_DOUBLE, MPI_SUM,
                     mField.get_comm());
       return array_sum;
     };
  
     const double *a;
     CHKERR VecGetArrayRead(f, &a);
  
     auto bc_mng = mField.getInterface<BcManager>();
     const auto problem_name = fe_ptr->problemPtr->getName();
     const auto nb_local_dofs = fe_ptr->problemPtr->nbLocDofsRow;
  
     std::array<double, LAST> total_reactions{0, 0, 0, 0, 0, 0};
  
     for (auto bc : bc_mng->getBcMapByBlockName()) {
       if (auto disp_bc = bc.second->dispBcPtr) {
  
         auto &bc_id = bc.first;
  
         auto regex_str = (boost::format("%s_(.*)") % problem_name).str();
         if (std::regex_match(bc_id, std::regex(regex_str))) {
  
           auto [field_name, block_name] =
               BcManager::extractStringFromBlockId(bc_id, problem_name);
  
           MOFEM_TAG_AND_LOG("WORLD", sevLevel, "Essential")
               << "EssentialPreProc<DisplacementCubitBcData>: " << problem_name
               << "_" << field_name << "_" << block_name;
           auto bit_number = mField.get_field_bit_number(field_name);
  
           FTensor::Tensor1<double, 3> t_off{0., 0., 0.};
           if (auto ext_disp_bc =
                   dynamic_cast<DisplacementCubitBcDataWithRotation const *>(
                       disp_bc.get())) {
             for (int a = 0; a != 3; ++a)
               t_off(a) = ext_disp_bc->rotOffset[a];
           }
  
           auto verts = bc.second->bcEnts.subset_by_type(MBVERTEX);
  
           FTensor::Index<'i', 3> i;
           FTensor::Index<'j', 3> j;
           FTensor::Index<'k', 3> k;
  
           auto get_coords_vec = [&]() {
             VectorDouble coords_vec(3 * verts.size());
             if (verts.size()) {
               CHKERR mField.get_moab().get_coords(verts, &*coords_vec.begin());
             }
             return coords_vec;
           };
  
           auto coords_vec = get_coords_vec();
           std::array<double, LAST> reactions{0, 0, 0, 0, 0, 0};
  
           for (auto pit = verts.const_pair_begin();
                pit != verts.const_pair_end(); ++pit) {
             auto [lo_uid, hi_uid] =
                 get_low_hi_uid_by_entities(bit_number, pit->first, pit->second);
             auto [lo, hi] = get_low_hi(lo_uid, hi_uid);
  
             for (; lo != hi; ++lo) {
               const auto loc_dof = (*lo)->getPetscLocalDofIdx();
               if (loc_dof < nb_local_dofs) {
  
                 const auto coeff = (*lo)->getDofCoeffIdx();
  
                 if (
  
                     ((disp_bc->data.flag1 || disp_bc->data.flag4) &&
                      coeff == 0) ||
                     ((disp_bc->data.flag2 || disp_bc->data.flag5) &&
                      coeff == 1) ||
                     ((disp_bc->data.flag3 || disp_bc->data.flag6) &&
                      coeff == 2)) {
  
                   const auto ent = (*lo)->getEnt();
                   reactions[coeff] += a[loc_dof];
  
                   auto force = [&]() {
                     FTensor::Tensor1<double, 3> t_force{0., 0., 0.};
                     t_force(coeff) = a[loc_dof];
                     return t_force;
                   };
  
                   auto coord = [&]() {
                     const auto idx = verts.index(ent);
                     FTensor::Tensor1<double, 3> t_coords{
                         coords_vec[3 * idx + X], coords_vec[3 * idx + Y],
                         coords_vec[3 * idx + Z]};
                     t_coords(i) -= t_off(i);
                     return t_coords;
                   };
  
                   auto moment = [&](auto &&t_force, auto &&t_coords) {
                     FTensor::Tensor1<double, 3> t_moment;
                     t_moment(i) =
                         (FTensor::levi_civita<double>(i, j, k) * t_coords(k)) *
                         t_force(j);
                     return t_moment;
                   };
  
                   auto t_moment = moment(force(), coord());
                   reactions[MX] += t_moment(X);
                   reactions[MY] += t_moment(Y);
                   reactions[MZ] += t_moment(Z);
                 }
               }
             }
           }
  
           FTensor::Tensor1<double, 3> t_force{reactions[X], reactions[Y],
                                               reactions[Z]};
           FTensor::Tensor1<double, 3> t_moment{reactions[MX], reactions[MY],
                                                reactions[MZ]};
           FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_total_force{
               &total_reactions[X], &total_reactions[Y], &total_reactions[Z]};
           FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_total_moment{
               &total_reactions[MX], &total_reactions[MY], &total_reactions[MZ]};
           t_total_force(i) += t_force(i);
           t_total_moment(i) +=
               t_moment(i) +
               (FTensor::levi_civita<double>(i, j, k) * t_off(k)) * t_force(j);
  
           auto mpi_reactions = mpi_array_reduce(reactions);
           if (printBlockName) {
             MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
                                 "Block %s Offset: %6.4e %6.4e %6.4e",
                                 block_name.c_str(), t_off(X), t_off(Y),
                                 t_off(Z));
             MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
                                 "Block %s Force: %6.4e %6.4e %6.4e",
                                 block_name.c_str(), mpi_reactions[X],
                                 mpi_reactions[Y], mpi_reactions[Z]);
             MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
                                 "Block %s Moment: %6.4e %6.4e %6.4e",
                                 block_name.c_str(), mpi_reactions[MX],
                                 mpi_reactions[MY], mpi_reactions[MZ]);
             if ((reactionPtr != nullptr) &&
                 (block_name.compare(reactionBlockName) == 0)) {
               (*reactionPtr)[X] = mpi_reactions[X];
               (*reactionPtr)[Y] = mpi_reactions[Y];
               (*reactionPtr)[Z] = mpi_reactions[Z];
             }
           } else {
             MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
                                 "Offset: %6.4e %6.4e %6.4e", t_off(X), t_off(Y),
                                 t_off(Z));
             MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
                                 "Force: %6.4e %6.4e %6.4e", mpi_reactions[X],
                                 mpi_reactions[Y], mpi_reactions[Z]);
             MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
                                 "Moment: %6.4e %6.4e %6.4e", mpi_reactions[MX],
                                 mpi_reactions[MY], mpi_reactions[MZ]);
           }
         }
       }
     }
  
     CHKERR VecRestoreArrayRead(f, &a);
  
     auto mpi_total_reactions = mpi_array_reduce(total_reactions);
     MOFEM_TAG_AND_LOG_C(
         "WORLD", sevLevel, "Essential", "Total force: %6.4e %6.4e %6.4e",
         mpi_total_reactions[X], mpi_total_reactions[Y], mpi_total_reactions[Z]);
     MOFEM_TAG_AND_LOG_C("WORLD", sevLevel, "Essential",
                         "Total moment: %6.4e %6.4e %6.4e",
                         mpi_total_reactions[MX], mpi_total_reactions[MY],
                         mpi_total_reactions[MZ]);
  
   } else {
     SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
             "Can not lock shared pointer");
   }
  
   MoFEMFunctionReturn(0);
 }