Element example

Description

This element is for a rigid link utilizing a Lagrange multiplier for enforcing the condition.

Code Listing

#include <Error_messager.h>

#include <Behavior.h>
#include <External_parameter.h>
#include <P_element.h>
#include <Global_matrix.h>
#include <Mesh.h>
#include <Utility_mesh.h>
#include <Node.h>
#include <Special_vector.h>
#include <GeomSpace.h>

#include <Print.h>

class RIGID_LINK : public D_ELEMENT {
  public :
    int nb_displ;
    virtual int dof_location(int)const;

    RIGID_LINK() { }
    virtual void setup_dofs(const char* key);

    virtual ~RIGID_LINK();
    virtual INTEGRATION_RESULT* internal_reaction(bool,VECTOR&,SMATRIX&,bool only_get_tg_matrix=FALSE);
};

INSTALL_ELEMENT(MECHANICAL_MESH,rigid_link,_2D,ST_LINE,RIGID_LINK)
INSTALL_ELEMENT(MECHANICAL_MESH,rigid_link,_3D,ST_LINE,RIGID_LINK)

int RIGID_LINK::dof_location(int idof)const
{ return (idof<nb_displ) ? DOF_NODE_ID : DOF_ELEMENT_ID;
}

void RIGID_LINK::setup_dofs(const char*)
{
  size_dofs(!node*geometry->space_dimension()+1);

  static DOF_TYPE u_name[3]={DOF::translate("U1",DOF_NODE_ID),
                             DOF::translate("U2",DOF_NODE_ID),
                             DOF::translate("U3",DOF_NODE_ID)};
  static DOF_TYPE RLF    = DOF::translate("RLF",DOF_ELEMENT_ID);

  int tot=0;
  for(int in=0;in<nb_node();in++) {
       for(int idime=0;idime<space_dimension();idime++) {
            dof_type[tot] = u_name[idime];
            map_dof_to_node[tot] = get_node(in);
            tot++;
       }
  }
  nb_displ = tot;

  dof_type[tot]=RLF;
  map_dof_to_node[tot] = NULL;
}

RIGID_LINK::~RIGID_LINK() {}

INTEGRATION_RESULT* RIGID_LINK::internal_reaction(bool ics,
       VECTOR& resi, SMATRIX& stiff,bool check)
{
  MATRIX elem_coord;
  VECTOR Dof(nb_dof()), dDof(nb_dof());
  int end = nb_dof()-1;
  get_elem_d_dof_tot(Dof);
  double force = Dof[end];

  VECTOR x1  = get_node(0)->get_coordw(END_OF_INCREMENT);
  VECTOR x2  = get_node(1)->get_coordw(END_OF_INCREMENT);

  VECTOR X1  = get_node(0)->get_coordw(BEGINNING_OF_PROBLEM);
  VECTOR X2  = get_node(1)->get_coordw(BEGINNING_OF_PROBLEM);

  VECTOR dX = X2-X1;
  VECTOR dx = x2-x1;
  double R0 = sqrt(dX|dX);
  double R1 = sqrt(dx|dx);
  VECTOR B  = dx*(1.0/R1);

  SMATRIX  dBdX  = B^B;
           dBdX.add_to_diagonal(1.0);
           dBdX *= 1.0/R1;
  SMATRIX mdBdX  = -1.0*dBdX;

  VECTOR subF1, subF2;
  int sd=space_dimension();

  subF1.reassign(sd,resi,0);  subF1= B*(-force);
  subF2.reassign(sd,resi,sd); subF2= B*force;
  resi[end] = R1 - R0;

  if (ics) {
     stiff = 0.0;
     MATRIX subS11(sd,sd,stiff,0,0);      subS11 = mdBdX*(-force);
     MATRIX subS12(sd,sd,stiff,0,sd);     subS12 =  dBdX*(-force);
     MATRIX subS21(sd,sd,stiff,sd,0);     subS21 = mdBdX*(force);
     MATRIX subS22(sd,sd,stiff,sd,sd);    subS22 =  dBdX*(force);

     for (int c=0;c<sd;c++) {
          stiff(c,   end) = -B[c];
          stiff(sd+c,end) =  B[c];

          stiff(end,c)    =  (x2[c]-x1[c])*(-1.0/R1);
          stiff(end,sd+c) =  (x2[c]-x1[c])*(1.0/R1);
     }

     stiff(end,end) = 0.0;
  }

  return NULL;
}