diff --git a/src/solidModels/Make/files b/src/solidModels/Make/files index ffe4cf9c2..c0b293a29 100644 --- a/src/solidModels/Make/files +++ b/src/solidModels/Make/files @@ -67,6 +67,8 @@ finiteVolume = finiteVolume $(finiteVolume)/gradSchemes/leastSquaresSolidInterfaceGrad/leastSquaresSolidInterfaceGrads.C $(finiteVolume)/gradSchemes/leastSquaresSolidInterfaceGrad/leastSquaresSolidInterfaceVectors.C +nonLinearGeometry/nonLinearGeometry.C + rheologyLaws = $(constitutiveModel)/rheologyLaws $(rheologyLaws)/rheologyLaw/rheologyLaw.C $(rheologyLaws)/rheologyLaw/newRheologyLaw.C diff --git a/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C b/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C index f32e6900d..07c40caa6 100644 --- a/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C +++ b/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C @@ -23,7 +23,6 @@ License Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA \*---------------------------------------------------------------------------*/ -//#define DEBUG Info<<"In file "<<__FILE__<<" at line "<<__LINE__<("cracked", dict, p.size()); cracked_ = Field("cracked", dict, p.size()); } - if(dict.found("curTractionN")) + if (dict.found("curTractionN")) { curTractionN_ = scalarField("curTractionN", dict, p.size()); } - if(dict.found("curTractionS")) + if (dict.found("curTractionS")) { curTractionS_ = scalarField("curTractionS", dict, p.size()); } - if(dict.found("oldTractionN")) + if (dict.found("oldTractionN")) { oldTractionN_ = scalarField("oldTractionN", dict, p.size()); } - if(dict.found("oldTractionS")) + if (dict.found("oldTractionS")) { oldTractionS_ = scalarField("oldTractionS", dict, p.size()); } - if(dict.found("deltaN")) + if (dict.found("deltaN")) { deltaN_ = scalarField("deltaN", dict, p.size()); } - if(dict.found("deltaS")) + if (dict.found("deltaS")) { deltaS_ = scalarField("deltaS", dict, p.size()); } - if(dict.found("oldDeltaN")) + if (dict.found("oldDeltaN")) { oldDeltaN_ = scalarField("oldDeltaN", dict, p.size()); } - if(dict.found("oldDeltaS")) + if (dict.found("oldDeltaS")) { oldDeltaS_ = scalarField("oldDeltaS", dict, p.size()); } - if(dict.found("curDeltaN")) + if (dict.found("curDeltaN")) { curDeltaN_ = scalarField("curDeltaN", dict, p.size()); } - if(dict.found("curDeltaS")) + if (dict.found("curDeltaS")) { curDeltaS_ = scalarField("curDeltaS", dict, p.size()); } - if(dict.found("unloadingDeltaEff")) + if (dict.found("unloadingDeltaEff")) { unloadingDeltaEff_ = scalarField("unloadingDeltaEff", dict, p.size()); } - if(dict.found("currentGI")) + if (dict.found("currentGI")) { currentGI_ = scalarField("currentGI", dict, p.size()); } - if(dict.found("currentGII")) + if (dict.found("currentGII")) { currentGII_ = scalarField("currentGII", dict, p.size()); } - if(dict.found("oldGI")) + if (dict.found("oldGI")) { oldGI_ = scalarField("oldGI", dict, p.size()); } - if(dict.found("oldGII")) + if (dict.found("oldGII")) { oldGII_ = scalarField("oldGII", dict, p.size()); } @@ -417,7 +419,7 @@ tmp solidCohesiveFvPatchVectorField::crackingAndDamage() const forAll(tcrackingAndDamage(), facei) { - if(cracked_[facei]) + if (cracked_[facei]) { tcrackingAndDamage()[facei] = 2.0; } @@ -471,7 +473,7 @@ bool solidCohesiveFvPatchVectorField::cracking() label sumCracked = 0; forAll(globalCracked, facei) { - if(globalCracked[facei] > 0.0) + if (globalCracked[facei] > 0.0) { sumCracked++; } @@ -689,7 +691,7 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() { // we force the penalty factor to be calculated here for the first time // as all processors must call this at the same time - if(contact_ && patch().size() > 0) + if (contact_ && patch().size() > 0) { // force calculation of penalty factor here penaltyFactor(); @@ -727,11 +729,11 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() const crackerFvMesh& crackerMesh = refCast(mesh); // global patch material field is needed for multimaterial cohesive laws - if(updateGlobalPatchMaterials_) + if (updateGlobalPatchMaterials_) { updateGlobalPatchMaterials_ = false; - if(mesh.objectRegistry::foundObject("materials")) + if (mesh.objectRegistry::foundObject("materials")) { scalarField localPatchMaterials = patch().lookupPatchField("materials").patchInternalField(); @@ -755,7 +757,7 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() // Patch displacement vectorField UPatch = *this; - if(fieldName_ == "DU") + if (fieldName_ == "DU") { UPatch += patch().lookupPatchField("U"); @@ -779,29 +781,18 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() globalIndex++; } - // calculate the actual traction instead of using what was previously set - // traction_ = - // tractionBoundaryGradient().traction - // ( - // patch().lookupPatchField("grad("+fieldName_+")"), - // fieldName_, - // patch(), - // orthotropic_, - // NamedEnum::names[nonLinear_] - // ); - //globalIndex = crackerMesh.localCrackStart(); for(label i = 0; i < patch().size(); i++) { // update deltas curDeltaN_[i] = n[i] & delta[i]; curDeltaS_[i] = mag( (I - sqr(n[i])) & delta[i]); // shearing - if(explicitSeparationDistance_) + if (explicitSeparationDistance_) { deltaN_[i] = oldDeltaN_[i]; deltaS_[i] = oldDeltaS_[i]; - if(mag(deltaN_[i]) < SMALL) + if (mag(deltaN_[i]) < SMALL) { deltaN_[i] = 2*SMALL; } @@ -824,9 +815,9 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() // update energies // stop calculating after cracking for convergence (because crack might jump in and out of damaged/failed // only update energies if there is loading - if( !cracked_[i] && (deltaEff > (unloadingDeltaEff_[i]-SMALL)) ) + if ( !cracked_[i] && (deltaEff > (unloadingDeltaEff_[i]-SMALL)) ) { - if((curTractionN_[i]+oldTractionN_[i]) > 0.0) // if the average normal stress is tensile + if ((curTractionN_[i]+oldTractionN_[i]) > 0.0) // if the average normal stress is tensile { // trapezoidal rule currentGI_[i] = oldGI_[i] + ((0.5*(curTractionN_[i]+oldTractionN_[i]))*(deltaN_[i]-oldDeltaN_[i])); @@ -856,16 +847,16 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() // loading // if the effective delta is greater than unloadingDeltaEff then there is loading // unloadingDeltaEff is the maximum previsouly reached deltaEff - if(deltaEff > (unloadingDeltaEff_[i]-SMALL) ) + if (deltaEff > (unloadingDeltaEff_[i]-SMALL) ) { // at the moment loading and unloading are the same // if total energy is dissipated, then fully crack face - //if( currentG > GIc[i] || cracked_[i] ) //Gc ) + //if ( currentG > GIc[i] || cracked_[i] ) //Gc ) // propagation - if( ((currentGI_[i]/GIc[i]) + (currentGII_[i]/GIIc[i])) >= 1 ) + if ( ((currentGI_[i]/GIc[i]) + (currentGII_[i]/GIIc[i])) >= 1 ) { //Pout << "GIc[i] is " << GIc[i] << ", curG is " << currentG << endl; - if(!cracked_[i]) Pout << "Face " << i << " is fully cracked" << endl; + if (!cracked_[i]) Pout << "Face " << i << " is fully cracked" << endl; cracked_[i] = true; @@ -880,7 +871,7 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() // To-do: we must calculate actual distances curNormalTraction = 0.0; curTangentialTraction = vector::zero; - if( contact_ && deltaN_[i] <= 0.0 ) + if ( contact_ && deltaN_[i] <= 0.0 ) { curNormalTraction = deltaN_[i]*penaltyFactor(); //Info << "penaltyFactor() is " << penaltyFactor() << endl; @@ -901,9 +892,9 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() } // damging face with positive normal delta - else if( deltaN_[i] > 0.0 ) + else if ( deltaN_[i] > 0.0 ) { - if(cracked_[i]) Pout << "Face " << i << " is un-cracked" << endl; + if (cracked_[i]) Pout << "Face " << i << " is un-cracked" << endl; cracked_[i] = false; @@ -937,7 +928,7 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() // damaging faces with negative normal delta else { - if(cracked_[i]) Pout << "Face " << i << " is un-cracked" << endl; + if (cracked_[i]) Pout << "Face " << i << " is un-cracked" << endl; cracked_[i] = false; //Pout << "Contact and shearing face " << i << endl; @@ -952,7 +943,7 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() // Simple penalty condition scalar penaltyFac = 0.0; - if(contact_) + if (contact_) { penaltyFac = penaltyFactor(); } @@ -997,15 +988,15 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() this->refGrad() = - tractionBoundaryGradient() - ( - traction_, - scalarField(traction_.size(), 0.0), - word(fieldName_), - patch(), - orthotropic_, - NamedEnum::names[nonLinear_] - )(); + tractionBoundaryGradient() + ( + traction_, + scalarField(traction_.size(), 0.0), + word(fieldName_), + patch(), + orthotropic_, + nonLinearGeometry::nonLinearNames_[nonLinear_] + )(); directionMixedFvPatchVectorField::updateCoeffs(); } @@ -1016,12 +1007,12 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() // approx penaltyFactor from mechanical properties // this can then be scaled using the penaltyScale // to-do: write equivalent for orthotropic - if(orthotropic_) - { - FatalError << "solidCohesiveFvPatchVectorField::calcPenaltyFactor()" - << " has yet to be written for orthotropic" - << exit(FatalError); - } + if (orthotropic_) + { + FatalError << "solidCohesiveFvPatchVectorField::calcPenaltyFactor()" + << " has yet to be written for orthotropic" + << exit(FatalError); + } const label patchID = patch().index(); const fvMesh& mesh = patch().boundaryMesh().mesh(); @@ -1080,22 +1071,11 @@ void solidCohesiveFvPatchVectorField::write(Ostream& os) const os.writeKeyword("curTimeIndex") << curTimeIndex_ << token::END_STATEMENT << nl; os.writeKeyword("contact") << contact_<< token::END_STATEMENT << nl; os.writeKeyword("penaltyScale") << penaltyScale_ << token::END_STATEMENT << nl; - os.writeKeyword("nonLinear") << nonLinearNames_[nonLinear_] << token::END_STATEMENT << nl; + os.writeKeyword("nonLinear") << nonLinearGeometry::nonLinearNames_[nonLinear_] << token::END_STATEMENT << nl; os.writeKeyword("orthotropic") << orthotropic_ << token::END_STATEMENT << nl; } -template<> -const char* Foam::NamedEnum::names[] = - { - "off", - "updatedLagrangian", - "totalLagrangian" - }; - -const Foam::NamedEnum -Foam::solidCohesiveFvPatchVectorField::nonLinearNames_; - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // makePatchTypeField(fvPatchVectorField, solidCohesiveFvPatchVectorField); diff --git a/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.H b/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.H index e34fe387d..731984dc1 100644 --- a/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.H +++ b/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.H @@ -29,7 +29,7 @@ Description Cohesive law fv patch field for arbitrary crack procedure, where the mode mixity is allowed to vary depending on how the crack opens. - + This procedure allows sigmaMax, tauMax, GIc, and GIIc to be independently specified. @@ -37,11 +37,11 @@ Description The dissipated fracture energy is constantly monitored, by integrating in time. - + Author Philip Cardiff UCD based on original Tukovic cohesive law boundary condition. - + SourceFiles solidCohesiveFvPatchVectorField.C @@ -52,6 +52,7 @@ SourceFiles #include "fvPatchFields.H" #include "directionMixedFvPatchFields.H" +#include "nonLinearGeometry.H" #include "cohesiveFvPatch.H" #include "Switch.H" #include "tractionBoundaryGradient.H" @@ -63,7 +64,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class solidCohesiveFvPatch Declaration + Class solidCohesiveFvPatch Declaration \*---------------------------------------------------------------------------*/ class solidCohesiveFvPatchVectorField @@ -115,7 +116,7 @@ class solidCohesiveFvPatchVectorField scalarField currentGII_; scalarField oldGII_; - + //- Current time index label curTimeIndex_; @@ -125,8 +126,8 @@ class solidCohesiveFvPatchVectorField scalar penaltyScale_; scalar frictionCoeff_; // Coulomb friction coefficient - // If yes, cohesive traction will be calculated - // using old separation distance + // If yes, cohesive traction will be calculated + // using old separation distance Switch explicitSeparationDistance_; scalarField curDeltaN_; scalarField curDeltaS_; @@ -135,18 +136,8 @@ class solidCohesiveFvPatchVectorField bool updateGlobalPatchMaterials_; scalarField globalPatchMaterials_; - //- non-linear solver options - enum nonLinearType - { - OFF, - UPDATED_LAGRANGIAN, - TOTAL_LAGRANGIAN - }; - - static const NamedEnum nonLinearNames_; - - //- if it is a non linear solver - nonLinearType nonLinear_; + //- Is it a non linear solver + nonLinearGeometry::nonLinearType nonLinear_; //- if it is an orthropic solver Switch orthotropic_; @@ -158,13 +149,13 @@ class solidCohesiveFvPatchVectorField void calcPenaltyFactor(); virtual scalar penaltyFactor() { - if(!penaltyFactorPtr_) - { - calcPenaltyFactor(); - } + if(!penaltyFactorPtr_) + { + calcPenaltyFactor(); + } - return *penaltyFactorPtr_; - } + return *penaltyFactorPtr_; + } public: diff --git a/src/solidModels/contactModels/frictionContactModels/frictionContactModel/frictionContactModel.H b/src/solidModels/contactModels/frictionContactModels/frictionContactModel/frictionContactModel.H index 515501b13..b3959a886 100644 --- a/src/solidModels/contactModels/frictionContactModels/frictionContactModel/frictionContactModel.H +++ b/src/solidModels/contactModels/frictionContactModels/frictionContactModel/frictionContactModel.H @@ -97,7 +97,7 @@ protected: //- Return reference to mesh const fvMesh& mesh() const { - return patch_.boundaryMesh().mesh(); + return patch_.boundaryMesh().mesh(); } @@ -116,14 +116,22 @@ public: dictionary, ( const word name, - const fvPatch& patch, + const fvPatch& patch, const dictionary& dict, - const label masterPatchID, - const label slavePatchID, - const label masterFaceZoneID, - const label slaveFaceZoneID + const label masterPatchID, + const label slavePatchID, + const label masterFaceZoneID, + const label slaveFaceZoneID ), - (name, patch, dict, masterPatchID, slavePatchID, masterFaceZoneID, slaveFaceZoneID) + ( + name, + patch, + dict, + masterPatchID, + slavePatchID, + masterFaceZoneID, + slaveFaceZoneID + ) ); @@ -135,10 +143,10 @@ public: const word& name, const fvPatch& patch, const dictionary& dict, - const label masterPatchID, - const label slavePatchID, - const label masterFaceZoneID, - const label slaveFaceZoneID + const label masterPatchID, + const label slavePatchID, + const label masterFaceZoneID, + const label slaveFaceZoneID ); @@ -150,10 +158,10 @@ public: const word& name, const fvPatch& patch, const dictionary& dict, - const label masterPatchID, - const label slavePatchID, - const label masterFaceZoneID, - const label slaveFaceZoneID + const label masterPatchID, + const label slavePatchID, + const label masterFaceZoneID, + const label slaveFaceZoneID ); @@ -173,18 +181,18 @@ public: //- Correct contatc model virtual void correct - ( - const vectorField& slavePressure, - const PrimitivePatch& masterFaceZonePatch, - const PrimitivePatch& slaveFaceZonePatch, - const intersection::algorithm alg, - const intersection::direction dir, - const word interpolationMethod, - const word fieldName, - const Switch orthotropic, - const word nonLinear, - const vectorField& slaveFaceNormals - ) = 0; + ( + const vectorField& slavePressure, + const PrimitivePatch& masterFaceZonePatch, + const PrimitivePatch& slaveFaceZonePatch, + const intersection::algorithm alg, + const intersection::direction dir, + const word interpolationMethod, + const word fieldName, + const Switch orthotropic, + const word nonLinear, + const vectorField& slaveFaceNormals + ) = 0; //- Return slave friction displacement virtual const vectorField& slaveDisp() const = 0; @@ -198,33 +206,33 @@ public: //- Return master patch ID virtual label masterPatchID() const { - return masterPatchID_; - } + return masterPatchID_; + } //- Return master patch ID virtual label slavePatchID() const { - return slavePatchID_; - } + return slavePatchID_; + } //- Return master face zone ID virtual label masterFaceZoneID() const { - return masterFaceZoneID_; - } + return masterFaceZoneID_; + } //- Return master face zone ID virtual label slaveFaceZoneID() const { - return slaveFaceZoneID_; - } + return slaveFaceZoneID_; + } //- Return stick slip faces field // virtual volScalarField& stickSlipFaces() virtual scalarField& stickSlipFaces() { - return stickSlipFaces_; - } + return stickSlipFaces_; + } //- Write model dictionary virtual void writeDict(Ostream& os) const {}; diff --git a/src/solidModels/fvPatchFields/fixedDisplacementOrSolidTraction/fixedDisplacementOrSolidTractionFvPatchVectorField.H b/src/solidModels/fvPatchFields/fixedDisplacementOrSolidTraction/fixedDisplacementOrSolidTractionFvPatchVectorField.H index 7c45f7833..a0edaa2bb 100644 --- a/src/solidModels/fvPatchFields/fixedDisplacementOrSolidTraction/fixedDisplacementOrSolidTractionFvPatchVectorField.H +++ b/src/solidModels/fvPatchFields/fixedDisplacementOrSolidTraction/fixedDisplacementOrSolidTractionFvPatchVectorField.H @@ -48,6 +48,7 @@ Author #include "fvPatchFields.H" #include "directionMixedFvPatchFields.H" +#include "nonLinearGeometry.H" #include "Switch.H" #include "interpolationTable.H" @@ -78,18 +79,8 @@ class fixedDisplacementOrSolidTractionFvPatchVectorField //- Displacement vectorField displacement_; - //- non-linear solver options - enum nonLinearType - { - OFF, - UPDATED_LAGRANGIAN, - TOTAL_LAGRANGIAN - }; - - static const NamedEnum nonLinearNames_; - - //- if it is a non linear solver - nonLinearType nonLinear_; + //- Is it a non linear solver + nonLinearGeometry::nonLinearType nonLinear_; //- if it is an orthropic solver Switch orthotropic_; diff --git a/src/solidModels/fvPatchFields/fixedDisplacementZeroShear/fixedDisplacementZeroShearFvPatchVectorField.C b/src/solidModels/fvPatchFields/fixedDisplacementZeroShear/fixedDisplacementZeroShearFvPatchVectorField.C index eb01266b6..bf1a8661d 100644 --- a/src/solidModels/fvPatchFields/fixedDisplacementZeroShear/fixedDisplacementZeroShearFvPatchVectorField.C +++ b/src/solidModels/fvPatchFields/fixedDisplacementZeroShear/fixedDisplacementZeroShearFvPatchVectorField.C @@ -43,7 +43,8 @@ namespace Foam // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -fixedDisplacementZeroShearFvPatchVectorField::fixedDisplacementZeroShearFvPatchVectorField +fixedDisplacementZeroShearFvPatchVectorField:: +fixedDisplacementZeroShearFvPatchVectorField ( const fvPatch& p, const DimensionedField& iF @@ -51,12 +52,13 @@ fixedDisplacementZeroShearFvPatchVectorField::fixedDisplacementZeroShearFvPatchV : directionMixedFvPatchVectorField(p, iF), fieldName_("undefined"), - nonLinear_(OFF), + nonLinear_(nonLinearGeometry::OFF), orthotropic_(false) {} -fixedDisplacementZeroShearFvPatchVectorField::fixedDisplacementZeroShearFvPatchVectorField +fixedDisplacementZeroShearFvPatchVectorField:: +fixedDisplacementZeroShearFvPatchVectorField ( const fixedDisplacementZeroShearFvPatchVectorField& ptf, const fvPatch& p, @@ -71,7 +73,8 @@ fixedDisplacementZeroShearFvPatchVectorField::fixedDisplacementZeroShearFvPatchV {} -fixedDisplacementZeroShearFvPatchVectorField::fixedDisplacementZeroShearFvPatchVectorField +fixedDisplacementZeroShearFvPatchVectorField:: +fixedDisplacementZeroShearFvPatchVectorField ( const fvPatch& p, const DimensionedField& iF, @@ -80,45 +83,48 @@ fixedDisplacementZeroShearFvPatchVectorField::fixedDisplacementZeroShearFvPatchV : directionMixedFvPatchVectorField(p, iF), fieldName_(dimensionedInternalField().name()), - nonLinear_(OFF), + nonLinear_(nonLinearGeometry::OFF), orthotropic_(false) { //- check if traction boundary is for non linear solver - if(dict.found("nonLinear")) - { - nonLinear_ = nonLinearNames_.read(dict.lookup("nonLinear"));; + if (dict.found("nonLinear")) + { + nonLinear_ = nonLinearGeometry::nonLinearNames_.read + ( + dict.lookup("nonLinear") + ); - if(nonLinear_ == UPDATED_LAGRANGIAN) - { - Info << "\tnonLinear set to updated Lagrangian" - << endl; - } - else if(nonLinear_ == TOTAL_LAGRANGIAN) - { - Info << "\tnonLinear set to total Lagrangian" - << endl; - } - } + if (nonLinear_ == nonLinearGeometry::UPDATED_LAGRANGIAN) + { + Info<< "\tnonLinear set to updated Lagrangian" + << endl; + } + else if (nonLinear_ == nonLinearGeometry::TOTAL_LAGRANGIAN) + { + Info<< "\tnonLinear set to total Lagrangian" + << endl; + } + } - if(dict.found("orthotropic")) - { - orthotropic_ = Switch(dict.lookup("orthotropic")); - Info << "\t\torthotropic set to " << orthotropic_ << endl; - } + if (dict.found("orthotropic")) + { + orthotropic_ = Switch(dict.lookup("orthotropic")); + Info << "\t\torthotropic set to " << orthotropic_ << endl; + } //- the leastSquares has zero non-orthogonal correction //- on the boundary //- so the gradient scheme should be extendedLeastSquares - // if(Foam::word(dimensionedInternalField().mesh()..schemesDict().gradScheme("grad(" + fieldName_ + ")")) != "extendedLeastSquares") + // if (Foam::word(dimensionedInternalField().mesh()..schemesDict().gradScheme("grad(" + fieldName_ + ")")) != "extendedLeastSquares") // { -// Warning << "The gradScheme for " << fieldName_ -// << " should be \"extendedLeastSquares 0\" for the boundary " -// << "non-orthogonal correction to be right" << endl; +// Warning << "The gradScheme for " << fieldName_ +// << " should be \"extendedLeastSquares 0\" for the boundary " +// << "non-orthogonal correction to be right" << endl; // } this->refGrad() = vector::zero; - - vectorField n = patch().nf(); + + vectorField n = patch().nf(); this->valueFraction() = sqr(n); if (dict.found("value")) @@ -128,18 +134,18 @@ fixedDisplacementZeroShearFvPatchVectorField::fixedDisplacementZeroShearFvPatchV else { FatalError << "value entry not found for patch " << patch().name() - << exit(FatalError); + << exit(FatalError); } this->refValue() = *this; Field normalValue = transform(valueFraction(), refValue()); - + Field gradValue = this->patchInternalField() + refGrad()/this->patch().deltaCoeffs(); - + Field transformGradValue = transform(I - valueFraction(), gradValue); - + Field::operator=(normalValue + transformGradValue); } @@ -198,14 +204,14 @@ void fixedDisplacementZeroShearFvPatchVectorField::updateCoeffs() //- or fix deformed normal //- I should add an option to choose which normal to fix - // if(nonLinear_ != OFF) + // if (nonLinear_ != OFF) // { // tensorField F = I + gradField; // tensorField Finv = inv(F); - // scalarField J = det(F); + // scalarField J = det(F); // vectorField nCurrent = Finv & n; // nCurrent /= mag(nCurrent); - // this->valueFraction() = sqr(nCurrent); + // this->valueFraction() = sqr(nCurrent); // } refGrad() = tractionBoundaryGradient() @@ -215,7 +221,7 @@ void fixedDisplacementZeroShearFvPatchVectorField::updateCoeffs() word(fieldName_), patch(), orthotropic_, - NamedEnum::names[nonLinear_] + nonLinearGeometry::nonLinearNames_[nonLinear_] )(); directionMixedFvPatchVectorField::updateCoeffs(); @@ -226,24 +232,20 @@ void fixedDisplacementZeroShearFvPatchVectorField::updateCoeffs() void fixedDisplacementZeroShearFvPatchVectorField::write(Ostream& os) const { directionMixedFvPatchVectorField::write(os); - os.writeKeyword("nonLinear") << nonLinearNames_[nonLinear_] << token::END_STATEMENT << nl; + os.writeKeyword("nonLinear") + << nonLinearGeometry::nonLinearNames_[nonLinear_] + << token::END_STATEMENT << nl; } -template<> -const char* Foam::NamedEnum::names[] = - { - "off", - "updatedLagrangian", - "totalLagrangian" - }; - -const Foam::NamedEnum -Foam::fixedDisplacementZeroShearFvPatchVectorField::nonLinearNames_; - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -makePatchTypeField(fvPatchVectorField, fixedDisplacementZeroShearFvPatchVectorField); +makePatchTypeField +( + fvPatchVectorField, + fixedDisplacementZeroShearFvPatchVectorField +); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/solidModels/fvPatchFields/fixedDisplacementZeroShearOrSolidTraction/fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField.H b/src/solidModels/fvPatchFields/fixedDisplacementZeroShearOrSolidTraction/fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField.H index aceea19fe..7eaffc093 100644 --- a/src/solidModels/fvPatchFields/fixedDisplacementZeroShearOrSolidTraction/fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField.H +++ b/src/solidModels/fvPatchFields/fixedDisplacementZeroShearOrSolidTraction/fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField.H @@ -48,6 +48,7 @@ Author #include "fvPatchFields.H" #include "directionMixedFvPatchFields.H" +#include "nonLinearGeometry.H" #include "Switch.H" #include "interpolationTable.H" @@ -81,20 +82,10 @@ class fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField //- Normal to fix for displacement phase vector fixedNormal_; - //- non-linear solver options - enum nonLinearType - { - OFF, - UPDATED_LAGRANGIAN, - TOTAL_LAGRANGIAN - }; - - static const NamedEnum nonLinearNames_; + //- Is it a non linear solver + nonLinearGeometry::nonLinearType nonLinear_; - //- if it is a non linear solver - nonLinearType nonLinear_; - - //- if it is an orthropic solver + //- Is it an orthropic solver Switch orthotropic_; //- The time series - 0 for traction, anything else for displacement. diff --git a/src/solidModels/fvPatchFields/fixedRotation/fixedRotationFvPatchVectorField.C b/src/solidModels/fvPatchFields/fixedRotation/fixedRotationFvPatchVectorField.C index 50d5a3c3b..d7b37236c 100644 --- a/src/solidModels/fvPatchFields/fixedRotation/fixedRotationFvPatchVectorField.C +++ b/src/solidModels/fvPatchFields/fixedRotation/fixedRotationFvPatchVectorField.C @@ -50,8 +50,8 @@ fixedRotationFvPatchVectorField::fixedRotationFvPatchVectorField rotationAngle_(0.0), rotationAxis_(vector::zero), rotationOrigin_(vector::zero), - origCf_(0,vector::zero), - nonLinear_(OFF) + origCf_(0, vector::zero), + nonLinear_(nonLinearGeometry::OFF) {} @@ -69,7 +69,7 @@ fixedRotationFvPatchVectorField::fixedRotationFvPatchVectorField rotationAxis_(ptf.rotationAxis_), rotationOrigin_(ptf.rotationOrigin_), origCf_(ptf.origCf_), - nonLinear_(OFF) + nonLinear_(nonLinearGeometry::OFF) {} @@ -86,33 +86,36 @@ fixedRotationFvPatchVectorField::fixedRotationFvPatchVectorField rotationAxis_(dict.lookup("rotationAxis")), rotationOrigin_(dict.lookup("rotationOrigin")), origCf_(patch().patch().faceCentres()), - nonLinear_(OFF) + nonLinear_(nonLinearGeometry::OFF) { //- check if traction boundary is for non linear solver - if(dict.found("nonLinear")) - { - nonLinear_ = nonLinearNames_.read(dict.lookup("nonLinear"));; + if (dict.found("nonLinear")) + { + nonLinear_ = nonLinearGeometry::nonLinearNames_.read + ( + dict.lookup("nonLinear") + ); - if(nonLinear_ == UPDATED_LAGRANGIAN) - { - Info << "\tnonLinear set to updated Lagrangian" - << endl; - } - else if(nonLinear_ == TOTAL_LAGRANGIAN) - { - Info << "\tnonLinear set to total Lagrangian" - << endl; - } - } + if (nonLinear_ == nonLinearGeometry::UPDATED_LAGRANGIAN) + { + Info<< "\tnonLinear set to updated Lagrangian" + << endl; + } + else if (nonLinear_ == nonLinearGeometry::TOTAL_LAGRANGIAN) + { + Info << "\tnonLinear set to total Lagrangian" + << endl; + } + } //- the leastSquares has zero non-orthogonal correction //- on the boundary //- so the gradient scheme should be extendedLeastSquares - if(Foam::word(dimensionedInternalField().mesh().schemesDict().gradScheme("grad(" + fieldName_ + ")")) != "extendedLeastSquares") + if (word(dimensionedInternalField().mesh().schemesDict().gradScheme("grad(" + fieldName_ + ")")) != "extendedLeastSquares") { - Warning << "The gradScheme for " << fieldName_ - << " should be \"extendedLeastSquares 0\" for the boundary " - << "non-orthogonal correction to be right" << endl; + Warning << "The gradScheme for " << fieldName_ + << " should be \"extendedLeastSquares 0\" for the boundary " + << "non-orthogonal correction to be right" << endl; } } @@ -155,7 +158,7 @@ snGrad() const { //- fixedValue snGrad with no correction // return (*this - patchInternalField())*this->patch().deltaCoeffs(); - + const fvPatchField& gradField = patch().lookupPatchField ( @@ -168,9 +171,9 @@ snGrad() const //- correction vector vectorField k = delta - n*(n&delta); - return + return ( - *this + *this - (patchInternalField() + (k&gradField.patchInternalField())) )*this->patch().deltaCoeffs(); } @@ -185,23 +188,23 @@ void fixedRotationFvPatchVectorField::updateCoeffs() tensor rotMat = RodriguesRotation(rotationAxis_, rotationAngle_); vectorField oldFaceCentres = dimensionedInternalField().mesh().C().boundaryField()[patch().index()]; - + vectorField newFaceCentres = (rotMat & (oldFaceCentres - rotationOrigin_)) + rotationOrigin_; vectorField disp = newFaceCentres - oldFaceCentres; - if(fieldName_ == "DU") + if (fieldName_ == "DU") { const fvPatchField& U = patch().lookupPatchField("U"); disp -= U; } - else if(fieldName_ != "U") + else if (fieldName_ != "U") { FatalError << "The displacement field should be U or DU" << exit(FatalError); } - + fvPatchField::operator== ( disp @@ -232,25 +235,21 @@ gradientBoundaryCoeffs() const void fixedRotationFvPatchVectorField::write(Ostream& os) const { fixedValueFvPatchVectorField::write(os); - os.writeKeyword("rotationAngle") << rotationAngle_ << token::END_STATEMENT << nl; - os.writeKeyword("rotationAxis") << rotationAxis_ << token::END_STATEMENT << nl; - os.writeKeyword("rotationOrigin") << rotationOrigin_ << token::END_STATEMENT << nl; - os.writeKeyword("nonLinear") << nonLinearNames_[nonLinear_] << token::END_STATEMENT << nl; + os.writeKeyword("rotationAngle") + << rotationAngle_ + << token::END_STATEMENT << nl; + os.writeKeyword("rotationAxis") + << rotationAxis_ + << token::END_STATEMENT << nl; + os.writeKeyword("rotationOrigin") + << rotationOrigin_ + << token::END_STATEMENT << nl; + os.writeKeyword("nonLinear") + << nonLinearGeometry::nonLinearNames_[nonLinear_] + << token::END_STATEMENT << nl; } -template<> -const char* Foam::NamedEnum::names[] = - { - "off", - "updatedLagrangian", - "totalLagrangian" - }; - -const Foam::NamedEnum -Foam::fixedRotationFvPatchVectorField::nonLinearNames_; - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // makePatchTypeField diff --git a/src/solidModels/fvPatchFields/fixedRotation/fixedRotationFvPatchVectorField.H b/src/solidModels/fvPatchFields/fixedRotation/fixedRotationFvPatchVectorField.H index 0563008ea..36ddb1410 100644 --- a/src/solidModels/fvPatchFields/fixedRotation/fixedRotationFvPatchVectorField.H +++ b/src/solidModels/fvPatchFields/fixedRotation/fixedRotationFvPatchVectorField.H @@ -42,6 +42,7 @@ Author #include "fvPatchFields.H" #include "fixedValueFvPatchFields.H" +#include "nonLinearGeometry.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -49,7 +50,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class fixedRotationFvPatch Declaration + Class fixedRotationFvPatch Declaration \*---------------------------------------------------------------------------*/ class fixedRotationFvPatchVectorField @@ -74,18 +75,9 @@ class fixedRotationFvPatchVectorField //- Initial mesh position const vectorField origCf_; - //- non-linear solver options - enum nonLinearType - { - OFF, - UPDATED_LAGRANGIAN, - TOTAL_LAGRANGIAN - }; - - static const NamedEnum nonLinearNames_; + //- Is it a non linear solver + nonLinearGeometry::nonLinearType nonLinear_; - //- if it is a non linear solver - nonLinearType nonLinear_; public: @@ -160,16 +152,17 @@ public: // Access functions - virtual nonLinearType& nonLinear() + const nonLinearGeometry::nonLinearType& nonLinear() const { return nonLinear_; } - virtual const NamedEnum& nonLinearNames() + nonLinearGeometry::nonLinearType& nonLinear() { - return nonLinearNames_; + return nonLinear_; } + // Evaluation functions //- Return patch-normal gradient diff --git a/src/solidModels/fvPatchFields/solidContact/solidContactFvPatchVectorField.C b/src/solidModels/fvPatchFields/solidContact/solidContactFvPatchVectorField.C index 6a4c7a38c..5cc82a276 100644 --- a/src/solidModels/fvPatchFields/solidContact/solidContactFvPatchVectorField.C +++ b/src/solidModels/fvPatchFields/solidContact/solidContactFvPatchVectorField.C @@ -26,7 +26,6 @@ Class solidContactFvPatchVectorField \*---------------------------------------------------------------------------*/ -//#define DEBUG Pout<<"file "<<__FILE__<<" line "<<__LINE__<, PrimitivePatch > - ( - *slaveFaceZonePatchPtr_, // from zone - *masterFaceZonePatchPtr_, // to zone - alg_, - dir_ - ); - } - else if(interpolationMethod_ == "ggi") // still created if master is rigid because it is used to interpolate normals - { - Info << "\tInterpolation of traction from slave to master: ggi" << endl; - slaveToMasterGgiInterpolatorPtr_ = - // new ggiZoneInterpolation - new GGIInterpolation< PrimitivePatch< face, List, pointField >, PrimitivePatch< face, List, pointField > > - ( - // masterFaceZonePatch, // master zone - // slaveFaceZonePatch, // slave zone - *masterFaceZonePatchPtr_, // master zone - *slaveFaceZonePatchPtr_, // slave zone - tensorField(0), - tensorField(0), - vectorField(0), - 0.0, - 0.0, - true, - ggiInterpolation::AABB - ); - } - else - { - SeriousError << "\n\nTraction interpolation method '" - << interpolationMethod_ << "' not found!\n" - << "\nValid interpolation methods are:\n" - << "ggi\n" - << "patchToPatch" - << endl - << exit(FatalError); - } - } + if (contactActive_) + { + const fvMesh& mesh = patch().boundaryMesh().mesh(); + label masterFaceZoneID = mesh.faceZones().findZoneID(masterFaceZoneName_); + label slaveFaceZoneID = mesh.faceZones().findZoneID(slaveFaceZoneName_); + if (masterFaceZoneID == -1) + { + FatalError << "face zone " << masterFaceZoneName_ << " not found!" + << exit(FatalError); + } + else if (slaveFaceZoneID == -1) + { + FatalError << "face zone " << masterFaceZoneName_ << " not found!" + << exit(FatalError); + } + + normalContactModelPtr_ = + normalContactModel::New + ( + dict.lookup("normalContactModel"), + patch(), + dict, + patch().index(), // master + shadowPatchID_, // slave + masterFaceZoneID_, // master face zone ID + slaveFaceZoneID_, // slave face zone ID + *masterFaceZonePatchPtr_, + *slaveFaceZonePatchPtr_ + ).ptr(); + frictionContactModelPtr_ = + frictionContactModel::New + ( + dict.lookup("frictionContactModel"), + patch(), + dict, + patch().index(), // master + shadowPatchID_, // slave + masterFaceZoneID_, // master face zone ID + slaveFaceZoneID_ // slave face zone ID + ).ptr(); + + //- create zoneToZone or ggiZone for interpolation of traction from + // slave to master only required if the contact is not rigid + if (interpolationMethod_ == "patchToPatch") + { + Info << "\tInterpolation of traction from slave to master: patchToPatch" << endl; + // slaveToMasterPatchToPatchInterpolatorPtr_ = + // new zoneToZoneInterpolation + // ( + // // slaveFaceZonePatch, // from zone + // // masterFaceZonePatch, // to zone + // *slaveFaceZonePatchPtr_, // from zone + // *masterFaceZonePatchPtr_, // to zone + // alg_, + // dir_ + // ); + slaveToMasterPatchToPatchInterpolatorPtr_ = + new PatchToPatchInterpolation + < + PrimitivePatch, + PrimitivePatch + > + ( + *slaveFaceZonePatchPtr_, // from zone + *masterFaceZonePatchPtr_, // to zone + alg_, + dir_ + ); + } + else if (interpolationMethod_ == "ggi") // still created if master is rigid because it is used to interpolate normals + { + Info << "\tInterpolation of traction from slave to master: ggi" << endl; + slaveToMasterGgiInterpolatorPtr_ = + // new ggiZoneInterpolation + new GGIInterpolation< PrimitivePatch< face, List, pointField >, PrimitivePatch< face, List, pointField > > + ( + // masterFaceZonePatch, // master zone + // slaveFaceZonePatch, // slave zone + *masterFaceZonePatchPtr_, // master zone + *slaveFaceZonePatchPtr_, // slave zone + tensorField(0), + tensorField(0), + vectorField(0), + 0.0, + 0.0, + true, + ggiInterpolation::AABB + ); + } else - { - Info << "\t**The contact is not active, contatc patches treated as traction-free**" << endl; - valueFraction() = symmTensor::zero; - } + { + SeriousError << "\n\nTraction interpolation method '" + << interpolationMethod_ << "' not found!\n" + << "\nValid interpolation methods are:\n" + << "ggi\n" + << "patchToPatch" + << endl + << exit(FatalError); + } + } + else + { + Info << "\t**The contact is not active, contatc patches treated as traction-free**" << endl; + valueFraction() = symmTensor::zero; + } } // end of if master - + if (dict.found("refValue")) { this->refValue() = vectorField("refValue", dict, p.size()); @@ -376,7 +390,7 @@ solidContactFvPatchVectorField::solidContactFvPatchVectorField { this->refValue() = vector::zero; } - + if (dict.found("refGradient")) { this->refGrad() = vectorField("refGradient", dict, p.size()); @@ -388,7 +402,7 @@ solidContactFvPatchVectorField::solidContactFvPatchVectorField if (dict.found("valueFraction")) { - this->valueFraction() = + this->valueFraction() = symmTensorField("valueFraction", dict, p.size()); } else @@ -508,210 +522,215 @@ void solidContactFvPatchVectorField::updateCoeffs() } //Info << "updating contact patch " << patch().name() << " ..." << flush; - if(contactActive_) + if (contactActive_) { - // for the first updateCoeffs, the slave needs to grab the conatct law pointers - if(!normalContactModelPtr_) - { - if(!master_) - { - // grab pointers - const volVectorField& Ufield = - this->db().objectRegistry::lookupObject(fieldName_); - - const solidContactFvPatchVectorField& Upatch = - refCast(Ufield.boundaryField()[shadowPatchID_]); - Info << "\tSlave contact patch " << patch().name() << " grabbing normalContactModel pointer from master" << endl; - normalContactModelPtr_ = Upatch.normalContactModelPtr(); - if(!normalContactModelPtr_) FatalError << "\nThe patch " << patch().name() << " has a NULL normalContactModel pointer!" << exit(FatalError); + // for the first updateCoeffs, the slave needs to grab the conatct law pointers + if (!normalContactModelPtr_) + { + if (!master_) + { + // grab pointers + const volVectorField& Ufield = + this->db().objectRegistry::lookupObject(fieldName_); - Info << "\tSlave contact patch " << patch().name() << " grabbing frictionContactModel pointer from master" << endl; - frictionContactModelPtr_ = Upatch.frictionContactModelPtr(); - if(!frictionContactModelPtr_) - FatalError << "\nThe patch " << patch().name() << " has a NULL frictionContactModel pointer!" << exit(FatalError); + const solidContactFvPatchVectorField& Upatch = + refCast(Ufield.boundaryField()[shadowPatchID_]); + Info << "\tSlave contact patch " << patch().name() << " grabbing normalContactModel pointer from master" << endl; + normalContactModelPtr_ = Upatch.normalContactModelPtr(); + if (!normalContactModelPtr_) FatalError << "\nThe patch " << patch().name() << " has a NULL normalContactModel pointer!" << exit(FatalError); - // lookup master's correction frequency - correctionFreq_ = Upatch.correctionFreq(); + Info << "\tSlave contact patch " << patch().name() << " grabbing frictionContactModel pointer from master" << endl; + frictionContactModelPtr_ = Upatch.frictionContactModelPtr(); + if (!frictionContactModelPtr_) + FatalError << "\nThe patch " << patch().name() << " has a NULL frictionContactModel pointer!" << exit(FatalError); - // get face zone patch pointers - masterFaceZonePatchPtr_ = Upatch.masterFaceZonePatchPtr(); - slaveFaceZonePatchPtr_ = Upatch.slaveFaceZonePatchPtr(); + // lookup master's correction frequency + correctionFreq_ = Upatch.correctionFreq(); - // get patchToPatch and GGI interpolator pointers - one of them is NULL - slaveToMasterGgiInterpolatorPtr_ = Upatch.slaveToMasterGgiInterpolatorPtr(); - slaveToMasterPatchToPatchInterpolatorPtr_ = Upatch.slaveToMasterPatchToPatchInterpolatorPtr(); - } - else - { - FatalErrorIn("solidContactFvPatchVectorField::updateCoeff()") - << "NULL contactLaw\ncontactLaw not created by master." - << abort(FatalError); - } - } + // get face zone patch pointers + masterFaceZonePatchPtr_ = Upatch.masterFaceZonePatchPtr(); + slaveFaceZonePatchPtr_ = Upatch.slaveFaceZonePatchPtr(); - // if it is a new time step then reset iCorr - if(curTimeIndex_ != this->db().time().timeIndex()) - { - curTimeIndex_ = this->db().time().timeIndex(); - iCorr_ = 0; + // get patchToPatch and GGI interpolator pointers - one of them is NULL + slaveToMasterGgiInterpolatorPtr_ = Upatch.slaveToMasterGgiInterpolatorPtr(); + slaveToMasterPatchToPatchInterpolatorPtr_ = Upatch.slaveToMasterPatchToPatchInterpolatorPtr(); + } + else + { + FatalErrorIn("solidContactFvPatchVectorField::updateCoeff()") + << "NULL contactLaw\ncontactLaw not created by master." + << abort(FatalError); + } + } - // update old face zone points - if(master_ && nonLinear_ == UPDATED_LAGRANGIAN) - { - oldMasterFaceZonePoints_ = masterFaceZonePatchPtr_->localPoints(); - oldSlaveFaceZonePoints_ = slaveFaceZonePatchPtr_->localPoints(); - } - } + // if it is a new time step then reset iCorr + if (curTimeIndex_ != this->db().time().timeIndex()) + { + curTimeIndex_ = this->db().time().timeIndex(); + iCorr_ = 0; - // the time taken to correct the contact may not be negligible - // so reduce the correctiion frequency can speed up the simulation - if(iCorr_++ % correctionFreq_ == 0 - || forceCorrection_ ) - { - forceCorrection_ = false; + // update old face zone points + if (master_ && nonLinear_ == nonLinearGeometry::UPDATED_LAGRANGIAN) + { + oldMasterFaceZonePoints_ = masterFaceZonePatchPtr_->localPoints(); + oldSlaveFaceZonePoints_ = slaveFaceZonePatchPtr_->localPoints(); + } + } - // master moves faceZone patches to the current deformed position - if(master_) - { - moveFaceZonePatches(); - } + // the time taken to correct the contact may not be negligible + // so reduce the correctiion frequency can speed up the simulation + if (iCorr_++ % correctionFreq_ == 0 + || forceCorrection_ ) + { + forceCorrection_ = false; - // just one of the patches updates the interpolators and - // corrects the contact laws - if(master_) //patchInChargeOfCorrection_) // hmmm the master needs to do the correction - { - // update interpolation classes - if(!rigidMaster_) - { - if(slaveToMasterGgiInterpolatorPtr_) - { - slaveToMasterGgiInterpolatorPtr_->movePoints(); - } - else if(slaveToMasterPatchToPatchInterpolatorPtr_) - { - slaveToMasterPatchToPatchInterpolatorPtr_->movePoints(); - } - else - { - FatalErrorIn("solidContactFvPatchVectorField::updateCoeff()") - << "Neither patchToPatch or GGI interpolator found!" - << abort(FatalError); - } - } + // master moves faceZone patches to the current deformed position + if (master_) + { + moveFaceZonePatches(); + } - // correct laws - // the normal model sets what face normal to use - // on the slave e.g. use actual face normals, or master normals - // interpolated to the slave, undeformed/deformed normals. - // the friction model then uses these face normals - vectorField slaveFaceNormals(patch().boundaryMesh().mesh().boundary()[shadowPatchID_].size(), vector::zero); - normalContactModelPtr_->correct - ( - *masterFaceZonePatchPtr_, - *slaveFaceZonePatchPtr_, - alg_, - dir_, - fieldName_, - orthotropic_, - NamedEnum::names[nonLinear_], - slaveFaceNormals, - slaveToMasterGgiInterpolatorPtr_ - ); + // just one of the patches updates the interpolators and + // corrects the contact laws + if (master_) //patchInChargeOfCorrection_) // hmmm the master needs to do the correction + { + // update interpolation classes + if (!rigidMaster_) + { + if (slaveToMasterGgiInterpolatorPtr_) + { + slaveToMasterGgiInterpolatorPtr_->movePoints(); + } + else if (slaveToMasterPatchToPatchInterpolatorPtr_) + { + slaveToMasterPatchToPatchInterpolatorPtr_->movePoints(); + } + else + { + FatalErrorIn("solidContactFvPatchVectorField::updateCoeff()") + << "Neither patchToPatch or GGI interpolator found!" + << abort(FatalError); + } + } - frictionContactModelPtr_->correct - ( - normalContactModelPtr_->slavePressure(), - *masterFaceZonePatchPtr_, - *slaveFaceZonePatchPtr_, - alg_, - dir_, - interpolationMethod_, - fieldName_, - orthotropic_, - NamedEnum::names[nonLinear_], - slaveFaceNormals - ); - } + // correct laws + // the normal model sets what face normal to use + // on the slave e.g. use actual face normals, or master normals + // interpolated to the slave, undeformed/deformed normals. + // the friction model then uses these face normals + vectorField slaveFaceNormals + ( + patch().boundaryMesh().mesh().boundary()[shadowPatchID_].size(), + vector::zero + ); - // set boundary conditions - if(!master_) - { - // set refValue, refGrad and valueFraction + normalContactModelPtr_->correct + ( + *masterFaceZonePatchPtr_, + *slaveFaceZonePatchPtr_, + alg_, + dir_, + fieldName_, + orthotropic_, + nonLinearGeometry::nonLinearNames_[nonLinear_], + slaveFaceNormals, + slaveToMasterGgiInterpolatorPtr_ + ); - // refValue - refValue() = normalContactModelPtr_->slaveDisp() + frictionContactModelPtr_->slaveDisp(); + frictionContactModelPtr_->correct + ( + normalContactModelPtr_->slavePressure(), + *masterFaceZonePatchPtr_, + *slaveFaceZonePatchPtr_, + alg_, + dir_, + interpolationMethod_, + fieldName_, + orthotropic_, + nonLinearGeometry::nonLinearNames_[nonLinear_], + slaveFaceNormals + ); + } - //refGrad - set traction - refGrad() = - tractionBoundaryGradient() - ( - frictionContactModelPtr_->slaveTraction()+normalContactModelPtr_->slavePressure(), - scalarField(patch().size(),0.0), - word(fieldName_), - patch(), - orthotropic_, - NamedEnum::names[nonLinear_] - )(); + // set boundary conditions + if (!master_) + { + // set refValue, refGrad and valueFraction - //valueFraction - valueFraction() = normalContactModelPtr_->slaveValueFrac() + frictionContactModelPtr_->slaveValueFrac(); - } - else - { - // master is always traction condition - // interpolate traction from slave - // Dirichlet-Neumann uses lagged traction - // penalty uses same as for slave traction - - if(rigidMaster_) - { - // set to master to traction free if it is rigid - refGrad() = - tractionBoundaryGradient() - ( - vectorField(patch().size(),vector::zero), - scalarField(patch().size(),0.0), - word(fieldName_), - patch(), - orthotropic_, - NamedEnum::names[nonLinear_] - )(); - } - else - { - refGrad() = - tractionBoundaryGradient() - ( - interpolateSlaveToMaster - ( - -frictionContactModelPtr_->slaveTraction() - -normalContactModelPtr_->slavePressure() - ), - scalarField(patch().size(),0.0), - word(fieldName_), - patch(), - orthotropic_, - NamedEnum::names[nonLinear_] - )(); - } - } - } // if correction freqeuncy + // refValue + refValue() = normalContactModelPtr_->slaveDisp() + frictionContactModelPtr_->slaveDisp(); + + //refGrad - set traction + refGrad() = + tractionBoundaryGradient() + ( + frictionContactModelPtr_->slaveTraction()+normalContactModelPtr_->slavePressure(), + scalarField(patch().size(),0.0), + word(fieldName_), + patch(), + orthotropic_, + nonLinearGeometry::nonLinearNames_[nonLinear_] + )(); + + //valueFraction + valueFraction() = normalContactModelPtr_->slaveValueFrac() + frictionContactModelPtr_->slaveValueFrac(); + } + else + { + // master is always traction condition + // interpolate traction from slave + // Dirichlet-Neumann uses lagged traction + // penalty uses same as for slave traction + + if (rigidMaster_) + { + // set to master to traction free if it is rigid + refGrad() = + tractionBoundaryGradient() + ( + vectorField(patch().size(),vector::zero), + scalarField(patch().size(),0.0), + word(fieldName_), + patch(), + orthotropic_, + nonLinearGeometry::nonLinearNames_[nonLinear_] + )(); + } + else + { + refGrad() = + tractionBoundaryGradient() + ( + interpolateSlaveToMaster + ( + -frictionContactModelPtr_->slaveTraction() + -normalContactModelPtr_->slavePressure() + ), + scalarField(patch().size(),0.0), + word(fieldName_), + patch(), + orthotropic_, + nonLinearGeometry::nonLinearNames_[nonLinear_] + )(); + } + } + } // if correction freqeuncy } // if contactActive // if the contact is not active then the patches behave as traction free else { - // set refGrad to traction free for master and slave - refGrad() = - tractionBoundaryGradient() - ( - vectorField(patch().size(),vector::zero), - scalarField(patch().size(),0.0), - word(fieldName_), - patch(), - orthotropic_, - NamedEnum::names[nonLinear_] - )(); + // set refGrad to traction free for master and slave + refGrad() = + tractionBoundaryGradient() + ( + vectorField(patch().size(),vector::zero), + scalarField(patch().size(),0.0), + word(fieldName_), + patch(), + orthotropic_, + nonLinearGeometry::nonLinearNames_[nonLinear_] + )(); } directionMixedFvPatchVectorField::updateCoeffs(); @@ -724,11 +743,11 @@ tmp solidContactFvPatchVectorField::interpolateSlaveToMaster const vectorField slaveField ) { - if(!master_) + if (!master_) { FatalError << "only the master may call the function " - "solidContactFvPatchVectorField::interpolateSlaveToMaster" - << exit(FatalError); + "solidContactFvPatchVectorField::interpolateSlaveToMaster" + << exit(FatalError); } const fvMesh& mesh = patch().boundaryMesh().mesh(); @@ -738,7 +757,7 @@ tmp solidContactFvPatchVectorField::interpolateSlaveToMaster // for now, the slave must be the slave const label slaveStart = mesh.boundaryMesh()[shadowPatchID_].start(); - + // global slave field vectorField globalSlaveField(slaveFaceZonePatchPtr_->size(), vector::zero); @@ -746,51 +765,51 @@ tmp solidContactFvPatchVectorField::interpolateSlaveToMaster forAll(slaveField, i) { globalSlaveField[mesh.faceZones()[slaveFaceZoneID_].whichFace(slaveStart + i)] = - slaveField[i]; + slaveField[i]; } //- exchange parallel data reduce(globalSlaveField, sumOp()); // sum because each face is only on one proc - + // select interpolator - patchToPatch or GGI vectorField globalMasterInterpField(masterFaceZonePatchPtr_->size(), vector::zero); - if(slaveToMasterPatchToPatchInterpolatorPtr_) + if (slaveToMasterPatchToPatchInterpolatorPtr_) { globalMasterInterpField = - slaveToMasterPatchToPatchInterpolatorPtr_->faceInterpolate - ( - globalSlaveField - ); + slaveToMasterPatchToPatchInterpolatorPtr_->faceInterpolate + ( + globalSlaveField + ); } - else if(slaveToMasterGgiInterpolatorPtr_) + else if (slaveToMasterGgiInterpolatorPtr_) { globalMasterInterpField = - slaveToMasterGgiInterpolatorPtr_->slaveToMaster - ( - globalSlaveField - ); + slaveToMasterGgiInterpolatorPtr_->slaveToMaster + ( + globalSlaveField + ); } else { FatalErrorIn("solidContactFvPatchVectorField::interpolateSlaveToMaster()") - << "interpolationMethod is not patchToPatch or GGI!" - << abort(FatalError); + << "interpolationMethod is not patchToPatch or GGI!" + << abort(FatalError); } // now put global back into local const label masterPatchStart = mesh.boundaryMesh()[patch().index()].start(); - + tmp tmasterInterpField(new vectorField(masterFaceZonePatchPtr_->size(), vector::zero)); vectorField& masterInterpField = tmasterInterpField(); forAll(masterInterpField, i) { masterInterpField[i] = - globalMasterInterpField - [ - mesh.faceZones()[masterFaceZoneID_].whichFace(masterPatchStart + i) - ]; + globalMasterInterpField + [ + mesh.faceZones()[masterFaceZoneID_].whichFace(masterPatchStart + i) + ]; } return tmasterInterpField; @@ -808,11 +827,11 @@ void solidContactFvPatchVectorField::moveFaceZonePatches() // to the vertices. And we move the vertices by these // interpolated displacements, so the global face zone patches // should be in the same deformed position on all procs. - if(!master_) + if (!master_) { FatalError << "Only the master may call the function " - "solidContactFvPatchVectorField::moveFaceZonePatches" - << exit(FatalError); + "solidContactFvPatchVectorField::moveFaceZonePatches" + << exit(FatalError); } // update face zone patch interpolators @@ -833,19 +852,19 @@ void solidContactFvPatchVectorField::moveFaceZonePatches() this->db().objectRegistry::lookupObject(fieldName_); vectorField localMasterU = dispField.boundaryField()[patch().index()]; vectorField localSlaveU = dispField.boundaryField()[shadowPatchID_]; - if(fieldName_ == "DU" + if (fieldName_ == "DU" && - nonLinear_ != UPDATED_LAGRANGIAN) + nonLinear_ != nonLinearGeometry::UPDATED_LAGRANGIAN) { const volVectorField& totalDispField = - this->db().objectRegistry::lookupObject("U"); + this->db().objectRegistry::lookupObject("U"); localMasterU += totalDispField.boundaryField()[patch().index()]; localSlaveU += totalDispField.boundaryField()[shadowPatchID_]; } - else if(fieldName_ != "U" && fieldName_ != "DU") + else if (fieldName_ != "U" && fieldName_ != "DU") { FatalError << "Displacement field must be U or DU!" - << exit(FatalError); + << exit(FatalError); } // add on initial position so that localSlaveU becomes current local @@ -857,18 +876,18 @@ void solidContactFvPatchVectorField::moveFaceZonePatches() forAll(localMasterU, i) { globalMasterU[mesh.faceZones()[masterFaceZoneID_].whichFace(masterPatchStart + i)] = - localMasterU[i]; + localMasterU[i]; } // put localSlaveU field into globalSlaveU forAll(localSlaveU, i) { globalSlaveU[mesh.faceZones()[slaveFaceZoneID_].whichFace(slavePatchStart + i)] = - localSlaveU[i]; + localSlaveU[i]; } - // if(Pstream::myProcNo() == 2) // && (facei == 20) ) + // if (Pstream::myProcNo() == 2) // && (facei == 20) ) // { - // Pout << nl << "localSlaveU[20] " << localSlaveU[20] << endl; + // Pout << nl << "localSlaveU[20] " << localSlaveU[20] << endl; // } //- exchange parallel data @@ -913,25 +932,25 @@ void solidContactFvPatchVectorField::moveFaceZonePatches() // forAll(mesh.boundary(), patchi) // { - // if(mesh.boundary()[patchi].type() == "symmetryPlane" - // || - // mesh.boundary()[patchi].type() == "empty") - // { - // Info << "patch " << mesh.boundary()[patchi].name() << " is a symmetryPlane or empty" - // << endl; + // if (mesh.boundary()[patchi].type() == "symmetryPlane" + // || + // mesh.boundary()[patchi].type() == "empty") + // { + // Info << "patch " << mesh.boundary()[patchi].name() << " is a symmetryPlane or empty" + // << endl; - // const labelList& meshPoints = mesh.boundary()[patchi].meshPoints(); - - // forAll(slaveBoundaryPointsGlobalIndex, pointi) - // { - // label pointGlobalID = slaveBoundaryPointsGlobalIndex[pointi]; - // if(mesh.boundary()[patchi].whichPoint(pointGlobalID) > -1) - // { - // // remove normal component - // globalMasterNewPoints[pointGlobalID] - // } - // } - // } + // const labelList& meshPoints = mesh.boundary()[patchi].meshPoints(); + + // forAll(slaveBoundaryPointsGlobalIndex, pointi) + // { + // label pointGlobalID = slaveBoundaryPointsGlobalIndex[pointi]; + // if (mesh.boundary()[patchi].whichPoint(pointGlobalID) > -1) + // { + // // remove normal component + // globalMasterNewPoints[pointGlobalID] + // } + // } + // } // } @@ -953,7 +972,7 @@ void solidContactFvPatchVectorField::moveFaceZonePatches() // masterFaceZonePatchPtr_->writeVTK // ( - // fileName("masterFaceZonePatch"+name(Pstream::myProcNo())), + // fileName("masterFaceZonePatch"+name(Pstream::myProcNo())), // masterFaceZonePatchPtr_->localFaces(), // masterFaceZonePatchPtr_->localPoints() // ); @@ -971,59 +990,59 @@ void solidContactFvPatchVectorField::moveFaceZonePatches() // slaveFaceZonePatchPtr_->localFaces(), // slaveFaceZonePatchPtr_->localPoints() // ); - - // if(mesh.time().value() > 0.00125) + + // if (mesh.time().value() > 0.00125) // { // masterFaceZonePatchPtr_->writeVTK - // ( - // fileName("masterFaceZonePatch"+name(Pstream::myProcNo())), - // masterFaceZonePatchPtr_->localFaces(), - // masterFaceZonePatchPtr_->localPoints() - // ); + // ( + // fileName("masterFaceZonePatch"+name(Pstream::myProcNo())), + // masterFaceZonePatchPtr_->localFaces(), + // masterFaceZonePatchPtr_->localPoints() + // ); // slaveFaceZonePatchPtr_->writeVTK - // ( - // fileName("slaveFaceZonePatch"+name(Pstream::myProcNo())), - // slaveFaceZonePatchPtr_->localFaces(), - // slaveFaceZonePatchPtr_->localPoints() - // ); + // ( + // fileName("slaveFaceZonePatch"+name(Pstream::myProcNo())), + // slaveFaceZonePatchPtr_->localFaces(), + // slaveFaceZonePatchPtr_->localPoints() + // ); // scalar a = 3; // reduce(a, sumOp()); // FatalError << "end" << exit(FatalError); // } - - // The patchToPatch, GGI and primitivePatch interpolators point to the + + // The patchToPatch, GGI and primitivePatch interpolators point to the // old primitive patches so I must delete these interpolators // it would be much nicer if I could just move the faceZonePatch points... - if(slaveToMasterPatchToPatchInterpolatorPtr_) + if (slaveToMasterPatchToPatchInterpolatorPtr_) { delete slaveToMasterPatchToPatchInterpolatorPtr_; slaveToMasterPatchToPatchInterpolatorPtr_ = - new PatchToPatchInterpolation, PrimitivePatch > - ( - *slaveFaceZonePatchPtr_, // from zone - *masterFaceZonePatchPtr_, // to zone - alg_, - dir_ - ); + new PatchToPatchInterpolation, PrimitivePatch > + ( + *slaveFaceZonePatchPtr_, // from zone + *masterFaceZonePatchPtr_, // to zone + alg_, + dir_ + ); } - else if(slaveToMasterGgiInterpolatorPtr_) + else if (slaveToMasterGgiInterpolatorPtr_) { delete slaveToMasterGgiInterpolatorPtr_; slaveToMasterGgiInterpolatorPtr_ = - new GGIInterpolation< PrimitivePatch< face, List, pointField >, PrimitivePatch< face, List, pointField > > - ( - *masterFaceZonePatchPtr_, // master zone - *slaveFaceZonePatchPtr_, // slave zone - tensorField(0), - tensorField(0), - vectorField(0), - 0.0, - 0.0, - true, - ggiInterpolation::AABB - ); + new GGIInterpolation< PrimitivePatch< face, List, pointField >, PrimitivePatch< face, List, pointField > > + ( + *masterFaceZonePatchPtr_, // master zone + *slaveFaceZonePatchPtr_, // slave zone + tensorField(0), + tensorField(0), + vectorField(0), + 0.0, + 0.0, + true, + ggiInterpolation::AABB + ); } // and primitive patch interpolators @@ -1045,64 +1064,64 @@ bool solidContactFvPatchVectorField::checkPatchAndFaceZones(const dictionary& di // check shadow patch word shadowName = dict.lookup("shadowPatch"); label shadowPatchID = patch().patch().boundaryMesh().findPatchID(shadowName); - if(shadowPatchID == -1) + if (shadowPatchID == -1) { FatalError << "shadowPatch " << shadowName << " not found for patch " - << patch().name() << exit(FatalError); + << patch().name() << exit(FatalError); } word curZoneName = patch().name()+"FaceZone"; word shadowZoneName = patch().boundaryMesh().mesh().boundary()[shadowPatchID].name() + "FaceZone"; label curPatchFaceZoneID = patch().boundaryMesh().mesh().faceZones().findZoneID(curZoneName); - if(curPatchFaceZoneID == -1) + if (curPatchFaceZoneID == -1) { FatalError << "faceZone " << curZoneName - << " not found and is required for the solidContact boundaries\n" - << "To create a faceZone from a patch, use the setSet and setsToZone utilities:\n" - << "\tsetSet\n" - << "\tfaceSet "<()) < 1) + if (returnReduce(patch().boundaryMesh().mesh().faceZones()[curPatchFaceZoneID].size(), sumOp