MultifieldPlasticity.hpp

Go to the documentation of this file.

/** \file MultifieldPlasticity
  * \brief Multifield plasticity
 
*/
 
/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
namespace MoFEM {
 
extern boost::weak_ptr<BlockMatContainer> block_mat_container_ptr;
 
template <typename EleOp>
struct OpBaseImpl<USER_ASSEMBLE, EleOp>
    : public OpBaseImpl<LAST_ASSEMBLE, EleOp> {
  using OpBaseImpl<LAST_ASSEMBLE, EleOp>::OpBaseImpl;
 
protected:
 
  MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data,
                          EntitiesFieldData::EntData &col_data,
                          const bool transpose) {
    MoFEMFunctionBegin;
 
    if (transpose) {
      this->locMatTranspose.resize(this->locMat.size2(), this->locMat.size1(),
                                   false);
      noalias(this->locMatTranspose) = trans(this->locMat);
    }
 
    if (auto bmc_ptr = block_mat_container_ptr.lock()) {
      // if (true) {
      // if (this->rowFieldName == "SIGMA" || this->colFieldName == "SIGMA") {
      //   // cout << "WE ARE ADDING THE BLOCK" << endl;
      //   // cout << this->rowFieldName << " " << this->colFieldName << endl;
 
      // }
      auto &bmc = *bmc_ptr;
      int *row_ind_ptr = &*row_data.getIndices().begin();
      int *col_ind_ptr = &*col_data.getIndices().begin();
      int row_nb_dofs = row_data.getIndices().size();
      int col_nb_dofs = col_data.getIndices().size();
      auto &K = this->locMat;
      auto &transK = this->locMatTranspose;
 
      auto set_block = [&](auto &row_name, auto &col_name, auto row_side,
                           auto col_side, auto row_type, auto col_type,
                           const auto &m, const auto &row_ind,
                           const auto &col_ind) {
        auto it = bmc.get<0>().find(boost::make_tuple(
            row_name, col_name, row_type, col_type, row_side, col_side));
        if (it != bmc.get<0>().end()) {
          it->setMat(m);
          it->setInd(row_ind, col_ind);
          it->setSetAtElement();
        } else
          bmc.insert(BlockMatData(row_name, col_name, row_type, col_type,
                                  row_side, col_side, m, row_ind, col_ind));
      };
      //FIXME: for contact fieldsplit
      if (/* true || */ this->rowFieldName != "SIGMA" &&
          this->colFieldName != "SIGMA") {
        if (auto e_ptr = row_data.getFieldEntities()[0]) {
          if (auto stored_data_ptr =
                  e_ptr->getSharedStoragePtr<EssentialBcStorage>()) {
            if (!this->onlyTranspose)
              set_block(this->rowFieldName, this->colFieldName, this->rowSide,
                        this->colSide, this->rowType, this->colType, K,
                        stored_data_ptr->entityIndices, col_data.getIndices());
            if (transpose)
              set_block(this->colFieldName, this->rowFieldName, this->colSide,
                        this->rowSide, this->colType, this->rowType, transK,
                        col_data.getIndices(), stored_data_ptr->entityIndices);
          } else {
            if (!this->onlyTranspose)
              set_block(this->rowFieldName, this->colFieldName, this->rowSide,
                        this->colSide, this->rowType, this->colType, K,
                        row_data.getIndices(), col_data.getIndices());
 
            if (transpose)
              set_block(this->colFieldName, this->rowFieldName, this->colSide,
                        this->rowSide, this->colType, this->rowType, transK,
                        col_data.getIndices(), row_data.getIndices());
          }
        } else
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "No field entities for assembly");
      }
      }
 
    // Assemble transpose
    if (transpose) {
      CHKERR MatSetValues<EssentialBcStorage>(
          this->getKSPB(), col_data, row_data,
          &*this->locMatTranspose.data().begin(), ADD_VALUES);
    }
 
    if (!this->onlyTranspose) {
      // assemble local matrix
      CHKERR MatSetValues<EssentialBcStorage>(
          this->getKSPB(), row_data, col_data, &*this->locMat.data().begin(),
          ADD_VALUES);
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data) {
    return VecSetValues<EssentialBcStorage>(
        this->getKSPf(), data, &*this->locF.data().begin(), ADD_VALUES);
  }
};
 
using EntData = EntitiesFieldData::EntData;
 
using DomainEle = VolumeElementForcesAndSourcesCore;
using DomainEleOp = DomainEle::UserDataOperator;
using DomainEleOpAssembly =
    FormsIntegrators<DomainEleOp>::Assembly<USER_ASSEMBLE>::OpBase;
 
using SkeletonEle = FaceElementForcesAndSourcesCore;
using SkeletonEleOp = SkeletonEle::UserDataOperator;
using SkeletonEleOpAssembly =
    FormsIntegrators<SkeletonEleOp>::Assembly<USER_ASSEMBLE>::OpBase;
 
using BoundaryEle = FaceElementForcesAndSourcesCore;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
using BoundaryEleOpAssembly =
    FormsIntegrators<BoundaryEleOp>::Assembly<USER_ASSEMBLE>::OpBase;
 
using DomainSideEle = VolumeElementForcesAndSourcesCoreOnSide;
using DomainSideEleOp = DomainSideEle::UserDataOperator;
using DomainSideEleOpAssembly =
    FormsIntegrators<DomainSideEleOp>::Assembly<USER_ASSEMBLE>::OpBase;
 
 
MoFEMErrorCode set_plastic_ep_dm(MoFEM::Interface &m_field,
                                 SmartPetscObj<DM> &dm,
                                 SmartPetscObj<DM> &dm_plastic);
MoFEMErrorCode set_plastic_tau_dm(MoFEM::Interface &m_field,
                                  SmartPetscObj<DM> &dm,
                                  SmartPetscObj<DM> &dm_plastic_tau);
MoFEMErrorCode set_contact_dm(MoFEM::Interface &m_field, SmartPetscObj<DM> &dm,
                              SmartPetscObj<DM> &dm_contact);
 
} // namespace MoFEM