applications/solvers/solidMechanics/elasticIncrAcpSolidFoam/updateCrack.H

   1 nFacesToBreak = 0;
   2 nCoupledFacesToBreak = 0;
   3 {
   4     // Check internal faces
   5
   6     // scalarField effTraction =
   7     //   cohesiveZone.internalField() *
   8     //   mag(traction.internalField());
   9     scalarField normalTraction =
  10         cohesiveZone.internalField()*
  11         ( n.internalField() & traction.internalField() );
  12
  13     // only consider tensile tractions
  14     normalTraction = max(normalTraction, scalar(0));
  15     scalarField shearTraction =
  16         cohesiveZone.internalField()*
  17         mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
  18
  19     // the traction fraction is monitored to decide which faces to break:
  20     // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
  21
  22     const surfaceScalarField sigmaMax = rheology.cohLaw().sigmaMax();
  23     const surfaceScalarField tauMax = rheology.cohLaw().tauMax();
  24
  25     const scalarField& sigmaMaxI = sigmaMax.internalField();
  26     const scalarField& tauMaxI = tauMax.internalField();
  27
  28     //scalarField effTractionFraction = effTraction/sigmaMax;
  29     // scalarField effTractionFraction =
  30     //   (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
  31     scalarField effTractionFraction(normalTraction.size(), 0.0);
  32
  33     if(cohesivePatchDUPtr)
  34     {
  35         effTractionFraction =
  36             (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
  37           + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
  38     }
  39     else
  40     {
  41         // solidCohesiveFixedModeMix only uses sigmaMax
  42         effTractionFraction =
  43             (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
  44           + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
  45     }
  46
  47     maxEffTractionFraction = gMax(effTractionFraction);
  48
  49     SLList<label> facesToBreakList;
  50     SLList<scalar> facesToBreakEffTractionFractionList;
  51
  52     forAll(effTractionFraction, faceI)
  53     {
  54         if (effTractionFraction[faceI] > 1.0)
  55         {
  56             facesToBreakList.insert(faceI);
  57             facesToBreakEffTractionFractionList.insert
  58             (
  59                 effTractionFraction[faceI]
  60             );
  61         }
  62     }
  63
  64     labelList facesToBreak(facesToBreakList);
  65     List<scalar> facesToBreakEffTractionFraction
  66     (
  67         facesToBreakEffTractionFractionList
  68     );
  69
  70     nFacesToBreak = facesToBreak.size();
  71
  72     // Break only one face per topo change
  73     if (nFacesToBreak > 1)
  74     {
  75         nFacesToBreak = 1;
  76     }
  77
  78     // philipc - select face with maximum effective traction fraction
  79     label faceToBreakIndex = -1;
  80     scalar faceToBreakEffTractionFraction = 0;
  81     forAll(facesToBreakEffTractionFraction, faceI)
  82     {
  83         if
  84         (
  85             facesToBreakEffTractionFraction[faceI]
  86           > faceToBreakEffTractionFraction
  87         )
  88         {
  89             faceToBreakEffTractionFraction =
  90                 facesToBreakEffTractionFraction[faceI];
  91             faceToBreakIndex = facesToBreak[faceI];
  92         }
  93     }
  94
  95     scalar gMaxEffTractionFraction =
  96         returnReduce(faceToBreakEffTractionFraction, maxOp<scalar>());
  97
  98     if (Pstream::parRun())
  99     {
 100         bool procHasFaceToBreak = false;
 101         if (nFacesToBreak > 0)
 102         {
 103             if
 104             (
 105                 mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction)
 106               < SMALL
 107             )
 108             {
 109                 // philipc - Maximum traction fraction is on this processor
 110                 procHasFaceToBreak = true;
 111             }
 112         }
 113
 114         // Check if maximum is present on more then one processors
 115         label procID = Pstream::nProcs();
 116         if (procHasFaceToBreak)
 117         {
 118             procID = Pstream::myProcNo();
 119         }
 120
 121         label minProcID =
 122             returnReduce<label>(procID, minOp<label>());
 123
 124         if (procID != minProcID)
 125         {
 126             nFacesToBreak = 0;
 127         }
 128     }
 129
 130
 131     // Check coupled (processor) patches
 132
 133     SLList<label> coupledFacesToBreakList;
 134     SLList<scalar> coupledFacesToBreakEffTractionFractionList;
 135     forAll(mesh.boundary(), patchI)
 136     {
 137         if (mesh.boundary()[patchI].coupled())
 138         {
 139             // scalarField pEffTraction =
 140             //   cohesiveZone.boundaryField()[patchI]*
 141             //   mag(traction.boundaryField()[patchI]);
 142             // scalarField pEffTractionFraction = pEffTraction/sigmaMax.boundaryField()[patchI];
 143
 144             scalarField pNormalTraction =
 145                 cohesiveZone.boundaryField()[patchI]*
 146                 ( n.boundaryField()[patchI] & traction.boundaryField()[patchI] );
 147
 148             // only consider tensile tractions
 149             pNormalTraction = max(pNormalTraction, scalar(0));
 150
 151             scalarField pShearTraction =
 152                 cohesiveZone.boundaryField()[patchI]*
 153                 mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );
 154
 155             // the traction fraction is monitored to decide which faces to break:
 156             // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
 157             const scalarField& pSigmaMax = sigmaMax.boundaryField()[patchI];
 158             const scalarField& pTauMax = tauMax.boundaryField()[patchI];
 159
 160             // scalarField pEffTractionFraction =
 161             //     (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
 162             scalarField pEffTractionFraction(pNormalTraction.size(), 0.0);
 163             if(cohesivePatchDUPtr)
 164             {
 165                 pEffTractionFraction =
 166                     (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
 167                   + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
 168             }
 169             else
 170             {
 171                 // solidCohesiveFixedModeMix only uses sigmaMax
 172                 pEffTractionFraction =
 173                     (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
 174                   + (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
 175             }
 176
 177             label start = mesh.boundaryMesh()[patchI].start();
 178
 179             forAll(pEffTractionFraction, faceI)
 180             {
 181                 if (pEffTractionFraction[faceI] > maxEffTractionFraction)
 182                 {
 183                     maxEffTractionFraction = pEffTractionFraction[faceI];
 184                 }
 185
 186                 if (pEffTractionFraction[faceI] > 1.0)
 187                 {
 188                     //Pout << "coupled face to break " << faceI << endl;
 189                     coupledFacesToBreakList.insert(start + faceI);
 190                     coupledFacesToBreakEffTractionFractionList.insert
 191                     (
 192                         pEffTractionFraction[faceI]
 193                     );
 194                 }
 195             }
 196         }
 197     }
 198
 199     labelList coupledFacesToBreak(coupledFacesToBreakList);
 200     List<scalar> coupledFacesToBreakEffTractionFraction
 201     (
 202         coupledFacesToBreakEffTractionFractionList
 203     );
 204
 205     nCoupledFacesToBreak = coupledFacesToBreak.size();
 206 //     Pout << "nCoupledFacesToBreak: " << nCoupledFacesToBreak << endl;
 207
 208     // Break only one face per topo change
 209     if (nCoupledFacesToBreak > 1)
 210     {
 211         nCoupledFacesToBreak = 1;
 212     }
 213
 214     // Select coupled face with maximum effective traction fraction
 215     label coupledFaceToBreakIndex = -1;
 216     scalar coupledFaceToBreakEffTractionFraction = 0;
 217     forAll(coupledFacesToBreakEffTractionFraction, faceI)
 218     {
 219         if
 220         (
 221             coupledFacesToBreakEffTractionFraction[faceI]
 222           > coupledFaceToBreakEffTractionFraction
 223         )
 224         {
 225             coupledFaceToBreakEffTractionFraction =
 226                 coupledFacesToBreakEffTractionFraction[faceI];
 227             coupledFaceToBreakIndex = coupledFacesToBreak[faceI];
 228         }
 229     }
 230
 231     scalar gMaxCoupledEffTractionFraction =
 232         returnReduce(coupledFaceToBreakEffTractionFraction, maxOp<scalar>());
 233
 234     if (Pstream::parRun())
 235     {
 236
 237         bool procHasCoupledFaceToBreak = false;
 238         if (nCoupledFacesToBreak > 0)
 239         {
 240             if
 241             (
 242                 mag(gMaxCoupledEffTractionFraction - coupledFaceToBreakEffTractionFraction)
 243               < SMALL
 244             )
 245             {
 246                 // Maximum traction fraction is on this processor
 247                 procHasCoupledFaceToBreak = true;
 248             }
 249         }
 250
 251         // Check if maximum is present on more then one processors
 252         label procID = Pstream::nProcs();
 253         if (procHasCoupledFaceToBreak)
 254         {
 255             procID = Pstream::myProcNo();
 256         }
 257
 258         label minProcID =
 259             returnReduce<label>(procID, minOp<label>());
 260
 261         if (procID != minProcID)
 262         {
 263             nCoupledFacesToBreak = 0;
 264         }
 265     }
 266
 267     if (gMaxCoupledEffTractionFraction > gMaxEffTractionFraction)
 268     {
 269         // Break coupled face
 270         nFacesToBreak = 0;
 271     }
 272     else
 273     {
 274         // Break internal face
 275         nCoupledFacesToBreak = 0;
 276     }
 277
 278     // Make sure that coupled faces are broken in pairs
 279     labelList ngbProc(Pstream::nProcs(), -1);
 280     labelList index(Pstream::nProcs(), -1);
 281     if (nCoupledFacesToBreak)
 282     {
 283         label patchID =
 284             mesh.boundaryMesh().whichPatch(coupledFaceToBreakIndex);
 285
 286         label start = mesh.boundaryMesh()[patchID].start();
 287         label localIndex = coupledFaceToBreakIndex - start;
 288
 289         const processorPolyPatch& procPatch =
 290             refCast<const processorPolyPatch>(mesh.boundaryMesh()[patchID]);
 291         label ngbProcNo = procPatch.neighbProcNo();
 292
 293         ngbProc[Pstream::myProcNo()] = ngbProcNo;
 294         index[Pstream::myProcNo()] = localIndex;
 295     }
 296
 297     if (returnReduce(nCoupledFacesToBreak, maxOp<label>()))
 298     {
 299         reduce(ngbProc, maxOp<labelList>());
 300         reduce(index, maxOp<labelList>());
 301
 302         for (label procI = 0; procI < Pstream::nProcs(); procI++)
 303         {
 304             if (procI != Pstream::myProcNo())
 305             {
 306                 if (ngbProc[procI] == Pstream::myProcNo())
 307                 {
 308                     forAll(mesh.boundaryMesh(), patchI)
 309                     {
 310                         if
 311                         (
 312                             mesh.boundaryMesh()[patchI].type()
 313                          == processorPolyPatch::typeName
 314                         )
 315                         {
 316                             const processorPolyPatch& procPatch =
 317                                 refCast<const processorPolyPatch>
 318                                 (
 319                                     mesh.boundaryMesh()[patchI]
 320                                 );
 321                             label ngbProcNo = procPatch.neighbProcNo();
 322
 323                             if (ngbProcNo == procI)
 324                             {
 325                                 label start =
 326                                     mesh.boundaryMesh()[patchI].start();
 327                                 coupledFaceToBreakIndex = start + index[procI];
 328                                 nCoupledFacesToBreak = 1;
 329                             }
 330                         }
 331                     }
 332                 }
 333             }
 334         }
 335     }
 336
 337
 338     vector faceToBreakTraction = vector::zero;
 339     vector faceToBreakNormal = vector::zero;
 340     scalar faceToBreakSigmaMax = 0.0;
 341     scalar faceToBreakTauMax = 0.0;
 342
 343     // Set faces to break
 344     if (nFacesToBreak > 0)
 345     {
 346         faceToBreakTraction = traction.internalField()[faceToBreakIndex];
 347         faceToBreakNormal = n.internalField()[faceToBreakIndex];
 348
 349         // Scale broken face traction
 350         // The scale factor is derived by solving the following eqn for alpha:
 351         // (alpha*tN/tNC)^2 + (alpha*tS/tSC)^2 = 1
 352         faceToBreakSigmaMax = sigmaMaxI[faceToBreakIndex];
 353         faceToBreakTauMax = tauMaxI[faceToBreakIndex];
 354         scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
 355         normalTrac = max(normalTrac, 0.0);
 356         scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
 357         scalar scaleFactor = 1;
 358         if(cohesivePatchDUPtr)
 359         {
 360             scaleFactor =
 361                 Foam::sqrt
 362                 (
 363                     1 /
 364                     (
 365                         (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
 366                       + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
 367                     )
 368                 );
 369         }
 370         else
 371         {
 372             // solidCohesiveFixedModeMix only uses sigmaMax
 373             scaleFactor =
 374                 Foam::sqrt
 375                 (
 376                     1 /
 377                     (
 378                         (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
 379                       + (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
 380                     )
 381                 );
 382         }
 383
 384         faceToBreakTraction *= scaleFactor;
 385
 386         topoChange = true;
 387     }
 388     else if (nCoupledFacesToBreak > 0)
 389     {
 390         label patchID =
 391             mesh.boundaryMesh().whichPatch(coupledFaceToBreakIndex);
 392         label start = mesh.boundaryMesh()[patchID].start();
 393         label localIndex = coupledFaceToBreakIndex - start;
 394
 395         faceToBreakTraction = traction.boundaryField()[patchID][localIndex];
 396         faceToBreakNormal = n.boundaryField()[patchID][localIndex];
 397
 398         // Scale broken face traction
 399         faceToBreakSigmaMax = sigmaMax.boundaryField()[patchID][localIndex];
 400         faceToBreakTauMax = tauMax.boundaryField()[patchID][localIndex];
 401         scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
 402         normalTrac = max(normalTrac, 0.0);
 403         scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
 404         scalar scaleFactor = 1;
 405         if(cohesivePatchDUPtr)
 406         {
 407             scaleFactor =
 408                 Foam::sqrt
 409                 (
 410                     1 /
 411                     (
 412                         (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
 413                       + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
 414                     )
 415                 );
 416         }
 417         else
 418         {
 419             // solidCohesiveFixedModeMix only uses sigmaMax
 420             scaleFactor =
 421                 Foam::sqrt
 422                 (
 423                     1 /
 424                     (
 425                         (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
 426                       + (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
 427                     )
 428                 );
 429         }
 430
 431         faceToBreakTraction *= scaleFactor;
 432
 433         topoChange = true;
 434     }
 435
 436     reduce(topoChange, orOp<bool>());
 437
 438     labelList faceToBreak(nFacesToBreak, faceToBreakIndex);
 439     boolList faceToBreakFlip(nFacesToBreak, false);
 440     labelList coupledFaceToBreak
 441     (
 442         nCoupledFacesToBreak,
 443         coupledFaceToBreakIndex
 444     );
 445
 446     reduce(nFacesToBreak, maxOp<label>());
 447     reduce(nCoupledFacesToBreak, maxOp<label>());
 448
 449     if (nFacesToBreak || nCoupledFacesToBreak)
 450     {
 451         Pout << "Internal face to break: " << faceToBreak << endl;
 452         Pout << "Coupled face to break: " << coupledFaceToBreak << endl;
 453
 454         mesh.setBreak(faceToBreak, faceToBreakFlip, coupledFaceToBreak);
 455         mesh.update();
 456
 457         const labelList& faceMap = mesh.topoChangeMap().faceMap();
 458         label start = mesh.boundaryMesh()[cohesivePatchID].start();
 459
 460         mu = rheology.mu();
 461         lambda = rheology.lambda();
 462         muf = fvc::interpolate(mu);
 463         lambdaf = fvc::interpolate(lambda);
 464
 465         // we need to modify propertiess after cracking otherwise momentum equation is wrong
 466         // but solidInterface seems to hold some information about old mesh
 467         // so we will delete it and make another
 468         // we could probably add a public clearout function
 469         // create new solidInterface
 470         if(rheology.solidInterfaceActive())
 471         {
 472             rheology.solInterface().clearOut();
 473             solidInterfacePtr->modifyProperties(muf, lambdaf);
 474         }
 475
 476         // All values on the new crack faces get set to zero
 477         // so we must manually correct them
 478         const vectorField DUpI =
 479             DU.boundaryField()[cohesivePatchID].patchInternalField();
 480         const vectorField oldDUpI =
 481             DU.oldTime().boundaryField()[cohesivePatchID].patchInternalField();
 482         const vectorField UpI =
 483             U.boundaryField()[cohesivePatchID].patchInternalField();
 484         const vectorField oldUpI =
 485             U.oldTime().boundaryField()[cohesivePatchID].patchInternalField();
 486         const symmTensorField sigmapI =
 487             sigma.boundaryField()[cohesivePatchID].patchInternalField();
 488         const scalarField muPI = mu.boundaryField()[cohesivePatchID].patchInternalField();
 489         const scalarField lambdaPI = lambda.boundaryField()[cohesivePatchID].patchInternalField();
 490
 491         // Global crack fields
 492         const vectorField globalDUpI = mesh.globalCrackField(DUpI);
 493         const vectorField globalOldDUpI = mesh.globalCrackField(oldDUpI);
 494         const vectorField globalUpI = mesh.globalCrackField(UpI);
 495         const vectorField globalOldUpI = mesh.globalCrackField(oldUpI);
 496         const symmTensorField globalsigmapI = mesh.globalCrackField(sigmapI);
 497         const scalarField globalMuPI = mesh.globalCrackField(muPI);
 498         const scalarField globalLambdaPI = mesh.globalCrackField(lambdaPI);
 499
 500         // cohesivePatchU.size()
 501         int cohesivePatchSize(cohesivePatchDUPtr ? cohesivePatchDUPtr->size() : cohesivePatchDUFixedModePtr->size());
 502
 503         // Initialise fields for new cohesive face
 504         const labelList& gcfa = mesh.globalCrackFaceAddressing();
 505         label globalIndex = mesh.localCrackStart();
 506         //for (label i=0; i<cohesivePatchDU.size(); i++)
 507         for (label i=0; i<cohesivePatchSize; i++)
 508         {
 509             label oldFaceIndex = faceMap[start+i];
 510
 511             // If new face
 512             if (oldFaceIndex == faceToBreakIndex)
 513             {
 514                 // set to average of old cell centres
 515                 // hmnnn it would be better to interpolate
 516                 // using weights... OK for now: future work
 517                 DU.boundaryField()[cohesivePatchID][i] =
 518                     0.5
 519                    *(
 520                        globalDUpI[globalIndex]
 521                      + globalDUpI[gcfa[globalIndex]]
 522                     );
 523                 DU.oldTime().boundaryField()[cohesivePatchID][i] =
 524                     0.5
 525                    *(
 526                        globalOldDUpI[globalIndex]
 527                      + globalOldDUpI[gcfa[globalIndex]]
 528                     );
 529                 U.boundaryField()[cohesivePatchID][i] =
 530                     0.5
 531                    *(
 532                        globalUpI[globalIndex]
 533                      + globalUpI[gcfa[globalIndex]]
 534                     );
 535                 U.oldTime().boundaryField()[cohesivePatchID][i] =
 536                     0.5
 537                    *(
 538                        globalOldUpI[globalIndex]
 539                      + globalOldUpI[gcfa[globalIndex]]
 540                     );
 541                 sigma.boundaryField()[cohesivePatchID][i] =
 542                     0.5
 543                    *(
 544                        globalsigmapI[globalIndex]
 545                      + globalsigmapI[gcfa[globalIndex]]
 546                     );
 547
 548                 // initialise mu and lambda on new faces
 549                 // set new face value to value of internal cell
 550                 muf.boundaryField()[cohesivePatchID][i] = globalMuPI[globalIndex];
 551                 lambdaf.boundaryField()[cohesivePatchID][i] = globalLambdaPI[globalIndex];
 552
 553                 globalIndex++;
 554             }
 555             else
 556             {
 557                 globalIndex++;
 558             }
 559         }
 560
 561         // we must calculate grad using interface
 562         // DU at the interface has not been calculated yet as interface.correct()
 563         // has not been called yet
 564         // not really a problem as gradDU is correct in second outer iteration
 565         // as long as this does not cause convergence problems for the first iterations.
 566         // we should be able to calculate the interface displacements without
 567         // having to call interface.correct()
 568         // todo: add calculateInterfaceDU() function
 569         // interface grad uses Gauss, we need least squares
 570         gradDU = fvc::grad(DU); // leastSquaresSolidInterface grad scheme
 571         //gradDU = solidInterfacePtr->grad(DU);
 572         //snGradDU = fvc::snGrad(DU);
 573
 574 #       include "calculateTraction.H"
 575         //if (nFacesToBreak || nCoupledFacesToBreak) mesh.write(); traction.write();
 576
 577         // Initialise initiation traction for new cohesive patch face
 578         // we also need to update the traction_ field in the crack boundary condition
 579         // this is because it cannot set itself during mapping.
 580         //for (label i=0; i<cohesivePatchDU.size(); i++)
 581         for (label i=0; i<cohesivePatchSize; i++)
 582         {
 583             label oldFaceIndex = faceMap[start+i];
 584
 585             // If new face
 586             if
 587             (
 588                 (oldFaceIndex == faceToBreakIndex)
 589              || (oldFaceIndex == coupledFaceToBreakIndex)
 590             )
 591             {
 592                 vector n0 =
 593                     mesh.Sf().boundaryField()[cohesivePatchID][i]
 594                    /mesh.magSf().boundaryField()[cohesivePatchID][i];
 595                 //vector n1 = -n0;
 596
 597                 if ((n0 & faceToBreakNormal) > SMALL)
 598                 {
 599                     traction.boundaryField()[cohesivePatchID][i] =
 600                         faceToBreakTraction;
 601
 602                     traction.oldTime().boundaryField()[cohesivePatchID][i] =
 603                         faceToBreakTraction;
 604
 605                     // this seems to slow convergence in some simple test cases...
 606                     // but surely it should be better update it
 607                     //cohesivePatchDU.traction()[i] = faceToBreakTraction;
 608                     if(cohesivePatchDUPtr)
 609                     {
 610                         cohesivePatchDUPtr->traction()[i] = faceToBreakTraction;
 611                     }
 612                     else
 613                     {
 614                         cohesivePatchDUFixedModePtr->traction()[i] = faceToBreakTraction;
 615                         cohesivePatchDUFixedModePtr->initiationTraction()[i] = faceToBreakTraction;
 616                     }
 617                 }
 618                 else
 619                 {
 620                     traction.boundaryField()[cohesivePatchID][i] =
 621                         -faceToBreakTraction;
 622                     traction.oldTime().boundaryField()[cohesivePatchID][i] =
 623                         -faceToBreakTraction;
 624
 625                     //cohesivePatchDU.traction()[i] = -faceToBreakTraction;
 626                     if(cohesivePatchDUPtr)
 627                     {
 628                         cohesivePatchDUPtr->traction()[i] = -faceToBreakTraction;
 629                     }
 630                     else
 631                     {
 632                         cohesivePatchDUFixedModePtr->traction()[i] =
 633                             -faceToBreakTraction;
 634                         cohesivePatchDUFixedModePtr->initiationTraction()[i] =
 635                             -faceToBreakTraction;
 636                     }
 637                 }
 638             }
 639         }
 640
 641         // hmmnn we only need a reference for very small groups of cells
 642         // turn off for now
 643         //#       include "updateReference.H"
 644     }
 645 }