v0.14.0
Public Member Functions | Protected Attributes | List of all members
MoFEM::EssentialPostProcLhs< MPCsType > Struct Reference

Specialization for MPCsType. More...

#include <src/boundary_conditions/EssentialMPCsData.hpp>

Collaboration diagram for MoFEM::EssentialPostProcLhs< MPCsType >:
[legend]

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)
 
virtual ~EssentialPostProcLhs ()=default
 
MoFEMErrorCode operator() ()
 

Protected Attributes

MoFEM::InterfacemField
 
boost::weak_ptr< FEMethodfePtr
 
double vDiag
 
SmartPetscObj< Mat > vLhs
 
SmartPetscObj< AO > vAO
 

Detailed Description

Specialization for MPCsType.

Template Parameters
Examples
nonlinear_elastic.cpp.

Definition at line 102 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) {}

◆ ~EssentialPostProcLhs()

Member Function Documentation

◆ operator()()

Definition at line 333 of file EssentialMPCsData.cpp.

333  {
334  MOFEM_LOG_CHANNEL("WORLD");
336 
337  if (auto fe_ptr = fePtr.lock()) {
338 
339  auto bc_mng = mField.getInterface<BcManager>();
340  auto is_mng = mField.getInterface<ISManager>();
341  const auto problem_name = fe_ptr->problemPtr->getName();
342  bool do_assembly = false;
343  for (auto bc : bc_mng->getBcMapByBlockName()) {
344  if (auto mpc_bc = bc.second->mpcPtr) {
345 
346  auto &bc_id = bc.first;
347 
348  auto regex_str = (boost::format("%s_(.*)") % problem_name).str();
349  if (std::regex_match(bc_id, std::regex(regex_str))) {
350 
351  auto [field_name, block_name] =
352  BcManager::extractStringFromBlockId(bc_id, problem_name);
353 
354  auto get_field_coeffs = [&](auto field_name) {
355  auto field_ptr = mField.get_field_structure(field_name);
356  return field_ptr->getNbOfCoeffs();
357  };
358 
359  const auto nb_field_coeffs = get_field_coeffs(field_name);
360  constexpr auto max_nb_dofs_per_node = 6;
361 
362  if (nb_field_coeffs > max_nb_dofs_per_node)
363  MOFEM_LOG("WORLD", Sev::error)
364  << "MultiPointConstraints PreProcLhs<MPCsType>: support only "
365  "up to 6 dofs per node.";
366  MOFEM_LOG("WORLD", Sev::noisy)
367  << "Apply MultiPointConstraints PreProc<MPCsType>: "
368  << problem_name << "_" << field_name << "_" << block_name;
369 
370  auto mpc_type = mpc_bc->mpcType;
371  switch (mpc_type) {
372  case MPC::COUPLING:
373  case MPC::TIE:
374  case MPC::RIGID_BODY:
375  break;
376  default:
377  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
378  "MPC type not implemented");
379  }
380 
381  // FIXME: there handle the vertices differently (block level)
382  auto master_verts = mpc_bc->mpcMasterEnts.subset_by_type(MBVERTEX);
383  auto slave_verts = mpc_bc->mpcSlaveEnts.subset_by_type(MBVERTEX);
384 
385  auto prb_name = fe_ptr->problemPtr->getName();
386  auto get_flag = [&](int idx) { return (&mpc_bc->data.flag1)[idx]; };
387  MOFEM_LOG("WORLD", Sev::verbose)
388  << "Size of master nodes: " << mpc_bc->mpcMasterEnts.size()
389  << " Size of slave nodes : " << mpc_bc->mpcSlaveEnts.size();
390 
391  if (mpc_bc->mpcMasterEnts.size() > mpc_bc->mpcSlaveEnts.size()) {
392  MOFEM_LOG("WORLD", Sev::error)
393  << "Size of master nodes < Size of slave nodes : "
394  << mpc_bc->mpcMasterEnts.size() << " > " << mpc_bc->mpcSlaveEnts.size();
395  // SETERRQ(PETSC_COMM_WORLD, MOFEM_OPERATION_UNSUCCESSFUL,
396  // "data inconsistency");
397  }
398 
399  auto add_is = [](auto is1, auto is2) {
400  IS is;
401  CHK_THROW_MESSAGE(ISExpand(is1, is2, &is), "is sum");
402  return SmartPetscObj<IS>(is);
403  };
404 
405  SmartPetscObj<IS> is_xyz_row_sum;
406  SmartPetscObj<IS> is_m_row[max_nb_dofs_per_node];
407  SmartPetscObj<IS> is_m_col[max_nb_dofs_per_node];
408  SmartPetscObj<IS> is_s_row[max_nb_dofs_per_node];
409  SmartPetscObj<IS> is_s_col[max_nb_dofs_per_node];
410 
411  for (int dd = 0; dd != nb_field_coeffs; dd++) {
412  if (get_flag(dd)) {
413  // SmartPetscObj<IS> is_xyz_m;
414  CHKERR is_mng->isCreateProblemFieldAndRank(
415  prb_name, ROW, field_name, dd, dd, is_m_row[dd],
416  &master_verts);
417  CHKERR is_mng->isCreateProblemFieldAndRank(
418  prb_name, COL, field_name, dd, dd, is_m_col[dd],
419  &master_verts);
420  CHKERR is_mng->isCreateProblemFieldAndRank(
421  prb_name, COL, field_name, dd, dd, is_s_col[dd],
422  &slave_verts);
423  CHKERR is_mng->isCreateProblemFieldAndRank(
424  prb_name, ROW, field_name, dd, dd, is_s_row[dd],
425  &slave_verts);
426  // ISView(is_s_row[dd], PETSC_VIEWER_STDOUT_WORLD);
427  // ISView(is_s_col[dd], PETSC_VIEWER_STDOUT_WORLD);
428 
429  if (!mpc_bc->isReprocitical) {
430  if (!is_xyz_row_sum) {
431  is_xyz_row_sum = is_s_row[dd];
432  } else {
433  is_xyz_row_sum = add_is(is_xyz_row_sum, is_s_row[dd]);
434  }
435  }
436  }
437  }
438 
439  // if (is_xyz_row_sum) {
440  SmartPetscObj<Mat> B =
441  vLhs ? vLhs : SmartPetscObj<Mat>(fe_ptr->B, true);
442  // The user is responsible for assembly if vLhs is provided
443  MatType typem;
444  CHKERR MatGetType(B, &typem);
445  CHKERR MatSetOption(B, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE);
446  PetscBool is_mat_baij = PETSC_FALSE;
447  PetscObjectTypeCompare((PetscObject)B, MATMPIBAIJ, &is_mat_baij);
448  if (is_mat_baij)
449  CHKERR MatSetOption(B, MAT_USE_HASH_TABLE, PETSC_TRUE);
450  if ((*fe_ptr->matAssembleSwitch) && !vLhs) {
451  if (*fe_ptr->matAssembleSwitch) {
452  CHKERR MatAssemblyBegin(B, MAT_FLUSH_ASSEMBLY);
453  CHKERR MatAssemblyEnd(B, MAT_FLUSH_ASSEMBLY);
454  *fe_ptr->matAssembleSwitch = false;
455  }
456  }
457 
458  if (vAO) {
459  MOFEM_LOG("WORLD", Sev::error) << "No support for AO yet";
460  }
461 
462  CHKERR MatZeroRowsIS(B, is_xyz_row_sum, 0, PETSC_NULL,
463  PETSC_NULL);
464  // }
465  auto set_mat_values = [&](auto row_is, auto col_is, double val, double perturb = 0) {
467  const int *row_index_ptr;
468  CHKERR ISGetIndices(row_is, &row_index_ptr);
469  const int *col_index_ptr;
470  CHKERR ISGetIndices(col_is, &col_index_ptr);
471  int row_size, col_size;
472  CHKERR ISGetLocalSize(row_is, &row_size);
473  CHKERR ISGetLocalSize(col_is, &col_size);
474 
475  MatrixDouble locMat(row_size, col_size);
476  fill(locMat.data().begin(), locMat.data().end(), 0.0);
477  for (int i = 0; i != row_size; i++)
478  for (int j = 0; j != col_size; j++)
479  if (i == j || col_size == 1)
480  locMat(i, j) = val;
481 
482  if (mpc_bc->isReprocitical) {
483  locMat *= perturb;
484  CHKERR ::MatSetValues(B, row_size, row_index_ptr, col_size,
485  col_index_ptr, &*locMat.data().begin(),
486  ADD_VALUES);
487  } else {
488  CHKERR ::MatSetValues(B, row_size, row_index_ptr, col_size,
489  col_index_ptr, &*locMat.data().begin(),
490  INSERT_VALUES);
491  }
492 
493  CHKERR ISRestoreIndices(row_is, &row_index_ptr);
494  CHKERR ISRestoreIndices(col_is, &col_index_ptr);
495 
497  };
498 
499  for (int dd = 0; dd != nb_field_coeffs; dd++) {
500  if (get_flag(dd)) {
501  do_assembly = true;
502  if (!mpc_bc->isReprocitical) {
503  CHKERR set_mat_values(is_s_row[dd], is_s_col[dd], vDiag);
504  CHKERR set_mat_values(is_s_row[dd], is_m_col[dd], -vDiag);
505  }
506  else
507  {
508  auto &pn = mpc_bc->mPenalty;
509  CHKERR set_mat_values(is_s_row[dd], is_s_col[dd], vDiag, pn);
510  CHKERR set_mat_values(is_s_row[dd], is_m_col[dd], -vDiag, pn);
511  CHKERR set_mat_values(is_m_row[dd], is_m_col[dd], vDiag, pn);
512  CHKERR set_mat_values(is_m_row[dd], is_s_col[dd], -vDiag, pn);
513  }
514  }
515  }
516  } // if regex
517  } // mpc loop
518  } // bc loop
519  if (do_assembly) {
520  SmartPetscObj<Mat> B = vLhs ? vLhs : SmartPetscObj<Mat>(fe_ptr->B, true);
521  CHKERR MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
522  CHKERR MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
523  }
524 
525  } else {
526  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
527  "Can not lock shared pointer");
528  } // if fe_ptr
530 }

Member Data Documentation

◆ fePtr

boost::weak_ptr<FEMethod> MoFEM::EssentialPostProcLhs< MPCsType >::fePtr
protected

Definition at line 115 of file EssentialMPCsData.hpp.

◆ mField

Definition at line 114 of file EssentialMPCsData.hpp.

◆ vAO

Definition at line 118 of file EssentialMPCsData.hpp.

◆ vDiag

Definition at line 116 of file EssentialMPCsData.hpp.

◆ vLhs

Definition at line 117 of file EssentialMPCsData.hpp.


The documentation for this struct was generated from the following files:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
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:1389
MoFEM::Field::getNbOfCoeffs
FieldCoefficientsNumber getNbOfCoeffs() const
Get number of field coefficients.
Definition: FieldMultiIndices.hpp:202
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1634
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::CoreInterface::get_field_structure
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
ROW
@ ROW
Definition: definitions.h:136
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::EssentialPostProcLhs< MPCsType >::fePtr
boost::weak_ptr< FEMethod > fePtr
Definition: EssentialMPCsData.hpp:115
MoFEM::EssentialPostProcLhs< MPCsType >::vDiag
double vDiag
Definition: EssentialMPCsData.hpp:116
MoFEM::MPC::RIGID_BODY
@ RIGID_BODY
COL
@ COL
Definition: definitions.h:136
MoFEM::MPC::TIE
@ TIE
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::EssentialPostProcLhs< MPCsType >::vAO
SmartPetscObj< AO > vAO
Definition: EssentialMPCsData.hpp:118
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::EssentialPostProcLhs< MPCsType >::vLhs
SmartPetscObj< Mat > vLhs
Definition: EssentialMPCsData.hpp:117
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::MPC::COUPLING
@ COUPLING
MoFEM::EssentialPostProcLhs< MPCsType >::mField
MoFEM::Interface & mField
Definition: EssentialMPCsData.hpp:114