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