1 | /*
|
---|
2 | * Project: MoleCuilder
|
---|
3 | * Description: creates and alters molecular systems
|
---|
4 | * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
|
---|
5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
---|
6 | */
|
---|
7 |
|
---|
8 | /*
|
---|
9 | * BaseShapes_impl.cpp
|
---|
10 | *
|
---|
11 | * Created on: Jun 18, 2010
|
---|
12 | * Author: crueger
|
---|
13 | */
|
---|
14 |
|
---|
15 | // include config.h
|
---|
16 | #ifdef HAVE_CONFIG_H
|
---|
17 | #include <config.h>
|
---|
18 | #endif
|
---|
19 |
|
---|
20 | #include "CodePatterns/MemDebug.hpp"
|
---|
21 |
|
---|
22 | #include "Shapes/BaseShapes.hpp"
|
---|
23 | #include "Shapes/BaseShapes_impl.hpp"
|
---|
24 | #include "Shapes/ShapeExceptions.hpp"
|
---|
25 | #include "Shapes/ShapeOps.hpp"
|
---|
26 |
|
---|
27 | #include "Helpers/defs.hpp"
|
---|
28 |
|
---|
29 | #include "CodePatterns/Assert.hpp"
|
---|
30 | #include "LinearAlgebra/Vector.hpp"
|
---|
31 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
32 | #include "LinearAlgebra/Line.hpp"
|
---|
33 | #include "LinearAlgebra/Plane.hpp"
|
---|
34 | #include "LinearAlgebra/LineSegment.hpp"
|
---|
35 | #include "LinearAlgebra/LineSegmentSet.hpp"
|
---|
36 |
|
---|
37 | #include <cmath>
|
---|
38 | #include <algorithm>
|
---|
39 |
|
---|
40 | // CYLINDER CODE
|
---|
41 | // ----------------------------------------------------------------------------
|
---|
42 | bool Cylinder_impl::isInside(const Vector &point) const {
|
---|
43 | return (Vector(point[0], point[1], 0.0).NormSquared() < 1.0+MYEPSILON) &&
|
---|
44 | (point[2] > -1.0-MYEPSILON) && (point[2] < 1.0+MYEPSILON);
|
---|
45 | }
|
---|
46 |
|
---|
47 | bool Cylinder_impl::isOnSurface(const Vector &point) const {
|
---|
48 | return fabs(Vector(point[0], point[1], 0.0).NormSquared()-1.0)<MYEPSILON &&
|
---|
49 | (point[2] > -1.0-MYEPSILON) && (point[2] < 1.0+MYEPSILON);
|
---|
50 |
|
---|
51 | }
|
---|
52 |
|
---|
53 | Vector Cylinder_impl::getNormal(const Vector &point) const throw(NotOnSurfaceException) {
|
---|
54 | if(!isOnSurface(point)){
|
---|
55 | throw NotOnSurfaceException() << ShapeVector(&point);
|
---|
56 | }
|
---|
57 |
|
---|
58 | if ((fabs(point[2]-1)<MYEPSILON) || (fabs(point[2])<MYEPSILON))
|
---|
59 | return Vector(0.0, 0.0, point[2]);
|
---|
60 | else
|
---|
61 | return Vector(point[0], point[1], 0.0);
|
---|
62 | }
|
---|
63 |
|
---|
64 | Vector Cylinder_impl::getCenter() const
|
---|
65 | {
|
---|
66 | return Vector(0.0, 0.0, 0.0);
|
---|
67 | }
|
---|
68 |
|
---|
69 | double Cylinder_impl::getRadius() const
|
---|
70 | {
|
---|
71 | return 1.0;
|
---|
72 | }
|
---|
73 |
|
---|
74 | double Cylinder_impl::getVolume() const
|
---|
75 | {
|
---|
76 | return M_PI*2.0; // pi r^2 h
|
---|
77 | }
|
---|
78 |
|
---|
79 | double Cylinder_impl::getSurfaceArea() const
|
---|
80 | {
|
---|
81 | return 2.0*M_PI*2.0; // 2 pi r h
|
---|
82 | }
|
---|
83 |
|
---|
84 | LineSegmentSet Cylinder_impl::getLineIntersections(const Line &line) const {
|
---|
85 | // TODO
|
---|
86 | ASSERT(0, "Cylinder_impl::getLineIntersections - not implemented.");
|
---|
87 | }
|
---|
88 |
|
---|
89 | std::string Cylinder_impl::toString() const
|
---|
90 | {
|
---|
91 | return "Cylinder()";
|
---|
92 | }
|
---|
93 |
|
---|
94 | enum ShapeType Cylinder_impl::getType() const
|
---|
95 | {
|
---|
96 | return CylinderType;
|
---|
97 | }
|
---|
98 |
|
---|
99 | std::vector<Vector> Cylinder_impl::getHomogeneousPointsOnSurface(const size_t N) const {
|
---|
100 | const double dnz = sqrt(N/M_PI);
|
---|
101 | const int nu = round(N/dnz);
|
---|
102 | const int nz = round(dnz);
|
---|
103 |
|
---|
104 | const double dphi = 2.0*M_PI/nu;
|
---|
105 | const double dz = 2.0/nz;
|
---|
106 |
|
---|
107 | std::vector<Vector> result;
|
---|
108 |
|
---|
109 | for(int useg=0; useg<nu; useg++)
|
---|
110 | for(int zseg=0; zseg<nu; zseg++)
|
---|
111 | result.push_back(Vector(cos(useg*dphi), sin(useg*dphi), zseg*dz-1.0));
|
---|
112 |
|
---|
113 | return result;
|
---|
114 | }
|
---|
115 |
|
---|
116 | std::vector<Vector> Cylinder_impl::getHomogeneousPointsInVolume(const size_t N) const {
|
---|
117 | // TODO
|
---|
118 | ASSERT(0, "Cylinder_impl::getHomogeneousPointsInVolume - not implemented.");
|
---|
119 | }
|
---|
120 |
|
---|
121 | Shape Cylinder() {
|
---|
122 | Shape::impl_ptr impl = Shape::impl_ptr(new Cylinder_impl());
|
---|
123 | return Shape(impl);
|
---|
124 | }
|
---|
125 |
|
---|
126 | Shape Cylinder(const Vector ¢er, const double xrot, const double yrot,
|
---|
127 | const double height, const double radius)
|
---|
128 | {
|
---|
129 | RealSpaceMatrix rot;
|
---|
130 | rot.setRotation(xrot, yrot, 0.0);
|
---|
131 |
|
---|
132 | return translate(
|
---|
133 | transform(
|
---|
134 | stretch(
|
---|
135 | Cylinder(),
|
---|
136 | Vector(radius, radius, height*0.5)),
|
---|
137 | rot),
|
---|
138 | center);
|
---|
139 | }
|
---|
140 | // ----------------------------------------------------------------------------
|
---|
141 |
|
---|
142 | bool Sphere_impl::isInside(const Vector &point) const{
|
---|
143 | return point.NormSquared()<=1.;
|
---|
144 | }
|
---|
145 |
|
---|
146 | bool Sphere_impl::isOnSurface(const Vector &point) const{
|
---|
147 | return fabs(point.NormSquared()-1.)<MYEPSILON;
|
---|
148 | }
|
---|
149 |
|
---|
150 | Vector Sphere_impl::getNormal(const Vector &point) const throw(NotOnSurfaceException){
|
---|
151 | if(!isOnSurface(point)){
|
---|
152 | throw NotOnSurfaceException() << ShapeVector(&point);
|
---|
153 | }
|
---|
154 | return point;
|
---|
155 | }
|
---|
156 |
|
---|
157 | Vector Sphere_impl::getCenter() const
|
---|
158 | {
|
---|
159 | return Vector(0.,0.,0.);
|
---|
160 | }
|
---|
161 |
|
---|
162 | double Sphere_impl::getRadius() const
|
---|
163 | {
|
---|
164 | return 1.;
|
---|
165 | }
|
---|
166 |
|
---|
167 | double Sphere_impl::getVolume() const
|
---|
168 | {
|
---|
169 | return (4./3.)*M_PI; // 4/3 pi r^3
|
---|
170 | }
|
---|
171 |
|
---|
172 | double Sphere_impl::getSurfaceArea() const
|
---|
173 | {
|
---|
174 | return 2.*M_PI; // 2 pi r^2
|
---|
175 | }
|
---|
176 |
|
---|
177 |
|
---|
178 | LineSegmentSet Sphere_impl::getLineIntersections(const Line &line) const{
|
---|
179 | LineSegmentSet res(line);
|
---|
180 | std::vector<Vector> intersections = line.getSphereIntersections();
|
---|
181 | if(intersections.size()==2){
|
---|
182 | res.insert(LineSegment(intersections[0],intersections[1]));
|
---|
183 | }
|
---|
184 | return res;
|
---|
185 | }
|
---|
186 |
|
---|
187 | std::string Sphere_impl::toString() const{
|
---|
188 | return "Sphere()";
|
---|
189 | }
|
---|
190 |
|
---|
191 | enum ShapeType Sphere_impl::getType() const
|
---|
192 | {
|
---|
193 | return SphereType;
|
---|
194 | }
|
---|
195 |
|
---|
196 | /**
|
---|
197 | * algorithm taken from http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere
|
---|
198 | * \param N number of points on surface
|
---|
199 | */
|
---|
200 | std::vector<Vector> Sphere_impl::getHomogeneousPointsOnSurface(const size_t N) const
|
---|
201 | {
|
---|
202 | std::vector<Vector> PointsOnSurface;
|
---|
203 | if (true) {
|
---|
204 | // Exactly N points but not symmetric.
|
---|
205 |
|
---|
206 | // This formula is derived by finding a curve on the sphere that spirals down from
|
---|
207 | // the north pole to the south pole keeping a constant distance between consecutive turns.
|
---|
208 | // The curve is then parametrized by arch length and evaluated in constant intervals.
|
---|
209 | double a = sqrt(N) * 2;
|
---|
210 | for (int i=0; i<N; i++){
|
---|
211 | double t0 = ((double)i + 0.5) / (double)N;
|
---|
212 | double t = (sqrt(t0) - sqrt(1.0 - t0) + 1.0) / 2.0 * M_PI;
|
---|
213 | Vector point;
|
---|
214 | point.Zero();
|
---|
215 | point[0] = sin(t) * sin(t * a);
|
---|
216 | point[1] = sin(t) * cos(t * a);
|
---|
217 | point[2] = cos(t);
|
---|
218 | PointsOnSurface.push_back(point);
|
---|
219 | }
|
---|
220 | ASSERT(PointsOnSurface.size() == N,
|
---|
221 | "Sphere_impl::getHomogeneousPointsOnSurface() did not create "
|
---|
222 | +::toString(N)+" but "+::toString(PointsOnSurface.size())+" points.");
|
---|
223 | } else {
|
---|
224 | // Symmetric but only approximately N points.
|
---|
225 | double a=4*M_PI/N;
|
---|
226 | double d= sqrt(a);
|
---|
227 | int Mtheta=int(M_PI/d);
|
---|
228 | double dtheta=M_PI/Mtheta;
|
---|
229 | double dphi=a/dtheta;
|
---|
230 | for (int m=0; m<Mtheta; m++)
|
---|
231 | {
|
---|
232 | double theta=M_PI*(m+0.5)/Mtheta;
|
---|
233 | int Mphi=int(2*M_PI*sin(theta)/dphi);
|
---|
234 | for (int n=0; n<Mphi;n++)
|
---|
235 | {
|
---|
236 | double phi= 2*M_PI*n/Mphi;
|
---|
237 | Vector point;
|
---|
238 | point.Zero();
|
---|
239 | point[0]=sin(theta)*cos(phi);
|
---|
240 | point[1]=sin(theta)*sin(phi);
|
---|
241 | point[2]=cos(theta);
|
---|
242 | PointsOnSurface.push_back(point);
|
---|
243 | }
|
---|
244 | }
|
---|
245 | }
|
---|
246 | return PointsOnSurface;
|
---|
247 | }
|
---|
248 |
|
---|
249 | std::vector<Vector> Sphere_impl::getHomogeneousPointsInVolume(const size_t N) const {
|
---|
250 | ASSERT(0,
|
---|
251 | "Sphere_impl::getHomogeneousPointsInVolume() - not implemented.");
|
---|
252 | return std::vector<Vector>();
|
---|
253 | }
|
---|
254 |
|
---|
255 | Shape Sphere(){
|
---|
256 | Shape::impl_ptr impl = Shape::impl_ptr(new Sphere_impl());
|
---|
257 | return Shape(impl);
|
---|
258 | }
|
---|
259 |
|
---|
260 | Shape Sphere(const Vector ¢er,double radius){
|
---|
261 | return translate(resize(Sphere(),radius),center);
|
---|
262 | }
|
---|
263 |
|
---|
264 | Shape Ellipsoid(const Vector ¢er, const Vector &radius){
|
---|
265 | return translate(stretch(Sphere(),radius),center);
|
---|
266 | }
|
---|
267 |
|
---|
268 | bool Cuboid_impl::isInside(const Vector &point) const{
|
---|
269 | return (point[0]>=0 && point[0]<=1) && (point[1]>=0 && point[1]<=1) && (point[2]>=0 && point[2]<=1);
|
---|
270 | }
|
---|
271 |
|
---|
272 | bool Cuboid_impl::isOnSurface(const Vector &point) const{
|
---|
273 | bool retVal = isInside(point);
|
---|
274 | // test all borders of the cuboid
|
---|
275 | // double fabs
|
---|
276 | retVal = retVal &&
|
---|
277 | (((fabs(point[0]-1.) < MYEPSILON) || (fabs(point[0]) < MYEPSILON)) ||
|
---|
278 | ((fabs(point[1]-1.) < MYEPSILON) || (fabs(point[1]) < MYEPSILON)) ||
|
---|
279 | ((fabs(point[2]-1.) < MYEPSILON) || (fabs(point[2]) < MYEPSILON)));
|
---|
280 | return retVal;
|
---|
281 | }
|
---|
282 |
|
---|
283 | Vector Cuboid_impl::getNormal(const Vector &point) const throw(NotOnSurfaceException){
|
---|
284 | if(!isOnSurface(point)){
|
---|
285 | throw NotOnSurfaceException() << ShapeVector(&point);
|
---|
286 | }
|
---|
287 | Vector res;
|
---|
288 | // figure out on which sides the Vector lies (maximum 3, when it is in a corner)
|
---|
289 | for(int i=NDIM;i--;){
|
---|
290 | if(fabs(fabs(point[i])-1)<MYEPSILON){
|
---|
291 | // add the scaled (-1/+1) Vector to the set of surface vectors
|
---|
292 | res[i] = point[i];
|
---|
293 | }
|
---|
294 | }
|
---|
295 | ASSERT(res.NormSquared()>=1 && res.NormSquared()<=3,"To many or to few sides found for this Vector");
|
---|
296 |
|
---|
297 | res.Normalize();
|
---|
298 | return res;
|
---|
299 | }
|
---|
300 |
|
---|
301 |
|
---|
302 | Vector Cuboid_impl::getCenter() const
|
---|
303 | {
|
---|
304 | return Vector(0.5,0.5,0.5);
|
---|
305 | }
|
---|
306 |
|
---|
307 | double Cuboid_impl::getRadius() const
|
---|
308 | {
|
---|
309 | return .5;
|
---|
310 | }
|
---|
311 |
|
---|
312 | double Cuboid_impl::getVolume() const
|
---|
313 | {
|
---|
314 | return 1.; // l^3
|
---|
315 | }
|
---|
316 |
|
---|
317 | double Cuboid_impl::getSurfaceArea() const
|
---|
318 | {
|
---|
319 | return 6.; // 6 * l^2
|
---|
320 | }
|
---|
321 |
|
---|
322 | LineSegmentSet Cuboid_impl::getLineIntersections(const Line &line) const{
|
---|
323 | LineSegmentSet res(line);
|
---|
324 | // get the intersection on each of the six faces
|
---|
325 | std::vector<Vector> intersections;
|
---|
326 | intersections.resize(2);
|
---|
327 | int c=0;
|
---|
328 | int x[2]={-1,+1};
|
---|
329 | for(int i=NDIM;i--;){
|
---|
330 | for(int p=0;p<2;++p){
|
---|
331 | if(c==2) goto end; // I know this sucks, but breaking two loops is stupid
|
---|
332 | Vector base;
|
---|
333 | base[i]=x[p];
|
---|
334 | // base now points to the surface and is normal to it at the same time
|
---|
335 | Plane p(base,base);
|
---|
336 | Vector intersection = p.GetIntersection(line);
|
---|
337 | if(isInside(intersection)){
|
---|
338 | // if we have a point on the edge it might already be contained
|
---|
339 | if(c==1 && intersections[0]==intersection)
|
---|
340 | continue;
|
---|
341 | intersections[c++]=intersection;
|
---|
342 | }
|
---|
343 | }
|
---|
344 | }
|
---|
345 | end:
|
---|
346 | if(c==2){
|
---|
347 | res.insert(LineSegment(intersections[0],intersections[1]));
|
---|
348 | }
|
---|
349 | return res;
|
---|
350 | }
|
---|
351 |
|
---|
352 | std::string Cuboid_impl::toString() const{
|
---|
353 | return "Cuboid()";
|
---|
354 | }
|
---|
355 |
|
---|
356 | enum ShapeType Cuboid_impl::getType() const
|
---|
357 | {
|
---|
358 | return CuboidType;
|
---|
359 | }
|
---|
360 |
|
---|
361 | /**
|
---|
362 | * \param N number of points on surface
|
---|
363 | */
|
---|
364 | std::vector<Vector> Cuboid_impl::getHomogeneousPointsOnSurface(const size_t N) const {
|
---|
365 | std::vector<Vector> PointsOnSurface;
|
---|
366 | ASSERT(false, "Cuboid_impl::getHomogeneousPointsOnSurface() not implemented yet");
|
---|
367 | return PointsOnSurface;
|
---|
368 | }
|
---|
369 |
|
---|
370 | std::vector<Vector> Cuboid_impl::getHomogeneousPointsInVolume(const size_t N) const {
|
---|
371 | ASSERT(0,
|
---|
372 | "Cuboid_impl::getHomogeneousPointsInVolume() - not implemented.");
|
---|
373 | return std::vector<Vector>();
|
---|
374 | }
|
---|
375 |
|
---|
376 | Shape Cuboid(){
|
---|
377 | Shape::impl_ptr impl = Shape::impl_ptr(new Cuboid_impl());
|
---|
378 | return Shape(impl);
|
---|
379 | }
|
---|
380 |
|
---|
381 | Shape Cuboid(const Vector &corner1, const Vector &corner2){
|
---|
382 | // make sure the two edges are upper left front and lower right back
|
---|
383 | Vector sortedC1;
|
---|
384 | Vector sortedC2;
|
---|
385 | for(int i=NDIM;i--;){
|
---|
386 | sortedC1[i] = std::min(corner1[i],corner2[i]);
|
---|
387 | sortedC2[i] = std::max(corner1[i],corner2[i]);
|
---|
388 | ASSERT(corner1[i]!=corner2[i],"Given points for cuboid edges did not define a valid space");
|
---|
389 | }
|
---|
390 | // get the middle point
|
---|
391 | Vector middle = (1./2.)*(sortedC1+sortedC2);
|
---|
392 | Vector factors = sortedC2-middle;
|
---|
393 | return translate(stretch(Cuboid(),factors),middle);
|
---|
394 | }
|
---|