| 754 | 
  | 
    } | 
| 755 | 
  | 
     | 
| 756 | 
  | 
#ifdef IS_MPI     | 
| 757 | 
< | 
    int worldRank = MPI::COMM_WORLD.Get_rank(); | 
| 758 | 
< | 
     | 
| 759 | 
< | 
    bool my_min_found = min_found; | 
| 760 | 
< | 
    bool my_max_found = max_found; | 
| 757 | 
> | 
    int worldRank; | 
| 758 | 
> | 
    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank); | 
| 759 | 
> | 
          | 
| 760 | 
> | 
    int my_min_found = min_found; | 
| 761 | 
> | 
    int my_max_found = max_found; | 
| 762 | 
  | 
 | 
| 763 | 
  | 
    // Even if we didn't find a minimum, did someone else? | 
| 764 | 
< | 
    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR); | 
| 764 | 
> | 
    MPI_Allreduce(&my_min_found, &min_found, 1, MPI_INT, MPI_LOR,  | 
| 765 | 
> | 
                  MPI_COMM_WORLD); | 
| 766 | 
  | 
    // Even if we didn't find a maximum, did someone else? | 
| 767 | 
< | 
    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR); | 
| 767 | 
> | 
    MPI_Allreduce(&my_max_found, &max_found, 1, MPI_INT, MPI_LOR,  | 
| 768 | 
> | 
                  MPI_COMM_WORLD); | 
| 769 | 
  | 
#endif | 
| 770 | 
  | 
 | 
| 771 | 
  | 
    if (max_found && min_found) { | 
| 784 | 
  | 
      min_vals.rank = worldRank;     | 
| 785 | 
  | 
       | 
| 786 | 
  | 
      // Who had the minimum? | 
| 787 | 
< | 
      MPI::COMM_WORLD.Allreduce(&min_vals, &min_vals,  | 
| 788 | 
< | 
                                1, MPI::REALTYPE_INT, MPI::MINLOC); | 
| 787 | 
> | 
      MPI_Allreduce(&min_vals, &min_vals,  | 
| 788 | 
> | 
                    1, MPI_REALTYPE_INT, MPI_MINLOC, MPI_COMM_WORLD); | 
| 789 | 
  | 
      min_val = min_vals.val; | 
| 790 | 
  | 
       | 
| 791 | 
  | 
      if (my_max_found) { | 
| 796 | 
  | 
      max_vals.rank = worldRank;     | 
| 797 | 
  | 
       | 
| 798 | 
  | 
      // Who had the maximum? | 
| 799 | 
< | 
      MPI::COMM_WORLD.Allreduce(&max_vals, &max_vals,  | 
| 800 | 
< | 
                                1, MPI::REALTYPE_INT, MPI::MAXLOC); | 
| 799 | 
> | 
      MPI_Allreduce(&max_vals, &max_vals,  | 
| 800 | 
> | 
                    1, MPI_REALTYPE_INT, MPI_MAXLOC, MPI_COMM_WORLD); | 
| 801 | 
  | 
      max_val = max_vals.val; | 
| 802 | 
  | 
#endif | 
| 803 | 
  | 
       | 
| 857 | 
  | 
           | 
| 858 | 
  | 
          Vector3d min_vel; | 
| 859 | 
  | 
          Vector3d max_vel = max_sd->getVel(); | 
| 860 | 
< | 
          MPI::Status status; | 
| 860 | 
> | 
          MPI_Status* status; | 
| 861 | 
  | 
 | 
| 862 | 
  | 
          // point-to-point swap of the velocity vector | 
| 863 | 
< | 
          MPI::COMM_WORLD.Sendrecv(max_vel.getArrayPointer(), 3, MPI::REALTYPE, | 
| 864 | 
< | 
                                   min_vals.rank, 0,  | 
| 865 | 
< | 
                                   min_vel.getArrayPointer(), 3, MPI::REALTYPE, | 
| 866 | 
< | 
                                   min_vals.rank, 0, status); | 
| 863 | 
> | 
          MPI_Sendrecv(max_vel.getArrayPointer(), 3, MPI_REALTYPE, | 
| 864 | 
> | 
                       min_vals.rank, 0,  | 
| 865 | 
> | 
                       min_vel.getArrayPointer(), 3, MPI_REALTYPE, | 
| 866 | 
> | 
                       min_vals.rank, 0, MPI_COMM_WORLD, status); | 
| 867 | 
  | 
           | 
| 868 | 
  | 
          switch(rnemdFluxType_) { | 
| 869 | 
  | 
          case rnemdKE : | 
| 874 | 
  | 
              Vector3d max_angMom = max_sd->getJ(); | 
| 875 | 
  | 
               | 
| 876 | 
  | 
              // point-to-point swap of the angular momentum vector | 
| 877 | 
< | 
              MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3,  | 
| 878 | 
< | 
                                       MPI::REALTYPE, min_vals.rank, 1,  | 
| 879 | 
< | 
                                       min_angMom.getArrayPointer(), 3,  | 
| 880 | 
< | 
                                       MPI::REALTYPE, min_vals.rank, 1,  | 
| 881 | 
< | 
                                       status); | 
| 877 | 
> | 
              MPI_Sendrecv(max_angMom.getArrayPointer(), 3,  | 
| 878 | 
> | 
                           MPI_REALTYPE, min_vals.rank, 1,  | 
| 879 | 
> | 
                           min_angMom.getArrayPointer(), 3,  | 
| 880 | 
> | 
                           MPI_REALTYPE, min_vals.rank, 1,  | 
| 881 | 
> | 
                           MPI_COMM_WORLD, status); | 
| 882 | 
  | 
               | 
| 883 | 
  | 
              max_sd->setJ(min_angMom); | 
| 884 | 
  | 
            } | 
| 903 | 
  | 
           | 
| 904 | 
  | 
          Vector3d max_vel; | 
| 905 | 
  | 
          Vector3d min_vel = min_sd->getVel(); | 
| 906 | 
< | 
          MPI::Status status; | 
| 906 | 
> | 
          MPI_Status* status; | 
| 907 | 
  | 
           | 
| 908 | 
  | 
          // point-to-point swap of the velocity vector | 
| 909 | 
< | 
          MPI::COMM_WORLD.Sendrecv(min_vel.getArrayPointer(), 3, MPI::REALTYPE, | 
| 910 | 
< | 
                                   max_vals.rank, 0,  | 
| 911 | 
< | 
                                   max_vel.getArrayPointer(), 3, MPI::REALTYPE, | 
| 912 | 
< | 
                                   max_vals.rank, 0, status); | 
| 909 | 
> | 
          MPI_Sendrecv(min_vel.getArrayPointer(), 3, MPI_REALTYPE, | 
| 910 | 
> | 
                       max_vals.rank, 0,  | 
| 911 | 
> | 
                       max_vel.getArrayPointer(), 3, MPI_REALTYPE, | 
| 912 | 
> | 
                       max_vals.rank, 0, MPI_COMM_WORLD, status); | 
| 913 | 
  | 
           | 
| 914 | 
  | 
          switch(rnemdFluxType_) { | 
| 915 | 
  | 
          case rnemdKE : | 
| 920 | 
  | 
              Vector3d max_angMom; | 
| 921 | 
  | 
               | 
| 922 | 
  | 
              // point-to-point swap of the angular momentum vector | 
| 923 | 
< | 
              MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3,  | 
| 924 | 
< | 
                                       MPI::REALTYPE, max_vals.rank, 1,  | 
| 925 | 
< | 
                                       max_angMom.getArrayPointer(), 3,  | 
| 926 | 
< | 
                                       MPI::REALTYPE, max_vals.rank, 1,  | 
| 927 | 
< | 
                                       status); | 
| 923 | 
> | 
              MPI_Sendrecv(min_angMom.getArrayPointer(), 3,  | 
| 924 | 
> | 
                           MPI_REALTYPE, max_vals.rank, 1,  | 
| 925 | 
> | 
                           max_angMom.getArrayPointer(), 3,  | 
| 926 | 
> | 
                           MPI_REALTYPE, max_vals.rank, 1,  | 
| 927 | 
> | 
                           MPI_COMM_WORLD, status); | 
| 928 | 
  | 
               | 
| 929 | 
  | 
              min_sd->setJ(max_angMom); | 
| 930 | 
  | 
            } | 
| 1093 | 
  | 
    Kcw *= 0.5; | 
| 1094 | 
  | 
 | 
| 1095 | 
  | 
#ifdef IS_MPI | 
| 1096 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM); | 
| 1097 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phy, 1, MPI::REALTYPE, MPI::SUM); | 
| 1098 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phz, 1, MPI::REALTYPE, MPI::SUM); | 
| 1099 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcx, 1, MPI::REALTYPE, MPI::SUM); | 
| 1100 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcy, 1, MPI::REALTYPE, MPI::SUM); | 
| 1101 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcz, 1, MPI::REALTYPE, MPI::SUM); | 
| 1096 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Phx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1097 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Phy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1098 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Phz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1099 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Pcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1100 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Pcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1101 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Pcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1102 | 
  | 
 | 
| 1103 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khx, 1, MPI::REALTYPE, MPI::SUM); | 
| 1104 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khy, 1, MPI::REALTYPE, MPI::SUM); | 
| 1105 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khz, 1, MPI::REALTYPE, MPI::SUM); | 
| 1106 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khw, 1, MPI::REALTYPE, MPI::SUM); | 
| 1103 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Khx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1104 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Khy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1105 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Khz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1106 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Khw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1107 | 
  | 
 | 
| 1108 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcx, 1, MPI::REALTYPE, MPI::SUM); | 
| 1109 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcy, 1, MPI::REALTYPE, MPI::SUM); | 
| 1110 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcz, 1, MPI::REALTYPE, MPI::SUM); | 
| 1111 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcw, 1, MPI::REALTYPE, MPI::SUM); | 
| 1108 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Kcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1109 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Kcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1110 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Kcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1111 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Kcw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1112 | 
  | 
#endif | 
| 1113 | 
  | 
 | 
| 1114 | 
  | 
    //solve coldBin coeff's first | 
| 1585 | 
  | 
    Kc *= 0.5; | 
| 1586 | 
  | 
     | 
| 1587 | 
  | 
#ifdef IS_MPI | 
| 1588 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM); | 
| 1589 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM); | 
| 1590 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lh[0], 3, MPI::REALTYPE, MPI::SUM); | 
| 1591 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lc[0], 3, MPI::REALTYPE, MPI::SUM); | 
| 1592 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mh, 1, MPI::REALTYPE, MPI::SUM); | 
| 1593 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kh, 1, MPI::REALTYPE, MPI::SUM); | 
| 1594 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mc, 1, MPI::REALTYPE, MPI::SUM); | 
| 1595 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kc, 1, MPI::REALTYPE, MPI::SUM); | 
| 1596 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ih.getArrayPointer(), 9,  | 
| 1597 | 
< | 
                              MPI::REALTYPE, MPI::SUM); | 
| 1598 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ic.getArrayPointer(), 9,  | 
| 1599 | 
< | 
                              MPI::REALTYPE, MPI::SUM); | 
| 1588 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Ph[0], 3, MPI_REALTYPE, MPI_SUM,  | 
| 1589 | 
> | 
                  MPI_COMM_WORLD); | 
| 1590 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Pc[0], 3, MPI_REALTYPE, MPI_SUM,  | 
| 1591 | 
> | 
                  MPI_COMM_WORLD); | 
| 1592 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Lh[0], 3, MPI_REALTYPE, MPI_SUM,  | 
| 1593 | 
> | 
                  MPI_COMM_WORLD); | 
| 1594 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Lc[0], 3, MPI_REALTYPE, MPI_SUM,  | 
| 1595 | 
> | 
                  MPI_COMM_WORLD); | 
| 1596 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Mh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1597 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Kh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1598 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Mc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1599 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, &Kc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1600 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, Ih.getArrayPointer(), 9,  | 
| 1601 | 
> | 
                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1602 | 
> | 
    MPI_Allreduce(MPI_IN_PLACE, Ic.getArrayPointer(), 9,  | 
| 1603 | 
> | 
                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 1604 | 
  | 
#endif | 
| 1605 | 
  | 
     | 
| 1606 | 
  | 
 | 
| 1881 | 
  | 
    RealType area = getDividingArea(); | 
| 1882 | 
  | 
    areaAccumulator_->add(area); | 
| 1883 | 
  | 
    Mat3x3d hmat = currentSnap_->getHmat(); | 
| 1884 | 
+ | 
    Vector3d u = angularMomentumFluxVector_; | 
| 1885 | 
+ | 
    u.normalize(); | 
| 1886 | 
+ | 
 | 
| 1887 | 
  | 
    seleMan_.setSelectionSet(evaluator_.evaluate()); | 
| 1888 | 
  | 
 | 
| 1889 | 
  | 
    int selei(0); | 
| 1890 | 
  | 
    StuntDouble* sd; | 
| 1891 | 
  | 
    int binNo; | 
| 1892 | 
+ | 
    RealType mass; | 
| 1893 | 
+ | 
    Vector3d vel;  | 
| 1894 | 
+ | 
    Vector3d rPos; | 
| 1895 | 
+ | 
    RealType KE; | 
| 1896 | 
+ | 
    Vector3d L; | 
| 1897 | 
+ | 
    Mat3x3d I; | 
| 1898 | 
+ | 
    RealType r2; | 
| 1899 | 
  | 
 | 
| 1900 | 
  | 
    vector<RealType> binMass(nBins_, 0.0); | 
| 1901 | 
< | 
    vector<RealType> binPx(nBins_, 0.0); | 
| 1902 | 
< | 
    vector<RealType> binPy(nBins_, 0.0); | 
| 1903 | 
< | 
    vector<RealType> binPz(nBins_, 0.0); | 
| 1904 | 
< | 
    vector<RealType> binOmegax(nBins_, 0.0); | 
| 1888 | 
< | 
    vector<RealType> binOmegay(nBins_, 0.0); | 
| 1889 | 
< | 
    vector<RealType> binOmegaz(nBins_, 0.0); | 
| 1901 | 
> | 
    vector<Vector3d> binP(nBins_, V3Zero); | 
| 1902 | 
> | 
    vector<RealType> binOmega(nBins_, 0.0); | 
| 1903 | 
> | 
    vector<Vector3d> binL(nBins_, V3Zero); | 
| 1904 | 
> | 
    vector<Mat3x3d>  binI(nBins_); | 
| 1905 | 
  | 
    vector<RealType> binKE(nBins_, 0.0); | 
| 1906 | 
  | 
    vector<int> binDOF(nBins_, 0); | 
| 1907 | 
  | 
    vector<int> binCount(nBins_, 0); | 
| 1941 | 
  | 
        binNo = int(rPos.length() / binWidth_); | 
| 1942 | 
  | 
      } | 
| 1943 | 
  | 
 | 
| 1944 | 
< | 
      RealType mass = sd->getMass(); | 
| 1945 | 
< | 
      Vector3d vel = sd->getVel(); | 
| 1946 | 
< | 
      Vector3d rPos = sd->getPos() - coordinateOrigin_; | 
| 1947 | 
< | 
      Vector3d aVel = cross(rPos, vel); | 
| 1948 | 
< | 
       | 
| 1944 | 
> | 
      mass = sd->getMass(); | 
| 1945 | 
> | 
      vel = sd->getVel(); | 
| 1946 | 
> | 
      rPos = sd->getPos() - coordinateOrigin_; | 
| 1947 | 
> | 
      KE = 0.5 * mass * vel.lengthSquare(); | 
| 1948 | 
> | 
      L = mass * cross(rPos, vel); | 
| 1949 | 
> | 
      I = outProduct(rPos, rPos) * mass; | 
| 1950 | 
> | 
      r2 = rPos.lengthSquare(); | 
| 1951 | 
> | 
      I(0, 0) += mass * r2; | 
| 1952 | 
> | 
      I(1, 1) += mass * r2; | 
| 1953 | 
> | 
      I(2, 2) += mass * r2; | 
| 1954 | 
> | 
 | 
| 1955 | 
> | 
      // Project the relative position onto a plane perpendicular to | 
| 1956 | 
> | 
      // the angularMomentumFluxVector: | 
| 1957 | 
> | 
      // Vector3d rProj = rPos - dot(rPos, u) * u; | 
| 1958 | 
> | 
      // Project the velocity onto a plane perpendicular to the | 
| 1959 | 
> | 
      // angularMomentumFluxVector: | 
| 1960 | 
> | 
      // Vector3d vProj = vel  - dot(vel, u) * u; | 
| 1961 | 
> | 
      // Compute angular velocity vector (should be nearly parallel to | 
| 1962 | 
> | 
      // angularMomentumFluxVector | 
| 1963 | 
> | 
      // Vector3d aVel = cross(rProj, vProj); | 
| 1964 | 
> | 
 | 
| 1965 | 
  | 
      if (binNo >= 0 && binNo < nBins_)  { | 
| 1966 | 
  | 
        binCount[binNo]++; | 
| 1967 | 
  | 
        binMass[binNo] += mass; | 
| 1968 | 
< | 
        binPx[binNo] += mass*vel.x(); | 
| 1969 | 
< | 
        binPy[binNo] += mass*vel.y(); | 
| 1970 | 
< | 
        binPz[binNo] += mass*vel.z(); | 
| 1971 | 
< | 
        binOmegax[binNo] += aVel.x(); | 
| 1941 | 
< | 
        binOmegay[binNo] += aVel.y(); | 
| 1942 | 
< | 
        binOmegaz[binNo] += aVel.z(); | 
| 1943 | 
< | 
        binKE[binNo] += 0.5 * (mass * vel.lengthSquare()); | 
| 1968 | 
> | 
        binP[binNo] += mass*vel; | 
| 1969 | 
> | 
        binKE[binNo] += KE; | 
| 1970 | 
> | 
        binI[binNo] += I; | 
| 1971 | 
> | 
        binL[binNo] += L; | 
| 1972 | 
  | 
        binDOF[binNo] += 3; | 
| 1973 | 
  | 
         | 
| 1974 | 
  | 
        if (sd->isDirectional()) { | 
| 1975 | 
  | 
          Vector3d angMom = sd->getJ(); | 
| 1976 | 
< | 
          Mat3x3d I = sd->getI(); | 
| 1976 | 
> | 
          Mat3x3d Ia = sd->getI(); | 
| 1977 | 
  | 
          if (sd->isLinear()) { | 
| 1978 | 
  | 
            int i = sd->linearAxis(); | 
| 1979 | 
  | 
            int j = (i + 1) % 3; | 
| 1980 | 
  | 
            int k = (i + 2) % 3; | 
| 1981 | 
< | 
            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +  | 
| 1982 | 
< | 
                                   angMom[k] * angMom[k] / I(k, k)); | 
| 1981 | 
> | 
            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +  | 
| 1982 | 
> | 
                                   angMom[k] * angMom[k] / Ia(k, k)); | 
| 1983 | 
  | 
            binDOF[binNo] += 2; | 
| 1984 | 
  | 
          } else { | 
| 1985 | 
< | 
            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) + | 
| 1986 | 
< | 
                                   angMom[1] * angMom[1] / I(1, 1) + | 
| 1987 | 
< | 
                                   angMom[2] * angMom[2] / I(2, 2)); | 
| 1985 | 
> | 
            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) + | 
| 1986 | 
> | 
                                   angMom[1] * angMom[1] / Ia(1, 1) + | 
| 1987 | 
> | 
                                   angMom[2] * angMom[2] / Ia(2, 2)); | 
| 1988 | 
  | 
            binDOF[binNo] += 3; | 
| 1989 | 
  | 
          } | 
| 1990 | 
  | 
        } | 
| 1992 | 
  | 
    } | 
| 1993 | 
  | 
     | 
| 1994 | 
  | 
#ifdef IS_MPI | 
| 1995 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0], | 
| 1996 | 
< | 
                              nBins_, MPI::INT, MPI::SUM); | 
| 1997 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binMass[0], | 
| 1998 | 
< | 
                              nBins_, MPI::REALTYPE, MPI::SUM); | 
| 1999 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPx[0], | 
| 2000 | 
< | 
                              nBins_, MPI::REALTYPE, MPI::SUM); | 
| 2001 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPy[0], | 
| 2002 | 
< | 
                              nBins_, MPI::REALTYPE, MPI::SUM); | 
| 2003 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPz[0], | 
| 2004 | 
< | 
                              nBins_, MPI::REALTYPE, MPI::SUM); | 
| 2005 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegax[0], | 
| 2006 | 
< | 
                              nBins_, MPI::REALTYPE, MPI::SUM); | 
| 2007 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegay[0], | 
| 2008 | 
< | 
                              nBins_, MPI::REALTYPE, MPI::SUM); | 
| 2009 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegaz[0], | 
| 2010 | 
< | 
                              nBins_, MPI::REALTYPE, MPI::SUM); | 
| 2011 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0], | 
| 2012 | 
< | 
                              nBins_, MPI::REALTYPE, MPI::SUM); | 
| 2013 | 
< | 
    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[0], | 
| 2014 | 
< | 
                              nBins_, MPI::INT, MPI::SUM); | 
| 1995 | 
> | 
 | 
| 1996 | 
> | 
    for (int i = 0; i < nBins_; i++) { | 
| 1997 | 
> | 
       | 
| 1998 | 
> | 
      MPI_Allreduce(MPI_IN_PLACE, &binCount[i], | 
| 1999 | 
> | 
                    1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 2000 | 
> | 
      MPI_Allreduce(MPI_IN_PLACE, &binMass[i], | 
| 2001 | 
> | 
                    1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 2002 | 
> | 
      MPI_Allreduce(MPI_IN_PLACE, &binP[i], | 
| 2003 | 
> | 
                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 2004 | 
> | 
      MPI_Allreduce(MPI_IN_PLACE, &binL[i], | 
| 2005 | 
> | 
                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 2006 | 
> | 
      MPI_Allreduce(MPI_IN_PLACE, &binI[i], | 
| 2007 | 
> | 
                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 2008 | 
> | 
      MPI_Allreduce(MPI_IN_PLACE, &binKE[i], | 
| 2009 | 
> | 
                    1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 2010 | 
> | 
      MPI_Allreduce(MPI_IN_PLACE, &binDOF[i], | 
| 2011 | 
> | 
                    1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 2012 | 
> | 
      //MPI_Allreduce(MPI_IN_PLACE, &binOmega[i], | 
| 2013 | 
> | 
      //                          1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 2014 | 
> | 
    } | 
| 2015 | 
> | 
     | 
| 2016 | 
  | 
#endif | 
| 2017 | 
  | 
 | 
| 2018 | 
< | 
    Vector3d vel; | 
| 1990 | 
< | 
    Vector3d aVel; | 
| 2018 | 
> | 
    Vector3d omega; | 
| 2019 | 
  | 
    RealType den; | 
| 2020 | 
  | 
    RealType temp; | 
| 2021 | 
  | 
    RealType z; | 
| 2032 | 
  | 
        den = binMass[i] * 3.0 * PhysicalConstants::densityConvert | 
| 2033 | 
  | 
          / (4.0 * M_PI * (pow(router,3) - pow(rinner,3))); | 
| 2034 | 
  | 
      } | 
| 2035 | 
< | 
      vel.x() = binPx[i] / binMass[i]; | 
| 2036 | 
< | 
      vel.y() = binPy[i] / binMass[i]; | 
| 2037 | 
< | 
      vel.z() = binPz[i] / binMass[i]; | 
| 2038 | 
< | 
      aVel.x() = binOmegax[i] / binCount[i]; | 
| 2039 | 
< | 
      aVel.y() = binOmegay[i] / binCount[i]; | 
| 2012 | 
< | 
      aVel.z() = binOmegaz[i] / binCount[i]; | 
| 2035 | 
> | 
      vel = binP[i] / binMass[i]; | 
| 2036 | 
> | 
 | 
| 2037 | 
> | 
      omega = binI[i].inverse() * binL[i]; | 
| 2038 | 
> | 
 | 
| 2039 | 
> | 
      // omega = binOmega[i] / binCount[i]; | 
| 2040 | 
  | 
 | 
| 2041 | 
  | 
      if (binCount[i] > 0) { | 
| 2042 | 
  | 
        // only add values if there are things to add | 
| 2059 | 
  | 
              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel); | 
| 2060 | 
  | 
              break; | 
| 2061 | 
  | 
            case ANGULARVELOCITY:   | 
| 2062 | 
< | 
              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel); | 
| 2062 | 
> | 
              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(omega); | 
| 2063 | 
  | 
              break; | 
| 2064 | 
  | 
            case DENSITY: | 
| 2065 | 
  | 
              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den); | 
| 2106 | 
  | 
     | 
| 2107 | 
  | 
#ifdef IS_MPI | 
| 2108 | 
  | 
    // If we're the root node, should we print out the results | 
| 2109 | 
< | 
    int worldRank = MPI::COMM_WORLD.Get_rank(); | 
| 2109 | 
> | 
    int worldRank; | 
| 2110 | 
> | 
    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank); | 
| 2111 | 
> | 
 | 
| 2112 | 
  | 
    if (worldRank == 0) { | 
| 2113 | 
  | 
#endif | 
| 2114 | 
  | 
      rnemdFile_.open(rnemdFileName_.c_str(), std::ios::out | std::ios::trunc ); |