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

Comparing trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp (file contents):
Revision 891 by tim, Wed Feb 22 20:35:16 2006 UTC vs.
Revision 981 by gezelter, Mon Jun 5 18:24:45 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"
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 < */
53 < bool HydrodynamicsModel::calcHydrodyanmicsProps(double eta) {
54 <    if (!createBeads(beads_)) {
55 <        std::cout << "can not create beads" << std::endl;
56 <        return false;
57 <    }
58 <    
59 <    int nbeads = beads_.size();
60 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
61 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
62 <    Mat3x3d I;
63 <    for (std::size_t i = 0; i < nbeads; ++i) {
64 <        for (std::size_t j = 0; j < nbeads; ++j) {
65 <            Mat3x3d Tij;
66 <            if (i != j ) {
67 <                Vector3d Rij = beads_[i].pos - beads_[j].pos;
68 <                double rij = Rij.length();
69 <                double rij2 = rij * rij;
70 <                double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;                
71 <                Mat3x3d tmpMat;
72 <                tmpMat = outProduct(beads_[i].pos, beads_[j].pos) / rij2;
73 <                double constant = 8.0 * NumericConstant::PI * eta * rij;
74 <                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
75 <            }else {
76 <                double constant = 1.0 / (6.0 * NumericConstant::PI * eta * beads_[i].radius);
77 <                Tij(0, 0) = constant;
78 <                Tij(1, 1) = constant;
79 <                Tij(2, 2) = constant;
80 <            }
81 <            B.setSubMatrix(i*3, j*3, Tij);
82 <        }
83 <    }
48 >  
49 >  bool HydrodynamicsModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
50 >    return false;
51 >  }
52 >  
53 >  void HydrodynamicsModel::writeHydroProps(std::ostream& os) {
54  
55 <    //invert B Matrix
56 <    invertMatrix(B, C);
55 >    Vector3d center;
56 >    Mat6x6d Xi, D;
57      
58 <    //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
89 <    std::vector<Mat3x3d> U;
90 <    for (int i = 0; i < nbeads; ++i) {
91 <        Mat3x3d currU;
92 <        currU.setupSkewMat(beads_[i].pos);
93 <        U.push_back(currU);
94 <    }
58 >    os << sd_->getType() << "\t";
59      
60 <    //calculate Xi matrix at arbitrary origin O
97 <    Mat3x3d Xitt;
98 <    Mat3x3d Xirr;
99 <    Mat3x3d Xitr;
100 <        
101 <    for (std::size_t i = 0; i < nbeads; ++i) {
102 <        for (std::size_t j = 0; j < nbeads; ++j) {
103 <            Mat3x3d Cij;
104 <            C.getSubMatrix(i*3, j*3, Cij);
105 <            
106 <            Xitt += Cij;
107 <            Xirr += U[i] * Cij;
108 <            Xitr += U[i] * Cij * U[j];            
109 <        }
110 <    }
60 >    //center of resistance
61  
62 <    //invert Xi to get Diffusion Tensor at arbitrary origin O
113 <    RectMatrix<double, 6, 6> Xi;    
114 <    RectMatrix<double, 6, 6> Do;
115 <    Xi.setSubMatrix(0, 0, Xitt);
116 <    Xi.setSubMatrix(0, 3, Xitr.transpose());
117 <    Xi.setSubMatrix(3, 0, Xitr);
118 <    Xi.setSubMatrix(3, 3, Xitt);
119 <    invertMatrix(Xi, Do);
62 >    center = cr_->getCOR();
63  
64 <    Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
65 <    Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
66 <    Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
67 <    Do.getSubMatrix(0, 0 , Dott);
125 <    Do.getSubMatrix(3, 0, Dotr);
126 <    Do.getSubMatrix(3, 3, Dorr);
64 >    os << center[0] <<  "\t" << center[1] <<  "\t" << center[2] <<  "\t";
65 >    
66 >    //resistance tensor at center of resistance
67 >    //translation
68  
69 <    //calculate center of diffusion
129 <    Mat3x3d tmpMat;
130 <    tmpMat(0, 0) = Dorr(1, 1) + Dorr(2, 2);
131 <    tmpMat(0, 1) = - Dorr(0, 1);
132 <    tmpMat(0, 2) = -Dorr(0, 2);
133 <    tmpMat(1, 0) = -Dorr(0, 1);
134 <    tmpMat(1, 1) = Dorr(0, 0)  + Dorr(2, 2);
135 <    tmpMat(1, 2) = -Dorr(1, 2);
136 <    tmpMat(2, 0) = -Dorr(0, 2);
137 <    tmpMat(2, 1) = -Dorr(1, 2);
138 <    tmpMat(2, 2) = Dorr(1, 1) + Dorr(0, 0);
69 >    Xi = cr_->getXi();
70  
71 <    Vector3d tmpVec;
72 <    tmpVec[0] = Dotr(1, 2) - Dotr(2, 1);
73 <    tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
74 <    tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
75 <        
76 <    Vector3d rod = tmpMat.inverse() * tmpVec;
71 >    os << Xi(0, 0) <<  "\t" << Xi(0, 1) <<  "\t" << Xi(0, 2) <<  "\t"
72 >       << Xi(1, 0) <<  "\t" << Xi(1, 1) <<  "\t" << Xi(1, 2) <<  "\t"
73 >       << Xi(2, 0) <<  "\t" << Xi(2, 1) <<  "\t" << Xi(2, 2) <<  "\t";
74 >    
75 >    //rotation-translation
76 >    os << Xi(0, 3) <<  "\t" << Xi(0, 4) <<  "\t" << Xi(0, 5) <<  "\t"
77 >       << Xi(1, 3) <<  "\t" << Xi(1, 4) <<  "\t" << Xi(1, 5) <<  "\t"
78 >       << Xi(2, 3) <<  "\t" << Xi(2, 4) <<  "\t" << Xi(2, 5) <<  "\t";
79 >    
80 >    //translation-rotation
81 >    os << Xi(3, 0) <<  "\t" << Xi(3, 1) <<  "\t" << Xi(3, 2) <<  "\t"
82 >       << Xi(4, 0) <<  "\t" << Xi(4, 1) <<  "\t" << Xi(4, 2) <<  "\t"
83 >       << Xi(5, 0) <<  "\t" << Xi(5, 1) <<  "\t" << Xi(5, 2) <<  "\t";
84 >    
85 >    //rotation
86 >    os << Xi(3, 3) <<  "\t" << Xi(3, 4) <<  "\t" << Xi(3, 5) <<  "\t"
87 >       << Xi(4, 3) <<  "\t" << Xi(4, 4) <<  "\t" << Xi(4, 5) <<  "\t"
88 >       << Xi(5, 3) <<  "\t" << Xi(5, 4) <<  "\t" << Xi(5, 5) <<  "\t";
89 >    
90 >    
91 >    //diffusion tensor at center of resistance
92 >    //translation
93  
94 <    //calculate Diffusion Tensor at center of diffusion
95 <    Mat3x3d Uod;
96 <    Uod.setupSkewMat(rod);
94 >    D = cr_->getD();
95 >
96 >    os << D(0, 0) <<  "\t" << D(0, 1) <<  "\t" << D(0, 2) <<  "\t"
97 >       << D(1, 0) <<  "\t" << D(1, 1) <<  "\t" << D(1, 2) <<  "\t"
98 >       << D(2, 0) <<  "\t" << D(2, 1) <<  "\t" << D(2, 2) <<  "\t";
99      
100 <    Mat3x3d Ddtt; //translational diffusion tensor at diffusion center
101 <    Mat3x3d Ddtr; //rotational diffusion tensor at diffusion center
102 <    Mat3x3d Ddrr; //translation-rotation couplingl diffusion tensor at diffusion tensor
100 >    //rotation-translation
101 >    os << D(0, 3) <<  "\t" << D(0, 4) <<  "\t" << D(0, 5) <<  "\t"
102 >       << D(1, 3) <<  "\t" << D(1, 4) <<  "\t" << D(1, 5) <<  "\t"
103 >       << D(2, 3) <<  "\t" << D(2, 4) <<  "\t" << D(2, 5) <<  "\t";
104      
105 <    Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr;
106 <    Ddrr = Dorr;
107 <    Ddtr = Dotr + Dorr * Uod;
105 >    //translation-rotation
106 >    os << D(3, 0) <<  "\t" << D(3, 1) <<  "\t" << D(3, 2) <<  "\t"
107 >       << D(4, 0) <<  "\t" << D(4, 1) <<  "\t" << D(4, 2) <<  "\t"
108 >       << D(5, 0) <<  "\t" << D(5, 1) <<  "\t" << D(5, 2) <<  "\t";
109      
110 <    props_.diffCenter = rod;
111 <    props_.transDiff = Ddtt;
112 <    props_.transRotDiff = Ddtr;
113 <    props_.rotDiff = Ddrr;
110 >    //rotation
111 >    os << D(3, 3) <<  "\t" << D(3, 4) <<  "\t" << D(3, 5) <<  "\t"
112 >       << D(4, 3) <<  "\t" << D(4, 4) <<  "\t" << D(4, 5) <<  "\t"
113 >       << D(5, 3) <<  "\t" << D(5, 4) <<  "\t" << D(5, 5) <<  "\t";
114 >    
115 >    //---------------------------------------------------------------------
116 >    
117 >    //center of diffusion
118  
119 <    return true;    
165 < }
119 >    center = cd_->getCOR();
120  
121 < void HydrodynamicsModel::writeBeads(std::ostream& os) {
121 >    os << center[0] <<  "\t" << center[1] <<  "\t" << center[2] <<  "\t";
122 >    
123 >    //resistance tensor at center of diffusion
124 >    //translation
125  
126 < }
126 >    Xi = cd_->getXi();
127  
128 < void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
128 >    os << Xi(0, 0) <<  "\t" << Xi(0, 1) <<  "\t" << Xi(0, 2) <<  "\t"
129 >       << Xi(1, 0) <<  "\t" << Xi(1, 1) <<  "\t" << Xi(1, 2) <<  "\t"
130 >       << Xi(2, 0) <<  "\t" << Xi(2, 1) <<  "\t" << Xi(2, 2) <<  "\t";
131 >    
132 >    //rotation-translation
133 >    os << Xi(0, 3) <<  "\t" << Xi(0, 4) <<  "\t" << Xi(0, 5) <<  "\t"
134 >       << Xi(1, 3) <<  "\t" << Xi(1, 4) <<  "\t" << Xi(1, 5) <<  "\t"
135 >       << Xi(2, 3) <<  "\t" << Xi(2, 4) <<  "\t" << Xi(2, 5) <<  "\t";
136 >    
137 >    //translation-rotation
138 >    os << Xi(3, 0) <<  "\t" << Xi(3, 1) <<  "\t" << Xi(3, 2) <<  "\t"
139 >       << Xi(4, 0) <<  "\t" << Xi(4, 1) <<  "\t" << Xi(4, 2) <<  "\t"
140 >       << Xi(5, 0) <<  "\t" << Xi(5, 1) <<  "\t" << Xi(5, 2) <<  "\t";
141 >    
142 >    //rotation
143 >    os << Xi(3, 3) <<  "\t" << Xi(3, 4) <<  "\t" << Xi(3, 5) <<  "\t"
144 >       << Xi(4, 3) <<  "\t" << Xi(4, 4) <<  "\t" << Xi(4, 5) <<  "\t"
145 >       << Xi(5, 3) <<  "\t" << Xi(5, 4) <<  "\t" << Xi(5, 5) <<  "\t";
146 >    
147 >    
148 >    //diffusion tensor at center of diffusion
149 >    //translation
150  
151 < }
151 >    D = cd_->getD();
152  
153 +    os << D(0, 0) <<  "\t" << D(0, 1) <<  "\t" << D(0, 2) <<  "\t"
154 +       << D(1, 0) <<  "\t" << D(1, 1) <<  "\t" << D(1, 2) <<  "\t"
155 +       << D(2, 0) <<  "\t" << D(2, 1) <<  "\t" << D(2, 2) <<  "\t";
156 +    
157 +    //rotation-translation
158 +    os << D(0, 3) <<  "\t" << D(0, 4) <<  "\t" << D(0, 5) <<  "\t"
159 +       << D(1, 3) <<  "\t" << D(1, 4) <<  "\t" << D(1, 5) <<  "\t"
160 +       << D(2, 3) <<  "\t" << D(2, 4) <<  "\t" << D(2, 5) <<  "\t";
161 +    
162 +    //translation-rotation
163 +    os << D(3, 0) <<  "\t" << D(3, 1) <<  "\t" << D(3, 2) <<  "\t"
164 +       << D(4, 0) <<  "\t" << D(4, 1) <<  "\t" << D(4, 2) <<  "\t"
165 +       << D(5, 0) <<  "\t" << D(5, 1) <<  "\t" << D(5, 2) <<  "\t";
166 +    
167 +    //rotation
168 +    os << D(3, 3) <<  "\t" << D(3, 4) <<  "\t" << D(3, 5) <<  "\t"
169 +       << D(4, 3) <<  "\t" << D(4, 4) <<  "\t" << D(4, 5) <<  "\t"
170 +       << D(5, 3) <<  "\t" << D(5, 4) <<  "\t" << D(5, 5) <<  "\n";    
171 +    
172 +  }
173 +  
174   }

Diff Legend

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