[kernel] correct implementation for two OSI
@@ -126,36 +126,59 @@ void MoreauJeanGOSI::fillDSLinks(Interaction &inter,
bool computeResidu = relation.requireResidu();
inter.initializeMemory(computeResidu,_steps);
- if (!(checkOSI(DSG.descriptor(ds1)) && checkOSI(DSG.descriptor(ds2))))
- {
- RuntimeException::selfThrow("MoreauJeanGOSI::fillDSLinks. The implementation is not correct for two different OSI for one interaction");
- }
+ // if (!(checkOSI(DSG.descriptor(ds1)) && checkOSI(DSG.descriptor(ds2))))
+ // {
+ // RuntimeException::selfThrow("MoreauJeanGOSI::fillDSLinks. 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;
+ unsigned int xfree =0;
if (relationType == Lagrangian)
{
- DSlink[LagrangianR::xfree].reset(new BlockVector());
- DSlink[LagrangianR::xfree]->insertPtr(workVds1[OneStepIntegrator::free]);
+ xfree = LagrangianR::xfree;
}
else if (relationType == NewtonEuler)
{
- DSlink[NewtonEulerR::xfree].reset(new BlockVector());
- DSlink[NewtonEulerR::xfree]->insertPtr(workVds1[OneStepIntegrator::free]);
+ xfree = NewtonEulerR::xfree;
}
-
if (ds1 != ds2)
{
- VectorOfVectors &workVds2 = *DSG.properties(DSG.descriptor(ds2)).workVectors;
- if (relationType == Lagrangian)
+ DEBUG_PRINT("ds1 != ds2\n");
+ if ((!DSlink[xfree]) || (DSlink[xfree]->numberOfBlocks() !=2 ))
+ DSlink[xfree].reset(new BlockVector(2));
+ }
+ else
{
- DSlink[LagrangianR::xfree]->insertPtr(workVds2[OneStepIntegrator::free]);
+ if ((!DSlink[xfree]) || (DSlink[xfree]->numberOfBlocks() !=1 ))
+ DSlink[xfree].reset(new BlockVector(1));
}
- else if (relationType == NewtonEuler)
+
+ 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;
+ DSlink[xfree]->setVectorPtr(0,workVds1[OneStepIntegrator::free]);
+ }
+ DEBUG_PRINTF("ds1->number() %i\n",ds1->number());
+ DEBUG_PRINTF("ds2->number() %i\n",ds2->number());
+
+
+ if (ds1 != ds2)
{
- DSlink[NewtonEulerR::xfree]->insertPtr(workVds2[OneStepIntegrator::free]);
+ 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;
+ DSlink[xfree]->setVectorPtr(1,workVds2[OneStepIntegrator::free]);
}
}
+
}
void MoreauJeanGOSI::initialize_nonsmooth_problems()
