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

Comparing trunk/OOPSE-4/src/applications/hydrodynamics/ApproximationModel.cpp (file contents):
Revision 2634 by tim, Fri Mar 17 23:20:35 2006 UTC vs.
Revision 2774 by tim, Thu May 25 17:02:00 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"
52   namespace oopse {
53   /**
54   * Reference:
# Line 56 | Line 57 | ApproximationModel::ApproximationModel(StuntDouble* sd
57   * Biophysical Journal, 75(6), 3044, 1999
58   */
59  
60 < ApproximationModel::ApproximationModel(StuntDouble* sd, SimInfo* info): HydrodynamicsModel(sd, info){
61 < /*
62 <    DynamicProperty::const_iterator iter;
63 <
63 <    iter = extraParams.find("Viscosity");
64 <    if (iter != extraParams.end()) {
65 <        boost::any param = iter->second;
66 <        viscosity = boost::any_cast<double>(param);
67 <    }else {
68 <        std::cout << "ApproximationModel Error\n" ;
69 <    }
70 <
71 <    iter = extraParams.find("Temperature");
72 <    if (iter != extraParams.end()) {
73 <        boost::any param = iter->second;
74 <        temperature = boost::any_cast<double>(param);
75 <    }else {
76 <        std::cout << "ApproximationModel Error\n" ;
77 <    }    
78 < */
79 < }
80 <
81 < bool ApproximationModel::calcHydroProps(Spheric* spheric, double viscosity, double temperature) {
82 <    return internalCalcHydroProps(static_cast<Shape*>(spheric), viscosity, temperature);
83 < }
84 <
85 < bool ApproximationModel::calcHydroProps(Ellipsoid* ellipsoid, double viscosity, double temperature) {
86 <    return internalCalcHydroProps(static_cast<Shape*>(ellipsoid), viscosity, temperature);
87 < }
88 < bool ApproximationModel::calcHydroProps(CompositeShape* compositeShape, double viscosity, double temperature) {
89 <    return internalCalcHydroProps(static_cast<Shape*>(compositeShape), viscosity, temperature);
90 < }
91 <
92 <
93 < bool ApproximationModel::internalCalcHydroProps(Shape* shape, double viscosity, double temperature) {
60 >  ApproximationModel::ApproximationModel(StuntDouble* sd, SimInfo* info): HydrodynamicsModel(sd, info){    
61 >  }
62 >  
63 >  void ApproximationModel::init() {
64      if (!createBeads(beads_)) {
65 <        std::cout << "can not create beads" << std::endl;
66 <        return false;
65 >      sprintf(painCave.errMsg, "ApproximationModel::init() : Can not create beads\n");
66 >      painCave.isFatal = 1;
67 >      simError();        
68      }
69 <
69 >    
70 >  }
71 >  
72 >  bool ApproximationModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
73 >    
74      bool ret = true;
75      HydroProps cr;
76      HydroProps cd;
77      calcHydroPropsAtCR(beads_, viscosity, temperature, cr);
78 <    calcHydroPropsAtCD(beads_, viscosity, temperature, cd);
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, 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;
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[i].radius*beads[i].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 <
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 <        }
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      
181 –
159      
160 +    
161      Mat3x3d tmp;
162      Mat3x3d tmpInv;
163      Vector3d tmpVec;
# Line 209 | 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 226 | 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 269 | Line 247 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
247  
248      return true;
249   }
250 <
251 < bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) {
252 <
250 >  
251 >  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) {
252 >    
253      int nbeads = beads.size();
254 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
255 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
254 >    DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
255 >    DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads);
256      Mat3x3d I;
257      I(0, 0) = 1.0;
258      I(1, 1) = 1.0;
259      I(2, 2) = 1.0;
260      
261      for (std::size_t i = 0; i < nbeads; ++i) {
262 <        for (std::size_t j = 0; j < nbeads; ++j) {
263 <            Mat3x3d Tij;
264 <            if (i != j ) {
265 <                Vector3d Rij = beads[i].pos - beads[j].pos;
266 <                double rij = Rij.length();
267 <                double rij2 = rij * rij;
268 <                double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[i].radius*beads[i].radius)) / rij2;                
269 <                Mat3x3d tmpMat;
270 <                tmpMat = outProduct(Rij, Rij) / rij2;
271 <                double constant = 8.0 * NumericConstant::PI * viscosity * rij;
272 <                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
273 <            }else {
274 <                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
275 <                Tij(0, 0) = constant;
276 <                Tij(1, 1) = constant;
277 <                Tij(2, 2) = constant;
278 <            }
279 <            B.setSubMatrix(i*3, j*3, Tij);
262 >      for (std::size_t j = 0; j < nbeads; ++j) {
263 >        Mat3x3d Tij;
264 >        if (i != j ) {
265 >          Vector3d Rij = beads[i].pos - beads[j].pos;
266 >          RealType rij = Rij.length();
267 >          RealType rij2 = rij * rij;
268 >          RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
269 >          Mat3x3d tmpMat;
270 >          tmpMat = outProduct(Rij, Rij) / rij2;
271 >          RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
272 >          RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0;
273 >          RealType tmp2 = 1.0 - sumSigma2OverRij2;
274 >          Tij = (tmp1 * I + tmp2 * tmpMat ) / constant;
275 >        }else {
276 >          RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
277 >          Tij(0, 0) = constant;
278 >          Tij(1, 1) = constant;
279 >          Tij(2, 2) = constant;
280          }
281 +        B.setSubMatrix(i*3, j*3, Tij);
282 +      }
283      }
284 <
284 >    
285      //invert B Matrix
286      invertMatrix(B, C);
287 <
287 >    
288      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
289      std::vector<Mat3x3d> U;
290      for (int i = 0; i < nbeads; ++i) {
291 <        Mat3x3d currU;
292 <        currU.setupSkewMat(beads[i].pos);
293 <        U.push_back(currU);
291 >      Mat3x3d currU;
292 >      currU.setupSkewMat(beads[i].pos);
293 >      U.push_back(currU);
294      }
295      
296      //calculate Xi matrix at arbitrary origin O
# Line 320 | Line 300 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
300  
301      //calculate the total volume
302  
303 <    double volume = 0.0;
303 >    RealType volume = 0.0;
304      for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) {
305 <        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
305 >      volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
306      }
307 <        
307 >    
308      for (std::size_t i = 0; i < nbeads; ++i) {
309 <        for (std::size_t j = 0; j < nbeads; ++j) {
310 <            Mat3x3d Cij;
311 <            C.getSubMatrix(i*3, j*3, Cij);
309 >      for (std::size_t j = 0; j < nbeads; ++j) {
310 >        Mat3x3d Cij;
311 >        C.getSubMatrix(i*3, j*3, Cij);
312              
313 <            Xitt += Cij;
314 <            Xitr += U[i] * Cij;
315 <            Xirr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;            
316 <        }
313 >        Xitt += Cij;
314 >        Xitr += U[i] * Cij;
315 >            //Xirr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;            
316 >        Xirr += -U[i] * Cij * U[j];
317 >      }
318      }
319 <
320 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
319 >    
320 >    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
321      Xitt *= convertConstant;
322      Xitr *= convertConstant;
323      Xirr *= convertConstant;
324 <
325 <    double kt = OOPSEConstant::kB * temperature;
326 <
324 >    
325 >    RealType kt = OOPSEConstant::kB * temperature;
326 >    
327      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
328      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
329      Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
330 <
330 >    
331      const static Mat3x3d zeroMat(0.0);
332      
333      Mat3x3d XittInv(0.0);
# Line 400 | Line 381 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
381      Ddrr = Dorr;
382      Ddtr = Dotr + Dorr * Uod;
383  
384 <    SquareMatrix<double, 6> Dd;
384 >    SquareMatrix<RealType, 6> Dd;
385      Dd.setSubMatrix(0, 0, Ddtt);
386      Dd.setSubMatrix(0, 3, Ddtr.transpose());
387      Dd.setSubMatrix(3, 0, Ddtr);
388      Dd.setSubMatrix(3, 3, Ddrr);    
389 <    SquareMatrix<double, 6> Xid;
389 >    SquareMatrix<RealType, 6> Xid;
390      Ddtt *= kt;
391      Ddtr *=kt;
392      Ddrr *= kt;
# Line 457 | Line 438 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
438      std::cout << Xidrr << std::endl;
439  
440      return true;
441 <      
442 < }
441 >    
442 >  }
443  
444 < /*
464 < void ApproximationModel::writeBeads(std::ostream& os) {
444 >  void ApproximationModel::writeBeads(std::ostream& os) {
445      std::vector<BeadParam>::iterator iter;
446      os << beads_.size() << std::endl;
447      os << "Generated by Hydro" << std::endl;
448      for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
449 <        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
449 >      os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
450      }
451 <
451 >    
452 >  }    
453   }
473 – */
474 –
475 –
476 – }

Diff Legend

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