contactTerms.cpp

//
// contact base term
//
// Description: contact term
//
//
// Author:  <Gauthier BECKER>, (C) 2011
//
// Copyright: See COPYING file that comes with this distribution
//
#include "contactTerms.h"
#include "unknownField.h"
template<> void rigidContactLinearTermBase<double>::get(MElement *ele, int npts, IntPt *GP,fullVector<double> &m) const
{
  const std::vector<int>& allComps = _spc->getConstRefToComp(); // all components of fields with contact
  nbcomp =  allComps.size(); 
  nbdofs = _spc->getNumKeys(ele); // total number of keys = keys in elements + number of rigid keys
  nbdofsGC = _spc->getNumKeysOfGravityCenter(); // number of rigid keys
  //
  m.resize(nbdofs);
  m.setAll(0.);
  
  nbdofsEle = nbdofs - nbdofsGC; 
  // get displacement of vertices (To know the current node positions)
  disp.clear(); R.clear();
  _spc->getKeys(ele,R);
  _ufield->get(R,disp);
  //
  verdisp[3] = disp[nbdofsEle];
  verdisp[4] = disp[nbdofsEle+1];
  verdisp[5] = disp[nbdofsEle+2];

  nbvertex = ele->getNumVertices();  
  for(int i=0;i<nbvertex;i++)
  {
    
    verdisp[0] = disp[i];
    verdisp[1] = disp[i+nbvertex];
    verdisp[2] = disp[i+2*nbvertex];
    
    for(int j=0;j<6;j++) 
      mvalue[j] =0.;
    
    ver = ele->getVertex(i);
    this->get(ver,verdisp,mvalue);
    // look if there is contact force
    double val = mvalue[0]*mvalue[0]+mvalue[1]*mvalue[1]+mvalue[2]*mvalue[2];
    if (val > 0){
      // in conact
      _allcontactNode->insert(ele,i,_lastNormalContact);
      for(int j=0;j<nbcomp;j++){
        m(i+j*nbvertex) += mvalue[j];
        m(nbdofsEle+j) += mvalue[j+3]; // count two times
      }
    }
  }
 // if(m.norm() != 0.)
//    m.print("contact");
}

template<> void defoDefoContactBilinearTermByPerturbation<double>::get(contactElement *cE, int npts, IntPt *GP, fullMatrix<double> &m) const
{
  int nbdofs = _spc->getNumKeys(cE);
  m.resize(nbdofs,nbdofs);
  m.setAll(0.); //here we need to fill slave node 1, x, slave node 2, x,slave node 3, x, ... slave node 1, y etc and then master node 1, x, master node 2, x, etc
  static std::vector<Dof> keys;
  keys.clear();
  _spc->getKeys(cE,keys);
  static fullVector<double> disp,Fplus,Fminus;
  disp.resize(nbdofs);
  disp.setAll(0.);
  _ufield->get(keys,disp); 
  defoDefoContactLinearTermBase<double> *defoDefoTerm = static_cast<defoDefoContactLinearTermBase<double> *> (_lterm);


  for(int j=0; j<nbdofs; j++)
  {
    disp(j)+=_pert;
    //defoDefoTerm->detectNodesInContact(cE,disp);  //we should do the detection before the perturbation if we put the function out
    defoDefoTerm->get(cE,npts,GP,Fplus,disp);
    disp(j)-=2.*_pert;
    //defoDefoTerm->detectNodesInContact(cE,disp);  //we should do the detection before the perturbation
    defoDefoTerm->get(cE,npts,GP,Fminus,disp);
    disp(j)+=_pert;
    for(int i=0; i<nbdofs; i++)
    {
      m(i,j)=(Fplus(i)-Fminus(i))/(2.*_pert); //
    }
  }
  defoDefoTerm->get(cE,npts,GP,Fminus,disp);// put back the right _nodesInContact
  //m.print("numerical defo defo stiffness matrix");
}