MortarContactProblem.hpp

Go to the documentation of this file.

/* 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/>. */
 
#ifndef __MORTAR_CONTACT_PROBLEM_HPP__
#define __MORTAR_CONTACT_PROBLEM_HPP__
 
struct MortarContactProblem : public SimpleContactProblem {
 
  using ContactOp = ContactPrismElementForcesAndSourcesCore::UserDataOperator;
 
  struct MortarContactPrismsData {
    Range pRisms; // All boundary surfaces
  };
 
  map<int, MortarContactPrismsData>
      setOfMortarContactPrism; ///< maps side set id with appropriate FluxData
 
  struct MortarContactElement
      : public SimpleContactProblem::ConvectMasterContactElement {
    MoFEM::Interface &mField;
    boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
        contactCommondataMultiIndex;
 
    bool newtonCotes;
    MortarContactElement(
        MoFEM::Interface &m_field,
        boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
            contact_commondata_multi_index,
        std::string spat_pos, std::string mat_pos, bool newton_cotes = false)
        : SimpleContactProblem::ConvectMasterContactElement(
              m_field, spat_pos, mat_pos, newton_cotes),
          mField(m_field),
          contactCommondataMultiIndex(contact_commondata_multi_index),
          newtonCotes(newton_cotes) {}
 
    int getRule(int order) { return -1; };
 
    // function to calculate area of triangle (copy from moab/CslamUtils.cpp)
    inline double area2D(double *a, double *b, double *c) {
      // (b-a)x(c-a) / 2
      return ((b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])) /
             2;
    }
    virtual MoFEMErrorCode setGaussPts(int order);
 
    // Destructor
    ~MortarContactElement() {}
  };
 
  struct CommonDataMortarContact
      : public SimpleContactProblem::CommonDataSimpleContact {
 
    boost::shared_ptr<VectorDouble> referenceGapPtr;
    boost::shared_ptr<MatrixDouble> referencePositionAtGaussPtsMasterPtr;
    boost::shared_ptr<MatrixDouble> referencePositionAtGaussPtsSlavePtr;
 
    bool addContactReferenceGap;
    CommonDataMortarContact(MoFEM::Interface &m_field)
        : SimpleContactProblem::CommonDataSimpleContact(m_field),
          mField(m_field) {
      addContactReferenceGap = false;
      referenceGapPtr = boost::make_shared<VectorDouble>();
      referencePositionAtGaussPtsMasterPtr = boost::make_shared<MatrixDouble>();
      referencePositionAtGaussPtsSlavePtr = boost::make_shared<MatrixDouble>();
    }
 
  private:
    MoFEM::Interface &mField;
  };
 
  struct LoadScale : public MethodForForceScaling {
 
    static double lAmbda;
 
    MoFEMErrorCode scaleNf(const FEMethod *fe, VectorDouble &nf) {
      MoFEMFunctionBegin;
      nf *= lAmbda;
      MoFEMFunctionReturn(0);
    }
  };
 
  MoFEMErrorCode addMortarContactElement(
      const string element_name, const string field_name,
      const string lagrange_field_name, Range &range_slave_master_prisms,
      string material_position_field_name = "MESH_NODE_POSITIONS",
      bool eigen_pos_flag = false,
      const string eigen_node_field_name = "EIGEN_SPATIAL_POSITIONS") {
    MoFEMFunctionBegin;
 
    CHKERR mField.add_finite_element(element_name, MF_ZERO);
 
    if (range_slave_master_prisms.size() > 0) {
 
      CHKERR mField.modify_finite_element_add_field_row(element_name,
                                                        lagrange_field_name);
 
      CHKERR mField.modify_finite_element_add_field_col(element_name,
                                                        field_name);
 
      CHKERR mField.modify_finite_element_add_field_col(element_name,
                                                        lagrange_field_name);
      CHKERR mField.modify_finite_element_add_field_row(element_name,
                                                        field_name);
 
      CHKERR mField.modify_finite_element_add_field_data(element_name,
                                                         lagrange_field_name);
 
      CHKERR mField.modify_finite_element_add_field_data(element_name,
                                                         field_name);
 
      if (eigen_pos_flag)
        CHKERR mField.modify_finite_element_add_field_data(
            element_name, eigen_node_field_name);
 
      CHKERR
      mField.modify_finite_element_add_field_data(element_name,
                                                  "MESH_NODE_POSITIONS");
 
      setOfMortarContactPrism[1].pRisms = range_slave_master_prisms;
 
      // Adding range_slave_master_prisms to Element element_name
      mField.add_ents_to_finite_element_by_type(range_slave_master_prisms,
                                                MBPRISM, element_name);
    }
 
    MoFEMFunctionReturn(0);
  }
 
  /**
   * @brief Function that adds field data for spatial positions and Lagrange
   * multipliers to rows and columns, provides access to field data and adds
   * prism entities to element.
   *
   * @param  element_name               String for the element name
   * @param  field_name                 String of field name for spatial
   * position
   * @param  lagrange_field_name         String of field name for Lagrange
   * multipliers
   * @param  mesh_node_field_name       String of field name for material
   * positions
   * @param  range_slave_master_prisms  Range for prism entities used to create
   * contact elements
   * @return                            Error code
   *
   */
  MoFEMErrorCode addMortarContactElementALE(const string element_name,
                                            const string field_name,
                                            const string mesh_node_field_name,
                                            const string lagrange_field_name,
                                            Range &range_slave_master_prisms) {
    MoFEMFunctionBegin;
 
    CHKERR mField.add_finite_element(element_name, MF_ZERO);
 
    if (range_slave_master_prisms.size() > 0) {
 
      CHKERR mField.modify_finite_element_add_field_row(element_name,
                                                        lagrange_field_name);
 
      CHKERR mField.modify_finite_element_add_field_col(element_name,
                                                        field_name);
      CHKERR mField.modify_finite_element_add_field_col(element_name,
                                                        mesh_node_field_name);
 
      CHKERR mField.modify_finite_element_add_field_col(element_name,
                                                        lagrange_field_name);
      CHKERR mField.modify_finite_element_add_field_row(element_name,
                                                        field_name);
      CHKERR mField.modify_finite_element_add_field_row(element_name,
                                                        mesh_node_field_name);
 
      CHKERR mField.modify_finite_element_add_field_data(element_name,
                                                         lagrange_field_name);
 
      CHKERR mField.modify_finite_element_add_field_data(element_name,
                                                         field_name);
 
      CHKERR mField.modify_finite_element_add_field_data(element_name,
                                                         mesh_node_field_name);
 
      setOfSimpleContactPrism[1].pRisms = range_slave_master_prisms;
 
      Range ents_to_add = range_slave_master_prisms;
      Range current_ents_with_fe;
      CHKERR mField.get_finite_element_entities_by_handle(element_name,
                                                          current_ents_with_fe);
      Range ents_to_remove;
      ents_to_remove = subtract(current_ents_with_fe, ents_to_add);
      CHKERR mField.remove_ents_from_finite_element(element_name,
                                                    ents_to_remove);
      CHKERR mField.add_ents_to_finite_element_by_type(ents_to_add, MBPRISM,
                                                       element_name);
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode addMortarContactFrictionElement(
      const string element_name, const string field_name,
      const string lagrange_field_name, const string prev_conv_field_name,
      const string tangent_lagrange_field_name,
      Range &range_slave_master_prisms,
      string material_position_field_name = "MESH_NODE_POSITIONS") {
    MoFEMFunctionBegin;
 
    CHKERR mField.add_finite_element(element_name, MF_ZERO);
 
    //============================================================================================================
    // C row as Lagrange_mul and col as DISPLACEMENT
 
    CHKERR mField.modify_finite_element_add_field_row(element_name,
                                                      lagrange_field_name);
 
    CHKERR mField.modify_finite_element_add_field_col(element_name, field_name);
 
    // CT col as Lagrange_mul and row as DISPLACEMENT
 
    CHKERR mField.modify_finite_element_add_field_col(element_name,
                                                      lagrange_field_name);
 
    CHKERR mField.modify_finite_element_add_field_row(element_name, field_name);
 
    // data
    CHKERR mField.modify_finite_element_add_field_data(element_name,
                                                       lagrange_field_name);
 
    CHKERR mField.modify_finite_element_add_field_data(element_name,
                                                       field_name);
 
    CHKERR
    mField.modify_finite_element_add_field_data(element_name,
                                                "MESH_NODE_POSITIONS");
 
    CHKERR
    mField.modify_finite_element_add_field_data(element_name,
                                                prev_conv_field_name);
 
    setOfMortarContactPrism[1].pRisms = range_slave_master_prisms;
    mField.add_ents_to_finite_element_by_type(range_slave_master_prisms,
                                              MBPRISM, element_name);
 
    MoFEMFunctionReturn(0);
  }
 
  double cnValue;
  boost::shared_ptr<double> cnValuePtr;
 
  MoFEM::Interface &mField;
  boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
      contactCommondataMultiIndex;
  bool newtonCotes;
 
  MortarContactProblem(
      MoFEM::Interface &m_field,
      boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
          contact_commondata_multi_index,
      boost::shared_ptr<double> cn_value_ptr, bool newton_cotes = false)
      : SimpleContactProblem(m_field, cn_value_ptr, newton_cotes),
        mField(m_field),
        contactCommondataMultiIndex(contact_commondata_multi_index),
        cnValuePtr(cn_value_ptr), newtonCotes(newton_cotes) {
    cnValue = *cnValuePtr.get();
  }
 
  /**
   * @brief Element used to integrate on slave surfaces. It convects integration
   * points on slaves, so that quantities like gap from master are evaluated at
   * correct points.
   *
   */
  struct MortarConvectMasterContactElement : public MortarContactElement {
 
    MortarConvectMasterContactElement(
        MoFEM::Interface &m_field,
        boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
            contact_commondata_multi_index,
        std::string spat_pos, std::string mat_pos, bool newton_cotes = false)
        : MortarContactElement(m_field, contact_commondata_multi_index,
                               spat_pos, mat_pos, newton_cotes),
          convectPtr(new ConvectSlaveIntegrationPts(this, spat_pos, mat_pos)) {}
 
    inline boost::shared_ptr<ConvectSlaveIntegrationPts> getConvectPtr() {
      return convectPtr;
    }
 
    int getRule(int order) { return -1; }
 
    MoFEMErrorCode setGaussPts(int order);
 
  protected:
    boost::shared_ptr<ConvectSlaveIntegrationPts> convectPtr;
  };
 
  /**
   * @brief Element used to integrate on master surfaces. It convects
   * integration points on slaves, so that quantities like Lagrange multiplier
   * from master are evaluated at correct points.
   *
   */
  struct MortarConvectSlaveContactElement
      : public MortarConvectMasterContactElement {
    using MortarConvectMasterContactElement::MortarConvectMasterContactElement;
 
    int getRule(int order) { return -1; }
 
    MoFEMErrorCode setGaussPts(int order);
  };
 
  MoFEMErrorCode setContactOperatorsRhs(
      boost::shared_ptr<MortarContactElement> fe_rhs_mortar_contact,
      boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
      string field_name, string lagrange_field_name, bool is_alm = false,
      bool is_eigen_pos_field = false,
      string eigen_pos_field_name = "EIGEN_SPATIAL_POSITIONS",
      bool use_reference_coordinates = false);
 
  MoFEMErrorCode setMasterForceOperatorsRhs(
      boost::shared_ptr<MortarContactElement> fe_rhs_mortar_contact,
      boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
      string field_name, string lagrange_field_name, bool is_alm = false,
      bool is_eigen_pos_field = false,
      string eigen_pos_field_name = "EIGEN_SPATIAL_POSITIONS",
      bool use_reference_coordinates = false);
 
  MoFEMErrorCode setContactOperatorsLhs(
      boost::shared_ptr<MortarContactElement> fe_lhs_mortar_contact,
      boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
      string field_name, string lagrange_field_name, bool is_alm = false,
      bool is_eigen_pos_field = false,
      string eigen_pos_field_name = "EIGEN_SPATIAL_POSITIONS",
      bool use_reference_coordinates = false);
 
  MoFEMErrorCode setContactOperatorsLhs(
      boost::shared_ptr<MortarConvectMasterContactElement>
          fe_lhs_simple_contact,
      boost::shared_ptr<CommonDataMortarContact> common_data_simple_contact,
      string field_name, string lagrange_field_name, bool is_alm = false,
      bool is_eigen_pos_field = false,
      string eigen_pos_field_name = "EIGEN_SPATIAL_POSITIONS",
      bool use_reference_coordinates = false);
 
  MoFEMErrorCode setMasterForceOperatorsLhs(
      boost::shared_ptr<MortarContactElement> fe_lhs_mortar_contact,
      boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
      string field_name, string lagrange_field_name, bool is_alm = false,
      bool is_eigen_pos_field = false,
      string eigen_pos_field_name = "EIGEN_SPATIAL_POSITIONS",
      bool use_reference_coordinates = false);
 
  MoFEMErrorCode setMasterForceOperatorsLhs(
      boost::shared_ptr<MortarConvectSlaveContactElement> fe_lhs_mortar_contact,
      boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
      string field_name, string lagrange_field_name, bool is_alm = false,
      bool is_eigen_pos_field = false,
      string eigen_pos_field_name = "EIGEN_SPATIAL_POSITIONS",
      bool use_reference_coordinates = false);
 
  MoFEMErrorCode setContactOperatorsForPostProc(
      boost::shared_ptr<MortarContactElement> fe_post_proc_mortar_contact,
      boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
      MoFEM::Interface &m_field, string field_name, string lagrange_field_name,
      moab::Interface &moab_out, bool alm_flag = false,
      bool is_eigen_pos_field = false,
      string eigen_pos_field_name = "EIGEN_SPATIAL_POSITIONS",
      bool is_displacements = false, Range post_proc_surface = Range(),
      double post_proc_gap_tol = 0.0);
 
  /**
   * @brief Function for the mortar contact element that sets the user data
   * post processing operators
   *
   * @param  fe_post_proc_simple_contact Pointer to the FE instance for post
   * processing
   * @param  common_data_simple_contact  Pointer to the common data for
   * mortar contact element
   * @param  field_name                  String of field name for spatial
   * positions
   * @param  mesh_node_field_name        String of field name for material
   * positions
   * @param  lagrange_field_name          String of field name for Lagrange
   * multipliers
   * @param  side_fe_name                String of 3D element adjacent to
   * the present contact element
   * @return                             Error code
   *
   */
  MoFEMErrorCode setContactOperatorsRhsALEMaterial(
      boost::shared_ptr<MortarContactElement> fe_rhs_simple_contact_ale,
      boost::shared_ptr<CommonDataMortarContact> common_data_simple_contact,
      const string field_name, const string mesh_node_field_name,
      const string lagrange_field_name, const string side_fe_name);
 
  /**
   * @brief Function for the mortar contact element that sets the user data
   * LHS-operators
   *
   * @param  fe_lhs_simple_contact_ale  Pointer to the FE instance for LHS
   * @param  common_data_simple_contact Pointer to the common data for
   * mortar contact element
   * @param  field_name                 String of field name for spatial
   * positions
   * @param  mesh_node_field_name       String of field name for material
   * positions
   * @param  lagrange_field_name         String of field name for Lagrange
   * multipliers
   * @param  side_fe_name               String of 3D element adjacent to the
   * present contact element
   * @return                            Error code
   *
   */
  MoFEMErrorCode setContactOperatorsLhsALEMaterial(
      boost::shared_ptr<MortarContactElement> fe_lhs_simple_contact_ale,
      boost::shared_ptr<CommonDataMortarContact> common_data_simple_contact,
      const string field_name, const string mesh_node_field_name,
      const string lagrange_field_name, const string side_fe_name);
 
  /**
   * @brief Function for the mortar contact element that sets the user data
   * LHS-operators
   *
   * @param  fe_lhs_simple_contact_ale  Pointer to the FE instance for LHS
   * @param  common_data_simple_contact Pointer to the common data for
   * mortar contact element
   * @param  field_name                 String of field name for spatial
   * positions
   * @param  mesh_node_field_name       String of field name for material
   * positions
   * @param  lagrange_field_name         String of field name for Lagrange
   * multipliers
   * @return                            Error code
   *
   */
  MoFEMErrorCode setContactOperatorsLhsALE(
      boost::shared_ptr<MortarContactElement> fe_lhs_simple_contact_ale,
      boost::shared_ptr<CommonDataMortarContact> common_data_simple_contact,
      const string field_name, const string mesh_node_field_name,
      const string lagrange_field_name);
 
}; // end of MortarContactProblem
#endif