---

stresses.old.c

Go to the documentation of this file.

#include "utilf.h"
#include "stresses.h"

   
void stresses(real2D S11, real2D S22, real2D S12, real2D cc, 
              real2D S11t, real2D S22t, 
#ifdef IMMERSED
              real ***isb_sxy,
#endif
              real sigma, real tau, 
              int nx, int ny)
{
  real sigmag25, h, sigmag, sigh, mag, dcx, dcy;
  int i, j;
  
  h = 1.0 / (nx - 2);
  sigh =.5 * tau * tau * sigma / h / h / h;
  for (i = 2; i <= nx; i++)
    for (j = 2; j <= ny; j++)
      {
#ifdef IMMERSED
        if (isb_sxy[i][j] != NULL) {
#ifdef CA90
          /* 90 degrees contact angle */
          dcx= 2. * (cc[i][j] - cc[i-1][j]);
          dcy= 0.0;
#else
          /* 180 degrees contact angle */
          dcx = 0.0;
          dcy = cc[i][j] + cc[i-1][j];
#endif
        }
        else {
          dcx= cc[i][j] + cc[i][j-1] - cc[i-1][j] - cc[i-1][j-1];
          dcy= cc[i][j] + cc[i-1][j] - cc[i][j-1] - cc[i-1][j-1];
        }
#else
        dcx= cc[i][j] + cc[i][j-1] - cc[i-1][j] - cc[i-1][j-1];
        dcy= cc[i][j] + cc[i-1][j] - cc[i][j-1] - cc[i-1][j-1];
#endif
        mag = dcx * dcx + dcy * dcy;
        if (mag > 0.0) {
          sigmag = sigh / sqrt(mag);
          sigmag25 = 0.25 * sigmag;
          S11t[i][j] = - sigmag25 * dcy * dcy;
          S22t[i][j] = - sigmag25 * dcx * dcx;
          S12[i][j] =  sigmag * dcx * dcy;
        }
        else
          S11t[i][j] = S22t[i][j] = S12[i][j] = 0.0;
      }
  for (i = 2; i <= nx-1; i++)
    for (j = 2; j <= ny-1; j++)
      {
        S11[i][j] = S11t[i][j] + S11t[i+1][j] 
          + S11t[i][j+1] + S11t[i+1][j+1];
        S22[i][j] = S22t[i][j] + S22t[i+1][j] + S22t[i][j+1] 
          + S22t[i+1][j+1];
      }
}

---