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 2597 by tim, Thu Feb 23 23:16:43 2006 UTC vs.
Revision 2768 by gezelter, Wed May 24 18:31:12 2006 UTC

# Line 38 | Line 38
38   * University of Notre Dame has been advised of the possibility of
39   * such damages.
40   */
41 + #include "applications/hydrodynamics/HydrodynamicsModel.hpp"
42 + #include "hydrodynamics/Shape.hpp"
43 + #include "hydrodynamics/Sphere.hpp"
44 + #include "hydrodynamics/Ellipsoid.hpp"
45 + #include "applications/hydrodynamics/CompositeShape.hpp"
46  
42 #include "applications/hydrodynamics/HydrodynamicsModel.hpp"
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:
50 < * Beatriz Carrasco and Jose Gracia de la Torre, Hydrodynamic Properties of Rigid Particles:
51 < * Comparison of Different Modeling and Computational Procedures.
52 < * Biophysical Journal, 75(6), 3044, 1999
53 < */
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;
79 <    }
48 >  
49 >  bool HydrodynamicsModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
50 >    return false;
51 >  }
52 >  
53 >  void HydrodynamicsModel::writeHydroProps(std::ostream& os) {
54      
81    int nbeads = beads_.size();
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;
55      
56 <    for (std::size_t i = 0; i < nbeads; ++i) {
90 <        for (std::size_t j = 0; j < nbeads; ++j) {
91 <            Mat3x3d Tij;
92 <            if (i != j ) {
93 <                Vector3d Rij = beads_[i].pos - beads_[j].pos;
94 <                double rij = Rij.length();
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;
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 * viscosity_ * beads_[i].radius);
103 <                Tij(0, 0) = constant;
104 <                Tij(1, 1) = constant;
105 <                Tij(2, 2) = constant;
106 <            }
107 <            B.setSubMatrix(i*3, j*3, Tij);
108 <        }
109 <    }
110 <
111 <    //invert B Matrix
112 <    invertMatrix(B, C);
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) {
116 <        Mat3x3d currU;
117 <        currU.setupSkewMat(beads_[i].pos);
118 <        U.push_back(currU);
119 <    }
56 >    os << sd_->getType() << "\t";
57      
58 <    //calculate Xi matrix at arbitrary origin O
59 <    Mat3x3d Xitt;
123 <    Mat3x3d Xirr;
124 <    Mat3x3d Xitr;
125 <        
126 <    for (std::size_t i = 0; i < nbeads; ++i) {
127 <        for (std::size_t j = 0; j < nbeads; ++j) {
128 <            Mat3x3d Cij;
129 <            C.getSubMatrix(i*3, j*3, Cij);
130 <            
131 <            Xitt += Cij;
132 <            Xirr += U[i] * Cij;
133 <            Xitr += U[i] * Cij * U[j];            
134 <        }
135 <    }
136 <
137 <    //invert Xi to get Diffusion Tensor at arbitrary origin O
138 <    RectMatrix<double, 6, 6> Xi;    
139 <    RectMatrix<double, 6, 6> Do;
140 <    Xi.setSubMatrix(0, 0, Xitt);
141 <    Xi.setSubMatrix(0, 3, Xitr.transpose());
142 <    Xi.setSubMatrix(3, 0, Xitr);
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
152 <
153 <    Mat3x3d XirrInv(0.0);
154 <    Mat3x3d XirrCopy;
155 <    XirrCopy = Xirr;
58 >    //center of resistance
59 >    os << cr_.center[0] <<  "\t" << cr_.center[1] <<  "\t" << cr_.center[2] <<  "\t";
60      
61 <    Mat3x3d XittInv(0.0);
62 <    Mat3x3d XittCopy;
63 <    XittCopy = Xitt;
64 <    invertMatrix(XittCopy, XittInv);
65 <
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();
61 >    //resistance tensor at center of resistance
62 >    //translation
63 >    os << cr_.Xi(0, 0) <<  "\t" << cr_.Xi(0, 1) <<  "\t" << cr_.Xi(0, 2) <<  "\t"
64 >       << cr_.Xi(1, 0) <<  "\t" << cr_.Xi(1, 1) <<  "\t" << cr_.Xi(1, 2) <<  "\t"
65 >       << cr_.Xi(2, 0) <<  "\t" << cr_.Xi(2, 1) <<  "\t" << cr_.Xi(2, 2) <<  "\t";
66      
67 <    if(!invertMatrix(tmp, tmpInv)) {
68 <        tmpInv = zeroMat;
69 <    }
70 <    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 <    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);
200 <
201 <    if(!invertMatrix(tmp, tmpInv)) {
202 <        tmpInv = zeroMat;
203 <    }
67 >    //rotation-translation
68 >    os << cr_.Xi(0, 3) <<  "\t" << cr_.Xi(0, 4) <<  "\t" << cr_.Xi(0, 5) <<  "\t"
69 >       << cr_.Xi(1, 3) <<  "\t" << cr_.Xi(1, 4) <<  "\t" << cr_.Xi(1, 5) <<  "\t"
70 >       << cr_.Xi(2, 3) <<  "\t" << cr_.Xi(2, 4) <<  "\t" << cr_.Xi(2, 5) <<  "\t";
71      
72 <    Vector3d rod = tmpInv * tmpVec;
73 <
74 <    //calculate Diffusion Tensor at center of diffusion
75 <    Mat3x3d Uod;
209 <    Uod.setupSkewMat(rod);
72 >    //translation-rotation
73 >    os << cr_.Xi(3, 0) <<  "\t" << cr_.Xi(3, 1) <<  "\t" << cr_.Xi(3, 2) <<  "\t"
74 >       << cr_.Xi(4, 0) <<  "\t" << cr_.Xi(4, 1) <<  "\t" << cr_.Xi(4, 2) <<  "\t"
75 >       << cr_.Xi(5, 0) <<  "\t" << cr_.Xi(5, 1) <<  "\t" << cr_.Xi(5, 2) <<  "\t";
76      
77 <    Mat3x3d Ddtt; //translational diffusion tensor at diffusion center
78 <    Mat3x3d Ddtr; //rotational diffusion tensor at diffusion center
79 <    Mat3x3d Ddrr; //translation-rotation couplingl diffusion tensor at diffusion tensor
77 >    //rotation
78 >    os << cr_.Xi(3, 3) <<  "\t" << cr_.Xi(3, 4) <<  "\t" << cr_.Xi(3, 5) <<  "\t"
79 >       << cr_.Xi(4, 3) <<  "\t" << cr_.Xi(4, 4) <<  "\t" << cr_.Xi(4, 5) <<  "\t"
80 >       << cr_.Xi(5, 3) <<  "\t" << cr_.Xi(5, 4) <<  "\t" << cr_.Xi(5, 5) <<  "\t";
81      
215    Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr;
216    Ddrr = Dorr;
217    Ddtr = Dotr + Dorr * Uod;
82      
83 <    props_.diffCenter = rod;
84 <    props_.transDiff = Ddtt;
85 <    props_.transRotDiff = Ddtr;
86 <    props_.rotDiff = Ddrr;
87 <
224 <    return true;    
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;
83 >    //diffusion tensor at center of resistance
84 >    //translation
85 >    os << cr_.D(0, 0) <<  "\t" << cr_.D(0, 1) <<  "\t" << cr_.D(0, 2) <<  "\t"
86 >       << cr_.D(1, 0) <<  "\t" << cr_.D(1, 1) <<  "\t" << cr_.D(1, 2) <<  "\t"
87 >       << cr_.D(2, 0) <<  "\t" << cr_.D(2, 1) <<  "\t" << cr_.D(2, 2) <<  "\t";
88      
89 <    /*
90 <    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\n"
91 <
92 <    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";
89 >    //rotation-translation
90 >    os << cr_.D(0, 3) <<  "\t" << cr_.D(0, 4) <<  "\t" << cr_.D(0, 5) <<  "\t"
91 >       << cr_.D(1, 3) <<  "\t" << cr_.D(1, 4) <<  "\t" << cr_.D(1, 5) <<  "\t"
92 >       << cr_.D(2, 3) <<  "\t" << cr_.D(2, 4) <<  "\t" << cr_.D(2, 5) <<  "\t";
93      
94 <    os << props_.transRotDiff(0, 0) << "\t" << props_.transRotDiff(0, 1) << "\t" << props_.transRotDiff(0, 2) << "\t"
95 <        << props_.transRotDiff(1, 0) << "\t" << props_.transRotDiff(1, 1) << "\t" << props_.transRotDiff(1, 2) << "\t"
96 <        << props_.transRotDiff(2, 0) << "\t" << props_.transRotDiff(2, 1) << "\t" << props_.transRotDiff(2, 2) << "\t"
97 <
98 <    os << props_.rotDiff(0, 0) << "\t" << props_.rotDiff(0, 1) << "\t" << props_.rotDiff(0, 2) << "\t"
99 <        << props_.rotDiff(1, 0) << "\t" << props_.rotDiff(1, 1) << "\t" << props_.rotDiff(1, 2) << "\t"
100 <        << props_.rotDiff(2, 0) << "\t" << props_.rotDiff(2, 1) << "\t" << props_.rotDiff(2, 2) << ";"
101 <        << std::endl;
102 <    */
94 >    //translation-rotation
95 >    os << cr_.D(3, 0) <<  "\t" << cr_.D(3, 1) <<  "\t" << cr_.D(3, 2) <<  "\t"
96 >       << cr_.D(4, 0) <<  "\t" << cr_.D(4, 1) <<  "\t" << cr_.D(4, 2) <<  "\t"
97 >       << cr_.D(5, 0) <<  "\t" << cr_.D(5, 1) <<  "\t" << cr_.D(5, 2) <<  "\t";
98 >    
99 >    //rotation
100 >    os << cr_.D(3, 3) <<  "\t" << cr_.D(3, 4) <<  "\t" << cr_.D(3, 5) <<  "\t"
101 >       << cr_.D(4, 3) <<  "\t" << cr_.D(4, 4) <<  "\t" << cr_.D(4, 5) <<  "\t"
102 >       << cr_.D(5, 3) <<  "\t" << cr_.D(5, 4) <<  "\t" << cr_.D(5, 5) <<  "\t";
103 >    
104 >    //---------------------------------------------------------------------
105 >    
106 >    //center of diffusion
107 >    os << cd_.center[0] <<  "\t" << cd_.center[1] <<  "\t" << cd_.center[2] <<  "\t";
108 >    
109 >    //resistance tensor at center of diffusion
110 >    //translation
111 >    os << cd_.Xi(0, 0) <<  "\t" << cd_.Xi(0, 1) <<  "\t" << cd_.Xi(0, 2) <<  "\t"
112 >       << cd_.Xi(1, 0) <<  "\t" << cd_.Xi(1, 1) <<  "\t" << cd_.Xi(1, 2) <<  "\t"
113 >       << cd_.Xi(2, 0) <<  "\t" << cd_.Xi(2, 1) <<  "\t" << cd_.Xi(2, 2) <<  "\t";
114 >    
115 >    //rotation-translation
116 >    os << cd_.Xi(0, 3) <<  "\t" << cd_.Xi(0, 4) <<  "\t" << cd_.Xi(0, 5) <<  "\t"
117 >       << cd_.Xi(1, 3) <<  "\t" << cd_.Xi(1, 4) <<  "\t" << cd_.Xi(1, 5) <<  "\t"
118 >       << cd_.Xi(2, 3) <<  "\t" << cd_.Xi(2, 4) <<  "\t" << cd_.Xi(2, 5) <<  "\t";
119 >    
120 >    //translation-rotation
121 >    os << cd_.Xi(3, 0) <<  "\t" << cd_.Xi(3, 1) <<  "\t" << cd_.Xi(3, 2) <<  "\t"
122 >       << cd_.Xi(4, 0) <<  "\t" << cd_.Xi(4, 1) <<  "\t" << cd_.Xi(4, 2) <<  "\t"
123 >       << cd_.Xi(5, 0) <<  "\t" << cd_.Xi(5, 1) <<  "\t" << cd_.Xi(5, 2) <<  "\t";
124 >    
125 >    //rotation
126 >    os << cd_.Xi(3, 3) <<  "\t" << cd_.Xi(3, 4) <<  "\t" << cd_.Xi(3, 5) <<  "\t"
127 >       << cd_.Xi(4, 3) <<  "\t" << cd_.Xi(4, 4) <<  "\t" << cd_.Xi(4, 5) <<  "\t"
128 >       << cd_.Xi(5, 3) <<  "\t" << cd_.Xi(5, 4) <<  "\t" << cd_.Xi(5, 5) <<  "\t";
129 >    
130 >    
131 >    //diffusion tensor at center of diffusion
132 >    //translation
133 >    os << cd_.D(0, 0) <<  "\t" << cd_.D(0, 1) <<  "\t" << cd_.D(0, 2) <<  "\t"
134 >       << cd_.D(1, 0) <<  "\t" << cd_.D(1, 1) <<  "\t" << cd_.D(1, 2) <<  "\t"
135 >       << cd_.D(2, 0) <<  "\t" << cd_.D(2, 1) <<  "\t" << cd_.D(2, 2) <<  "\t";
136 >    
137 >    //rotation-translation
138 >    os << cd_.D(0, 3) <<  "\t" << cd_.D(0, 4) <<  "\t" << cd_.D(0, 5) <<  "\t"
139 >       << cd_.D(1, 3) <<  "\t" << cd_.D(1, 4) <<  "\t" << cd_.D(1, 5) <<  "\t"
140 >       << cd_.D(2, 3) <<  "\t" << cd_.D(2, 4) <<  "\t" << cd_.D(2, 5) <<  "\t";
141 >    
142 >    //translation-rotation
143 >    os << cd_.D(3, 0) <<  "\t" << cd_.D(3, 1) <<  "\t" << cd_.D(3, 2) <<  "\t"
144 >       << cd_.D(4, 0) <<  "\t" << cd_.D(4, 1) <<  "\t" << cd_.D(4, 2) <<  "\t"
145 >       << cd_.D(5, 0) <<  "\t" << cd_.D(5, 1) <<  "\t" << cd_.D(5, 2) <<  "\t";
146 >    
147 >    //rotation
148 >    os << cd_.D(3, 3) <<  "\t" << cd_.D(3, 4) <<  "\t" << cd_.D(3, 5) <<  "\t"
149 >       << cd_.D(4, 3) <<  "\t" << cd_.D(4, 4) <<  "\t" << cd_.D(4, 5) <<  "\t"
150 >       << cd_.D(5, 3) <<  "\t" << cd_.D(5, 4) <<  "\t" << cd_.D(5, 5) <<  "\n";    
151 >    
152 >  }
153 >  
154   }
272
273 }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines