40 |
|
*/ |
41 |
|
|
42 |
|
#include "applications/hydrodynamics/RoughShell.hpp" |
43 |
< |
|
43 |
> |
#include "applications/hydrodynamics/ShapeBuilder.hpp" |
44 |
> |
#include "brains/SimInfo.hpp" |
45 |
|
namespace oopse { |
46 |
< |
|
47 |
< |
RoughShell::RoughShell(StuntDouble* sd, const DynamicProperty& extraParams) : HydrodynamicsModel(sd, extraParams), sdShape_(sd){ |
48 |
< |
DynamicProperty::const_iterator iter = extraParams.find("Sigma"); |
49 |
< |
if (iter != extraParams.end()) { |
50 |
< |
boost::any param = iter->second; |
51 |
< |
sigma_ = boost::any_cast<double>(param); |
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 |
– |
} |