Unverified Commit 08b135ce authored by ckadding's avatar ckadding Committed by GitHub
Browse files

Apply uniform LAMMPS formatting

parent a6ba5508

src/USER-MISC/compute_pressure_cylinder.cpp

+38 −61
Original line number Diff line number Diff line
@@ -152,16 +152,14 @@ void ComputePressureCyl::init() 


  double phi;



  for (int iphi=0;iphi<nphi;iphi++)

  {

  for (int iphi = 0; iphi < nphi; iphi++) {

    phi=((double)iphi)*MY_PI/180.0;

    tangent[iphi]=tan(phi);

    ephi_x[iphi]=-sin(phi);

    ephi_y[iphi]=cos(phi);

  }



  for (int iq=0;iq<nbins;iq++)

  {

  for (int iq = 0; iq < nbins; iq++) {

    R[iq]=((double)iq+0.5)*bin_width;

    Rinv[iq]=1.0/R[iq];

    R2[iq]=R[iq]*R[iq];
@@ -173,8 +171,8 @@ void ComputePressureCyl::init() 


  invVbin[0]=1.0/((zhi-zlo)*MY_PI*R2kin[0]);

  PzAinv[0]=1.0/(MY_PI*R2kin[0]*((double)nzbins));

  for (int jq=1;jq<nbins;jq++)

  {

   

  for (int jq = 1; jq < nbins; jq++) {

    invVbin[jq]=1.0/((zhi-zlo)*MY_PI*(R2kin[jq]-R2kin[jq-1]));

    PzAinv[jq]=1.0/(MY_PI*(R2kin[jq]-R2kin[jq-1])*((double)nzbins));

  }
@@ -213,8 +211,7 @@ void ComputePressureCyl::compute_array() 
  int ibin;



  // clear pressures

  for (ibin=0;ibin<nbins;ibin++)

  {

  for (ibin = 0; ibin < nbins; ibin++) {

    density_temp[ibin]=0.0;

    density_all[ibin]=0.0;

    Pr_temp[ibin]=0.0;
@@ -254,8 +251,7 @@ void ComputePressureCyl::compute_array() 


  // calculate number density (by radius)

  double temp_R2;

  for (i=0;i<nlocal;i++) if (x[i][2]<zhi && x[i][2]>zlo)

  {

  for (i = 0; i < nlocal; i++) if ((x[i][2] < zhi) && (x[i][2] > zlo)) {

    temp_R2=x[i][0]*x[i][0]+x[i][1]*x[i][1];

    if (temp_R2 > R2kin[nbins-1]) continue; // outside of Rmax
@@ -289,8 +285,7 @@ void ComputePressureCyl::compute_array() 
  double sqrtD;

  double lower_z,upper_z;



  for (ii = 0; ii < inum; ii++)

  {

  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    if (!(mask[i] & groupbit)) continue;
@@ -304,8 +299,7 @@ void ComputePressureCyl::compute_array() 


    r1=x[i][0]*x[i][0]+x[i][1]*x[i][1];



    for (jj = 0; jj < jnum; jj++)

    {

    for (jj = 0; jj < jnum; jj++) {

      j = jlist[jj];

      factor_lj = special_lj[sbmask(j)];

      factor_coul = special_coul[sbmask(j)];
@@ -315,22 +309,17 @@ void ComputePressureCyl::compute_array() 


      // itag = jtag is possible for long cutoffs that include images of self

      // do calculation only on appropriate processor

      if (newton_pair == 0 && j >= nlocal)

      {

      if (newton_pair == 0 && j >= nlocal) {

        jtag = tag[j];

        if (itag > jtag)

        {

        if (itag > jtag) {

          if ((itag+jtag) % 2 == 0) continue;

        }

        else if (itag < jtag)

        {

        else if (itag < jtag) {

          if ((itag+jtag) % 2 == 1) continue;

        }

        else

        {

        else {

          if (x[j][2] < ztmp) continue;

          if (x[j][2] == ztmp)

          {

          if (x[j][2] == ztmp) {

            if (x[j][1] < ytmp) continue;

            if (x[j][1] == ytmp && x[j][0] < xtmp) continue;

          }
@@ -344,8 +333,7 @@ void ComputePressureCyl::compute_array() 
      r2=x[j][0]*x[j][0]+x[j][1]*x[j][1];



      // ri is smaller of r1 and r2

      if (r2<r1)

      {

      if (r2 < r1) {

        risq=r2;

        rjsq=r1;

        xi=x[j][0];
@@ -355,8 +343,7 @@ void ComputePressureCyl::compute_array() 
        dy=x[i][1]-x[j][1];

        dz=x[i][2]-x[j][2];

      }

      else

      {

      else {

        risq=r1;

        rjsq=r2;

        xi=x[i][0];
@@ -377,8 +364,7 @@ void ComputePressureCyl::compute_array() 
      B=2.0*(xi*dx+yi*dy);



      // normal pressure contribution P_rhorho

      for (ibin=0;ibin<nbins;ibin++)

      {

      for (ibin = 0; ibin < nbins; ibin++) {

        // completely inside of R

        if (rjsq < R2[ibin]) continue;
@@ -392,10 +378,9 @@ void ComputePressureCyl::compute_array() 
        alpha1=0.5*(-B+sqrtD)/A;

        alpha2=0.5*(-B-sqrtD)/A;



        if (alpha1>0.0 && alpha1<1.0)

        {

        if ((alpha1 > 0.0) && (alpha1 < 1.0)) {

          zR=zi+alpha1*dz;

          if (zR<zhi && zR>zlo)

          if ((zR < zhi) && (zR > zlo))

          {

            xR=xi+alpha1*dx;

            yR=yi+alpha1*dy;
@@ -404,11 +389,9 @@ void ComputePressureCyl::compute_array() 
            Pr_temp[ibin]+=fn;

          }

        }

        if (alpha2>0.0 && alpha2<1.0)

        {

        if ((alpha2 > 0.0) && (alpha2 < 1.0)) {

          zR=zi+alpha2*dz;

          if (zR<zhi && zR>zlo)

          {

          if ((zR < zhi) && (zR > zlo)) {

            xR=xi+alpha2*dx;

            yR=yi+alpha2*dy;

            fn=fpair*fabs(xR*dx+yR*dy);
@@ -419,16 +402,15 @@ void ComputePressureCyl::compute_array() 
      }



      // azimuthal pressure contribution (P_phiphi)

      for (int iphi=0;iphi<nphi;iphi++)

      {

      for (int iphi = 0; iphi < nphi; iphi++) {

        alpha=(yi-xi*tangent[iphi])/(dx*tangent[iphi]-dy);



        // no intersection with phi surface

        if (alpha>=1.0 || alpha<=0.0) continue;

        if ((alpha >= 1.0) || (alpha <= 0.0)) continue;



        // no contribution (outside of averaging region)

        zL=zi+alpha*dz;

        if (zL>zhi || zL<zlo) continue;

        if ((zL > zhi) || (zL < zlo)) continue;



        xL=xi+alpha*dx;

        yL=yi+alpha*dy;
@@ -441,19 +423,17 @@ void ComputePressureCyl::compute_array() 
        ftphi=fabs(dx*ephi_x[iphi]+dy*ephi_y[iphi])*fpair;



        // add to appropriate bin

        for (ibin=0;ibin<nbins;ibin++) if (L2<R2kin[ibin])

        {

        for (ibin = 0; ibin < nbins; ibin++) if (L2 < R2kin[ibin]) {

          Pphi_temp[ibin]+=ftphi;

          break;

        }

      }



      // z pressure contribution (P_zz)

      for (int zbin=0;zbin<nzbins;zbin++)

      {

      for (int zbin = 0; zbin < nzbins; zbin++) {

        // check if interaction contributes

        if (x[i][2]>binz[zbin] && x[j][2]>binz[zbin]) continue;

        if (x[i][2]<binz[zbin] && x[j][2]<binz[zbin]) continue;

        if ((x[i][2] > binz[zbin]) && (x[j][2] > binz[zbin])) continue;

        if ((x[i][2] < binz[zbin]) && (x[j][2] < binz[zbin])) continue;



        alpha=(binz[zbin]-zi)/dz;
@@ -467,8 +447,7 @@ void ComputePressureCyl::compute_array() 
        ftz=fabs(dz)*fpair;



        // add to appropriate bin

        for (ibin=0;ibin<nbins;ibin++) if (L2<R2kin[ibin])

        {

        for (ibin = 0; ibin < nbins; ibin++) if (L2 < R2kin[ibin]) {

          Pz_temp[ibin]+=ftz;

          break;

        }
@@ -477,8 +456,7 @@ void ComputePressureCyl::compute_array() 
  }



  // calculate pressure (force over area)

  for (ibin=0;ibin<nbins;ibin++)

  {

  for (ibin = 0; ibin < nbins; ibin++) {

    Pr_temp[ibin]*=PrAinv[ibin]*Rinv[ibin];

    Pphi_temp[ibin]*=PphiAinv;

    Pz_temp[ibin]*=PzAinv[ibin];
@@ -490,8 +468,7 @@ void ComputePressureCyl::compute_array() 
  MPI_Allreduce(Pz_temp,Pz_all,nbins,MPI_DOUBLE,MPI_SUM,world);



  // populate array

  for (ibin=0;ibin<nbins;ibin++)

  {

  for (ibin = 0; ibin < nbins; ibin++) {

    array[ibin][0]=R[ibin];

    array[ibin][2]=Pr_all[ibin]*nktv2p;

    array[ibin][3]=Pphi_all[ibin]*nktv2p;

CRA Git | Maintained and supported by SUSTech CRA and CCSE