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

Comparing trunk/OOPSE-4/src/applications/hydrodynamics/HydrodynamicsModel.cpp (file contents):
Revision 2596 by tim, Wed Feb 22 20:35:16 2006 UTC vs.
Revision 2597 by tim, Thu Feb 23 23:16:43 2006 UTC

# Line 43 | Line 43 | namespace oopse {
43   #include "math/LU.hpp"
44   #include "math/DynamicRectMatrix.hpp"
45   #include "math/SquareMatrix3.hpp"
46 + #include "utils/OOPSEConstant.hpp"
47   namespace oopse {
48   /**
49   * Reference:
# Line 50 | Line 51 | bool HydrodynamicsModel::calcHydrodyanmicsProps(double
51   * Comparison of Different Modeling and Computational Procedures.
52   * Biophysical Journal, 75(6), 3044, 1999
53   */
54 < bool HydrodynamicsModel::calcHydrodyanmicsProps(double eta) {
54 >
55 > HydrodynamicsModel::HydrodynamicsModel(StuntDouble* sd, const DynamicProperty& extraParams) : sd_(sd){
56 >    DynamicProperty::const_iterator iter;
57 >
58 >    iter = extraParams.find("Viscosity");
59 >    if (iter != extraParams.end()) {
60 >        boost::any param = iter->second;
61 >        viscosity_ = boost::any_cast<double>(param);
62 >    }else {
63 >        std::cout << "HydrodynamicsModel Error\n" ;
64 >    }
65 >
66 >    iter = extraParams.find("Temperature");
67 >    if (iter != extraParams.end()) {
68 >        boost::any param = iter->second;
69 >        temperature_ = boost::any_cast<double>(param);
70 >    }else {
71 >        std::cout << "HydrodynamicsModel Error\n" ;
72 >    }    
73 > }
74 >
75 > bool HydrodynamicsModel::calcHydrodyanmicsProps() {
76      if (!createBeads(beads_)) {
77          std::cout << "can not create beads" << std::endl;
78          return false;
# Line 60 | Line 82 | bool HydrodynamicsModel::calcHydrodyanmicsProps(double
82      DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
83      DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
84      Mat3x3d I;
85 +    I(0, 0) = 1.0;
86 +    I(1, 1) = 1.0;
87 +    I(2, 2) = 1.0;
88 +    
89      for (std::size_t i = 0; i < nbeads; ++i) {
90          for (std::size_t j = 0; j < nbeads; ++j) {
91              Mat3x3d Tij;
# Line 70 | Line 96 | bool HydrodynamicsModel::calcHydrodyanmicsProps(double
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;
99 <                double constant = 8.0 * NumericConstant::PI * eta * rij;
99 >                double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
100                  Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
101              }else {
102 <                double constant = 1.0 / (6.0 * NumericConstant::PI * eta * beads_[i].radius);
102 >                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity_ * beads_[i].radius);
103                  Tij(0, 0) = constant;
104                  Tij(1, 1) = constant;
105                  Tij(2, 2) = constant;
# Line 84 | Line 110 | bool HydrodynamicsModel::calcHydrodyanmicsProps(double
110  
111      //invert B Matrix
112      invertMatrix(B, C);
87 –    
113      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
114      std::vector<Mat3x3d> U;
115      for (int i = 0; i < nbeads; ++i) {
# Line 115 | Line 140 | bool HydrodynamicsModel::calcHydrodyanmicsProps(double
140      Xi.setSubMatrix(0, 0, Xitt);
141      Xi.setSubMatrix(0, 3, Xitr.transpose());
142      Xi.setSubMatrix(3, 0, Xitr);
143 <    Xi.setSubMatrix(3, 3, Xitt);
144 <    invertMatrix(Xi, Do);
143 >    Xi.setSubMatrix(3, 3, Xirr);
144 >    //invertMatrix(Xi, Do);
145 >    double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
146 >    //Do *= kt;    
147  
148 +
149      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
150      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
151      Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
124 –    Do.getSubMatrix(0, 0 , Dott);
125 –    Do.getSubMatrix(3, 0, Dotr);
126 –    Do.getSubMatrix(3, 3, Dorr);
152  
153 +    Mat3x3d XirrInv(0.0);
154 +    Mat3x3d XirrCopy;
155 +    XirrCopy = Xirr;
156 +    
157 +    Mat3x3d XittInv(0.0);
158 +    Mat3x3d XittCopy;
159 +    XittCopy = Xitt;
160 +    invertMatrix(XittCopy, XittInv);
161 +
162 +    Mat3x3d tmp;
163 +    Mat3x3d tmpInv;
164 +    tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
165 +
166 +    const static Mat3x3d zeroMat(0.0);
167 +    if (!invertMatrix(tmp, tmpInv)) {
168 +        tmpInv = zeroMat;
169 +    }
170 +
171 +    Dott = kt * tmpInv;
172 +    Dotr = -kt*XirrInv * Xitr * tmpInv;
173 +
174 +    tmp = Xirr - Xitr * XittInv * Xitr.transpose();
175 +    
176 +    if(!invertMatrix(tmp, tmpInv)) {
177 +        tmpInv = zeroMat;
178 +    }
179 +    Dorr = kt * tmpInv;
180 +
181 +    //Do.getSubMatrix(0, 0 , Dott);
182 +    //Do.getSubMatrix(3, 0, Dotr);
183 +    //Do.getSubMatrix(3, 3, Dorr);
184 +
185      //calculate center of diffusion
186 <    Mat3x3d tmpMat;
187 <    tmpMat(0, 0) = Dorr(1, 1) + Dorr(2, 2);
188 <    tmpMat(0, 1) = - Dorr(0, 1);
189 <    tmpMat(0, 2) = -Dorr(0, 2);
190 <    tmpMat(1, 0) = -Dorr(0, 1);
191 <    tmpMat(1, 1) = Dorr(0, 0)  + Dorr(2, 2);
192 <    tmpMat(1, 2) = -Dorr(1, 2);
193 <    tmpMat(2, 0) = -Dorr(0, 2);
194 <    tmpMat(2, 1) = -Dorr(1, 2);
138 <    tmpMat(2, 2) = Dorr(1, 1) + Dorr(0, 0);
186 >    tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
187 >    tmp(0, 1) = - Dorr(0, 1);
188 >    tmp(0, 2) = -Dorr(0, 2);
189 >    tmp(1, 0) = -Dorr(0, 1);
190 >    tmp(1, 1) = Dorr(0, 0)  + Dorr(2, 2);
191 >    tmp(1, 2) = -Dorr(1, 2);
192 >    tmp(2, 0) = -Dorr(0, 2);
193 >    tmp(2, 1) = -Dorr(1, 2);
194 >    tmp(2, 2) = Dorr(1, 1) + Dorr(0, 0);
195  
196      Vector3d tmpVec;
197      tmpVec[0] = Dotr(1, 2) - Dotr(2, 1);
198      tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
199      tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
144 –        
145 –    Vector3d rod = tmpMat.inverse() * tmpVec;
200  
201 +    if(!invertMatrix(tmp, tmpInv)) {
202 +        tmpInv = zeroMat;
203 +    }
204 +    
205 +    Vector3d rod = tmpInv * tmpVec;
206 +
207      //calculate Diffusion Tensor at center of diffusion
208      Mat3x3d Uod;
209      Uod.setupSkewMat(rod);
# Line 165 | Line 225 | void HydrodynamicsModel::writeBeads(std::ostream& os)
225   }
226  
227   void HydrodynamicsModel::writeBeads(std::ostream& os) {
228 +    std::vector<BeadParam>::iterator iter;
229 +    os << beads_.size() << std::endl;
230 +    os << "Generated by Hydro" << std::endl;
231 +    for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
232 +        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
233 +    }
234  
235   }
236  
237   void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
238 +    os << "//viscosity = " << viscosity_ << std::endl;
239 +    os << "//temperature = " << temperature_<< std::endl;
240 +    std::vector<BeadParam>::iterator iter;
241 +    os << sd_->getType() << "\n";
242  
243 +    os << "//diffusion center" << std::endl;
244 +    os << props_.diffCenter << std::endl;
245 +
246 +    os << "//translational diffusion tensor" << std::endl;
247 +    os << props_.transDiff << std::endl;
248 +
249 +    os << "//translational diffusion tensor" << std::endl;
250 +    os << props_.transRotDiff << std::endl;
251 +
252 +    os << "//rotational diffusion tensor" << std::endl;
253 +    os << props_.rotDiff << std::endl;
254 +    
255 +    /*
256 +    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\n"
257 +
258 +    os << props_.transDiff(0, 0) << "\t" << props_.transDiff(0, 1) << "\t" << props_.transDiff(0, 2) << "\t"
259 +        << props_.transDiff(1, 0) << "\t" << props_.transDiff(1, 1) << "\t" << props_.transDiff(1, 2) << "\t"
260 +        << props_.transDiff(2, 0) << "\t" << props_.transDiff(2, 1) << "\t" << props_.transDiff(2, 2) << "\n";
261 +    
262 +    os << props_.transRotDiff(0, 0) << "\t" << props_.transRotDiff(0, 1) << "\t" << props_.transRotDiff(0, 2) << "\t"
263 +        << props_.transRotDiff(1, 0) << "\t" << props_.transRotDiff(1, 1) << "\t" << props_.transRotDiff(1, 2) << "\t"
264 +        << props_.transRotDiff(2, 0) << "\t" << props_.transRotDiff(2, 1) << "\t" << props_.transRotDiff(2, 2) << "\t"
265 +
266 +    os << props_.rotDiff(0, 0) << "\t" << props_.rotDiff(0, 1) << "\t" << props_.rotDiff(0, 2) << "\t"
267 +        << props_.rotDiff(1, 0) << "\t" << props_.rotDiff(1, 1) << "\t" << props_.rotDiff(1, 2) << "\t"
268 +        << props_.rotDiff(2, 0) << "\t" << props_.rotDiff(2, 1) << "\t" << props_.rotDiff(2, 2) << ";"
269 +        << std::endl;
270 +    */
271   }
272  
273   }

Diff Legend

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