applications/solvers/solidMechanics/elasticIncrSolidFoam/writeForceDisplacement.H

   1 //- write force displacement to file
   2
   3 label leftPatchID = mesh.boundaryMesh().findPatchID("leftClamp");
   4 if(leftPatchID == -1)
   5   {
   6     FatalError << "Cannot find patch left for calculating force" << endl;
   7   }
   8
   9 //- calculate force in x direction on leftClamp patch
  10 scalar leftForce = gSum
  11 (
  12     vector(1, 0, 0) &
  13     (mesh.boundary()[leftPatchID].Sf() & sigma.boundaryField()[leftPatchID])
  14 );
  15
  16 //- patchIntegrate utility integrates it this way but this is worng because the sigma tensor should
  17 //- be dotted with the surface normal to give the actual traction/force
  18 //- you cannot just take the component of the sigma tensor
  19 //scalar leftForcePatchIntegrateMethod = gSum
  20 //(
  21 //    mesh.magSf().boundaryField()[leftPatchID]*
  22 //    sigma.boundaryField()[leftPatchID].component(symmTensor::XY)
  23 //);
  24
  25 vector gaugeU1 = vector::zero;
  26 vector gaugeU2 = vector::zero;
  27 if(gaugeFaceID1 != -1)
  28   {
  29     gaugeU1 = U.boundaryField()[gaugeFacePatchID1][gaugeFaceID1];
  30   }
  31 if(gaugeFaceID2 != -1)
  32   {
  33     gaugeU2 = U.boundaryField()[gaugeFacePatchID2][gaugeFaceID2];
  34   }
  35
  36 //- reduce across procs
  37 reduce(gaugeU1, sumOp<vector>());
  38 reduce(gaugeU2, sumOp<vector>());
  39
  40 Pout << "gaugeU1 is " << gaugeU1 << nl
  41      << "gaugeU2 is " << gaugeU2 << endl;
  42
  43 scalar gaugeDisp = mag(gaugeU1 - gaugeU2);
  44
  45 //- write to file
  46 if(Pstream::master())
  47   {
  48     OFstream& forceDispFile = *filePtr;
  49     forceDispFile << 1000*gaugeDisp << "\t" << -1*leftForce << endl;
  50   }