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 977 by tim, Thu May 25 16:27:27 2006 UTC vs.
Revision 1177 by xsun, Tue Aug 14 17:40:33 2007 UTC

# Line 72 | Line 72 | namespace oopse {
72    bool ApproximationModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
73      
74      bool ret = true;
75 <    HydroProps cr;
76 <    HydroProps cd;
75 >    HydroProp* cr = new HydroProp();
76 >    HydroProp* cd = new HydroProp();
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 –    
81      return true;    
82    }
83    
84 <  bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) {
84 >  bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProp* cr) {
85      
86      int nbeads = beads.size();
87      DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
# Line 103 | Line 102 | namespace oopse {
102                Mat3x3d tmpMat;
103                tmpMat = outProduct(Rij, Rij) / rij2;
104                RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
105 <              Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
105 >              RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0;
106 >              RealType tmp2 = 1.0 - sumSigma2OverRij2;
107 >              Tij = (tmp1 * I + tmp2 * tmpMat ) / constant;
108              }else {
109                RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
110                Tij(0, 0) = constant;
# Line 211 | Line 212 | namespace oopse {
212      Xirtr *= OOPSEConstant::kb * temperature;
213      Xirrr *= OOPSEConstant::kb * temperature;
214      
215 +    Mat6x6d Xi, D;
216  
217 <    cr.center = ror;
218 <    cr.Xi.setSubMatrix(0, 0, Xirtt);
219 <    cr.Xi.setSubMatrix(0, 3, Xirtr);
220 <    cr.Xi.setSubMatrix(3, 0, Xirtr);
221 <    cr.Xi.setSubMatrix(3, 3, Xirrr);
222 <    cr.D.setSubMatrix(0, 0, Drtt);
223 <    cr.D.setSubMatrix(0, 3, Drrt);
224 <    cr.D.setSubMatrix(3, 0, Drtr);
225 <    cr.D.setSubMatrix(3, 3, Drrr);    
217 >    cr->setCOR(ror);
218 >
219 >    Xi.setSubMatrix(0, 0, Xirtt);
220 >    Xi.setSubMatrix(0, 3, Xirtr);
221 >    Xi.setSubMatrix(3, 0, Xirtr);
222 >    Xi.setSubMatrix(3, 3, Xirrr);
223 >
224 >    cr->setXi(Xi);
225 >
226 >    D.setSubMatrix(0, 0, Drtt);
227 >    D.setSubMatrix(0, 3, Drrt);
228 >    D.setSubMatrix(3, 0, Drtr);
229 >    D.setSubMatrix(3, 3, Drrr);    
230 >
231 >    cr->setD(D);
232      
233      std::cout << "-----------------------------------------\n";
234      std::cout << "center of resistance :" << std::endl;
# Line 246 | Line 254 | namespace oopse {
254      return true;
255   }
256    
257 <  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) {
257 >  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProp* cd) {
258      
259      int nbeads = beads.size();
260      DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
# Line 267 | Line 275 | namespace oopse {
275            Mat3x3d tmpMat;
276            tmpMat = outProduct(Rij, Rij) / rij2;
277            RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
278 <          Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
278 >          RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0;
279 >          RealType tmp2 = 1.0 - sumSigma2OverRij2;
280 >          Tij = (tmp1 * I + tmp2 * tmpMat ) / constant;
281          }else {
282            RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
283            Tij(0, 0) = constant;
# Line 394 | Line 404 | namespace oopse {
404      //Xid /= OOPSEConstant::energyConvert;
405      Xid *= OOPSEConstant::kb * temperature;
406  
407 <    cr.center = rod;
398 <    cr.D.setSubMatrix(0, 0, Ddtt);
399 <    cr.D.setSubMatrix(0, 3, Ddtr);
400 <    cr.D.setSubMatrix(3, 0, Ddtr);
401 <    cr.D.setSubMatrix(3, 3, Ddrr);
402 <    cr.Xi = Xid;
407 >    Mat6x6d Xi, D;
408  
409 +    cd->setCOR(rod);
410 +
411 +    cd->setXi(Xid);
412 +
413 +    D.setSubMatrix(0, 0, Ddtt);
414 +    D.setSubMatrix(0, 3, Ddtr);
415 +    D.setSubMatrix(3, 0, Ddtr);
416 +    D.setSubMatrix(3, 3, Ddrr);
417 +
418 +    cd->setD(D);
419 +
420      std::cout << "viscosity = " << viscosity << std::endl;
421      std::cout << "temperature = " << temperature << std::endl;
422      std::cout << "center of diffusion :" << std::endl;

Diff Legend

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