#include <src/boundary_conditions/EssentialMPCsData.hpp>

Collaboration diagram for MoFEM::EssentialPostProcLhs< MPCsType >:

Public Member Functions
	EssentialPostProcLhs (MoFEM::Interface &m_field, boost::shared_ptr< FEMethod > fe_ptr, double diag=1, SmartPetscObj< Mat > lhs=nullptr, SmartPetscObj< AO > ao=nullptr)

MoFEMErrorCode	operator() ()

Protected Attributes
MoFEM::Interface &	mField

boost::weak_ptr< FEMethod >	fePtr

double	vDiag

SmartPetscObj< Mat >	vLhs

SmartPetscObj< AO >	vAO

Detailed Description

Specialization for MPCsType.

Template Parameters

Examples: nonlinear_elastic.cpp.

Definition at line 99 of file EssentialMPCsData.hpp.

Constructor & Destructor Documentation

◆ EssentialPostProcLhs()

MoFEM::EssentialPostProcLhs< MPCsType >::EssentialPostProcLhs	(	MoFEM::Interface &	m_field,
		boost::shared_ptr< FEMethod >	fe_ptr,
		double	diag = `1`,
		SmartPetscObj< Mat >	lhs = `nullptr`,
		SmartPetscObj< AO >	ao = `nullptr`
	)

Definition at line 328 of file EssentialMPCsData.cpp.

331 : mField(m_field), fePtr(fe_ptr), vDiag(diag), vLhs(lhs), vAO(ao) {}

Member Function Documentation

◆ operator()()

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

Definition at line 333 of file EssentialMPCsData.cpp.

                                                           {
   MOFEM_LOG_CHANNEL("WORLD");
   MoFEMFunctionBegin;
  
   if (auto fe_ptr = fePtr.lock()) {
  
     auto bc_mng = mField.getInterface<BcManager>();
     auto is_mng = mField.getInterface<ISManager>();
     const auto problem_name = fe_ptr->problemPtr->getName();
     bool do_assembly = false;
     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.";
           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");
           }
  
           // FIXME: there handle the vertices differently (block level)
           auto master_verts = mpc_bc->mpcMasterEnts.subset_by_type(MBVERTEX);
           auto slave_verts = mpc_bc->mpcSlaveEnts.subset_by_type(MBVERTEX);
  
           auto prb_name = fe_ptr->problemPtr->getName();
           auto get_flag = [&](int idx) { return (&mpc_bc->data.flag1)[idx]; };
           MOFEM_LOG("WORLD", Sev::verbose)
               << "Size of master nodes: " << mpc_bc->mpcMasterEnts.size()
               << " Size of slave nodes : " << mpc_bc->mpcSlaveEnts.size();
  
           if (mpc_bc->mpcMasterEnts.size() > mpc_bc->mpcSlaveEnts.size()) {
             MOFEM_LOG("WORLD", Sev::error)
                 << "Size of master nodes < Size of slave nodes : "
                 << mpc_bc->mpcMasterEnts.size() << " > " << mpc_bc->mpcSlaveEnts.size();
             // SETERRQ(PETSC_COMM_WORLD, MOFEM_OPERATION_UNSUCCESSFUL,
             //         "data inconsistency");
           }
  
           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_xyz_row_sum;
           SmartPetscObj<IS> is_m_row[max_nb_dofs_per_node];
           SmartPetscObj<IS> is_m_col[max_nb_dofs_per_node];
           SmartPetscObj<IS> is_s_row[max_nb_dofs_per_node];
           SmartPetscObj<IS> is_s_col[max_nb_dofs_per_node];
  
           for (int dd = 0; dd != nb_field_coeffs; dd++) {
             if (get_flag(dd)) {
               // SmartPetscObj<IS> is_xyz_m;
               CHKERR is_mng->isCreateProblemFieldAndRank(
                   prb_name, ROW, field_name, dd, dd, is_m_row[dd],
                   &master_verts);
               CHKERR is_mng->isCreateProblemFieldAndRank(
                   prb_name, COL, field_name, dd, dd, is_m_col[dd],
                   &master_verts);
               CHKERR is_mng->isCreateProblemFieldAndRank(
                   prb_name, COL, field_name, dd, dd, is_s_col[dd],
                   &slave_verts);
               CHKERR is_mng->isCreateProblemFieldAndRank(
                   prb_name, ROW, field_name, dd, dd, is_s_row[dd],
                   &slave_verts);
               // ISView(is_s_row[dd], PETSC_VIEWER_STDOUT_WORLD);
               // ISView(is_s_col[dd], PETSC_VIEWER_STDOUT_WORLD);
  
               if (!mpc_bc->isReprocitical) {
                 if (!is_xyz_row_sum) {
                   is_xyz_row_sum = is_s_row[dd];
                 } else {
                   is_xyz_row_sum = add_is(is_xyz_row_sum, is_s_row[dd]);
                 }
               }
             }
           }
  
           // if (is_xyz_row_sum) {
             SmartPetscObj<Mat> B =
                 vLhs ? vLhs : SmartPetscObj<Mat>(fe_ptr->B, true);
             // The user is responsible for assembly if vLhs is provided
             MatType typem;
             CHKERR MatGetType(B, &typem);
             CHKERR MatSetOption(B, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE);
             if (typem == MATMPIBAIJ)
               CHKERR MatSetOption(B, MAT_USE_HASH_TABLE, PETSC_TRUE);
  
             if ((*fe_ptr->matAssembleSwitch) && !vLhs) {
               if (*fe_ptr->matAssembleSwitch) {
                 CHKERR MatAssemblyBegin(B, MAT_FLUSH_ASSEMBLY);
                 CHKERR MatAssemblyEnd(B, MAT_FLUSH_ASSEMBLY);
                 *fe_ptr->matAssembleSwitch = false;
               }
             }
  
             if (vAO) {
               MOFEM_LOG("WORLD", Sev::error) << "No support for AO yet";
             }
  
             CHKERR MatZeroRowsIS(B, is_xyz_row_sum, 0, PETSC_NULL,
                                  PETSC_NULL);
           // }
             auto set_mat_values = [&](auto row_is, auto col_is, double val, double perturb = 0) {
               MoFEMFunctionBeginHot;
               const int *row_index_ptr;
               CHKERR ISGetIndices(row_is, &row_index_ptr);
               const int *col_index_ptr;
               CHKERR ISGetIndices(col_is, &col_index_ptr);
               int row_size, col_size;
               CHKERR ISGetLocalSize(row_is, &row_size);
               CHKERR ISGetLocalSize(col_is, &col_size);
  
               MatrixDouble locMat(row_size, col_size);
               fill(locMat.data().begin(), locMat.data().end(), 0.0);
               for (int i = 0; i != row_size; i++)
                 for (int j = 0; j != col_size; j++)
                   if (i == j || col_size == 1)
                     locMat(i, j) = val;
  
               if (mpc_bc->isReprocitical) {
                 locMat *= perturb;
                 CHKERR ::MatSetValues(B, row_size, row_index_ptr, col_size,
                                       col_index_ptr, &*locMat.data().begin(),
                                       ADD_VALUES);
               } else {
                 CHKERR ::MatSetValues(B, row_size, row_index_ptr, col_size,
                                       col_index_ptr, &*locMat.data().begin(),
                                       INSERT_VALUES);
               }
  
               CHKERR ISRestoreIndices(row_is, &row_index_ptr);
               CHKERR ISRestoreIndices(col_is, &col_index_ptr);
  
               MoFEMFunctionReturnHot(0);
             };
  
             for (int dd = 0; dd != nb_field_coeffs; dd++) {
               if (get_flag(dd)) {
                 do_assembly = true;
                 if (!mpc_bc->isReprocitical) {
                   CHKERR set_mat_values(is_s_row[dd], is_s_col[dd], vDiag);
                   CHKERR set_mat_values(is_s_row[dd], is_m_col[dd], -vDiag);
                 } 
                 else
                 {
                   auto &pn = mpc_bc->mPenalty;
                   CHKERR set_mat_values(is_s_row[dd], is_s_col[dd], vDiag, pn);
                   CHKERR set_mat_values(is_s_row[dd], is_m_col[dd], -vDiag, pn);
                   CHKERR set_mat_values(is_m_row[dd], is_m_col[dd], vDiag, pn);
                   CHKERR set_mat_values(is_m_row[dd], is_s_col[dd], -vDiag, pn);
                 }
               }
             }
         } // if regex
       } // mpc loop
     }   // bc loop
     if (do_assembly) {
       SmartPetscObj<Mat> B = vLhs ? vLhs : SmartPetscObj<Mat>(fe_ptr->B, true);
       CHKERR MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
       CHKERR MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
     }
  
   } else {
     SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
             "Can not lock shared pointer");
   } // if fe_ptr
   MoFEMFunctionReturn(0);
 }