Specialization for DisplacementCubitBcData. More...

#include <src/boundary_conditions/EssentialMPCsData.hpp>

Collaboration diagram for MoFEM::EssentialPostProcRhs< MPCsType >:

Public Member Functions
	EssentialPostProcRhs (MoFEM::Interface &m_field, boost::shared_ptr< FEMethod > fe_ptr, double diag=1, SmartPetscObj< Vec > rhs=nullptr)

MoFEMErrorCode	operator() ()

Protected Attributes
MoFEM::Interface &	mField

boost::weak_ptr< FEMethod >	fePtr

double	vDiag

SmartPetscObj< Vec >	vRhs

Detailed Description

Specialization for DisplacementCubitBcData.

Template Parameters

Examples: nonlinear_elastic.cpp.

Definition at line 80 of file EssentialMPCsData.hpp.

Constructor & Destructor Documentation

◆ EssentialPostProcRhs()

MoFEM::EssentialPostProcRhs< MPCsType >::EssentialPostProcRhs	(	MoFEM::Interface &	m_field,
		boost::shared_ptr< FEMethod >	fe_ptr,
		double	diag = `1`,
		SmartPetscObj< Vec >	rhs = `nullptr`
	)

Definition at line 531 of file EssentialMPCsData.cpp.

534 : mField(m_field), fePtr(fe_ptr), vDiag(diag), vRhs(rhs) {}

Member Function Documentation

◆ operator()()

MoFEMErrorCode MoFEM::EssentialPostProcRhs< MPCsType >::operator() ( )

Definition at line 536 of file EssentialMPCsData.cpp.

                                                           {
   MOFEM_LOG_CHANNEL("WORLD");
   MoFEMFunctionBegin;
  
   if (auto fe_method_ptr = fePtr.lock()) {
  
     auto bc_mng = mField.getInterface<BcManager>();
     auto is_mng = mField.getInterface<ISManager>();
     auto vec_mng = mField.getInterface<VecManager>();
     const auto problem_name = fe_method_ptr->problemPtr->getName();
     // get connectivity from edge and assign master, slave ents
     for (auto bc : bc_mng->getBcMapByBlockName()) {
       if (auto mpc_bc = bc.second->mpcPtr) {
  
         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);
  
           auto get_field_coeffs = [&](auto field_name) {
             auto field_ptr = mField.get_field_structure(field_name);
             return field_ptr->getNbOfCoeffs();
           };
  
           const auto nb_field_coeffs = get_field_coeffs(field_name);
           constexpr auto max_nb_dofs_per_node = 6;
  
           if (nb_field_coeffs > max_nb_dofs_per_node)
             MOFEM_LOG("WORLD", Sev::error)
                 << "MultiPointConstraints PreProcLhs<MPCsType>: support only "
                    "up to 6 dofs per node for now.";
  
           MOFEM_LOG("WORLD", Sev::noisy)
               << "Apply MultiPointConstraints PreProc<MPCsType>: "
               << problem_name << "_" << field_name << "_" << block_name;
  
           auto mpc_type = mpc_bc->mpcType;
           switch (mpc_type) {
           case MPC::COUPLING:
           case MPC::TIE:
           case MPC::RIGID_BODY:
             break;
           default:
             SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                     "MPC type not implemented");
           }
  
           auto master_verts = mpc_bc->mpcMasterEnts.subset_by_type(MBVERTEX);
           auto slave_verts = mpc_bc->mpcSlaveEnts.subset_by_type(MBVERTEX);
  
           auto prb_name = fe_method_ptr->problemPtr->getName();
  
           auto get_flag = [&](int idx) { return (&mpc_bc->data.flag1)[idx]; };
  
           auto add_is = [](auto is1, auto is2) {
             IS is;
             CHK_THROW_MESSAGE(ISExpand(is1, is2, &is), "is sum");
             return SmartPetscObj<IS>(is);
           };
  
           SmartPetscObj<IS> is_m[max_nb_dofs_per_node];
           SmartPetscObj<IS> is_s[max_nb_dofs_per_node];
  
           for (int dd = 0; dd != nb_field_coeffs; dd++) {
             if (get_flag(dd)) {
         
               CHKERR is_mng->isCreateProblemFieldAndRankLocal(
                   prb_name, ROW, field_name, dd, dd, is_m[dd], &master_verts);
               CHKERR is_mng->isCreateProblemFieldAndRankLocal(
                   prb_name, ROW, field_name, dd, dd, is_s[dd], &slave_verts);
   
             }
           }
  
           {
  
             if (auto fe_ptr = fePtr.lock()) {
               auto snes_ctx = fe_ptr->snes_ctx;
               auto ts_ctx = fe_ptr->ts_ctx;
               const bool is_nonlinear = snes_ctx != FEMethod::CTX_SNESNONE ||
                                         ts_ctx != FEMethod::CTX_TSNONE;
               // is_nonlinear = is_nonlinear || mpc_bc->isReprocitical;
               const int contrb = mpc_bc->isReprocitical ? 1 : 0;
               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 vec_set_values = [&](auto is_xyz_m, auto is_xyz_s, double val) {
                 MoFEMFunctionBeginHot;
                 const int *m_index_ptr;
                 CHKERR ISGetIndices(is_xyz_m, &m_index_ptr);
                 const int *s_index_ptr;
                 CHKERR ISGetIndices(is_xyz_s, &s_index_ptr);
                 int size_m;
                 CHKERR ISGetLocalSize(is_xyz_m, &size_m);
                 int size_s;
                 CHKERR ISGetLocalSize(is_xyz_s, &size_s);
  
                 double *f;
                 CHKERR VecGetArray(F, &f);
                 // if (size_m > size_s)
                 //   MOFEM_LOG("WORLD", Sev::error)
                 //       << "Size of master IS > Size of slave IS : " << size_m
                 //       << " > " << size_s;
                 // if (size_m == 0)
                 //   MOFEM_LOG("WORLD", Sev::error)
                 //       << "Size of master IS is " << size_m;
                 
                   if (is_nonlinear) {
                     auto x = fe_ptr->x;
                     auto tmp_x = vectorDuplicate(F);
  
                     CHKERR vec_mng->setLocalGhostVector(problem_name, ROW,
                                                         tmp_x, INSERT_VALUES,
                                                         SCATTER_FORWARD);
                     const double *v;
                     const double *u;
  
                     CHKERR VecGetArrayRead(tmp_x, &u);
                     CHKERR VecGetArrayRead(x, &v);
                     if (size_m && size_s) {
                       for (auto i = 0; i != size_s; ++i) {
                           auto m_idx = size_m == 1 ? 0 : i;
                           f[s_index_ptr[i]] = val * (v[s_index_ptr[i]] - v[m_index_ptr[m_idx]]) + f[s_index_ptr[i]] * contrb ;
                       }
                     }
                     CHKERR VecRestoreArrayRead(x, &v);
                     CHKERR VecRestoreArrayRead(tmp_x, &u);
                   } else {
                     if (size_m && size_s) 
                     for (auto i = 0; i != size_s; ++i) {
                       f[s_index_ptr[i]] = 0;
                     }
                   }
                 
                 CHKERR VecRestoreArray(F, &f);
                 CHKERR ISRestoreIndices(is_xyz_m, &m_index_ptr);
                 CHKERR ISRestoreIndices(is_xyz_s, &s_index_ptr);
  
                 MoFEMFunctionReturnHot(0);
               };
  
               for (int dd = 0; dd != nb_field_coeffs; dd++)
                 if (get_flag(dd)) {
  
                   if (!mpc_bc->isReprocitical) {
                     CHKERR vec_set_values(is_m[dd], is_s[dd], vDiag);
                   } else {
                     auto &pn = mpc_bc->mPenalty;
                     CHKERR vec_set_values(is_m[dd], is_s[dd], pn);
                     CHKERR vec_set_values(is_s[dd], is_m[dd], pn);
                   }
                 }
             };
  
             // User is responsible for assembly if vLhs is provided
  
             // ISView(is_xyz_m, PETSC_VIEWER_STDOUT_WORLD);
             // CHKERR MatZeroRowsColumnsIS(B, is_xyz_m, vDiag, PETSC_NULL,
             //                             PETSC_NULL);
           } 
         }
       }
     }
   } else {
     SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
             "Can not lock shared pointer");
   }
  
   MoFEMFunctionReturn(0);
 }