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 972 by gezelter, Wed May 24 18:31:12 2006 UTC vs.
Revision 977 by tim, Thu May 25 16:27:27 2006 UTC

# Line 69 | Line 69 | namespace oopse {
69      
70    }
71    
72 <  bool ApproximationModel::calcHydroProps(Shape* shape, double viscosity, double temperature) {
72 >  bool ApproximationModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
73      
74      bool ret = true;
75      HydroProps cr;
# Line 82 | Line 82 | namespace oopse {
82      return true;    
83    }
84    
85 <  bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) {
85 >  bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& 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;
# Line 97 | Line 97 | namespace oopse {
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;                
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 <              double constant = 8.0 * NumericConstant::PI * viscosity * rij;
105 >              RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
106                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
107              }else {
108 <              double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
108 >              RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
109                Tij(0, 0) = constant;
110                Tij(1, 1) = constant;
111                Tij(2, 2) = constant;
# Line 132 | Line 132 | namespace oopse {
132      
133      //calculate the total volume
134      
135 <    double volume = 0.0;
135 >    RealType 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      }
# Line 149 | Line 149 | namespace oopse {
149        }
150      }
151      
152 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
152 >    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
153      Xiott *= convertConstant;
154      Xiotr *= convertConstant;
155      Xiorr *= convertConstant;
# Line 185 | Line 185 | namespace oopse {
185      Xirrr = Xiorr - Uor * Xiott * Uor + Xiotr * Uor - Uor * Xiotr.transpose();
186      
187  
188 <    SquareMatrix<double,6> Xir6x6;
189 <    SquareMatrix<double,6> Dr6x6;
188 >    SquareMatrix<RealType,6> Xir6x6;
189 >    SquareMatrix<RealType,6> Dr6x6;
190  
191      Xir6x6.setSubMatrix(0, 0, Xirtt);
192      Xir6x6.setSubMatrix(0, 3, Xirtr.transpose());
# Line 202 | Line 202 | namespace oopse {
202      Dr6x6.getSubMatrix(0, 3, Drrt);
203      Dr6x6.getSubMatrix(3, 0, Drtr);
204      Dr6x6.getSubMatrix(3, 3, Drrr);
205 <    double kt = OOPSEConstant::kB * temperature ;
205 >    RealType kt = OOPSEConstant::kB * temperature ;
206      Drtt *= kt;
207      Drrt *= kt;
208      Drtr *= kt;
# Line 246 | Line 246 | namespace oopse {
246      return true;
247   }
248    
249 <  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) {
249 >  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) {
250      
251      int nbeads = beads.size();
252 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
253 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
252 >    DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
253 >    DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads);
254      Mat3x3d I;
255      I(0, 0) = 1.0;
256      I(1, 1) = 1.0;
# Line 261 | Line 261 | namespace oopse {
261          Mat3x3d Tij;
262          if (i != j ) {
263            Vector3d Rij = beads[i].pos - beads[j].pos;
264 <          double rij = Rij.length();
265 <          double rij2 = rij * rij;
266 <          double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
264 >          RealType rij = Rij.length();
265 >          RealType rij2 = rij * rij;
266 >          RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
267            Mat3x3d tmpMat;
268            tmpMat = outProduct(Rij, Rij) / rij2;
269 <          double constant = 8.0 * NumericConstant::PI * viscosity * rij;
269 >          RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
270            Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
271          }else {
272 <          double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
272 >          RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
273            Tij(0, 0) = constant;
274            Tij(1, 1) = constant;
275            Tij(2, 2) = constant;
# Line 296 | Line 296 | namespace oopse {
296  
297      //calculate the total volume
298  
299 <    double volume = 0.0;
299 >    RealType volume = 0.0;
300      for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) {
301        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
302      }
# Line 313 | Line 313 | namespace oopse {
313        }
314      }
315      
316 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
316 >    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
317      Xitt *= convertConstant;
318      Xitr *= convertConstant;
319      Xirr *= convertConstant;
320      
321 <    double kt = OOPSEConstant::kB * temperature;
321 >    RealType kt = OOPSEConstant::kB * temperature;
322      
323      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
324      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
# Line 377 | Line 377 | namespace oopse {
377      Ddrr = Dorr;
378      Ddtr = Dotr + Dorr * Uod;
379  
380 <    SquareMatrix<double, 6> Dd;
380 >    SquareMatrix<RealType, 6> Dd;
381      Dd.setSubMatrix(0, 0, Ddtt);
382      Dd.setSubMatrix(0, 3, Ddtr.transpose());
383      Dd.setSubMatrix(3, 0, Ddtr);
384      Dd.setSubMatrix(3, 3, Ddrr);    
385 <    SquareMatrix<double, 6> Xid;
385 >    SquareMatrix<RealType, 6> Xid;
386      Ddtt *= kt;
387      Ddtr *=kt;
388      Ddrr *= kt;

Diff Legend

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