ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp
(Generate patch)

Comparing trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp (file contents):
Revision 892 by tim, Thu Feb 23 23:16:43 2006 UTC vs.
Revision 893 by tim, Fri Feb 24 21:17:05 2006 UTC

# Line 95 | Line 95 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
95                  double rij2 = rij * rij;
96                  double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;                
97                  Mat3x3d tmpMat;
98 <                tmpMat = outProduct(beads_[i].pos, beads_[j].pos) / rij2;
98 >                tmpMat = outProduct(Rij, Rij) / rij2;
99                  double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
100                  Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
101              }else {
# Line 105 | Line 105 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
105                  Tij(2, 2) = constant;
106              }
107              B.setSubMatrix(i*3, j*3, Tij);
108 +            std::cout << Tij << std::endl;
109          }
110      }
111  
112 +    std::cout << "B=\n"
113 +                  << B << std::endl;
114      //invert B Matrix
115      invertMatrix(B, C);
116 +
117 +    std::cout << "C=\n"
118 +                  << C << std::endl;
119 +
120      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
121      std::vector<Mat3x3d> U;
122      for (int i = 0; i < nbeads; ++i) {
# Line 122 | Line 129 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
129      Mat3x3d Xitt;
130      Mat3x3d Xirr;
131      Mat3x3d Xitr;
132 +
133 +    //calculate the total volume
134 +
135 +    double volume = 0.0;
136 +    for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
137 +        volume = 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
138 +    }
139          
140      for (std::size_t i = 0; i < nbeads; ++i) {
141          for (std::size_t j = 0; j < nbeads; ++j) {
# Line 129 | Line 143 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
143              C.getSubMatrix(i*3, j*3, Cij);
144              
145              Xitt += Cij;
146 <            Xirr += U[i] * Cij;
147 <            Xitr += U[i] * Cij * U[j];            
146 >            Xitr += U[i] * Cij;
147 >            Xirr += -U[i] * Cij * U[j];            
148 >            //Xirr += -U[i] * Cij * U[j] + (0.166*6 * viscosity_ * volume) * I;            
149          }
150      }
151  
# Line 150 | Line 165 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
165      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
166      Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
167  
168 <    Mat3x3d XirrInv(0.0);
154 <    Mat3x3d XirrCopy;
155 <    XirrCopy = Xirr;
168 >    const static Mat3x3d zeroMat(0.0);
169      
170      Mat3x3d XittInv(0.0);
171 <    Mat3x3d XittCopy;
172 <    XittCopy = Xitt;
173 <    invertMatrix(XittCopy, XittInv);
171 >    XittInv = Xitt.inverse();
172 >    
173 >    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
174 >    //( a11 a12 0)
175 >    //( a21 a22 0)
176 >    //( 0    0    0)
177 >    Mat3x3d XirrInv;
178 >    XirrInv = Xirr.inverse();
179  
180      Mat3x3d tmp;
181      Mat3x3d tmpInv;
182      tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
183 +    tmpInv = tmp.inverse();
184  
166 –    const static Mat3x3d zeroMat(0.0);
167 –    if (!invertMatrix(tmp, tmpInv)) {
168 –        tmpInv = zeroMat;
169 –    }
170 –
185      Dott = kt * tmpInv;
186 <    Dotr = -kt*XirrInv * Xitr * tmpInv;
186 >    Dotr = -kt*XirrInv * Xitr * tmpInv* 1.0E8;
187  
188 <    tmp = Xirr - Xitr * XittInv * Xitr.transpose();
188 >    tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
189 >    tmpInv = tmp.inverse();
190      
191 <    if(!invertMatrix(tmp, tmpInv)) {
177 <        tmpInv = zeroMat;
178 <    }
179 <    Dorr = kt * tmpInv;
191 >    Dorr = kt * tmpInv*1.0E16;
192  
193      //Do.getSubMatrix(0, 0 , Dott);
194      //Do.getSubMatrix(3, 0, Dotr);
# Line 198 | Line 210 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
210      tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
211      tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
212  
213 <    if(!invertMatrix(tmp, tmpInv)) {
202 <        tmpInv = zeroMat;
203 <    }
213 >    tmpInv = tmp.inverse();
214      
215      Vector3d rod = tmpInv * tmpVec;
216  
# Line 246 | Line 256 | void HydrodynamicsModel::writeDiffCenterAndDiffTensor(
256      os << "//translational diffusion tensor" << std::endl;
257      os << props_.transDiff << std::endl;
258  
259 <    os << "//translational diffusion tensor" << std::endl;
259 >    os << "//translation-rotation coupling diffusion tensor" << std::endl;
260      os << props_.transRotDiff << std::endl;
261  
262      os << "//rotational diffusion tensor" << std::endl;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines