MPC example

Description

This code imposes that a node set remains straight (mpc5).

Code Listing

class MPC5 : public RELATIONSHIP {
    NNSET* nset;
    ARRAY<double> eta;
    ARRAY<NODE*> nodes;
    NODE *n1, *n2;
  public :
    MPC5();
    virtual void initialize(ASCII_FILE&,MESH&);
    void set_relationship(MESH&);
};

// DECLARE_OBJECT(RELATIONSHIP,MPC5,mpc5);

MPC5::MPC5() {}

void MPC5::initialize(ASCII_FILE& file,MESH& mesh)
{ if(!mesh.areYouA("MECHANICAL_MESH"))
     INPUT_ERROR("MPC5 requires a mechanical mesh");
  nset=(NNSET*)mesh.find_nset(file.getSTRING());
  NNSET& ns = *nset;
  int dim = nset->get_dimension();

  if(!ns<2) INPUT_ERROR("The node set size is too small for: "+ns.get_name());
  VECTOR nn(dim), n1n2(dim);

  // compute coefficients
  n1=(NODE*)ns[0]; n2=(NODE*)ns[1];
  make_VECTOR(*n1,*n2,n1n2); double dist = norm(n1n2);
  eta.resize(!ns-2); nodes.resize(!ns-2);
  int count=0;
  for(int in1=2;in1<!ns;in1++) {
     //if(ns[in1]==n1 || ns[in1]==n2) continue;
     make_VECTOR(*n1,*ns[in1],nn); // d1=norm(nn);
     nodes[count]=(NODE*)ns[in1]; eta[count]=(n1n2|nn)/dist/dist;
     count++;
  }
}

void MPC5::set_relationship(MESH&)
{ static DOF_TYPE Displ_name[3]={DOF::translate("U1",DOF_NODE_ID),
                                 DOF::translate("U2",DOF_NODE_ID),
                                 DOF::translate("U3",DOF_NODE_ID)};
  ARRAY<DOF*> master_dof(2);
  ARRAY<double> coef(2); coef[0]=1.;
  DOF *slave;
  DOF_TYPE type;

  for(int idim=0;idim<nset->get_dimension();idim++) {
      type=Displ_name[idim];
      master_dof[0]=n1->get_dof(type);
      master_dof[1]=n2->get_dof(type);
      for(int in=0;in<!nodes;in++) {
          coef[0]=1.-eta[in]; coef[1]=eta[in];
          slave = nodes[in]->get_dof(type);
          slave->set_mpc(master_dof,coef,0.);
      }
  }
}