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

Comparing trunk/src/applications/hydrodynamics/AnalyticalModel.cpp (file contents):
Revision 909 by tim, Tue Mar 21 00:26:55 2006 UTC vs.
Revision 1208 by xsun, Wed Jan 16 20:19:28 2008 UTC

# Line 39 | Line 39
39   * such damages.
40   */
41   #include "applications/hydrodynamics/AnalyticalModel.hpp"
42 < #include "applications/hydrodynamics/Spheric.hpp"
43 < #include "applications/hydrodynamics/Ellipsoid.hpp"
42 > #include "hydrodynamics/Sphere.hpp"
43 > #include "hydrodynamics/Ellipsoid.hpp"
44   #include "applications/hydrodynamics/CompositeShape.hpp"
45   #include "math/LU.hpp"
46   namespace oopse {
47 < bool AnalyticalModel::calcHydroProps(Spheric* spheric, double viscosity, double temperature) {
48 <
49 <    double radius = spheric->getRadius();
50 <    HydroProps props;
51 <    props.center =V3Zero;
52 <    double Xitt  = 6.0 * NumericConstant::PI * viscosity * radius;
53 <    double Xirr = 8.0 * NumericConstant::PI * viscosity * radius * radius * radius;
54 <    props.Xi(0, 0) = Xitt;
55 <    props.Xi(1, 1) = Xitt;
56 <    props.Xi(2, 2) = Xitt;
57 <    props.Xi(3, 3) = Xirr;
58 <    props.Xi(4, 4) = Xirr;
59 <    props.Xi(5, 5) = Xirr;
47 >  
48 >  bool AnalyticalModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
49      
50 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
51 <    props.Xi *= convertConstant;
52 <    Mat6x6d XiCopy = props.Xi;
53 <    invertMatrix(XiCopy, props.D);
54 <    double kt = OOPSEConstant::kB * temperature;
55 <    props.D *= kt;
56 <    props.Xi *= OOPSEConstant::kb * temperature;
57 <
58 <    setCR(props);
59 <    setCD(props);
60 <
61 <    return true;
50 >    HydroProp* props;
51 >    Sphere* sphere = dynamic_cast<Sphere*>(shape);        
52 >    if (sphere != NULL) {
53 >      props = sphere->getHydroProp(viscosity, temperature);
54 >      setCR(props);
55 >      setCD(props);
56 >      return true;
57 >    } else {
58 >      Ellipsoid* ellipsoid = dynamic_cast<Ellipsoid*>(shape);        
59 >      if (ellipsoid != NULL) {
60 >        props = ellipsoid->getHydroProp(viscosity, temperature);
61 >        setCR(props);
62 >        setCD(props);
63 >        return true;
64 >      } else {
65 >        CompositeShape* composite = dynamic_cast<CompositeShape*>(shape);
66 >        if (composite != NULL) {
67 > //           props = composite->getHydroProp(viscosity, temperature);
68 > //           setCR(props);
69 > //           setCD(props);
70 > //           return true;
71 >          return false;
72 >        } else {
73 >          sprintf( painCave.errMsg,
74 >                   "Could not figure out what kind of shape this is!\n");
75 >          painCave.severity = OOPSE_ERROR;
76 >          painCave.isFatal = 1;
77 >          simError();    
78 >          return false;
79 >        }
80 >      }
81 >    }          
82 >  }
83      
84 +  void AnalyticalModel::writeBeads(std::ostream& os) {
85 +    os << "1\n";
86 +    os << "Generated by Hydro\n";
87 +    Vector3d pos = sd_->getPos();
88 +    os << sd_->getType() << "\t" << pos[0] << "\t" << pos[1] << "\t" << pos[2] << std::endl;
89 +  }  
90   }
75 –
76 – /**
77 – * Reference:
78 – * (2) F. Perrin , J. Phys. Radium, [7] 5, 497-511, 1934
79 – * (3) F. Perrin, J. Phys. Radium, [7] 7, 1-11, 1936
80 – */        
81 – bool AnalyticalModel::calcHydroProps(Ellipsoid* ellipsoid, double viscosity, double temperature) {
82 –
83 –    double rMajor = ellipsoid->getRMajor();
84 –    double rMinor = ellipsoid->getRMinor();
85 –
86 –    double a = rMinor;
87 –    double b = rMajor;
88 –    double a2 = a * a;
89 –    double b2 = b* b;
90 –    
91 –    double p = a /b;
92 –    double S;
93 –    if (p > 1.0) { //prolate
94 –        S = 2.0/sqrt(a2 - b2) * log((a + sqrt(a2-b2))/b);
95 –    } else { //oblate
96 –        S = 2.0/sqrt(b2 - a2) * atan(sqrt(b2-a2)/a);
97 –    }
98 –
99 –    double P = 1.0/(a2 - b2) * (S - 2.0/a);
100 –    double Q = 0.5/(a2-b2) * (2.0*a/b2 - S);
101 –
102 –    double transMinor = 16.0 * NumericConstant::PI * viscosity * (a2 - b2) /((2.0*a2-b2)*S -2.0*a);
103 –    double transMajor = 32.0 * NumericConstant::PI * viscosity * (a2 - b2) /((2.0*a2-3.0*b2)*S +2.0*a);
104 –    double rotMinor = 32.0/3.0 * NumericConstant::PI * viscosity *(a2 - b2) * b2 /(2.0*a -b2*S);
105 –    double rotMajor = 32.0/3.0 * NumericConstant::PI * viscosity *(a2*a2 - b2*b2)/((2.0*a2-b2)*S-2.0*a);
106 –    
107 –        
108 –    HydroProps props;
109 –
110 –    props.Xi(0,0) = transMajor;
111 –    props.Xi(1,1) = transMajor;
112 –    props.Xi(2,2) = transMinor;
113 –    props.Xi(3,3) = rotMajor;
114 –    props.Xi(4,4) = rotMajor;
115 –    props.Xi(5,5) = rotMinor;
116 –    
117 –    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
118 –    props.Xi *= convertConstant;    
119 –    
120 –    Mat6x6d XiCopy = props.Xi;
121 –    invertMatrix(XiCopy, props.D);
122 –    double kt = OOPSEConstant::kB * temperature;
123 –    props.D *= kt;
124 –    props.Xi *= OOPSEConstant::kb * temperature;
125 –    
126 –    setCR(props);
127 –    setCD(props);
128 –
129 –    return true;
130 – }
131 –
132 – bool AnalyticalModel::calcHydroProps(CompositeShape* compositexShape, double viscosity, double temperature) {
133 –    return false;
134 – }
135 –        
136 –
137 –
138 – }

Diff Legend

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