applications/solvers/solidMechanics/elasticPlasticNonLinULSolidFoam/writeFields.H

   1 if (runTime.outputTime())
   2   {
   3     volScalarField epsilonEq
   4       (
   5        IOobject
   6        (
   7         "epsilonEq",
   8         runTime.timeName(),
   9         mesh,
  10         IOobject::NO_READ,
  11         IOobject::AUTO_WRITE
  12         ),
  13        sqrt((2.0/3.0)*magSqr(dev(epsilon)))
  14        );
  15
  16     Info<< "Max epsilonEq = " << max(epsilonEq).value()
  17         << endl;
  18
  19     volScalarField sigmaEq
  20       (
  21        IOobject
  22        (
  23         "sigmaEq",
  24         runTime.timeName(),
  25         mesh,
  26         IOobject::NO_READ,
  27         IOobject::AUTO_WRITE
  28         ),
  29        sqrt((3.0/2.0)*magSqr(dev(sigma)))
  30        );
  31
  32     Info<< "Max sigmaEq = " << max(sigmaEq).value()
  33         << endl;
  34
  35     pointMesh pMesh(mesh);
  36     pointScalarField contactPointGap
  37       (
  38        IOobject
  39        (
  40         "contactPointGap",
  41         runTime.timeName(),
  42         mesh,
  43         IOobject::NO_READ,
  44         IOobject::AUTO_WRITE
  45         ),
  46        pMesh,
  47        dimensionedScalar("zero", dimless, 0.0)
  48        );
  49
  50     forAll(mesh.boundary(), patchi)
  51       {
  52         if(DU.boundaryField()[patchi].type() == solidContactFvPatchVectorField::typeName)
  53           {
  54             const solidContactFvPatchVectorField& DUpatch =
  55               refCast<const solidContactFvPatchVectorField>
  56               (DU.boundaryField()[patchi]);
  57
  58                 if(!DUpatch.master())
  59                   {
  60                     const labelList& meshPoints = mesh.boundaryMesh()[patchi].meshPoints();
  61                     const scalarField gap = DUpatch.normalContactModelPtr()->slaveContactPointGap();
  62                     forAll(meshPoints, pointi)
  63                       {
  64                         contactPointGap[meshPoints[pointi]] = gap[pointi];
  65                       }
  66                   }
  67           }
  68       }
  69
  70
  71     //- total force
  72     /*    forAll(mesh.boundary(), patchi)
  73       {
  74         Info << "Patch " << mesh.boundary()[patchi].name() << endl;
  75         vectorField totalForce = mesh.Sf().boundaryField()[patchi] & sigma.boundaryField()[patchi];
  76
  77         vector force = sum( totalForce );
  78         Info << "\ttotal force is " << force << " N" << endl;
  79
  80         tensorField F = I + gradDU.boundaryField()[patchi];
  81         tensorField Finv = inv(F);
  82         //vectorField nCurrent = Finv & n.boundaryField()[patchi];
  83         //nCurrent /= mag(nCurrent);
  84         //scalar normalForce = sum( nCurrent & totalForce );
  85         scalar normalForce = sum( n.boundaryField()[patchi] & totalForce );
  86         Info << "\tnormal force is " << normalForce << " N" << endl;
  87         //scalar shearForce = mag(sum( (I - sqr(nCurrent)) & totalForce ));
  88         scalar shearForce = mag(sum( (I - sqr(n.boundaryField()[patchi])) & totalForce ));
  89         Info << "\tshear force is " << shearForce << " N" << endl;
  90
  91         // if(mesh.boundary()[patchi].type() != "empty")
  92         //   {
  93         //     vector Sf0 = Sf.boundaryField()[patchi][0];
  94         //     symmTensor sigma0 = sigma.boundaryField()[patchi][0];
  95         //     Info << "sigmab[0] is " << sigma0 << nl
  96         //       << "Sfb is " << Sf0 << nl
  97         //       << "force is " << (Sf.boundaryField()[patchi][0]&sigma.boundaryField()[patchi][0]) << nl
  98         //       << "Sfx*sigmaxx " << (Sf0[vector::X]*sigma0[symmTensor::XX]) << nl
  99         //       << "Sfy*sigmaxy " << (Sf0[vector::Y]*sigma0[symmTensor::XY]) << nl
 100         //       << "Sfx*sigmayx " << (Sf0[vector::X]*sigma0[symmTensor::XY]) << nl
 101         //       << "Sfy*sigmayy " << (Sf0[vector::Y]*sigma0[symmTensor::YY]) << nl
 102         //       << endl;
 103             //vector SfTL(-0.000137451, 0.00383599, -4.76878e-20);
 104             // vector SfTL = Finv[0] & vector(0,0.004,0);
 105             // Info << "SfTLx*sigmaxx " << (SfTL[vector::X]*sigma0[symmTensor::XX]) << nl
 106             //   << "SfTLy*sigmaxy " << (SfTL[vector::Y]*sigma0[symmTensor::XY]) << nl
 107             //   << "SfTLx*sigmayx " << (SfTL[vector::X]*sigma0[symmTensor::XY]) << nl
 108             //   << "SfTLy*sigmayy " << (SfTL[vector::Y]*sigma0[symmTensor::YY]) << nl
 109             //   << endl;
 110         //        }
 111         }*/
 112
 113     runTime.write();
 114   }