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

Comparing trunk/src/applications/hydrodynamics/RoughShell.cpp (file contents):
Revision 898 by tim, Wed Mar 15 15:51:44 2006 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42   #include "applications/hydrodynamics/RoughShell.hpp"
43 <
44 < namespace oopse {
45 <
46 < RoughShell::RoughShell(StuntDouble* sd, const DynamicProperty& extraParams) : HydrodynamicsModel(sd, extraParams), sdShape_(sd){
47 <    DynamicProperty::const_iterator iter = extraParams.find("Sigma");
48 <    if (iter != extraParams.end()) {
49 <        boost::any param = iter->second;
50 <        sigma_ = boost::any_cast<double>(param);
43 > #include "applications/hydrodynamics/ShapeBuilder.hpp"
44 > #include "brains/SimInfo.hpp"
45 > namespace OpenMD {
46 >  
47 >  RoughShell::RoughShell(StuntDouble* sd, SimInfo* info) : ApproximationModel(sd, info){
48 >    shape_=ShapeBuilder::createShape(sd);
49 >    Globals* simParams = info->getSimParams();
50 >    if (simParams->haveBeadSize()) {
51 >      sigma_ = simParams->getBeadSize();
52      }else {
53 <        std::cout << "RoughShell Model Error\n" ;
53 >      
54      }
55 < }
56 <
57 < struct BeadLattice {
55 >  }
56 >  
57 >  struct BeadLattice {
58      Vector3d origin;
59 <    double radius;
59 >    RealType radius;
60      bool interior;
61 < };
62 <
63 < struct ExteriorFunctor : public std::unary_function<BeadLattice, bool>{
64 <
61 >  };
62 >  
63 >  struct ExteriorFunctor : public std::unary_function<BeadLattice, bool>{
64 >    
65      bool operator() (const BeadLattice& bead) {
66 <        return !bead.interior;
66 >      return !bead.interior;
67      }
68 <
69 < };
70 <
71 < struct InteriorFunctor  : public std::unary_function<BeadLattice, bool>{
72 <
68 >    
69 >  };
70 >  
71 >  struct InteriorFunctor  : public std::unary_function<BeadLattice, bool>{
72 >    
73      bool operator() (const BeadLattice& bead) {
74 <        return bead.interior;
74 >      return bead.interior;
75      }
76 <
77 < };
78 < bool RoughShell::createBeads(std::vector<BeadParam>& beads) {
79 <    std::pair<Vector3d, Vector3d> boxBoundary = sdShape_.getBox();
80 <    double len = boxBoundary.second[0] - boxBoundary.first[0];
81 <    int numLattices = static_cast<int>(len/sigma_) + 1;
76 >    
77 >  };
78 >  bool RoughShell::createBeads(std::vector<BeadParam>& beads) {
79 >    std::pair<Vector3d, Vector3d> boxBoundary = shape_->getBoundingBox();
80 >    RealType firstMin = std::min(std::min(boxBoundary.first[0], boxBoundary.first[1]), boxBoundary.first[2]);
81 >    RealType secondMax = std::max(std::max(boxBoundary.second[0], boxBoundary.second[1]), boxBoundary.second[2]);
82 >    RealType len = secondMax - firstMin;
83 >    int numLattices = static_cast<int>(len/sigma_) + 2;
84      Grid3D<BeadLattice>  grid(numLattices, numLattices, numLattices);
85 <
85 >    
86      //fill beads
87      for (int i = 0; i < numLattices; ++i) {
88 <        for (int j = 0; j < numLattices; ++j) {
89 <            for (int k = 0; k < numLattices; ++k) {
90 <                BeadLattice& currentBead = grid(i, j, k);
91 <                currentBead.origin = Vector3d(i*sigma_ + boxBoundary.first[0], j *sigma_ + boxBoundary.first[1], k*sigma_+ boxBoundary.first[2]);
92 <                currentBead.radius = sigma_;
93 <                currentBead.interior = sdShape_.isInterior(grid(i, j, k).origin);                
91 <            }
88 >      for (int j = 0; j < numLattices; ++j) {
89 >        for (int k = 0; k < numLattices; ++k) {
90 >          BeadLattice& currentBead = grid(i, j, k);
91 >          currentBead.origin = Vector3d((i-1)*sigma_ + boxBoundary.first[0], (j-1) *sigma_ + boxBoundary.first[1], (k-1)*sigma_+ boxBoundary.first[2]);
92 >          currentBead.radius = sigma_ / 2.0;
93 >          currentBead.interior = shape_->isInterior(grid(i, j, k).origin);                
94          }
95 +      }
96      }
97 <
97 >    
98      //remove embedded beads
99      for (int i = 0; i < numLattices; ++i) {
100 <        for (int j = 0; j < numLattices; ++j) {
101 <            for (int k = 0; k < numLattices; ++k) {
102 <                 std::vector<BeadLattice> neighborCells = grid.getAllNeighbors(i, j, k);
103 <                 //if one of its neighbor cells is exterior, current cell is on the surface
104 <
105 <                 std::vector<BeadLattice>::iterator ei = std::find_if(neighborCells.begin(), neighborCells.end(), ExteriorFunctor());                
106 <                 std::vector<BeadLattice>::iterator ii = std::find_if(neighborCells.begin(), neighborCells.end(), InteriorFunctor());                
107 <                
108 <                  if (ei != neighborCells.end() && ii != neighborCells.end()) {
109 <                      BeadParam surfaceBead;
110 <                      surfaceBead.atomName = "Bead";
111 <                      surfaceBead.pos = grid(i, j, k).origin;
112 <                      surfaceBead.radius = grid(i, j, k).radius;
110 <                      beads.push_back(surfaceBead);
111 <                  }
112 <
100 >      for (int j = 0; j < numLattices; ++j) {
101 >        for (int k = 0; k < numLattices; ++k) {
102 >          std::vector<BeadLattice> neighborCells = grid.getAllNeighbors(i, j, k);
103 >          //if one of its neighbor cells is exterior, current cell is on the surface
104 >          
105 >          if (grid(i, j, k).interior){
106 >            
107 >            bool allNeighBorIsInterior = true;
108 >            for (std::vector<BeadLattice>::iterator l = neighborCells.begin(); l != neighborCells.end(); ++l) {
109 >              if (!l->interior) {
110 >                allNeighBorIsInterior = false;
111 >                break;
112 >              }
113              }
114 +            
115 +            if (allNeighBorIsInterior)
116 +              continue;
117 +            
118 +            BeadParam surfaceBead;
119 +            surfaceBead.atomName = "H";
120 +            surfaceBead.pos = grid(i, j, k).origin;
121 +            surfaceBead.radius = grid(i, j, k).radius;
122 +            beads.push_back(surfaceBead);                    
123 +            
124 +          }
125          }
126 +      }
127      }
128 <    
128 >    
129      return true;
130 +  }
131 +  
132   }
119 –
120 – }

Diff Legend

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