C++ Interaction::linkToDSVariables方法代码示例

void NewtonEulerR::computeOutput(double time, Interaction& inter, unsigned int derivativeNumber)
{

  DEBUG_BEGIN("NewtonEulerR::computeOutput(...)\n");
  DEBUG_PRINTF("with time = %f and derivativeNumber = %i starts\n", time, derivativeNumber);

  VectorOfBlockVectors& DSlink = inter.linkToDSVariables();
  SiconosVector& y = *inter.y(derivativeNumber);
  BlockVector& q = *DSlink[NewtonEulerR::q0];


  if (derivativeNumber == 0)
  {
    computeh(time, q, y);
  }
  else
  {
    /* \warning  V.A. 15/04/2016
     * We decide finally not to update the Jacobian there. To be discussed
     */
    // computeJachq(time, inter, DSlink[NewtonEulerR::q0]);
    // computeJachqT(inter, DSlink[NewtonEulerR::q0]);

    if (derivativeNumber == 1)
    {
      assert(_jachqT);
      assert(DSlink[NewtonEulerR::velocity]);
      DEBUG_EXPR(_jachqT->display();); DEBUG_EXPR((*DSlink[NewtonEulerR::velocity]).display(););

void FirstOrderLinearR::initialize(Interaction& inter)
{

  FirstOrderR::initialize(inter);

  // get interesting size
  unsigned int sizeY = inter.dimension();
  unsigned int sizeX = inter.getSizeOfDS();

  VectorOfBlockVectors& DSlink = inter.linkToDSVariables();
  unsigned int sizeZ = DSlink[FirstOrderR::z]->size();
  VectorOfSMatrices& relationMat = inter.relationMatrices();
  VectorOfVectors & relationVec= inter.relationVectors();

  if (!_C && _pluginJachx->fPtr)
    relationMat[FirstOrderR::mat_C].reset(new SimpleMatrix(sizeY, sizeX));
  if (!_D && _pluginJachlambda->fPtr)
    relationMat[FirstOrderR::mat_D].reset(new SimpleMatrix(sizeY, sizeY));
  if (!_B && _pluginJacglambda->fPtr)
    relationMat[FirstOrderR::mat_B].reset(new SimpleMatrix(sizeX, sizeY));
  if (!_F && _pluginf->fPtr)
    relationMat[FirstOrderR::mat_F].reset(new SimpleMatrix(sizeY, sizeZ));
  if (!_e && _plugine->fPtr)
    relationVec[FirstOrderR::e].reset(new SiconosVector(sizeY));

  checkSize(inter);
}

void FirstOrderLinearR::computeInput(double time, Interaction& inter, unsigned int level)
{

  
  SiconosVector& lambda = *inter.lambda(level);
  VectorOfBlockVectors& DSlink = inter.linkToDSVariables();
  BlockVector& z = *DSlink[FirstOrderR::z];
  SP::SiconosVector z_vec(new SiconosVector(z));
  computeg(time, lambda, *z_vec, *DSlink[FirstOrderR::r]);
  *DSlink[FirstOrderR::z] = *z_vec;
}

void FirstOrderLinearR::checkSize(Interaction& inter)
{

  VectorOfBlockVectors& DSlink = inter.linkToDSVariables();

  // get interesting size
  unsigned int sizeY = inter.dimension();
  unsigned int sizeX = inter.getSizeOfDS();
  unsigned int sizeZ = DSlink[FirstOrderR::z]->size();

  // Check if various operators sizes are consistent.
  // Reference: interaction.

  if (_C)
  {
    if (_C->size(0) == 0)
      _C->resize(sizeX, sizeY);
    else
      assert((_C->size(0) == sizeY && _C->size(1) == sizeX) && "FirstOrderLinearR::initialize , inconsistent size between C and Interaction.");
  }
  if (_B)
  {
    if (_B->size(0) == 0)
    _B->resize(sizeY, sizeX);
  else
    assert((_B->size(1) == sizeY && _B->size(0) == sizeX) && "FirstOrderLinearR::initialize , inconsistent size between B and interaction.");
  }
  // C and B are the minimum inputs. The others may remain null.

  if (_D)
  {
    if (_D->size(0) == 0)
      _D->resize(sizeY, sizeY);
    else
      assert((_D->size(0) == sizeY || _D->size(1) == sizeY) && "FirstOrderLinearR::initialize , inconsistent size between C and D.");
  }

  if (_F)
  {
    if (_F->size(0) == 0)
      _F->resize(sizeY, sizeZ);
    else
      assert(((_F->size(0) == sizeY) && (_F->size(1) == sizeZ)) && "FirstOrderLinearR::initialize , inconsistent size between C and F.");
  }

  if (_e)
  {
    if (_e->size() == 0)
      _e->resize(sizeY);
    else
      assert(_e->size() == sizeY && "FirstOrderLinearR::initialize , inconsistent size between C and e.");
  }
}

void NewtonEulerR::initialize(Interaction& inter)
{


  DEBUG_BEGIN("NewtonEulerR::initialize(Interaction& inter)\n");

  unsigned int ySize = inter.dimension();
  unsigned int xSize = inter.getSizeOfDS();
  unsigned int qSize = 7 * (xSize / 6);

  if (!_jachq)
    _jachq.reset(new SimpleMatrix(ySize, qSize));
  else
  {
    if (_jachq->size(0) == 0)
    {
      // if the matrix dim are null
      _jachq->resize(ySize, qSize);
    }
    else
    {
      assert((_jachq->size(1) == qSize && _jachq->size(0) == ySize) ||
             (printf("NewtonEuler::initializeWorkVectorsAndMatrices _jachq->size(1) = %d ,_qsize = %d , _jachq->size(0) = %d ,_ysize =%d \n", _jachq->size(1), qSize, _jachq->size(0), ySize) && false) ||
             ("NewtonEuler::initializeWorkVectorsAndMatrices inconsistent sizes between _jachq matrix and the interaction." && false));
    }
  }

  DEBUG_EXPR(_jachq->display());

  if (! _jachqT)
    _jachqT.reset(new SimpleMatrix(ySize, xSize));

  if (! _T)
  {
    _T.reset(new SimpleMatrix(7, 6));
    _T->zero();
    _T->setValue(0, 0, 1.0);
    _T->setValue(1, 1, 1.0);
    _T->setValue(2, 2, 1.0);
  }
  DEBUG_EXPR(_jachqT->display());
  VectorOfBlockVectors& DSlink = inter.linkToDSVariables();
  if (!_contactForce)
  {
    _contactForce.reset(new SiconosVector(DSlink[NewtonEulerR::p1]->size()));
    _contactForce->zero();
  }
  DEBUG_END("NewtonEulerR::initialize(Interaction& inter)\n");
}

void FirstOrderLinearR::computeOutput(double time, Interaction& inter, unsigned int level)
{
  DEBUG_BEGIN("FirstOrderLinearR::computeOutput \n");
  VectorOfBlockVectors& DSlink = inter.linkToDSVariables();
  BlockVector& z = *DSlink[FirstOrderR::z];
  BlockVector& x = *DSlink[FirstOrderR::x];

  SP::SiconosVector z_vec(new SiconosVector(z));
  SiconosVector& y = *inter.y(level);
  SiconosVector& lambda = *inter.lambda(level);

  computeh(time, x, lambda, *z_vec, y);

  *DSlink[FirstOrderR::z] = *z_vec;

  DEBUG_END("FirstOrderLinearR::computeOutput \n");
}

void D1MinusLinearOSI::initializeWorkVectorsForInteraction(Interaction &inter,
				     InteractionProperties& interProp,
				     DynamicalSystemsGraph & DSG)
{

  DEBUG_BEGIN("D1MinusLinearOSI::initializeWorkVectorsForInteraction(Interaction &inter, InteractionProperties& interProp, DynamicalSystemsGraph & DSG)\n");
  SP::DynamicalSystem ds1= interProp.source;
  SP::DynamicalSystem ds2= interProp.target;
  assert(ds1);
  assert(ds2);
  DEBUG_PRINTF("interaction number %i\n", inter.number());
  VectorOfBlockVectors& DSlink = inter.linkToDSVariables();

  if (!interProp.workVectors)
  {
    interProp.workVectors.reset(new VectorOfVectors);
    interProp.workVectors->resize(D1MinusLinearOSI::WORK_INTERACTION_LENGTH);
  }
  if (!interProp.workBlockVectors)
  {
    interProp.workBlockVectors.reset(new VectorOfBlockVectors);
    interProp.workBlockVectors->resize(D1MinusLinearOSI::BLOCK_WORK_LENGTH);
  }

  VectorOfVectors& inter_work = *interProp.workVectors;
  VectorOfBlockVectors& inter_work_block = *interProp.workBlockVectors;

  Relation &relation =  *inter.relation();
  RELATION::TYPES relationType = relation.getType();

  inter_work[D1MinusLinearOSI::OSNSP_RHS].reset(new SiconosVector(inter.dimension()));

  // Check if interations levels (i.e. y and lambda sizes) are compliant with the current osi.
  _check_and_update_interaction_levels(inter);
  // Initialize/allocate memory buffers in interaction.
  inter.initializeMemory(_steps);

  if (!(checkOSI(DSG.descriptor(ds1)) && checkOSI(DSG.descriptor(ds2))))
  {
    std::cout << "checkOSI(DSG.descriptor(ds1)): "
              << std::boolalpha
              << checkOSI(DSG.descriptor(ds1)) << std::endl;
    std::cout << "checkOSI(DSG.descriptor(ds2)): "
              << std::boolalpha
              << checkOSI(DSG.descriptor(ds2)) << std::endl;


    RuntimeException::selfThrow("D1MinusLinearOSI::initializeWorkVectorsForInteraction. The implementation is not correct for two different OSI for one interaction");
  }


  /* allocate and set work vectors for the osi */
  unsigned int xfree = D1MinusLinearOSI::xfree;
  DEBUG_PRINTF("ds1->number() %i\n",ds1->number());
  DEBUG_PRINTF("ds2->number() %i\n",ds2->number());

  if (ds1 != ds2)
  {
    DEBUG_PRINT("ds1 != ds2\n");
    if ((!inter_work_block[xfree]) || (inter_work_block[xfree]->numberOfBlocks() !=2 ))
      inter_work_block[xfree].reset(new BlockVector(2));
  }
  else
  {
    if ((!inter_work_block[xfree]) || (inter_work_block[xfree]->numberOfBlocks() !=1 ))
      inter_work_block[xfree].reset(new BlockVector(1));
  }

  if(checkOSI(DSG.descriptor(ds1)))
  {
    DEBUG_PRINTF("ds1->number() %i is taken into account\n", ds1->number());
    assert(DSG.properties(DSG.descriptor(ds1)).workVectors);
    VectorOfVectors &workVds1 = *DSG.properties(DSG.descriptor(ds1)).workVectors;
    inter_work_block[xfree]->setVectorPtr(0,workVds1[D1MinusLinearOSI::FREE]);
  }
  if (ds1 != ds2)
  {
    DEBUG_PRINT("ds1 != ds2\n");
    if(checkOSI(DSG.descriptor(ds2)))
    {
      DEBUG_PRINTF("ds2->number() %i is taken into account\n",ds2->number());
      assert(DSG.properties(DSG.descriptor(ds2)).workVectors);
      VectorOfVectors &workVds2 = *DSG.properties(DSG.descriptor(ds2)).workVectors;
      inter_work_block[xfree]->setVectorPtr(1,workVds2[D1MinusLinearOSI::FREE]);
    }
  }

  DEBUG_EXPR(inter_work_block[xfree]->display(););

void SchatzmanPaoliOSI::initializeWorkVectorsForInteraction(Interaction &inter,
				      InteractionProperties& interProp,
				      DynamicalSystemsGraph & DSG)
{
  SP::DynamicalSystem ds1= interProp.source;
  SP::DynamicalSystem ds2= interProp.target;
  assert(ds1);
  assert(ds2);

  VectorOfBlockVectors& DSlink = inter.linkToDSVariables();

  if (!interProp.workVectors)
  {
    interProp.workVectors.reset(new VectorOfVectors);
    interProp.workVectors->resize(SchatzmanPaoliOSI::WORK_INTERACTION_LENGTH);
  }

  if (!interProp.workBlockVectors)
  {
    interProp.workBlockVectors.reset(new VectorOfBlockVectors);
    interProp.workBlockVectors->resize(SchatzmanPaoliOSI::BLOCK_WORK_LENGTH);
  }

  VectorOfVectors& inter_work = *interProp.workVectors;
  VectorOfBlockVectors& inter_work_block = *interProp.workBlockVectors;

  Relation &relation =  *inter.relation();

  inter_work[SchatzmanPaoliOSI::OSNSP_RHS].reset(new SiconosVector(inter.dimension()));


  RELATION::TYPES relationType = relation.getType();

  // Check if interations levels (i.e. y and lambda sizes) are compliant with the current osi.
  _check_and_update_interaction_levels(inter);
  // Initialize/allocate memory buffers in interaction.
  inter.initializeMemory(_steps);

  if (!(checkOSI(DSG.descriptor(ds1)) && checkOSI(DSG.descriptor(ds2))))
  {
    RuntimeException::selfThrow("SchatzmanPaoliOSI::initializeWorkVectorsForInteraction. The implementation is not correct for two different OSI for one interaction");
  }

  /* allocate and set work vectors for the osi */
  VectorOfVectors &workVds1 = *DSG.properties(DSG.descriptor(ds1)).workVectors;
  if (relationType == Lagrangian)
  {
    LagrangianDS& lds = *std11::static_pointer_cast<LagrangianDS> (ds1);
    DSlink[LagrangianR::p0].reset(new BlockVector());
    DSlink[LagrangianR::p0]->insertPtr(lds.p(0));

    inter_work_block[SchatzmanPaoliOSI::xfree].reset(new BlockVector());
    inter_work_block[SchatzmanPaoliOSI::xfree]->insertPtr(workVds1[SchatzmanPaoliOSI::FREE]);

  }
  else if (relationType == NewtonEuler)
  {
    inter_work_block[SchatzmanPaoliOSI::xfree].reset(new BlockVector());
    inter_work_block[SchatzmanPaoliOSI::xfree]->insertPtr(workVds1[SchatzmanPaoliOSI::FREE]);
  }

  if (ds1 != ds2)
  {
    VectorOfVectors &workVds2 = *DSG.properties(DSG.descriptor(ds2)).workVectors;
    if (relationType == Lagrangian)
    {
      inter_work_block[SchatzmanPaoliOSI::xfree]->insertPtr(workVds2[SchatzmanPaoliOSI::FREE]);
      LagrangianDS& lds = *std11::static_pointer_cast<LagrangianDS> (ds2);
      DSlink[LagrangianR::p0]->insertPtr(lds.p(0));
    }
    else if (relationType == NewtonEuler)
    {
      inter_work_block[SchatzmanPaoliOSI::xfree]->insertPtr(workVds2[SchatzmanPaoliOSI::FREE]);
    }
  }
}