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

Comparing trunk/src/applications/hydrodynamics/ApproximationModel.cpp (file contents):
Revision 954 by tim, Wed May 10 01:44:48 2006 UTC vs.
Revision 981 by gezelter, Mon Jun 5 18:24:45 2006 UTC

# Line 44 | Line 44
44   #include "math/DynamicRectMatrix.hpp"
45   #include "math/SquareMatrix3.hpp"
46   #include "utils/OOPSEConstant.hpp"
47 < #include "applications/hydrodynamics/Spheric.hpp"
48 < #include "applications/hydrodynamics/Ellipsoid.hpp"
47 > #include "hydrodynamics/Sphere.hpp"
48 > #include "hydrodynamics/Ellipsoid.hpp"
49   #include "applications/hydrodynamics/CompositeShape.hpp"
50   #include "math/LU.hpp"
51   #include "utils/simError.h"
# Line 57 | Line 57 | namespace oopse {
57   * Biophysical Journal, 75(6), 3044, 1999
58   */
59  
60 < ApproximationModel::ApproximationModel(StuntDouble* sd, SimInfo* info): HydrodynamicsModel(sd, info){
61 <
62 < }
63 <
64 < bool ApproximationModel::calcHydroProps(Spheric* spheric, double viscosity, double temperature) {
65 <    return internalCalcHydroProps(static_cast<Shape*>(spheric), viscosity, temperature);
66 < }
67 <
68 < bool ApproximationModel::calcHydroProps(Ellipsoid* ellipsoid, double viscosity, double temperature) {
69 <    return internalCalcHydroProps(static_cast<Shape*>(ellipsoid), viscosity, temperature);
70 < }
71 < bool ApproximationModel::calcHydroProps(CompositeShape* compositeShape, double viscosity, double temperature) {
72 <    return internalCalcHydroProps(static_cast<Shape*>(compositeShape), viscosity, temperature);
73 < }
74 <
75 < void ApproximationModel::init() {
60 >  ApproximationModel::ApproximationModel(StuntDouble* sd, SimInfo* info): HydrodynamicsModel(sd, info){    
61 >  }
62 >  
63 >  void ApproximationModel::init() {
64      if (!createBeads(beads_)) {
65        sprintf(painCave.errMsg, "ApproximationModel::init() : Can not create beads\n");
66        painCave.isFatal = 1;
67        simError();        
68      }
69 <
70 < }
71 <
72 < bool ApproximationModel::internalCalcHydroProps(Shape* shape, double viscosity, double temperature) {
73 <
69 >    
70 >  }
71 >  
72 >  bool ApproximationModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
73 >    
74      bool ret = true;
75 <    HydroProps cr;
76 <    HydroProps cd;
75 >    HydroProp* cr;
76 >    HydroProp* cd;
77      calcHydroPropsAtCR(beads_, viscosity, temperature, cr);
78      //calcHydroPropsAtCD(beads_, viscosity, temperature, cd);
79      setCR(cr);
80      setCD(cd);
81      
82      return true;    
83 < }
84 <
85 < bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) {
86 <
83 >  }
84 >  
85 >  bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProp* cr) {
86 >    
87      int nbeads = beads.size();
88 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
89 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
88 >    DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
89 >    DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads);
90      Mat3x3d I;
91      I(0, 0) = 1.0;
92      I(1, 1) = 1.0;
93      I(2, 2) = 1.0;
94      
95      for (std::size_t i = 0; i < nbeads; ++i) {
96 <        for (std::size_t j = 0; j < nbeads; ++j) {
97 <            Mat3x3d Tij;
96 >      for (std::size_t j = 0; j < nbeads; ++j) {
97 >        Mat3x3d Tij;
98              if (i != j ) {
99 <                Vector3d Rij = beads[i].pos - beads[j].pos;
100 <                double rij = Rij.length();
101 <                double rij2 = rij * rij;
102 <                double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
103 <                Mat3x3d tmpMat;
104 <                tmpMat = outProduct(Rij, Rij) / rij2;
105 <                double constant = 8.0 * NumericConstant::PI * viscosity * rij;
106 <                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
99 >              Vector3d Rij = beads[i].pos - beads[j].pos;
100 >              RealType rij = Rij.length();
101 >              RealType rij2 = rij * rij;
102 >              RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
103 >              Mat3x3d tmpMat;
104 >              tmpMat = outProduct(Rij, Rij) / rij2;
105 >              RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
106 >              RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0;
107 >              RealType tmp2 = 1.0 - sumSigma2OverRij2;
108 >              Tij = (tmp1 * I + tmp2 * tmpMat ) / constant;
109              }else {
110 <                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
111 <                Tij(0, 0) = constant;
112 <                Tij(1, 1) = constant;
113 <                Tij(2, 2) = constant;
110 >              RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
111 >              Tij(0, 0) = constant;
112 >              Tij(1, 1) = constant;
113 >              Tij(2, 2) = constant;
114              }
115              B.setSubMatrix(i*3, j*3, Tij);
116 <        }
116 >      }
117      }
118 <
118 >    
119      //invert B Matrix
120      invertMatrix(B, C);
121      
122      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
123      std::vector<Mat3x3d> U;
124      for (int i = 0; i < nbeads; ++i) {
125 <        Mat3x3d currU;
126 <        currU.setupSkewMat(beads[i].pos);
127 <        U.push_back(currU);
125 >      Mat3x3d currU;
126 >      currU.setupSkewMat(beads[i].pos);
127 >      U.push_back(currU);
128      }
129      
130      //calculate Xi matrix at arbitrary origin O
131      Mat3x3d Xiott;
132      Mat3x3d Xiorr;
133      Mat3x3d Xiotr;
134 <
134 >    
135      //calculate the total volume
136 <
137 <    double volume = 0.0;
136 >    
137 >    RealType volume = 0.0;
138      for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) {
139 <        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
139 >      volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
140      }
141 <        
141 >    
142      for (std::size_t i = 0; i < nbeads; ++i) {
143 <        for (std::size_t j = 0; j < nbeads; ++j) {
144 <            Mat3x3d Cij;
145 <            C.getSubMatrix(i*3, j*3, Cij);
146 <            
147 <            Xiott += Cij;
148 <            Xiotr += U[i] * Cij;
149 <            //Xiorr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;    
150 <            Xiorr += -U[i] * Cij * U[j];
151 <        }
143 >      for (std::size_t j = 0; j < nbeads; ++j) {
144 >        Mat3x3d Cij;
145 >        C.getSubMatrix(i*3, j*3, Cij);
146 >        
147 >        Xiott += Cij;
148 >        Xiotr += U[i] * Cij;
149 >        //Xiorr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;    
150 >        Xiorr += -U[i] * Cij * U[j];
151 >      }
152      }
153 <
154 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
153 >    
154 >    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
155      Xiott *= convertConstant;
156      Xiotr *= convertConstant;
157      Xiorr *= convertConstant;
158      
169 –
159      
160 +    
161      Mat3x3d tmp;
162      Mat3x3d tmpInv;
163      Vector3d tmpVec;
# Line 197 | Line 187 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
187      Xirrr = Xiorr - Uor * Xiott * Uor + Xiotr * Uor - Uor * Xiotr.transpose();
188      
189  
190 <    SquareMatrix<double,6> Xir6x6;
191 <    SquareMatrix<double,6> Dr6x6;
190 >    SquareMatrix<RealType,6> Xir6x6;
191 >    SquareMatrix<RealType,6> Dr6x6;
192  
193      Xir6x6.setSubMatrix(0, 0, Xirtt);
194      Xir6x6.setSubMatrix(0, 3, Xirtr.transpose());
# Line 214 | Line 204 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
204      Dr6x6.getSubMatrix(0, 3, Drrt);
205      Dr6x6.getSubMatrix(3, 0, Drtr);
206      Dr6x6.getSubMatrix(3, 3, Drrr);
207 <    double kt = OOPSEConstant::kB * temperature ;
207 >    RealType kt = OOPSEConstant::kB * temperature ;
208      Drtt *= kt;
209      Drrt *= kt;
210      Drtr *= kt;
# Line 223 | Line 213 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
213      Xirtr *= OOPSEConstant::kb * temperature;
214      Xirrr *= OOPSEConstant::kb * temperature;
215      
216 +    Mat6x6d Xi, D;
217  
218 <    cr.center = ror;
219 <    cr.Xi.setSubMatrix(0, 0, Xirtt);
220 <    cr.Xi.setSubMatrix(0, 3, Xirtr);
221 <    cr.Xi.setSubMatrix(3, 0, Xirtr);
222 <    cr.Xi.setSubMatrix(3, 3, Xirrr);
223 <    cr.D.setSubMatrix(0, 0, Drtt);
224 <    cr.D.setSubMatrix(0, 3, Drrt);
225 <    cr.D.setSubMatrix(3, 0, Drtr);
226 <    cr.D.setSubMatrix(3, 3, Drrr);    
218 >    cr->setCOR(ror);
219 >
220 >    Xi.setSubMatrix(0, 0, Xirtt);
221 >    Xi.setSubMatrix(0, 3, Xirtr);
222 >    Xi.setSubMatrix(3, 0, Xirtr);
223 >    Xi.setSubMatrix(3, 3, Xirrr);
224 >
225 >    cr->setXi(Xi);
226 >
227 >    D.setSubMatrix(0, 0, Drtt);
228 >    D.setSubMatrix(0, 3, Drrt);
229 >    D.setSubMatrix(3, 0, Drtr);
230 >    D.setSubMatrix(3, 3, Drrr);    
231 >
232 >    cr->setD(D);
233      
234      std::cout << "-----------------------------------------\n";
235      std::cout << "center of resistance :" << std::endl;
# Line 257 | Line 254 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
254  
255      return true;
256   }
257 <
258 < bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) {
259 <
257 >  
258 >  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProp* cr) {
259 >    
260      int nbeads = beads.size();
261 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
262 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
261 >    DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
262 >    DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads);
263      Mat3x3d I;
264      I(0, 0) = 1.0;
265      I(1, 1) = 1.0;
266      I(2, 2) = 1.0;
267      
268      for (std::size_t i = 0; i < nbeads; ++i) {
269 <        for (std::size_t j = 0; j < nbeads; ++j) {
270 <            Mat3x3d Tij;
271 <            if (i != j ) {
272 <                Vector3d Rij = beads[i].pos - beads[j].pos;
273 <                double rij = Rij.length();
274 <                double rij2 = rij * rij;
275 <                double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
276 <                Mat3x3d tmpMat;
277 <                tmpMat = outProduct(Rij, Rij) / rij2;
278 <                double constant = 8.0 * NumericConstant::PI * viscosity * rij;
279 <                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
280 <            }else {
281 <                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
282 <                Tij(0, 0) = constant;
283 <                Tij(1, 1) = constant;
284 <                Tij(2, 2) = constant;
285 <            }
286 <            B.setSubMatrix(i*3, j*3, Tij);
269 >      for (std::size_t j = 0; j < nbeads; ++j) {
270 >        Mat3x3d Tij;
271 >        if (i != j ) {
272 >          Vector3d Rij = beads[i].pos - beads[j].pos;
273 >          RealType rij = Rij.length();
274 >          RealType rij2 = rij * rij;
275 >          RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
276 >          Mat3x3d tmpMat;
277 >          tmpMat = outProduct(Rij, Rij) / rij2;
278 >          RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
279 >          RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0;
280 >          RealType tmp2 = 1.0 - sumSigma2OverRij2;
281 >          Tij = (tmp1 * I + tmp2 * tmpMat ) / constant;
282 >        }else {
283 >          RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
284 >          Tij(0, 0) = constant;
285 >          Tij(1, 1) = constant;
286 >          Tij(2, 2) = constant;
287          }
288 +        B.setSubMatrix(i*3, j*3, Tij);
289 +      }
290      }
291 <
291 >    
292      //invert B Matrix
293      invertMatrix(B, C);
294 <
294 >    
295      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
296      std::vector<Mat3x3d> U;
297      for (int i = 0; i < nbeads; ++i) {
298 <        Mat3x3d currU;
299 <        currU.setupSkewMat(beads[i].pos);
300 <        U.push_back(currU);
298 >      Mat3x3d currU;
299 >      currU.setupSkewMat(beads[i].pos);
300 >      U.push_back(currU);
301      }
302      
303      //calculate Xi matrix at arbitrary origin O
# Line 308 | Line 307 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
307  
308      //calculate the total volume
309  
310 <    double volume = 0.0;
310 >    RealType volume = 0.0;
311      for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) {
312 <        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
312 >      volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
313      }
314 <        
314 >    
315      for (std::size_t i = 0; i < nbeads; ++i) {
316 <        for (std::size_t j = 0; j < nbeads; ++j) {
317 <            Mat3x3d Cij;
318 <            C.getSubMatrix(i*3, j*3, Cij);
316 >      for (std::size_t j = 0; j < nbeads; ++j) {
317 >        Mat3x3d Cij;
318 >        C.getSubMatrix(i*3, j*3, Cij);
319              
320 <            Xitt += Cij;
321 <            Xitr += U[i] * Cij;
320 >        Xitt += Cij;
321 >        Xitr += U[i] * Cij;
322              //Xirr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;            
323 <            Xirr += -U[i] * Cij * U[j];
324 <        }
323 >        Xirr += -U[i] * Cij * U[j];
324 >      }
325      }
326 <
327 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
326 >    
327 >    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
328      Xitt *= convertConstant;
329      Xitr *= convertConstant;
330      Xirr *= convertConstant;
331 <
332 <    double kt = OOPSEConstant::kB * temperature;
333 <
331 >    
332 >    RealType kt = OOPSEConstant::kB * temperature;
333 >    
334      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
335      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
336      Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
337 <
337 >    
338      const static Mat3x3d zeroMat(0.0);
339      
340      Mat3x3d XittInv(0.0);
# Line 389 | Line 388 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
388      Ddrr = Dorr;
389      Ddtr = Dotr + Dorr * Uod;
390  
391 <    SquareMatrix<double, 6> Dd;
391 >    SquareMatrix<RealType, 6> Dd;
392      Dd.setSubMatrix(0, 0, Ddtt);
393      Dd.setSubMatrix(0, 3, Ddtr.transpose());
394      Dd.setSubMatrix(3, 0, Ddtr);
395      Dd.setSubMatrix(3, 3, Ddrr);    
396 <    SquareMatrix<double, 6> Xid;
396 >    SquareMatrix<RealType, 6> Xid;
397      Ddtt *= kt;
398      Ddtr *=kt;
399      Ddrr *= kt;
# Line 406 | Line 405 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
405      //Xid /= OOPSEConstant::energyConvert;
406      Xid *= OOPSEConstant::kb * temperature;
407  
408 <    cr.center = rod;
410 <    cr.D.setSubMatrix(0, 0, Ddtt);
411 <    cr.D.setSubMatrix(0, 3, Ddtr);
412 <    cr.D.setSubMatrix(3, 0, Ddtr);
413 <    cr.D.setSubMatrix(3, 3, Ddrr);
414 <    cr.Xi = Xid;
408 >    Mat6x6d Xi, D;
409  
410 +    cr->setCOR(rod);
411 +
412 +    cr->setXi(Xid);
413 +
414 +    D.setSubMatrix(0, 0, Ddtt);
415 +    D.setSubMatrix(0, 3, Ddtr);
416 +    D.setSubMatrix(3, 0, Ddtr);
417 +    D.setSubMatrix(3, 3, Ddrr);
418 +
419 +    cr->setD(D);
420 +
421      std::cout << "viscosity = " << viscosity << std::endl;
422      std::cout << "temperature = " << temperature << std::endl;
423      std::cout << "center of diffusion :" << std::endl;
# Line 446 | Line 451 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
451      std::cout << Xidrr << std::endl;
452  
453      return true;
454 <      
455 < }
454 >    
455 >  }
456  
457 <
453 < void ApproximationModel::writeBeads(std::ostream& os) {
457 >  void ApproximationModel::writeBeads(std::ostream& os) {
458      std::vector<BeadParam>::iterator iter;
459      os << beads_.size() << std::endl;
460      os << "Generated by Hydro" << std::endl;
461      for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
462 <        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
462 >      os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
463      }
464 <
464 >    
465 >  }    
466   }
462 –
463 –
464 –
465 – }

Diff Legend

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