src/solidModels/constitutiveModel/cohesiveLaws/multiMaterialCohesiveLaw/multiMaterialCohesiveLaw.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 "multiMaterialCohesiveLaw.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "zeroGradientFvPatchFields.H"
  29 #include "fvc.H"
  30 #include "crackerFvMesh.H"
  31 #include "multiMaterial.H"
  32 #include "constitutiveModel.H"
  33 #include "cohesiveFvPatch.H"
  34 #include "cohesiveLaw.H"
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40     defineTypeNameAndDebug(multiMaterialCohesiveLaw, 0);
  41     addToRunTimeSelectionTable
  42     (
  43         cohesiveLaw,
  44         multiMaterialCohesiveLaw,
  45         dictionary
  46     );
  47 }
  48
  49
  50 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  51
  52 Foam::tmp<Foam::scalarField> Foam::multiMaterialCohesiveLaw::indicator
  53 (
  54     const label i
  55 ) const
  56 {
  57     const scalarField& mat = materials_.internalField();
  58
  59     tmp<scalarField> tresult(new scalarField(mat.size(), 0.0));
  60     scalarField& result = tresult();
  61
  62     forAll (mat, matI)
  63     {
  64         if (mat[matI] > i - SMALL && mat[matI] < i + 1 - SMALL)
  65         {
  66             result[matI] = 1.0;
  67         }
  68     }
  69
  70     return tresult;
  71 }
  72
  73
  74 Foam::scalar Foam::multiMaterialCohesiveLaw::indicator
  75 (
  76     const label index,
  77     const label cellID
  78 ) const
  79 {
  80     const scalar mat = materials_.internalField()[cellID];
  81     scalar result = 0.0;
  82
  83     if (mat > index - SMALL && mat < index + 1 - SMALL)
  84       {
  85     result = 1.0;
  86       }
  87
  88     return result;
  89 }
  90
  91 Foam::scalar Foam::multiMaterialCohesiveLaw::interfaceID
  92 (
  93     const label mat1,
  94     const label mat2
  95 ) const
  96 {
  97   word mat1name = (*this)[mat1].name();
  98   word mat2name = (*this)[mat2].name();
  99
 100   word interfaceName("interface_"+mat1name+"_"+mat2name);
 101   label interfaceLawID = -1;
 102   forAll(interfaceCohesiveLaws_, lawI)
 103   {
 104       if (interfaceCohesiveLaws_[lawI].name() == interfaceName)
 105       {
 106           interfaceLawID = lawI;
 107           break;
 108       }
 109   }
 110   if (interfaceLawID == -1)
 111   {
 112       // flip name
 113       interfaceName = word("interface_"+mat2name+"_"+mat1name);
 114       forAll(interfaceCohesiveLaws_, lawI)
 115       {
 116           if (interfaceCohesiveLaws_[lawI].name() == interfaceName)
 117           {
 118               interfaceLawID = lawI;
 119               break;
 120           }
 121       }
 122       if (interfaceLawID == -1)
 123       {
 124           FatalError
 125               << "Cannot find cohesive interfaceLaw "
 126               << word("interface_"+mat1name+"_"+mat2name) << " or "
 127               << interfaceName << nl
 128               << "One of these should be defined!"
 129               << exit(FatalError);
 130       }
 131   }
 132
 133   return interfaceLawID;
 134 }
 135 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 136
 137 // Construct from dictionary
 138 Foam::multiMaterialCohesiveLaw::multiMaterialCohesiveLaw
 139 (
 140     const word& name,
 141     const volSymmTensorField& sigma,
 142     const dictionary& dict
 143 )
 144 :
 145     cohesiveLaw(name, sigma, dict),
 146     PtrList<cohesiveLaw>(),
 147     materials_
 148     (
 149         sigma.mesh().objectRegistry::lookupObject<volScalarField>("materials")
 150     ),
 151     // (
 152     //  IOobject
 153     //  (
 154     //   "materials",
 155     //   mesh().time().timeName(),
 156     //   mesh(),
 157     //   IOobject::MUST_READ,
 158     //   IOobject::AUTO_WRITE
 159     //   ),
 160     //  mesh()
 161     //  ),
 162     interfaceCohesiveLaws_() //,
 163     // materialsSurf_
 164     // (
 165     //  IOobject
 166     //  (
 167     //   "materialsSurf",
 168     //   mesh().time().timeName(),
 169     //   mesh(),
 170     //   IOobject::NO_READ,
 171     //   IOobject::NO_WRITE
 172     //   ),
 173     //  mesh(),
 174     //  dimensionedScalar("zero", dimless, 0.0)
 175     //  )
 176 {
 177     PtrList<cohesiveLaw>& laws = *this;
 178
 179     PtrList<entry> lawEntries(dict.lookup("laws"));
 180     laws.setSize(lawEntries.size());
 181
 182     forAll (laws, lawI)
 183     {
 184         laws.set
 185         (
 186             lawI,
 187             cohesiveLaw::New
 188             (
 189                 lawEntries[lawI].keyword(),
 190                 sigma,
 191                 lawEntries[lawI].dict()
 192             )
 193         );
 194     }
 195
 196     if
 197     (
 198         min(materials_).value() < 0
 199      || max(materials_).value() > laws.size() + SMALL
 200     )
 201     {
 202         FatalErrorIn
 203         (
 204             "multiMaterialCohesiveLaw::multiMaterialCohesiveLaw\n"
 205             "(\n"
 206             "    const word& name,\n"
 207             "    const volSymmTensorField& sigma,\n"
 208             "    const dictionary& dict\n"
 209             ")"
 210         )   << "Invalid definition of material indicator field.  "
 211             << "Number of materials: " << laws.size()
 212             << " max index: " << max(materials_)
 213             << ".  Should be " << laws.size() - 1
 214             << abort(FatalError);
 215     }
 216
 217     // cohesive laws must be the same size as rheology laws
 218     const constitutiveModel& conModel =
 219         mesh().objectRegistry::lookupObject<constitutiveModel>
 220         ("rheologyProperties");
 221     if (conModel.law().type() == "multiMaterial")
 222       {
 223     const multiMaterial& mulMatLaw =
 224         refCast<const multiMaterial>(conModel.law());
 225
 226     if (laws.size() != mulMatLaw.size())
 227       {
 228         FatalError
 229             << "There should be the same number of cohesive laws "
 230             << "as rheology laws" << nl
 231             << "Currently there are " << mulMatLaw.size()
 232             << " rheology laws "
 233             << "and " << laws.size() << " cohesive laws!"
 234             << exit(FatalError);
 235       }
 236       }
 237
 238     // set interfaceCohesiveLaws_
 239     PtrList<entry> interfaceLawEntries(dict.lookup("interfaceLaws"));
 240     // if (interfaceLawEntries.size() != int(laws.size()*(laws.size()-1)/2))
 241     if ( mag(interfaceLawEntries.size() - (laws.size()*(laws.size()-1)/2))
 242         > SMALL )
 243       {
 244     // number of interfaces is a trianular number of number of materials
 245     // ((n)*(n-1)/2)
 246           FatalError
 247               << "There are " << interfaceLawEntries.size()
 248               << " interface cohesive"
 249               << " laws defined, but there should be "
 250               << (laws.size()*(laws.size()-1)/2)
 251               << "," << nl
 252               << "as there are " << laws.size()
 253               << " materials in cohesive laws" << nl
 254               << exit(FatalError);
 255       }
 256     interfaceCohesiveLaws_.setSize(interfaceLawEntries.size());
 257     forAll (interfaceCohesiveLaws_, lawI)
 258       {
 259         interfaceCohesiveLaws_.set
 260       (
 261        lawI,
 262        cohesiveLaw::New
 263        (
 264         interfaceLawEntries[lawI].keyword(),
 265         sigma,
 266         interfaceLawEntries[lawI].dict()
 267             )
 268        );
 269       }
 270
 271
 272     // Set materialsSurf
 273     // materialsSurf_ = fvc::interpolate(materials_);
 274     // forAll(mesh().boundary(), patchi)
 275     //   {
 276     //      materialsSurf_.boundaryField()[patchi] =
 277     //        materials_.boundaryField()[patchi].patchInternalField();
 278     //   }
 279     // // Fix interface values
 280     // const labelList owner = mesh().owner();
 281     // const labelList neighbour = mesh().neighbour();
 282     // forAll (materialsSurf_.internalField(), faceI)
 283     //   {
 284     //      // round value to integer and check difference
 285     //      // if it is small then the face is not on a multi-material
 286     //      // interface
 287     //      scalar matIDscalar = materialsSurf_.internalField()[faceI];
 288     //      label matID = int(matIDscalar);
 289     //      if (mag(matIDscalar - matID) > SMALL)
 290     //        {
 291     //          // find which interface it is on
 292     //          const label own = owner[faceI];
 293     //          const label nei = neighbour[faceI];
 294
 295     //          materialsSurf_.internalField()[faceI] =
 296     //            interfaceID(materials_[own], materials_[nei]) + laws.size();
 297     //        }
 298     //   }
 299
 300     // philipc
 301     // processor boundaries are meant to hold the patchNeighbourField
 302     // but the proc boundaries are interpolated by decomposePar
 303     // so we must correct them. Now the proc boundaries hold the
 304     // patchNiehgbourField
 305     //materials_.correctBoundaryConditions(); // done by rheologyLaw
 306 }
 307
 308
 309 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 310
 311 Foam::multiMaterialCohesiveLaw::~multiMaterialCohesiveLaw()
 312 {}
 313
 314
 315 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 316
 317 Foam::tmp<Foam::volScalarField>
 318 Foam::multiMaterialCohesiveLaw::materials() const
 319 {
 320     return materials_;
 321 }
 322
 323 Foam::tmp<Foam::surfaceScalarField>
 324 Foam::multiMaterialCohesiveLaw::sigmaMax() const
 325 {
 326     tmp<surfaceScalarField> tresult
 327     (
 328         new surfaceScalarField
 329         (
 330             IOobject
 331             (
 332             "sigmaMax",
 333             mesh().time().timeName(),
 334             mesh(),
 335             IOobject::NO_READ,
 336             IOobject::NO_WRITE
 337         ),
 338             mesh(),
 339         dimensionedScalar("zero", dimForce/dimArea, 0)
 340     )
 341     );
 342     surfaceScalarField& result = tresult();
 343
 344     // Accumulate data for all fields
 345     const PtrList<cohesiveLaw>& laws = *this;
 346     volScalarField indic
 347       (
 348        IOobject
 349        (
 350     "indic",
 351     mesh().time().timeName(),
 352     mesh(),
 353     IOobject::NO_READ,
 354     IOobject::NO_WRITE
 355     ),
 356        mesh(),
 357        dimensionedScalar("zero", dimless, 0),
 358        zeroGradientFvPatchScalarField::typeName
 359        );
 360     forAll (laws, lawI)
 361       {
 362     // interface fields will not be correct but we fix them after
 363     indic.internalField() = indicator(lawI)();
 364     surfaceScalarField indicatorSurf = fvc::interpolate(indic);
 365     // fix boundary fields
 366     forAll(mesh().boundary(), patchi)
 367       {
 368         indicatorSurf.boundaryField()[patchi] =
 369           indic.boundaryField()[patchi].patchInternalField();
 370       }
 371         result += indicatorSurf*laws[lawI].sigmaMax()();
 372
 373       }
 374
 375     // fix interfaces
 376     // forAll surfaces check if surface is a material interface
 377     // material indicator should read non integer
 378     // Get the two materials it is an interface of
 379     // Look up value of sigmaMaxf in dictionary
 380     // Overwrite existing value on surface with new value
 381
 382     const labelList& owner = mesh().owner();
 383     const labelList& neighbour = mesh().neighbour();
 384
 385     forAll(result.internalField(), faceI)
 386     {
 387       label ownerMat = label(materials_[owner[faceI]]);
 388       label neighbourMat = label(materials_[neighbour[faceI]]);
 389
 390         if (ownerMat != neighbourMat)
 391         {
 392       result.internalField()[faceI] = interfaceSigmaMax(ownerMat, neighbourMat);
 393     }
 394     }
 395
 396     forAll(mesh().boundary(), patchI)
 397     {
 398
 399         if (mesh().boundary()[patchI].type() == cohesiveFvPatch::typeName)
 400         {
 401       const crackerFvMesh& crackerMesh = refCast<const crackerFvMesh>(mesh());
 402       // we must use owner values
 403       scalarField localPatchMaterials =
 404           materials_.boundaryField()[patchI].patchInternalField();
 405       scalarField globalPatchMaterials =
 406           crackerMesh.globalCrackField(localPatchMaterials);
 407       // crackerMesh.globalCrackField(materials_.boundaryField()[patchI]);
 408       const labelList& gcfa = crackerMesh.globalCrackFaceAddressing();
 409       label globalIndex = crackerMesh.localCrackStart();
 410
 411       forAll(mesh().boundaryMesh()[patchI], facei)
 412         {
 413           label ownerMat = label(globalPatchMaterials[globalIndex]);
 414           label neighbourMat = label(globalPatchMaterials[gcfa[globalIndex]]);
 415
 416           if (ownerMat != neighbourMat)
 417                 {
 418                     result.boundaryField()[patchI][facei] =
 419                         interfaceSigmaMax(ownerMat, neighbourMat);
 420                 }
 421           globalIndex++;
 422             }
 423         }
 424
 425     // philipc
 426     else if (mesh().boundary()[patchI].coupled())
 427       {
 428         const scalarField ownerMatField =
 429             materials_.boundaryField()[patchI].patchInternalField();
 430         const scalarField neighbourMatField =
 431             materials_.boundaryField()[patchI].patchNeighbourField();
 432         //const labelList& faceCells = mesh().boundary()[patchI].faceCells();
 433
 434         forAll(mesh().boundaryMesh()[patchI], facei)
 435           {
 436         label ownerMat = label(ownerMatField[facei]);
 437         label neighbourMat = label(neighbourMatField[facei]);
 438
 439         if (ownerMat != neighbourMat)
 440           {
 441               // result.boundaryField()[patchI][facei] = sigmaMax();
 442               result.boundaryField()[patchI][facei] =
 443                   interfaceSigmaMax(ownerMat, neighbourMat);
 444           }
 445           }
 446       }
 447     }
 448
 449     return tresult;
 450 }
 451
 452
 453 Foam::tmp<Foam::surfaceScalarField>
 454 Foam::multiMaterialCohesiveLaw::tauMax() const
 455 {
 456     tmp<surfaceScalarField> tresult
 457     (
 458         new surfaceScalarField
 459         (
 460             IOobject
 461             (
 462             "tauMax",
 463             mesh().time().timeName(),
 464             mesh(),
 465             IOobject::NO_READ,
 466             IOobject::NO_WRITE
 467         ),
 468             mesh(),
 469         dimensionedScalar("zero", dimForce/dimArea, 0)
 470     )
 471     );
 472     surfaceScalarField& result = tresult();
 473
 474     const PtrList<cohesiveLaw>& laws = *this;
 475     volScalarField indic
 476       (
 477        IOobject
 478        (
 479     "indic",
 480     mesh().time().timeName(),
 481     mesh(),
 482     IOobject::NO_READ,
 483     IOobject::NO_WRITE
 484     ),
 485        mesh(),
 486        dimensionedScalar("zero", dimless, 0),
 487        zeroGradientFvPatchScalarField::typeName
 488        );
 489     forAll (laws, lawI)
 490       {
 491     // interface fields will not be correct but we fix them after
 492     indic.internalField() = indicator(lawI)();
 493     surfaceScalarField indicatorSurf = fvc::interpolate(indic);
 494     // fix boundary fields
 495     forAll(mesh().boundary(), patchi)
 496       {
 497         indicatorSurf.boundaryField()[patchi] =
 498           indic.boundaryField()[patchi].patchInternalField();
 499       }
 500         result += indicatorSurf*laws[lawI].tauMax()();
 501       }
 502
 503     // forAll surfaces check if surface is a material interface
 504     // material indicator should read non integer
 505     // Get the two materials it is an interface of
 506     // Look up value of tauMaxf in dictionary
 507     // Overwrite existing value on surface with new value
 508
 509     const labelList& owner = mesh().owner();
 510     const labelList& neighbour = mesh().neighbour();
 511
 512     forAll(result.internalField(), faceI)
 513     {
 514       label ownerMat = label(materials_[owner[faceI]]);
 515       label neighbourMat = label(materials_[neighbour[faceI]]);
 516
 517         if (ownerMat != neighbourMat)
 518         {
 519       result.internalField()[faceI] = interfaceTauMax(ownerMat, neighbourMat);
 520     }
 521     }
 522
 523     forAll(mesh().boundary(), patchI)
 524     {
 525
 526         if (mesh().boundary()[patchI].type() == cohesiveFvPatch::typeName)
 527         {
 528       //label size = (mesh().boundary()[patchI].size())/2;
 529       const crackerFvMesh& crackerMesh = refCast<const crackerFvMesh>(mesh());
 530       // we must use owner values
 531       scalarField localPatchMaterials =
 532           materials_.boundaryField()[patchI].patchInternalField();
 533       scalarField globalPatchMaterials =
 534           crackerMesh.globalCrackField(localPatchMaterials);
 535       // crackerMesh.globalCrackField(materials_.boundaryField()[patchI]);
 536       const labelList& gcfa = crackerMesh.globalCrackFaceAddressing();
 537       label globalIndex = crackerMesh.localCrackStart();
 538
 539       forAll(mesh().boundaryMesh()[patchI], facei)
 540         {
 541           label ownerMat = label(globalPatchMaterials[globalIndex]);
 542           label neighbourMat = label(globalPatchMaterials[gcfa[globalIndex]]);
 543
 544           if (ownerMat != neighbourMat)
 545         {
 546           result.boundaryField()[patchI][facei] =
 547               interfaceTauMax(ownerMat, neighbourMat);
 548                 }
 549           globalIndex++;
 550             }
 551         }
 552
 553     // philipc
 554     else if (mesh().boundary()[patchI].coupled())
 555       {
 556         const scalarField ownerMatField =
 557             materials_.boundaryField()[patchI].patchInternalField();
 558         const scalarField neighbourMatField =
 559             materials_.boundaryField()[patchI].patchNeighbourField();
 560         //const labelList& faceCells = mesh().boundary()[patchI].faceCells();
 561
 562         forAll(mesh().boundaryMesh()[patchI], facei)
 563           {
 564         label ownerMat = label(ownerMatField[facei]);
 565         label neighbourMat = label(neighbourMatField[facei]);
 566
 567         if (ownerMat != neighbourMat)
 568           {
 569             result.boundaryField()[patchI][facei] =
 570                 interfaceTauMax(ownerMat, neighbourMat);
 571           }
 572           }
 573       }
 574     }
 575
 576     return tresult;
 577 }
 578
 579 Foam::tmp<Foam::surfaceScalarField> Foam::multiMaterialCohesiveLaw::GIc() const
 580 {
 581     tmp<surfaceScalarField> tresult
 582     (
 583         new surfaceScalarField
 584         (
 585             IOobject
 586             (
 587                 "GIc",
 588                 mesh().time().timeName(),
 589                 mesh(),
 590                 IOobject::NO_READ,
 591                 IOobject::NO_WRITE
 592             ),
 593             mesh(),
 594         dimensionedScalar("zero", dimForce*dimLength/dimArea, 0)
 595     )
 596     );
 597     surfaceScalarField& result = tresult();
 598
 599     const PtrList<cohesiveLaw>& laws = *this;
 600     volScalarField indic
 601       (
 602        IOobject
 603        (
 604     "indic",
 605     mesh().time().timeName(),
 606     mesh(),
 607     IOobject::NO_READ,
 608     IOobject::NO_WRITE
 609     ),
 610        mesh(),
 611        dimensionedScalar("zero", dimless, 0),
 612        zeroGradientFvPatchScalarField::typeName
 613        );
 614     forAll (laws, lawI)
 615       {
 616     // interface fields will not be correct but we fix them after
 617     indic.internalField() = indicator(lawI)();
 618     surfaceScalarField indicatorSurf = fvc::interpolate(indic);
 619     // fix boundary fields
 620     forAll(mesh().boundary(), patchi)
 621       {
 622         indicatorSurf.boundaryField()[patchi] =
 623           indic.boundaryField()[patchi].patchInternalField();
 624       }
 625         result += indicatorSurf*laws[lawI].GIc()();
 626       }
 627
 628     const labelList& owner = mesh().owner();
 629     const labelList& neighbour = mesh().neighbour();
 630
 631     forAll(result.internalField(), faceI)
 632     {
 633       label ownerMat = label(materials_[owner[faceI]]);
 634       label neighbourMat = label(materials_[neighbour[faceI]]);
 635
 636         if (ownerMat != neighbourMat)
 637         {
 638         result.internalField()[faceI] = interfaceGIc(ownerMat, neighbourMat);
 639     }
 640     }
 641
 642     forAll(mesh().boundary(), patchI)
 643     {
 644
 645         if (mesh().boundary()[patchI].type() == cohesiveFvPatch::typeName)
 646         {
 647       //label size = (mesh().boundary()[patchI].size())/2;
 648       // const labelList& fCells = mesh().boundary()[patchI].faceCells();
 649       scalarField localPatchMaterials =
 650           materials_.boundaryField()[patchI].patchInternalField();
 651       const crackerFvMesh& crackerMesh = refCast<const crackerFvMesh>(mesh());
 652       scalarField globalPatchMaterials =
 653           crackerMesh.globalCrackField(localPatchMaterials);
 654       const labelList& gcfa = crackerMesh.globalCrackFaceAddressing();
 655       label globalIndex = crackerMesh.localCrackStart();
 656
 657       forAll(mesh().boundaryMesh()[patchI], facei)
 658         {
 659           label ownerMat = label(globalPatchMaterials[globalIndex]);
 660           label neighbourMat = label(globalPatchMaterials[gcfa[globalIndex]]);
 661
 662                 if (ownerMat != neighbourMat)
 663                 {
 664                     result.boundaryField()[patchI][facei] =
 665                         interfaceGIc(ownerMat, neighbourMat);
 666                 }
 667         globalIndex++;
 668             }
 669         }
 670     else if (mesh().boundary()[patchI].coupled())
 671       {
 672         const scalarField ownerMatField =
 673             materials_.boundaryField()[patchI].internalField();
 674         const scalarField neighbourMatField =
 675             materials_.boundaryField()[patchI].patchNeighbourField();
 676
 677         forAll(mesh().boundaryMesh()[patchI], facei)
 678         {
 679             label ownerMat = label(ownerMatField[facei]);
 680             label neighbourMat = label(neighbourMatField[facei]);
 681
 682             if (ownerMat != neighbourMat)
 683             {
 684                 //result.boundaryField()[patchI][facei] = iterGIc();
 685                 result.boundaryField()[patchI][facei] =
 686                     interfaceGIc(ownerMat, neighbourMat);
 687             }
 688         }
 689       }
 690     }
 691
 692     return tresult;
 693 }
 694
 695 Foam::tmp<Foam::surfaceScalarField> Foam::multiMaterialCohesiveLaw::GIIc() const
 696 {
 697     tmp<surfaceScalarField> tresult
 698     (
 699         new surfaceScalarField
 700         (
 701             IOobject
 702             (
 703                 "GIIc",
 704                 mesh().time().timeName(),
 705                 mesh(),
 706                 IOobject::NO_READ,
 707                 IOobject::NO_WRITE
 708             ),
 709             mesh(),
 710         dimensionedScalar("zero", dimForce*dimLength/dimArea, 0)
 711     )
 712     );
 713     surfaceScalarField& result = tresult();
 714
 715     const PtrList<cohesiveLaw>& laws = *this;
 716     volScalarField indic
 717       (
 718        IOobject
 719        (
 720     "indic",
 721     mesh().time().timeName(),
 722     mesh(),
 723     IOobject::NO_READ,
 724     IOobject::NO_WRITE
 725     ),
 726        mesh(),
 727        dimensionedScalar("zero", dimless, 0),
 728        zeroGradientFvPatchScalarField::typeName
 729        );
 730     forAll (laws, lawI)
 731       {
 732     // interface fields will not be correct but we fix them after
 733     indic.internalField() = indicator(lawI)();
 734     surfaceScalarField indicatorSurf = fvc::interpolate(indic);
 735     // fix boundary fields
 736     forAll(mesh().boundary(), patchi)
 737       {
 738         indicatorSurf.boundaryField()[patchi] =
 739           indic.boundaryField()[patchi].patchInternalField();
 740       }
 741         result += indicatorSurf*laws[lawI].GIIc()();
 742       }
 743
 744     const labelList& owner = mesh().owner();
 745     const labelList& neighbour = mesh().neighbour();
 746
 747     forAll(result.internalField(), faceI)
 748     {
 749       label ownerMat = label(materials_[owner[faceI]]);
 750       label neighbourMat = label(materials_[neighbour[faceI]]);
 751
 752         if (ownerMat != neighbourMat)
 753         {
 754       result.internalField()[faceI] = interfaceGIIc(ownerMat, neighbourMat);
 755     }
 756     }
 757
 758     forAll(mesh().boundary(), patchI)
 759     {
 760
 761         if (mesh().boundary()[patchI].type() == cohesiveFvPatch::typeName)
 762         {
 763       //label size = (mesh().boundary()[patchI].size())/2;
 764       // const labelList& fCells = mesh().boundary()[patchI].faceCells();
 765       scalarField localPatchMaterials =
 766           materials_.boundaryField()[patchI].patchInternalField();
 767       const crackerFvMesh& crackerMesh = refCast<const crackerFvMesh>(mesh());
 768       scalarField globalPatchMaterials =
 769           crackerMesh.globalCrackField(localPatchMaterials);
 770       const labelList& gcfa = crackerMesh.globalCrackFaceAddressing();
 771       label globalIndex = crackerMesh.localCrackStart();
 772
 773       forAll(mesh().boundaryMesh()[patchI], facei)
 774         {
 775           label ownerMat = label(globalPatchMaterials[globalIndex]);
 776           label neighbourMat = label(globalPatchMaterials[gcfa[globalIndex]]);
 777
 778           if (ownerMat != neighbourMat)
 779           {
 780               result.boundaryField()[patchI][facei] =
 781                   interfaceGIIc(ownerMat, neighbourMat);
 782           }
 783           globalIndex++;
 784         }
 785         }
 786
 787     else if (mesh().boundary()[patchI].coupled())
 788       {
 789         const scalarField ownerMatField =
 790             materials_.boundaryField()[patchI].internalField();
 791         const scalarField neighbourMatField =
 792             materials_.boundaryField()[patchI].patchNeighbourField();
 793
 794         forAll(mesh().boundaryMesh()[patchI], facei)
 795           {
 796                 label ownerMat = label(ownerMatField[facei]);
 797                 label neighbourMat = label(neighbourMatField[facei]);
 798
 799                 if (ownerMat != neighbourMat)
 800           {
 801               result.boundaryField()[patchI][facei] =
 802                   interfaceGIIc(ownerMat, neighbourMat);
 803           }
 804           }
 805       }
 806
 807     }
 808
 809     return tresult;
 810 }
 811
 812 Foam::scalar Foam::multiMaterialCohesiveLaw::interfaceSigmaMax
 813 (
 814     double mat1,
 815     double mat2
 816 ) const
 817 {
 818   label interfaceLawID = interfaceID(mat1, mat2);
 819
 820   return interfaceCohesiveLaws_[interfaceLawID].sigmaMaxValue();
 821 }
 822
 823 Foam::scalar Foam::multiMaterialCohesiveLaw::
 824 interfaceTauMax
 825 (
 826     double mat1,
 827     double mat2
 828 ) const
 829 {
 830   label interfaceLawID = interfaceID(mat1, mat2);
 831
 832   return interfaceCohesiveLaws_[interfaceLawID].tauMaxValue();
 833 }
 834
 835 Foam::scalar Foam::multiMaterialCohesiveLaw::interfaceGIc
 836 (
 837     double mat1,
 838     double mat2
 839 ) const
 840 {
 841   label interfaceLawID = interfaceID(mat1, mat2);
 842
 843   return interfaceCohesiveLaws_[interfaceLawID].GIcValue();
 844 }
 845
 846 Foam::scalar Foam::multiMaterialCohesiveLaw::interfaceGIIc
 847 (
 848     double mat1,
 849     double mat2
 850 ) const
 851 {
 852   label interfaceLawID = interfaceID(mat1, mat2);
 853
 854   return interfaceCohesiveLaws_[interfaceLawID].GIIcValue();
 855 }
 856
 857
 858 void Foam::multiMaterialCohesiveLaw::damageTractions
 859 (
 860  scalar& tN,
 861  scalar& tS,
 862  const scalar deltaN,
 863  const scalar deltaS,
 864  const scalar GI,
 865  const scalar GII,
 866  const label faceID,
 867  const scalarField& globalPatchMaterials
 868  ) const
 869 {
 870   // Find out which cohesive law does the face belong to
 871   const crackerFvMesh& crackerMesh = refCast<const crackerFvMesh>(mesh());
 872   const labelList& gcfa = crackerMesh.globalCrackFaceAddressing();
 873   label ownerMat = label(globalPatchMaterials[faceID]);
 874   label neighbourMat = label(globalPatchMaterials[gcfa[faceID]]);
 875
 876   if (ownerMat != neighbourMat)
 877     {
 878       label matID = interfaceID(ownerMat, neighbourMat);
 879
 880       // face is on multi-material interface
 881       interfaceCohesiveLaws_[matID].damageTractions
 882           (tN, tS, deltaN, deltaS, GI, GII, faceID, globalPatchMaterials);
 883     }
 884   else
 885     {
 886       // face is within one material
 887       // call material law function
 888       label matID = ownerMat;
 889       (*this)[matID].damageTractions
 890           (tN, tS, deltaN, deltaS, GI, GII, faceID, globalPatchMaterials);
 891     }
 892 }
 893
 894
 895 Foam::tmp<Foam::surfaceVectorField>
 896 Foam::multiMaterialCohesiveLaw::interfaceTraction
 897 (
 898  surfaceVectorField n,
 899  volVectorField U,
 900  volTensorField gradU,
 901  volScalarField mu,
 902  volScalarField lambda
 903  ) const
 904 {
 905   tmp<surfaceVectorField> tresult
 906     (
 907         new surfaceVectorField
 908         (
 909             IOobject
 910             (
 911             "traction",
 912             mesh().time().timeName(),
 913             mesh(),
 914             IOobject::NO_READ,
 915             IOobject::NO_WRITE
 916         ),
 917             mesh(),
 918         dimensionedVector("zero", dimForce/dimArea, vector::zero)
 919     )
 920     );
 921     surfaceVectorField& result = tresult();
 922
 923     surfaceScalarField mu1 = fvc::interpolate(mu);
 924     surfaceScalarField lambda1 = fvc::interpolate(lambda);
 925
 926     surfaceTensorField sGradU =
 927         ((I - n*n)&fvc::interpolate(gradU));
 928
 929     vectorField UI = U.internalField();
 930
 931     // result =
 932     //   (2*mu1 + lambda1)*fvc::snGrad(U)
 933     //   - (mu1 + lambda1)*(n&sGradU)
 934     //   + mu1*(sGradU&n)
 935     //   + lambda1*tr(sGradU)*n;
 936     result =
 937         2*mu1
 938         *(n&fvc::interpolate(symm(gradU)))
 939         + lambda1*n*fvc::interpolate(tr(gradU));
 940
 941     //    const labelList& owner = mesh().owner();
 942     //const labelList& neighbour = mesh().neighbour();
 943
 944     /*        forAll(result.internalField(), faceI)
 945     {
 946             label ownerMat = materials_[owner[faceI]];
 947             label neighbourMat = materials_[neighbour[faceI]];
 948
 949             if (ownerMat != neighbourMat)
 950             {
 951         vector ownN = n.internalField()[faceI];
 952         vector ngbN = -n.internalField()[faceI];
 953
 954         tensor ownSGradU =
 955         ((I-ownN*ownN)&gradU.internalField()[owner[faceI]]);
 956         tensor ngbSGradU =
 957         ((I-ngbN*ngbN)&gradU.internalField()[neighbour[faceI]]);
 958
 959         scalar ownTrSGradUt = tr(ownSGradU&(I-ownN*ownN));
 960         scalar ngbTrSGradUt = tr(ngbSGradU&(I-ngbN*ngbN));
 961
 962         vector ownSGradUn = (ownSGradU&ownN);
 963         vector ngbSGradUn = (ngbSGradU&ngbN);
 964
 965         scalar ownMu = mu.internalField()[owner[faceI]];
 966         scalar ngbMu = mu.internalField()[neighbour[faceI]];
 967
 968         scalar ownLambda = lambda.internalField()[owner[faceI]];
 969         scalar ngbLambda = lambda.internalField()[neighbour[faceI]];
 970
 971         vector ownUt = ((I-ownN*ownN)&UI[owner[faceI]]);
 972         vector ngbUt = ((I-ngbN*ngbN)&UI[neighbour[faceI]]);
 973
 974         vector ownUn = ownN*(ownN&U.internalField()[owner[faceI]]);
 975         vector ngbUn = ngbN*(ngbN&U.internalField()[neighbour[faceI]]);
 976
 977         scalar ownDn = mesh().weights().internalField()[faceI]
 978                  * (1.0/mesh().deltaCoeffs().internalField()[faceI]);
 979         scalar ngbDn =
 980         (1.0/mesh().deltaCoeffs().internalField()[faceI]) - ownDn;
 981
 982         scalar own2ML = (2*ownMu + ownLambda);
 983         scalar ngb2ML = (2*ngbMu + ngbLambda);
 984
 985         scalar face2ML = 1.0/((1.0-mesh().weights().internalField()[faceI])
 986                    / own2ML + mesh().weights().internalField()[faceI]/ngb2ML);
 987         scalar faceMu = 1.0/((1.0-mesh().weights().internalField()[faceI])
 988                   / ownMu + mesh().weights().internalField()[faceI]/ngbMu);
 989
 990         vector tractionN =
 991                     face2ML*(ownUn - ngbUn)/(ownDn+ngbDn)
 992                   + ((own2ML*ngbDn*ngbLambda*ngbN*ngbTrSGradUt)
 993                   + (ngb2ML*ownDn*ownLambda*ngbN*ownTrSGradUt))
 994                   / ((own2ML*ngbDn) + (ngb2ML*ownDn));
 995
 996         vector tractionT =
 997                 faceMu*(ownUt - ngbUt)/(ownDn+ngbDn)
 998                   + ((ownMu*ngbMu*ngbDn*ngbSGradUn)
 999                   + (ownMu*ngbMu*ownDn*ownSGradUn))
1000                   / ((ownMu*ngbDn) + (ngbMu*ownDn));
1001
1002         //      Info << tractionN << " *** " << tractionT <<endl;
1003
1004         result[faceI] = tractionN + tractionT;
1005         //vector traction = tractionN + tractionT;
1006
1007         //vector ownDelta =
1008         mesh().C()[neighbour[faceI]] - mesh().C()[owner[faceI]];
1009         }
1010         }
1011     */
1012     // philipc
1013     // cracks and bi-material interfaces may be on processor boundaries
1014     // so we must correct tractions there too
1015     forAll(mesh().boundary(), patchI)
1016       {
1017         /*
1018  //if (mesh().boundaryMesh()[patchI].type() == processorPolyPatch::typeName)
1019         if (mesh().boundary()[patchI].coupled())
1020           {
1021         const scalarField ownerMatField =
1022         materials_.boundaryField()[patchI].patchInternalField();
1023         const scalarField neighbourMatField =
1024         materials_.boundaryField()[patchI].patchNeighbourField();
1025         const labelList& faceCells = mesh().boundary()[patchI].faceCells();
1026         const tensorField ngbGradUField =
1027         gradU.boundaryField()[patchI].patchNeighbourField();
1028         const scalarField ngbMuField =
1029         mu.boundaryField()[patchI].patchNeighbourField();
1030         const scalarField ngbLambdaField =
1031         lambda.boundaryField()[patchI].patchNeighbourField();
1032         const vectorField ngbUField =
1033         U.boundaryField()[patchI].patchNeighbourField();
1034         const scalarField weights = mesh().weights().boundaryField()[patchI];
1035
1036         forAll(mesh().boundaryMesh()[patchI], faceI)
1037           {
1038             label faceCelli = faceCells[faceI];
1039             label ownerMat = ownerMatField[faceI];
1040             label neighbourMat = neighbourMatField[faceI];
1041
1042             if (ownerMat != neighbourMat)
1043               {
1044             const vector& ownN = n.boundaryField()[patchI][faceI];
1045             const vector& ngbN = -n.boundaryField()[patchI][faceI];
1046             const tensor& ngbGradU = ngbGradUField[faceI];
1047             const vector& ngbU = ngbUField[faceI];
1048
1049             //tensor ownSGradU =
1050 //            ((I-ownN*ownN)&gradU.boundaryField()[patchI][owner[faceI]]);
1051             tensor ownSGradU = ((I-ownN*ownN)&gradU.internalField()[faceCelli]);
1052             //tensor ngbSGradU =
1053             // ((I-ngbN*ngbN)&gradU.boundaryField()[patchI][neighbour[faceI]]);
1054             tensor ngbSGradU = ((I-ngbN*ngbN)&ngbGradU);
1055
1056             scalar ownTrSGradUt = tr(ownSGradU&(I-ownN*ownN));
1057             scalar ngbTrSGradUt = tr(ngbSGradU&(I-ngbN*ngbN));
1058
1059             vector ownSGradUn = (ownSGradU&ownN);
1060             vector ngbSGradUn = (ngbSGradU&ngbN);
1061
1062             scalar ownMu = mu.internalField()[faceCelli];
1063             //scalar ngbMu = mu.boundaryField()[patchI][neighbour[faceI]];
1064             scalar ngbMu = ngbMuField[faceI];
1065
1066             scalar ownLambda = lambda.internalField()[faceCelli];
1067             //scalar ngbLambda =
1068             //lambda.boundaryField()[patchI][neighbour[faceI]];
1069             scalar ngbLambda = ngbLambdaField[faceI];
1070
1071             vector ownUt = ((I-ownN*ownN)&UI[owner[faceI]]);
1072             //vector ngbUt = ((I-ngbN*ngbN)&UI[neighbour[faceI]]);
1073             vector ngbUt = ((I-ngbN*ngbN)&ngbU);
1074
1075             vector ownUn = ownN*(ownN&U.internalField()[faceCelli]);
1076             //vector ngbUn =
1077             //ngbN*(ngbN&U.boundaryField()[patchI][neighbour[faceI]]);
1078             vector ngbUn = ngbN*(ngbN&ngbU);
1079
1080             //scalar ownDn = mesh().weights().boundaryField()[patchI][faceI]
1081             scalar ownDn = weights[faceI]
1082               * (1.0/mesh().deltaCoeffs().boundaryField()[patchI][faceI]);
1083             scalar ngbDn =
1084             (1.0/mesh().deltaCoeffs().boundaryField()[patchI][faceI]) - ownDn;
1085
1086             scalar own2ML = (2*ownMu + ownLambda);
1087             scalar ngb2ML = (2*ngbMu + ngbLambda);
1088
1089             // scalar face2ML =
1090             // 1.0/((1.0-mesh().weights().boundaryField()[patchI][faceI])
1091             // / own2ML
1092             //+ mesh().weights().boundaryField()[patchI][faceI]/ngb2ML);
1093             scalar face2ML = 1.0/((1.0-weights[faceI])
1094                           / own2ML + weights[faceI]/ngb2ML);
1095             // scalar faceMu =
1096             // 1.0/((1.0-mesh().weights().boundaryField()[patchI][faceI])
1097             // / ownMu + mesh().weights().boundaryField()[patchI][faceI]/ngbMu);
1098             scalar faceMu = 1.0/((1.0-weights[faceI])
1099                          / ownMu + weights[faceI]/ngbMu);
1100
1101             vector tractionN =
1102               face2ML*(ownUn - ngbUn)/(ownDn+ngbDn)
1103               + ((own2ML*ngbDn*ngbLambda*ngbN*ngbTrSGradUt)
1104               + (ngb2ML*ownDn*ownLambda*ngbN*ownTrSGradUt))
1105               / ((own2ML*ngbDn) + (ngb2ML*ownDn));
1106
1107             vector tractionT =
1108               faceMu*(ownUt - ngbUt)/(ownDn+ngbDn)
1109               + ((ownMu*ngbMu*ngbDn*ngbSGradUn)
1110               + (ownMu*ngbMu*ownDn*ownSGradUn))
1111               / ((ownMu*ngbDn) + (ngbMu*ownDn));
1112
1113             result.boundaryField()[patchI][faceI] = tractionN + tractionT;
1114               }
1115           }
1116           }
1117
1118           else*/
1119           if (mesh().boundary()[patchI].type() == cohesiveFvPatch::typeName)
1120           {
1121               // I think mu and lambda are wrong on new crack faces
1122               // so use internal cell values
1123               const scalarField muPatch =
1124                   mu.boundaryField()[patchI].patchInternalField();
1125               const scalarField lambdaPatch =
1126                   lambda.boundaryField()[patchI].patchInternalField();
1127               result.boundaryField()[patchI] =
1128                   (2*muPatch*n.boundaryField()[patchI]
1129                    &symm(gradU.boundaryField()[patchI]))
1130                   + (lambdaPatch*tr(gradU.boundaryField()[patchI])
1131                      *n.boundaryField()[patchI]);
1132           }
1133
1134       }
1135
1136      return tresult;
1137 }
1138
1139 void Foam::multiMaterialCohesiveLaw::correct()
1140 {
1141     PtrList<cohesiveLaw>& laws = *this;
1142
1143     forAll (laws, lawI)
1144     {
1145         laws[lawI].correct();
1146     }
1147 }
1148
1149
1150 // ************************************************************************* //