C++ PenetrationInfo类代码示例

void
ContactMaster::addPoints()
{
  _point_to_info.clear();

  std::map<dof_id_type, PenetrationInfo *>::iterator
    it  = _penetration_locator._penetration_info.begin(),
    end = _penetration_locator._penetration_info.end();

  for (; it!=end; ++it)
  {
    PenetrationInfo * pinfo = it->second;

    // Skip this pinfo if there are no DOFs on this node.
    if ( ! pinfo || pinfo->_node->n_comp(_sys.number(), _vars(_component)) < 1 )
      continue;

    if ( pinfo->isCaptured() )
    {
      addPoint(pinfo->_elem, pinfo->_closest_point);
      _point_to_info[pinfo->_closest_point] = pinfo;
      computeContactForce(pinfo);
    }
  }
}

Real
MechanicalContactConstraint::computeQpResidual(Moose::ConstraintType type)
{
  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
  computeContactForce(pinfo);
  Real resid = pinfo->_contact_force(_component);

  switch(type)
  {
    case Moose::Slave:
      if (_formulation == CF_DEFAULT)
      {
        //Real distance = (*_current_node)(_component) - pinfo->_closest_point(_component);
        //Real pen_force = _penalty * distance;
        RealVectorValue distance_vec(*_current_node - pinfo->_closest_point);
        RealVectorValue pen_force(_penalty * distance_vec);

        if (_model == CM_FRICTIONLESS || _model == CM_EXPERIMENTAL)
          resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;

        else if (_model == CM_GLUED || _model == CM_TIED || _model == CM_COULOMB)
          resid += pen_force(_component);

      }
      return _test_slave[_i][_qp] * resid;
    case Moose::Master:
      return _test_master[_i][_qp] * -resid;
  }

  return 0;
}

Real
SlaveConstraint::computeQpResidual()
{
  PenetrationInfo * pinfo = _point_to_info[_current_point];
  const Node * node = pinfo->_node;

  Real resid = pinfo->_contact_force(_component);

  const Real area = nodalArea(*pinfo);

  if (_formulation == CF_DEFAULT)
  {
    RealVectorValue distance_vec(_mesh.node(node->id()) - pinfo->_closest_point);
    RealVectorValue pen_force(_penalty * distance_vec);
    if (_normalize_penalty)
      pen_force *= area;

    if (_model == CM_FRICTIONLESS)
      resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;

    else if (_model == CM_GLUED || _model == CM_COULOMB)
      resid += pen_force(_component);

  }

  return _test[_i][_qp] * resid;
}

Real
ContactMaster::computeQpResidual()
{
  PenetrationInfo * pinfo = _point_to_info[_current_point];
  Real resid = -pinfo->_contact_force(_component);
  return _test[_i][_qp] * resid;
}

Real
GluedContactConstraint::computeQpResidual(Moose::ConstraintType type)
{
  switch (type)
  {
    case Moose::Slave:
    {
      PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
      Real distance = (*_current_node)(_component) - pinfo->_closest_point(_component);
      Real pen_force = _penalty * distance;

      Real resid = pen_force;
      pinfo->_contact_force(_component) = resid;
      pinfo->_mech_status=PenetrationInfo::MS_STICKING;

      return _test_slave[_i][_qp] * resid;
    }
    case Moose::Master:
    {
      PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];

      long int dof_number = _current_node->dof_number(0, _vars(_component), 0);
      Real resid = _residual_copy(dof_number);

      pinfo->_contact_force(_component) = -resid;
      pinfo->_mech_status=PenetrationInfo::MS_STICKING;

      return _test_master[_i][_qp] * resid;
    }
  }

  return 0;
}

Real
MultiDContactConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
{
  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];

  double slave_jac = 0;
  switch (type)
  {
  case Moose::SlaveSlave:
    switch (_model)
    {
    case CM_FRICTIONLESS:

      slave_jac = pinfo->_normal(_component) * pinfo->_normal(_component) * ( _penalty*_phi_slave[_j][_qp] - (*_jacobian)(_current_node->dof_number(0, _var.number(), 0), _connected_dof_indices[_j]) );
      break;

    case CM_GLUED:
/*
      resid = pen_force(_component)
        - res_vec(_component)
        ;
*/
      break;

    default:
      mooseError("Invalid or unavailable contact model");
      break;
    }
    return _test_slave[_i][_qp] * slave_jac;
  case Moose::SlaveMaster:
    switch (_model)
    {
    case CM_FRICTIONLESS:

      slave_jac = pinfo->_normal(_component) * pinfo->_normal(_component) * ( -_penalty*_phi_master[_j][_qp] );
      break;

    case CM_GLUED:
/*
      resid = pen_force(_component)
        - res_vec(_component)
        ;
*/
      break;

    default:
      mooseError("Invalid or unavailable contact model");
      break;
    }
    return _test_slave[_i][_qp] * slave_jac;
  case Moose::MasterSlave:
    slave_jac = (*_jacobian)(_current_node->dof_number(0, _var.number(), 0), _connected_dof_indices[_j]);
    return slave_jac*_test_master[_i][_qp];
  case Moose::MasterMaster:
    return 0;
  }
  return 0;
}

Real
OneDContactConstraint::computeQpSlaveValue()
{
  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
  Moose::err<<std::endl
           <<"Popping out node: "<<_current_node->id()<<std::endl
           <<"Closest Point x: "<<pinfo->_closest_point(0)<<std::endl
           <<"Current Node x: "<<(*_current_node)(0)<<std::endl
           <<"Current Value: "<<_u_slave[_qp]<<std::endl<<std::endl;

  return pinfo->_closest_point(0) - ((*_current_node)(0) - _u_slave[_qp]);
}

Real
MultiDContactConstraint::computeQpSlaveValue()
{
  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
/*
  Moose::err<<std::endl
           <<"Popping out node: "<<_current_node->id()<<std::endl
           <<"Closest Point "<<_component<<": "<<pinfo->_closest_point(_component)<<std::endl
           <<"Current Node "<<_component<<": "<<(*_current_node)(_component)<<std::endl
           <<"Current Value: "<<_u_slave[_qp]<<std::endl
           <<"New Value: "<<pinfo->_closest_point(_component) - ((*_current_node)(_component) - _u_slave[_qp])<<std::endl
           <<"Change: "<<_u_slave[_qp] - (pinfo->_closest_point(_component) - ((*_current_node)(_component) - _u_slave[_qp]))<<std::endl<<std::endl;
*/
  return  pinfo->_closest_point(_component) - ((*_current_node)(_component) - _u_slave[_qp]);
}

Real
OneDContactConstraint::computeQpResidual(Moose::ConstraintType type)
{
  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];

  switch (type)
  {
  case Moose::Slave:
    // return (_u_slave[_qp] - _u_master[_qp])*_test_slave[_i][_qp];
    return ((*_current_node)(0) - pinfo->_closest_point(0)) * _test_slave[_i][_qp];
  case Moose::Master:
    double slave_resid = _residual_copy(_current_node->dof_number(0, _var.number(), 0));
    return slave_resid*_test_master[_i][_qp];
  }
  return 0;
}

void
SlaveConstraint::addPoints()
{
  _point_to_info.clear();

  std::map<dof_id_type, PenetrationInfo *>::iterator
      it = _penetration_locator._penetration_info.begin(),
      end = _penetration_locator._penetration_info.end();

  const auto & node_to_elem_map = _mesh.nodeToElemMap();
  for (; it != end; ++it)
  {
    PenetrationInfo * pinfo = it->second;

    // Skip this pinfo if there are no DOFs on this node.
    if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
      continue;

    dof_id_type slave_node_num = it->first;
    const Node * node = pinfo->_node;

    if (pinfo->isCaptured() && node->processor_id() == processor_id())
    {
      // Find an element that is connected to this node that and that is also on this processor
      auto node_to_elem_pair = node_to_elem_map.find(slave_node_num);
      mooseAssert(node_to_elem_pair != node_to_elem_map.end(), "Missing node in node to elem map");
      const std::vector<dof_id_type> & connected_elems = node_to_elem_pair->second;

      Elem * elem = NULL;

      for (unsigned int i = 0; i < connected_elems.size() && !elem; ++i)
      {
        Elem * cur_elem = _mesh.elemPtr(connected_elems[i]);
        if (cur_elem->processor_id() == processor_id())
          elem = cur_elem;
      }

      mooseAssert(elem,
                  "Couldn't find an element on this processor that is attached to the slave node!");

      addPoint(elem, *node);
      _point_to_info[*node] = pinfo;
    }
  }
}

bool
MechanicalContactConstraint::shouldApply()
{
  bool in_contact = false;

  std::map<dof_id_type, PenetrationInfo *>::iterator found =
    _penetration_locator._penetration_info.find(_current_node->id());
  if ( found != _penetration_locator._penetration_info.end() )
  {
    PenetrationInfo * pinfo = found->second;
    if ( pinfo != NULL && pinfo->isCaptured() )
    {
      in_contact = true;

      // This does the contact force once per constraint, rather than once per quad point and for
      // both master and slave cases.
      computeContactForce(pinfo);
    }
  }

  return in_contact;
}

unsigned int
FrictionalContactProblem::numLocalFrictionalConstraints()
{
  GeometricSearchData & displaced_geom_search_data = getDisplacedProblem()->geomSearchData();
  std::map<std::pair<unsigned int, unsigned int>, PenetrationLocator *> * penetration_locators = &displaced_geom_search_data._penetration_locators;

  unsigned int num_constraints(0);

  for (std::map<std::pair<unsigned int, unsigned int>, PenetrationLocator *>::iterator plit = penetration_locators->begin();
      plit != penetration_locators->end();
      ++plit)
  {
    PenetrationLocator & pen_loc = *plit->second;

    bool frictional_contact_this_interaction = false;

    std::map<std::pair<int,int>,InteractionParams>::iterator ipit;
    std::pair<int,int> ms_pair(pen_loc._master_boundary,pen_loc._slave_boundary);
    ipit = _interaction_params.find(ms_pair);
    if (ipit != _interaction_params.end())
      frictional_contact_this_interaction = true;

    if (frictional_contact_this_interaction)
    {
      std::vector<dof_id_type> & slave_nodes = pen_loc._nearest_node._slave_nodes;

      for (unsigned int i=0; i<slave_nodes.size(); i++)
      {
        dof_id_type slave_node_num = slave_nodes[i];

        PenetrationInfo * pinfo = pen_loc._penetration_info[slave_node_num];
        if (pinfo)
          if (pinfo->isCaptured())
            ++num_constraints;
      }
    }
  }
  return num_constraints;
}

void
ContactMaster::updateContactSet(bool beginning_of_step)
{
  std::map<dof_id_type, PenetrationInfo *>::iterator
    it  = _penetration_locator._penetration_info.begin(),
    end = _penetration_locator._penetration_info.end();
  for (; it!=end; ++it)
  {
    const dof_id_type slave_node_num = it->first;
    PenetrationInfo * pinfo = it->second;

    // Skip this pinfo if there are no DOFs on this node.
    if ( ! pinfo || pinfo->_node->n_comp(_sys.number(), _vars(_component)) < 1 )
      continue;

    if (beginning_of_step)
    {
      pinfo->_locked_this_step = 0;
      pinfo->_starting_elem = it->second->_elem;
      pinfo->_starting_side_num = it->second->_side_num;
      pinfo->_starting_closest_point_ref = it->second->_closest_point_ref;
      pinfo->_contact_force_old = pinfo->_contact_force;
      pinfo->_accumulated_slip_old = pinfo->_accumulated_slip;
      pinfo->_frictional_energy_old = pinfo->_frictional_energy;
    }

    const Real contact_pressure = -(pinfo->_normal * pinfo->_contact_force) / nodalArea(*pinfo);
    const Real distance = pinfo->_normal * (pinfo->_closest_point - _mesh.node(slave_node_num));

    // Capture
    if ( ! pinfo->isCaptured() && MooseUtils::absoluteFuzzyGreaterEqual(distance, 0, _capture_tolerance))
    {
      pinfo->capture();

      // Increment the lock count every time the node comes back into contact from not being in contact.
      if (_formulation == CF_KINEMATIC)
        ++pinfo->_locked_this_step;
    }
    // Release
    else if (_model != CM_GLUED &&
             pinfo->isCaptured() &&
             _tension_release >= 0 &&
             -contact_pressure >= _tension_release &&
             pinfo->_locked_this_step < 2)
    {
      pinfo->release();
      pinfo->_contact_force.zero();
    }

    if (_formulation == CF_AUGMENTED_LAGRANGE && pinfo->isCaptured())
      pinfo->_lagrange_multiplier -= getPenalty(*pinfo) * distance;
  }
}

Real
PenetrationAux::computeValue()
{
  const Node * current_node = NULL;

  if (_nodal)
    current_node = _current_node;
  else
    current_node = _mesh.getQuadratureNode(_current_elem, _current_side, _qp);

  PenetrationInfo * pinfo = _penetration_locator._penetration_info[current_node->id()];

  Real retVal(NotPenetrated);

  if (pinfo)
    switch (_quantity)
    {
      case PA_DISTANCE:
        retVal = pinfo->_distance;
        break;
      case PA_TANG_DISTANCE:
        retVal = pinfo->_tangential_distance;
        break;
      case PA_NORMAL_X:
        retVal = pinfo->_normal(0);
        break;
      case PA_NORMAL_Y:
        retVal = pinfo->_normal(1);
        break;
      case PA_NORMAL_Z:
        retVal = pinfo->_normal(2);
        break;
      case PA_CLOSEST_POINT_X:
        retVal = pinfo->_closest_point(0);
        break;
      case PA_CLOSEST_POINT_Y:
        retVal = pinfo->_closest_point(1);
        break;
      case PA_CLOSEST_POINT_Z:
        retVal = pinfo->_closest_point(2);
        break;
      case PA_ELEM_ID:
        retVal = static_cast<Real>(pinfo->_elem->id() + 1);
        break;
      case PA_SIDE:
        retVal = static_cast<Real>(pinfo->_side_num);
        break;
      case PA_INCREMENTAL_SLIP_MAG:
        retVal = pinfo->isCaptured() ? pinfo->_incremental_slip.norm() : 0;
        break;
      case PA_INCREMENTAL_SLIP_X:
        retVal = pinfo->isCaptured() ? pinfo->_incremental_slip(0) : 0;
        break;
      case PA_INCREMENTAL_SLIP_Y:
        retVal = pinfo->isCaptured() ? pinfo->_incremental_slip(1) : 0;
        break;
      case PA_INCREMENTAL_SLIP_Z:
        retVal = pinfo->isCaptured() ? pinfo->_incremental_slip(2) : 0;
        break;
      case PA_ACCUMULATED_SLIP:
        retVal = pinfo->_accumulated_slip;
        break;
      case PA_FORCE_X:
        retVal = pinfo->_contact_force(0);
        break;
      case PA_FORCE_Y:
        retVal = pinfo->_contact_force(1);
        break;
      case PA_FORCE_Z:
        retVal = pinfo->_contact_force(2);
        break;
      case PA_NORMAL_FORCE_MAG:
        retVal = -pinfo->_contact_force * pinfo->_normal;
        break;
      case PA_NORMAL_FORCE_X:
        retVal = (pinfo->_contact_force * pinfo->_normal) * pinfo->_normal(0);
        break;
      case PA_NORMAL_FORCE_Y:
        retVal = (pinfo->_contact_force * pinfo->_normal) * pinfo->_normal(1);
        break;
      case PA_NORMAL_FORCE_Z:
        retVal = (pinfo->_contact_force * pinfo->_normal) * pinfo->_normal(2);
        break;
      case PA_TANGENTIAL_FORCE_MAG:
      {
        RealVectorValue contact_force_normal((pinfo->_contact_force * pinfo->_normal) *
                                             pinfo->_normal);
        RealVectorValue contact_force_tangential(pinfo->_contact_force - contact_force_normal);
        retVal = contact_force_tangential.norm();
        break;
      }
      case PA_TANGENTIAL_FORCE_X:
        retVal =
            pinfo->_contact_force(0) - (pinfo->_contact_force * pinfo->_normal) * pinfo->_normal(0);
        break;
      case PA_TANGENTIAL_FORCE_Y:
        retVal =
            pinfo->_contact_force(1) - (pinfo->_contact_force * pinfo->_normal) * pinfo->_normal(1);
        break;
      case PA_TANGENTIAL_FORCE_Z:
//.........这里部分代码省略.........

Real
MultiDContactConstraint::computeQpResidual(Moose::ConstraintType type)
{
  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
  const Node * node = pinfo->_node;

  RealVectorValue res_vec;
  // Build up residual vector
  for (unsigned int i=0; i<_mesh_dimension; ++i)
  {
    dof_id_type dof_number = node->dof_number(0, _vars(i), 0);
    res_vec(i) = _residual_copy(dof_number);
  }

  const RealVectorValue distance_vec(_mesh.node(node->id()) - pinfo->_closest_point);
  const RealVectorValue pen_force(_penalty * distance_vec);
  Real resid = 0;

  switch (type)
  {
  case Moose::Slave:
    switch (_model)
    {
    case CM_FRICTIONLESS:

      resid = pinfo->_normal(_component) * (pinfo->_normal * ( pen_force - res_vec ));
      break;

    case CM_GLUED:

      resid = pen_force(_component)
        - res_vec(_component)
        ;

      break;

    default:
      mooseError("Invalid or unavailable contact model");
      break;
    }
    return _test_slave[_i][_qp] * resid;
  case Moose::Master:
    switch (_model)
    {
    case CM_FRICTIONLESS:

      resid = pinfo->_normal(_component) * (pinfo->_normal * res_vec);
      break;

    case CM_GLUED:
      resid = res_vec(_component);
      break;

    default:
      mooseError("Invalid or unavailable contact model");
      break;
    }
    return _test_master[_i][_qp] * resid;
  }
  return 0;
}