src/solidModels/arbitraryCrack/solidCohesiveFixedModeMix/solidCohesiveFixedModeMixFvPatchVectorField.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "solidCohesiveFixedModeMixFvPatchVectorField.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "volFields.H"
  29 #include "constitutiveModel.H"
  30 #include "regionSplit.H"
  31 #include "crackerFvMesh.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  39
  40 solidCohesiveFixedModeMixFvPatchVectorField::
  41 solidCohesiveFixedModeMixFvPatchVectorField
  42 (
  43     const fvPatch& p,
  44     const DimensionedField<vector, volMesh>& iF
  45 )
  46 :
  47     directionMixedFvPatchVectorField(p, iF),
  48     fieldName_("undefined"),
  49     relaxationFactor_(1.0),
  50     traction_(p.size(), vector::zero),
  51     initiationTraction_(p.size(), vector::zero),
  52     separationDistance_(p.size(), vector::zero),
  53     oldSeparationDistance_(p.size(), vector::zero),
  54     unloadingSeparationDistance_(p.size(), vector::zero),
  55     minUnloadingSeparationDistance_(0.0),
  56     beta_(0.0),
  57     explicitSeparationDistance_(false),
  58     contact_(false),
  59     contactConstant_(0),
  60     curTimeIndex_(-1)
  61 {}
  62
  63
  64 solidCohesiveFixedModeMixFvPatchVectorField::
  65 solidCohesiveFixedModeMixFvPatchVectorField
  66 (
  67     const fvPatch& p,
  68     const DimensionedField<vector, volMesh>& iF,
  69     const dictionary& dict
  70 )
  71 :
  72     directionMixedFvPatchVectorField(p, iF),
  73     fieldName_(dimensionedInternalField().name()),
  74     relaxationFactor_(readScalar(dict.lookup("relaxationFactor"))),
  75     traction_(p.size(), vector::zero),
  76     initiationTraction_(p.size(), vector::zero),
  77     separationDistance_(p.size(), vector::zero),
  78     oldSeparationDistance_(p.size(), vector::zero),
  79     unloadingSeparationDistance_(p.size(), vector::zero),
  80     minUnloadingSeparationDistance_
  81     (
  82         readScalar(dict.lookup("minUnloadingSeparationDistance"))
  83     ),
  84     beta_(readScalar(dict.lookup("beta"))),
  85     explicitSeparationDistance_(dict.lookup("explicitSeparationDistance")),
  86     contact_(dict.lookup("contact")),
  87     contactConstant_(readScalar(dict.lookup("contactConstant"))),
  88     curTimeIndex_(-1)
  89 {
  90   Info << "Creating solidCohesiveFixedModeMix boundary condition for "
  91        << patch().name() << " patch" << nl
  92        << "NOTE: this method only uses sigmaMax and GIc from "
  93        << "cohesiveProperties; tauMax and GIIc are ignored!"
  94        << endl;
  95
  96     if (minUnloadingSeparationDistance_ < 0.001)
  97     {
  98         minUnloadingSeparationDistance_ = 0.001;
  99     }
 100
 101     if (minUnloadingSeparationDistance_ > 1.0)
 102     {
 103         minUnloadingSeparationDistance_ = 1.0;
 104     }
 105
 106     Info << "minUnloadingSeparationDistance: "
 107         << minUnloadingSeparationDistance_ << endl;
 108
 109 //     if (beta_ < (1 - SMALL))
 110 //     {
 111 //         beta_ = 0;
 112 //     }
 113 //     else
 114 //     {
 115 //         beta_ = 1;
 116 //     }
 117
 118     Info << "beta: " << beta_ << endl;
 119
 120     if (dict.found("refValue"))
 121     {
 122         this->refValue() = vectorField("refValue", dict, p.size());
 123     }
 124     else
 125     {
 126         this->refValue() = vector::zero;
 127     }
 128
 129     if (dict.found("refGradient"))
 130     {
 131         this->refGrad() = vectorField("refGradient", dict, p.size());
 132     }
 133     else
 134     {
 135         this->refGrad() = vector::zero;
 136     }
 137
 138     if (dict.found("valueFraction"))
 139     {
 140         this->valueFraction() =
 141             symmTensorField("valueFraction", dict, p.size());
 142     }
 143     else
 144     {
 145         this->valueFraction() = symmTensor::zero;
 146     }
 147
 148     if (dict.found("value"))
 149     {
 150         vectorField::operator=(vectorField("value", dict, p.size()));
 151     }
 152     else
 153     {
 154         vectorField normalValue = transform(valueFraction(), refValue());
 155
 156         vectorField gradValue =
 157             this->patchInternalField() + refGrad()/this->patch().deltaCoeffs();
 158
 159         vectorField transformGradValue =
 160             transform(I - valueFraction(), gradValue);
 161
 162         vectorField::operator=(normalValue + transformGradValue);
 163     }
 164
 165     if (dict.found("traction"))
 166     {
 167         traction_ =
 168             vectorField("traction", dict, p.size());
 169     }
 170
 171     if (dict.found("initiationTraction"))
 172     {
 173         initiationTraction_ =
 174             vectorField("initiationTraction", dict, p.size());
 175     }
 176
 177     if (dict.found("separationDistance"))
 178     {
 179         separationDistance_ =
 180             vectorField("separationDistance", dict, p.size());
 181     }
 182
 183     if (dict.found("oldSeparationDistance"))
 184     {
 185         separationDistance_ =
 186             vectorField("oldSeparationDistance", dict, p.size());
 187     }
 188
 189     if (dict.found("unloadingSeparationDistance"))
 190     {
 191         unloadingSeparationDistance_ =
 192             vectorField("unloadingSeparationDistance", dict, p.size());
 193     }
 194 }
 195
 196
 197 solidCohesiveFixedModeMixFvPatchVectorField::
 198 solidCohesiveFixedModeMixFvPatchVectorField
 199 (
 200     const solidCohesiveFixedModeMixFvPatchVectorField& cpf
 201 )
 202 :
 203     directionMixedFvPatchVectorField(cpf),
 204     fieldName_(cpf.fieldName_),
 205     relaxationFactor_(cpf.relaxationFactor_),
 206     traction_(cpf.traction_),
 207     initiationTraction_(cpf.initiationTraction_),
 208     separationDistance_(cpf.separationDistance_),
 209     oldSeparationDistance_(cpf.oldSeparationDistance_),
 210     unloadingSeparationDistance_(cpf.unloadingSeparationDistance_),
 211     minUnloadingSeparationDistance_(cpf.minUnloadingSeparationDistance_),
 212     beta_(cpf.beta_),
 213     explicitSeparationDistance_(cpf.explicitSeparationDistance_),
 214     contact_(cpf.contact_),
 215     contactConstant_(cpf.contactConstant_),
 216     curTimeIndex_(-1)
 217 {}
 218
 219
 220 solidCohesiveFixedModeMixFvPatchVectorField::
 221 solidCohesiveFixedModeMixFvPatchVectorField
 222 (
 223     const solidCohesiveFixedModeMixFvPatchVectorField& cpf,
 224     const fvPatch& p,
 225     const DimensionedField<vector, volMesh>& iF,
 226     const fvPatchFieldMapper& mapper
 227 )
 228 :
 229     directionMixedFvPatchVectorField(cpf, p, iF, mapper),
 230     fieldName_(cpf.fieldName_),
 231     relaxationFactor_(cpf.relaxationFactor_),
 232     traction_(cpf.traction_, mapper),
 233     initiationTraction_(cpf.initiationTraction_, mapper),
 234     separationDistance_(cpf.separationDistance_, mapper),
 235     oldSeparationDistance_(cpf.oldSeparationDistance_, mapper),
 236     unloadingSeparationDistance_(cpf.unloadingSeparationDistance_, mapper),
 237     minUnloadingSeparationDistance_(cpf.minUnloadingSeparationDistance_),
 238     beta_(cpf.beta_),
 239     explicitSeparationDistance_(cpf.explicitSeparationDistance_),
 240     contact_(cpf.contact_),
 241     contactConstant_(cpf.contactConstant_),
 242     curTimeIndex_(-1)
 243 {}
 244
 245
 246 solidCohesiveFixedModeMixFvPatchVectorField::
 247 solidCohesiveFixedModeMixFvPatchVectorField
 248 (
 249     const solidCohesiveFixedModeMixFvPatchVectorField& cpf,
 250     const DimensionedField<vector, volMesh>& iF
 251 )
 252 :
 253     directionMixedFvPatchVectorField(cpf, iF),
 254     fieldName_(cpf.fieldName_),
 255     relaxationFactor_(cpf.relaxationFactor_),
 256     traction_(cpf.traction_),
 257     initiationTraction_(cpf.initiationTraction_),
 258     separationDistance_(cpf.separationDistance_),
 259     oldSeparationDistance_(cpf.oldSeparationDistance_),
 260     unloadingSeparationDistance_(cpf.unloadingSeparationDistance_),
 261     minUnloadingSeparationDistance_(cpf.minUnloadingSeparationDistance_),
 262     beta_(cpf.beta_),
 263     explicitSeparationDistance_(cpf.explicitSeparationDistance_),
 264     contact_(cpf.contact_),
 265     contactConstant_(cpf.contactConstant_),
 266     curTimeIndex_(-1)
 267 {}
 268
 269
 270 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 271
 272
 273 tmp<scalarField>
 274 solidCohesiveFixedModeMixFvPatchVectorField::relativeSeparationDistance() const
 275 {
 276     tmp<scalarField> tRelativeSeparationDistance
 277     (
 278         new scalarField(size(), 0.0)
 279     );
 280
 281     vectorField n = patch().nf();
 282
 283     const constitutiveModel& rheology =
 284         this->db().objectRegistry::lookupObject<constitutiveModel>
 285         ("rheologyProperties");
 286     label patchID = patch().index();
 287     const scalarField sigmaMax =
 288         rheology.cohLaw().sigmaMax()().boundaryField()[patchID];
 289     const scalarField GIc = rheology.cohLaw().GIc()().boundaryField()[patchID];
 290     const scalarField deltaC = GIc/sigmaMax;
 291     int numCrackedFaces = 0;
 292
 293     for (label i=0; i<size(); i++)
 294     {
 295         vector curSepDist =
 296             separationDistance_[i];
 297
 298         if (mag(n[i]&curSepDist) < 0)
 299         {
 300             curSepDist -= n[i]*(n[i]&curSepDist);
 301         }
 302
 303         tRelativeSeparationDistance()[i] =
 304             // mag(curSepDist)/law().deltaC().value();
 305             mag(curSepDist)/deltaC[i];
 306
 307     if (tRelativeSeparationDistance()[i] > 1.0) numCrackedFaces++;
 308     }
 309
 310     Info<< "Relative separation distance, max: "
 311         << gMax(tRelativeSeparationDistance())
 312      << ", avr: " << gAverage(tRelativeSeparationDistance())
 313      << ", min: " << gMin(tRelativeSeparationDistance())
 314      << ", numCrackedFaces: " << numCrackedFaces << endl;
 315
 316     return tRelativeSeparationDistance;
 317 }
 318
 319
 320 void solidCohesiveFixedModeMixFvPatchVectorField::autoMap
 321 (
 322     const fvPatchFieldMapper& m
 323 )
 324 {
 325     // if (cohesiveLawPtr_ == NULL)
 326     // {
 327     //     FatalErrorIn("cohesiveFvPatchVectorField::autoMap")
 328     //         << "NULL cohesive law"
 329     //             << abort(FatalError);
 330     // }
 331
 332     directionMixedFvPatchVectorField::autoMap(m);
 333
 334     traction_.autoMap(m);
 335     initiationTraction_.autoMap(m);
 336     separationDistance_.autoMap(m);
 337     oldSeparationDistance_.autoMap(m);
 338     unloadingSeparationDistance_.autoMap(m);
 339
 340     // Must use primitive mesh data because fvMesh data is packed in fields
 341     // and cannot be accessed during mapping.  HJ, 12/Dec/2008
 342     vectorField n = patch().patch().faceNormals();
 343
 344     const labelList& addressing = m.directAddressing();
 345
 346     label nNewFaces = m.size() - m.sizeBeforeMapping();
 347
 348     if ( (patch().size()==1) && (nNewFaces == 1) )
 349     {
 350         label i=0;
 351         this->valueFraction()[i] = symmTensor::zero;
 352         traction_[i] = vector::zero; // set in solver
 353         initiationTraction_[i] = vector::zero; // set in solver
 354         separationDistance_[i] = vector::zero;
 355         oldSeparationDistance_[i] = vector::zero;
 356         unloadingSeparationDistance_[i] = vector::zero;
 357     }
 358     else if ( (patch().size()==2) && (nNewFaces == 1) )
 359     {
 360         label i=1;
 361         this->valueFraction()[i] = symmTensor::zero;
 362         traction_[i] = vector::zero; //law().sigmaMax().value()*n[i];
 363         initiationTraction_[i] = vector::zero; //law().sigmaMax().value()*n[i];
 364         separationDistance_[i] = vector::zero;
 365         oldSeparationDistance_[i] = vector::zero;
 366         unloadingSeparationDistance_[i] = vector::zero;
 367     }
 368     else if ( (patch().size()==2) && (nNewFaces == 2) )
 369     {
 370         label i=0;
 371         this->valueFraction()[i] = symmTensor::zero;
 372         traction_[i] = vector::zero; //law().sigmaMax().value()*n[i];
 373         initiationTraction_[i] = vector::zero; //law().sigmaMax().value()*n[i];
 374         separationDistance_[i] = vector::zero;
 375         oldSeparationDistance_[i] = vector::zero;
 376         unloadingSeparationDistance_[i] = vector::zero;
 377         i=1;
 378         this->valueFraction()[i] = symmTensor::zero;
 379         traction_[i] = vector::zero; //law().sigmaMax().value()*n[i];
 380         initiationTraction_[i] = vector::zero; //law().sigmaMax().value()*n[i];
 381         separationDistance_[i] = vector::zero;
 382         oldSeparationDistance_[i] = vector::zero;
 383         unloadingSeparationDistance_[i] = vector::zero;
 384     }
 385     else
 386     {
 387         for (label i = 1; i < patch().size(); i++)
 388         {
 389             if (addressing[i] == 0)
 390             {
 391           //Info << "correcting automap for face " << i << " to zero" << endl;
 392                 this->valueFraction()[i] = symmTensor::zero;
 393                 traction_[i] = vector::zero; //law().sigmaMax().value()*n[i];
 394                 initiationTraction_[i] = vector::zero;
 395                 separationDistance_[i] = vector::zero;
 396                 oldSeparationDistance_[i] = vector::zero;
 397                 unloadingSeparationDistance_[i] = vector::zero;
 398             }
 399         }
 400     }
 401
 402 //     label sizeByTwo = patch().size()/2;
 403 //     for (label i = 0; i < sizeByTwo; i++)
 404 //     {
 405 //         if (addressing[i] == addressing[sizeByTwo + i])
 406 //         {
 407 //             this->valueFraction()[i] = symmTensor::zero;
 408 //             this->valueFraction()[sizeByTwo + i] = symmTensor::zero;
 409
 410 //             traction_[i] = law().sigmaMax().value()*n[i];
 411 //             traction_[sizeByTwo + i] =
 412 //                 law().sigmaMax().value()*n[sizeByTwo + i];
 413
 414 //             initiationTraction_[i] = law().sigmaMax().value()*n[i];
 415 //             initiationTraction_[sizeByTwo + i] =
 416 //                 law().sigmaMax().value()*n[sizeByTwo + i];
 417
 418 //             separationDistance_[i] = vector::zero;
 419 //             separationDistance_[sizeByTwo + i] = vector::zero;
 420
 421 //             oldSeparationDistance_[i] = vector::zero;
 422 //             oldSeparationDistance_[sizeByTwo + i] = vector::zero;
 423
 424 //             unloadingSeparationDistance_[i] = vector::zero;
 425 //             unloadingSeparationDistance_[sizeByTwo + i] = vector::zero;
 426 //         }
 427 //     }
 428 }
 429
 430
 431 // Reverse-map the given fvPatchField onto this fvPatchField
 432 void solidCohesiveFixedModeMixFvPatchVectorField::rmap
 433 (
 434     const fvPatchVectorField& ptf,
 435     const labelList& addr
 436 )
 437 {
 438     directionMixedFvPatchVectorField::rmap(ptf, addr);
 439
 440     const solidCohesiveFixedModeMixFvPatchVectorField& dmptf =
 441         refCast<const solidCohesiveFixedModeMixFvPatchVectorField>(ptf);
 442
 443     // No need to grab the cohesive zone pointer more than once
 444     // if (!cohesiveLawPtr_)
 445     // {
 446     //     cohesiveLawPtr_ = dmptf.cohesiveLawPtr_->clone().ptr();
 447     // }
 448
 449     relaxationFactor_ = dmptf.relaxationFactor_;
 450     traction_.rmap(dmptf.traction_, addr);
 451     initiationTraction_.rmap(dmptf.initiationTraction_, addr);
 452     separationDistance_.rmap(dmptf.separationDistance_, addr);
 453     oldSeparationDistance_.rmap(dmptf.oldSeparationDistance_, addr);
 454     unloadingSeparationDistance_.rmap
 455     (
 456         dmptf.unloadingSeparationDistance_,
 457         addr
 458     );
 459     minUnloadingSeparationDistance_ = dmptf.minUnloadingSeparationDistance_;
 460     beta_ = dmptf.beta_;
 461     contact_ = dmptf.contact_;
 462     contactConstant_ = dmptf.contactConstant_;
 463 }
 464
 465
 466 // Update the coefficients associated with the patch field
 467 void solidCohesiveFixedModeMixFvPatchVectorField::updateCoeffs()
 468 {
 469     if (updated())
 470     {
 471         return;
 472     }
 473
 474     vectorField n = patch().nf();
 475
 476     // lookup cohesive law from rheology
 477     // Note: this method only uses sigmaMax and GIc
 478     // tauMax and GIIc are ignored
 479     const constitutiveModel& rheology =
 480         this->db().objectRegistry::lookupObject<constitutiveModel>
 481         ("rheologyProperties");
 482     label patchID = patch().index();
 483     const scalarField sigmaMax =
 484         rheology.cohLaw().sigmaMax()().boundaryField()[patchID];
 485     //const scalarField tauMax =
 486     //    rheology.cohLaw().tauMax()().boundaryField()[patchID];
 487     const scalarField GIc =
 488         rheology.cohLaw().GIc()().boundaryField()[patchID];
 489     //const scalarField GIIc =
 490     //    rheology.cohLaw().GIIc()().boundaryField()[patchID];
 491     const scalarField deltaC = GIc/sigmaMax;
 492
 493     if (curTimeIndex_ != this->db().time().timeIndex())
 494     {
 495         forAll(unloadingSeparationDistance_, faceI)
 496         {
 497             vector curSepDist =
 498                 separationDistance_[faceI];
 499
 500             if (explicitSeparationDistance_)
 501             {
 502                 curSepDist = oldSeparationDistance_[faceI];
 503             }
 504
 505             scalar curNormalSepDist = (n[faceI]&curSepDist);
 506
 507             if (curNormalSepDist < 0)
 508             {
 509                 curSepDist -= n[faceI]*(n[faceI]&curSepDist);
 510                 curNormalSepDist = 0;
 511             }
 512
 513             vector curTangentialSepDist =
 514                 ((I - n[faceI]*n[faceI])&curSepDist);
 515
 516             if (beta_ < SMALL) // Mode I
 517             {
 518                 if
 519                 (
 520                     curNormalSepDist
 521                     > minUnloadingSeparationDistance_*deltaC[faceI]
 522                 )
 523                 {
 524                     if
 525                     (
 526                         curNormalSepDist
 527                       > mag(unloadingSeparationDistance_[faceI])
 528                     )
 529                     {
 530                         unloadingSeparationDistance_[faceI] =
 531                             curNormalSepDist*n[faceI];
 532                     }
 533                 }
 534             }
 535             else // Mixed mode I/II
 536             {
 537                 scalar curEffSepDist =
 538                     sqrt
 539                     (
 540                         sqr(curNormalSepDist)
 541                       + sqr(beta_)*magSqr(curTangentialSepDist)
 542                     );
 543
 544                 if
 545                 (
 546                     curEffSepDist
 547                     > minUnloadingSeparationDistance_*deltaC[faceI]
 548                 )
 549                 {
 550                     if
 551                     (
 552                         curEffSepDist
 553                       > mag(unloadingSeparationDistance_[faceI])
 554                     )
 555                     {
 556                         unloadingSeparationDistance_[faceI] =
 557                             curNormalSepDist*n[faceI]
 558                           + beta_*curTangentialSepDist;
 559                     }
 560                 }
 561             }
 562         }
 563
 564         oldSeparationDistance_ = separationDistance_;
 565         curTimeIndex_ = this->db().time().timeIndex();
 566     }
 567
 568     // Get face cells regions
 569     const unallocLabelList& faceCells = patch().faceCells();
 570
 571     const fvMesh& mesh = patch().boundaryMesh().mesh();
 572
 573     if (!isA<crackerFvMesh>(mesh))
 574     {
 575         FatalErrorIn
 576         (
 577             "void solidCohesiveFixedModeMixFvPatchVectorField::"
 578             "updateCoeffs() const"
 579         )   << "Mesh should be of type: " << crackerFvMesh::typeName
 580             << abort(FatalError);
 581     }
 582
 583     const crackerFvMesh& crackerMesh = refCast<const crackerFvMesh>(mesh);
 584
 585     const regionSplit& regions = crackerMesh.regions();
 586
 587     labelField faceCellRegion(size(), -1);
 588     forAll(faceCellRegion, faceI)
 589     {
 590         faceCellRegion[faceI] = regions[faceCells[faceI]];
 591     }
 592     labelField globalFaceCellRegion =
 593         crackerMesh.globalCrackField(faceCellRegion);
 594
 595     // Looking up rheology
 596     const fvPatchField<scalar>& mu =
 597       patch().lookupPatchField<volScalarField, scalar>("mu");
 598     const fvPatchField<scalar>& lambda =
 599       patch().lookupPatchField<volScalarField, scalar>("lambda");
 600
 601     const fvPatchField<tensor>& gradU =
 602         patch().lookupPatchField<volTensorField, tensor>("grad(U)");
 603
 604     // Patch displacement
 605     const vectorField& UPatch = *this;
 606
 607     // Global displacement
 608     vectorField globalUPatch = crackerMesh.globalCrackField(UPatch);
 609
 610     // Update separation distance
 611     vectorField newSeparationDistance(patch().size(), vector::zero);
 612     const labelList& gcfa = crackerMesh.globalCrackFaceAddressing();
 613     label globalIndex = crackerMesh.localCrackStart();
 614     for (label i = 0; i < patch().size(); i++)
 615     {
 616         newSeparationDistance[i] =
 617             globalUPatch[gcfa[globalIndex]]
 618           - globalUPatch[globalIndex];
 619
 620         globalIndex++;
 621     }
 622
 623 //     label sizeByTwo = patch().size()/2;
 624 //     for (label i = 0; i < sizeByTwo; i++)
 625 //     {
 626 //         newSeparationDistance[i] = UPatch[sizeByTwo + i] - UPatch[i];
 627 //         newSeparationDistance[sizeByTwo + i] = -newSeparationDistance[i];
 628 //     }
 629
 630 //     newSeparationDistance =
 631 //         beta_*((I-n*n)&newSeparationDistance)
 632 //       + (n&newSeparationDistance)*n;
 633
 634     if (explicitSeparationDistance_)
 635     {
 636         separationDistance_ = newSeparationDistance;
 637     }
 638     else
 639     {
 640         separationDistance_ =
 641             separationDistance_
 642           + relaxationFactor_
 643            *(newSeparationDistance - separationDistance_);
 644     }
 645
 646 //     for (label i=0; i<sizeByTwo; i++)
 647     globalIndex = crackerMesh.localCrackStart();
 648     for (label i = 0; i < patch().size(); i++)
 649     {
 650         vector curSepDist = separationDistance_[i];
 651
 652         if (explicitSeparationDistance_)
 653         {
 654             curSepDist = oldSeparationDistance_[i];
 655         }
 656
 657         scalar curNormalSepDist = (n[i]&curSepDist);
 658         scalar curNegNormalSepDist = 0;
 659
 660         if (curNormalSepDist < 0)
 661         {
 662             curSepDist -= n[i]*(n[i]&curSepDist);
 663             curNegNormalSepDist = curNormalSepDist;
 664             curNormalSepDist = 0;
 665         }
 666
 667         vector curTangentialSepDist =
 668             ((I - n[i]*n[i])&curSepDist);
 669
 670         scalar curNormalTraction = 0;
 671         vector curTangentialTraction = vector::zero;
 672
 673 //         if (faceCellRegion[i] == faceCellRegion[sizeByTwo+i])
 674         if
 675         (
 676             globalFaceCellRegion[globalIndex]
 677          == globalFaceCellRegion[gcfa[globalIndex]]
 678         )
 679         {
 680             if (beta_ > SMALL)
 681             {
 682                 // Mode II or Mixed mode I/II
 683
 684                 scalar curEffSepDist =
 685                     sqrt
 686                     (
 687                         sqr(curNormalSepDist)
 688                       + sqr(beta_)*mag(curTangentialSepDist)
 689                     );
 690
 691                 if
 692                 (
 693                     curEffSepDist
 694                   > (mag(unloadingSeparationDistance_[i]) - SMALL)
 695                 )
 696                 {
 697                     // Loading
 698
 699             if (curEffSepDist > deltaC[i])
 700               {
 701             curNormalTraction = 0.0;
 702             curTangentialTraction = vector::zero;
 703               }
 704             else
 705               {
 706             curNormalTraction = (n[i]&initiationTraction_[i]);
 707             curTangentialTraction =
 708               ((I - n[i]*n[i])&initiationTraction_[i]);
 709               }
 710
 711             // scale traction based on cohesive law
 712             // Not needed for Dugdale as traction stays at initiation value
 713                     // if (curNormalTraction >= 0)
 714                     // {
 715                     //     // Tension
 716
 717                     //     curNormalTraction *=
 718             //        /sigmaMax[i];
 719                     //        //  law().traction(curEffSepDist)
 720                     //        // /law().sigmaMax().value();
 721
 722                     //     curTangentialTraction *=
 723             //        /sigmaMax[i];
 724                     //        //  law().traction(curEffSepDist)
 725                     //        // /law().sigmaMax().value();
 726                     // }
 727                     // else
 728                     // {
 729                     //     // Compression
 730
 731                     //     curTangentialTraction *=
 732                     //         law().traction(curEffSepDist)
 733                     //        /law().sigmaMax().value();
 734                     // }
 735                 }
 736                 else
 737                 {
 738                     // Unloading
 739
 740                     scalar unloadingNormalTraction =
 741                         (n[i]&initiationTraction_[i]);
 742                     vector unloadingTangentialTraction =
 743                         ((I - n[i]*n[i])&initiationTraction_[i]);
 744
 745                     if (unloadingNormalTraction >= 0)
 746                     {
 747               // scaling not needed for Dugdale
 748                         // unloadingNormalTraction *=
 749                         //     law().traction
 750                         //     (
 751                         //         mag(unloadingSeparationDistance_[i])
 752                         //     )
 753                         //    /law().sigmaMax().value();
 754
 755                         // unloadingTangentialTraction *=
 756                         //     law().traction
 757                         //     (
 758                         //         mag(unloadingSeparationDistance_[i])
 759                         //     )
 760                         //    /law().sigmaMax().value();
 761
 762                         curNormalTraction =
 763                             unloadingNormalTraction
 764                            *(
 765                                 curEffSepDist
 766                                /mag(unloadingSeparationDistance_[i])
 767                             );
 768
 769                         curTangentialTraction =
 770                             unloadingTangentialTraction
 771                            *(
 772                                 curEffSepDist
 773                                /mag(unloadingSeparationDistance_[i])
 774                             );
 775                     }
 776                     else
 777                     {
 778               // scaling not needed for Dugdale
 779                         // unloadingTangentialTraction *=
 780                         //     law().traction
 781                         //     (
 782                         //         mag(unloadingSeparationDistance_[i])
 783                         //     )
 784                         //    /law().sigmaMax().value();
 785
 786                         curTangentialTraction =
 787                             unloadingTangentialTraction
 788                            *(
 789                                 curEffSepDist
 790                                /mag(unloadingSeparationDistance_[i])
 791                             );
 792
 793                         curNormalTraction = unloadingNormalTraction;
 794                     }
 795                 }
 796             }
 797             else
 798             {
 799                 // Mode I
 800
 801                 if
 802                 (
 803                     curNormalSepDist
 804                   > (mag(unloadingSeparationDistance_[i]) - SMALL)
 805                 )
 806                 {
 807                     //Loading
 808
 809             if (curNormalSepDist > deltaC[i])
 810               {
 811             curNormalTraction = 0.0;
 812             curTangentialTraction = vector::zero;
 813               }
 814             else
 815               {
 816             curNormalTraction =
 817               (n[i]&initiationTraction_[i]);
 818
 819             curTangentialTraction =
 820               ((I - n[i]*n[i])&initiationTraction_[i]);
 821               }
 822
 823               // scaling not needed for Dugdale
 824                     // curNormalTraction *=
 825                     //     law().traction(curNormalSepDist)
 826                     //    /law().sigmaMax().value();
 827
 828               // scaling not needed for Dugdale
 829                     // curTangentialTraction *=
 830                     //     law().traction(curNormalSepDist)
 831                     //    /law().sigmaMax().value();
 832                 }
 833                 else
 834                 {
 835                     // Unloading
 836
 837                     scalar unloadingNormalTraction =
 838               (n[i]&initiationTraction_[i]);
 839               // scaling not needed for Dugdale
 840               //*law().traction(mag(unloadingSeparationDistance_[i]))
 841               ///law().sigmaMax().value();
 842
 843                     vector unloadingTangentialTraction =
 844               ((I - n[i]*n[i])&initiationTraction_[i]);
 845               // scaling not needed for Dugdale
 846                        // *law().traction(mag(unloadingSeparationDistance_[i]))
 847                        // /law().sigmaMax().value();
 848
 849                     curNormalTraction =
 850                         unloadingNormalTraction
 851                        *(
 852                            curNormalSepDist
 853                           /mag(unloadingSeparationDistance_[i])
 854                         );
 855
 856 //                     if (curNegNormalSepDist < 0)
 857 //                     {
 858 //                         curNormalTraction =
 859 //                             unloadingNormalTraction
 860 //                            *(
 861 //                                 curNegNormalSepDist
 862 //                                /mag(unloadingSeparationDistance_[i])
 863 //                             );
 864 //                     }
 865
 866                     curTangentialTraction =
 867                         unloadingTangentialTraction
 868                        *(
 869                            curNormalSepDist
 870                           /mag(unloadingSeparationDistance_[i])
 871                         );
 872                 }
 873             }
 874         }
 875
 876         // Correct direction of tangential traction
 877         // to be aligned with tangential separation distance
 878         if (mag(curTangentialSepDist) > SMALL)
 879         {
 880             curTangentialTraction =
 881                 mag(curTangentialTraction)
 882                *curTangentialSepDist
 883                /mag(curTangentialSepDist);
 884         }
 885
 886         // Contact
 887         if (contact_ && (curNegNormalSepDist < 0))
 888         {
 889       scalar m =
 890         (1.0/(1.0-contactConstant_))
 891         *sigmaMax[i]/deltaC[i];
 892         //      *law().sigmaMax().value()
 893       ///law().deltaC().value();
 894
 895             curNormalTraction = m*curNegNormalSepDist;
 896
 897             // Info << "Contact: " << curNegNormalSepDist/law().deltaC().value()
 898             //     << ", " << curNormalTraction << endl;
 899             //Info << "Contact: " << curNegNormalSepDist/deltaC[i]
 900         //  << ", " << curNormalTraction << endl;
 901         }
 902
 903         traction_[i] = curNormalTraction*n[i] + curTangentialTraction;
 904 //         traction_[sizeByTwo + i] = -traction_[i];
 905     }
 906
 907     this->refGrad() =
 908     (
 909         traction_
 910       - (n & (mu*gradU.T() - (mu + lambda)*gradU))
 911       - n*lambda*tr(gradU)
 912     )/(2.0*mu + lambda);
 913
 914     directionMixedFvPatchVectorField::updateCoeffs();
 915 }
 916
 917
 918 // Write
 919 void solidCohesiveFixedModeMixFvPatchVectorField::write(Ostream& os) const
 920 {
 921     directionMixedFvPatchVectorField::write(os);
 922     traction_.writeEntry("traction", os);
 923     initiationTraction_.writeEntry("initiationTraction", os);
 924     separationDistance_.writeEntry("separationDistance", os);
 925     oldSeparationDistance_.writeEntry("oldSeparationDistance", os);
 926     unloadingSeparationDistance_.writeEntry("unloadingSeparationDistance", os);
 927     //os.writeKeyword("cohesiveLaw") << law().type()
 928     //  << token::END_STATEMENT << nl;
 929     os.writeKeyword("relaxationFactor") << relaxationFactor_
 930         << token::END_STATEMENT << nl;
 931     os.writeKeyword("beta") << beta_
 932         << token::END_STATEMENT << nl;
 933     os.writeKeyword("explicitSeparationDistance")
 934         << explicitSeparationDistance_
 935             << token::END_STATEMENT << nl;
 936     os.writeKeyword("contact") << contact_
 937         << token::END_STATEMENT << nl;
 938     os.writeKeyword("contactConstant") << contactConstant_
 939         << token::END_STATEMENT << nl;
 940     os.writeKeyword("minUnloadingSeparationDistance")
 941         << minUnloadingSeparationDistance_
 942             << token::END_STATEMENT << nl;
 943 }
 944
 945
 946 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 947
 948 makePatchTypeField
 949 (
 950     fvPatchVectorField,
 951     solidCohesiveFixedModeMixFvPatchVectorField
 952 );
 953
 954 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 955
 956 } // End namespace Foam
 957
 958 // ************************************************************************* //