[bcf653] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
[0aa122] | 4 | * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
|
---|
[94d5ac6] | 5 | *
|
---|
| 6 | *
|
---|
| 7 | * This file is part of MoleCuilder.
|
---|
| 8 | *
|
---|
| 9 | * MoleCuilder is free software: you can redistribute it and/or modify
|
---|
| 10 | * it under the terms of the GNU General Public License as published by
|
---|
| 11 | * the Free Software Foundation, either version 2 of the License, or
|
---|
| 12 | * (at your option) any later version.
|
---|
| 13 | *
|
---|
| 14 | * MoleCuilder is distributed in the hope that it will be useful,
|
---|
| 15 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 16 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 17 | * GNU General Public License for more details.
|
---|
| 18 | *
|
---|
| 19 | * You should have received a copy of the GNU General Public License
|
---|
| 20 | * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
|
---|
[bcf653] | 21 | */
|
---|
| 22 |
|
---|
[e38447] | 23 | /*
|
---|
| 24 | * BaseShapes_impl.cpp
|
---|
| 25 | *
|
---|
| 26 | * Created on: Jun 18, 2010
|
---|
| 27 | * Author: crueger
|
---|
| 28 | */
|
---|
| 29 |
|
---|
[bf3817] | 30 | // include config.h
|
---|
| 31 | #ifdef HAVE_CONFIG_H
|
---|
| 32 | #include <config.h>
|
---|
| 33 | #endif
|
---|
| 34 |
|
---|
[ad011c] | 35 | #include "CodePatterns/MemDebug.hpp"
|
---|
[bbbad5] | 36 |
|
---|
[e38447] | 37 | #include "Shapes/BaseShapes.hpp"
|
---|
| 38 | #include "Shapes/BaseShapes_impl.hpp"
|
---|
[b94634] | 39 | #include "Shapes/ShapeExceptions.hpp"
|
---|
[f3526d] | 40 | #include "Shapes/ShapeOps.hpp"
|
---|
[e38447] | 41 |
|
---|
[e4fe8d] | 42 | #include "Helpers/defs.hpp"
|
---|
[5de9da] | 43 |
|
---|
[ad011c] | 44 | #include "CodePatterns/Assert.hpp"
|
---|
[57f243] | 45 | #include "LinearAlgebra/Vector.hpp"
|
---|
[0eb8f4] | 46 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
[6c438f] | 47 | #include "LinearAlgebra/Line.hpp"
|
---|
| 48 | #include "LinearAlgebra/Plane.hpp"
|
---|
| 49 | #include "LinearAlgebra/LineSegment.hpp"
|
---|
| 50 | #include "LinearAlgebra/LineSegmentSet.hpp"
|
---|
[c6f395] | 51 |
|
---|
[5de9da] | 52 | #include <cmath>
|
---|
[d76a7c] | 53 | #include <algorithm>
|
---|
[e38447] | 54 |
|
---|
[0eb8f4] | 55 | // CYLINDER CODE
|
---|
| 56 | // ----------------------------------------------------------------------------
|
---|
| 57 | bool Cylinder_impl::isInside(const Vector &point) const {
|
---|
| 58 | return (Vector(point[0], point[1], 0.0).NormSquared() < 1.0+MYEPSILON) &&
|
---|
| 59 | (point[2] > -1.0-MYEPSILON) && (point[2] < 1.0+MYEPSILON);
|
---|
| 60 | }
|
---|
| 61 |
|
---|
| 62 | bool Cylinder_impl::isOnSurface(const Vector &point) const {
|
---|
| 63 | return fabs(Vector(point[0], point[1], 0.0).NormSquared()-1.0)<MYEPSILON &&
|
---|
| 64 | (point[2] > -1.0-MYEPSILON) && (point[2] < 1.0+MYEPSILON);
|
---|
| 65 |
|
---|
| 66 | }
|
---|
| 67 |
|
---|
| 68 | Vector Cylinder_impl::getNormal(const Vector &point) const throw(NotOnSurfaceException) {
|
---|
| 69 | if(!isOnSurface(point)){
|
---|
| 70 | throw NotOnSurfaceException() << ShapeVector(&point);
|
---|
| 71 | }
|
---|
| 72 |
|
---|
| 73 | if ((fabs(point[2]-1)<MYEPSILON) || (fabs(point[2])<MYEPSILON))
|
---|
| 74 | return Vector(0.0, 0.0, point[2]);
|
---|
| 75 | else
|
---|
| 76 | return Vector(point[0], point[1], 0.0);
|
---|
| 77 | }
|
---|
| 78 |
|
---|
| 79 | Vector Cylinder_impl::getCenter() const
|
---|
| 80 | {
|
---|
| 81 | return Vector(0.0, 0.0, 0.0);
|
---|
| 82 | }
|
---|
| 83 |
|
---|
| 84 | double Cylinder_impl::getRadius() const
|
---|
| 85 | {
|
---|
| 86 | return 1.0;
|
---|
| 87 | }
|
---|
| 88 |
|
---|
| 89 | double Cylinder_impl::getVolume() const
|
---|
| 90 | {
|
---|
| 91 | return M_PI*2.0; // pi r^2 h
|
---|
| 92 | }
|
---|
| 93 |
|
---|
| 94 | double Cylinder_impl::getSurfaceArea() const
|
---|
| 95 | {
|
---|
| 96 | return 2.0*M_PI*2.0; // 2 pi r h
|
---|
| 97 | }
|
---|
| 98 |
|
---|
| 99 | LineSegmentSet Cylinder_impl::getLineIntersections(const Line &line) const {
|
---|
[f4a863] | 100 | const Vector origin = line.getOrigin();
|
---|
| 101 | const Vector direction = line.getDirection();
|
---|
| 102 |
|
---|
| 103 | const Vector e(direction[0], direction[1], 0.0);
|
---|
| 104 | const Vector f(origin[0], origin[1], 0.0);
|
---|
| 105 | const double A = e.ScalarProduct(e);
|
---|
| 106 | const double B = 2.0*e.ScalarProduct(f);
|
---|
| 107 | const double C = f.ScalarProduct(f) - 1.0;
|
---|
| 108 |
|
---|
| 109 | std::vector<double> solutions;
|
---|
| 110 |
|
---|
| 111 | // Common routine to solve quadratic quations, anywhere?
|
---|
| 112 | const double neg_p_half = -B/(2.0*A);
|
---|
| 113 | const double q = C/A;
|
---|
| 114 | const double radicant = neg_p_half*neg_p_half-q;
|
---|
| 115 |
|
---|
| 116 | if (radicant > 0.0) {
|
---|
| 117 | const double root = sqrt(radicant);
|
---|
| 118 | solutions.push_back(neg_p_half+root);
|
---|
| 119 | const double sln2 = neg_p_half-root;
|
---|
| 120 | if (sln2 != solutions.back())
|
---|
| 121 | solutions.push_back(sln2);
|
---|
| 122 | }
|
---|
| 123 |
|
---|
| 124 | // Now get parameter for intersection with z-Planes.
|
---|
| 125 | const double origin_z = origin[2];
|
---|
| 126 | const double dir_z = direction[2];
|
---|
| 127 |
|
---|
| 128 | if (dir_z != 0.0) {
|
---|
| 129 | solutions.push_back((-1.0-origin_z)/dir_z);
|
---|
| 130 | solutions.push_back((1.0-origin_z)/dir_z);
|
---|
| 131 | }
|
---|
| 132 |
|
---|
| 133 | // Calculate actual vectors from obtained parameters and check,
|
---|
| 134 | // if they are actual intersections.
|
---|
| 135 | std::vector<Vector> intersections;
|
---|
| 136 |
|
---|
| 137 | for(unsigned int i=0; i<solutions.size(); i++) {
|
---|
| 138 | const Vector check_me(origin + direction*solutions[i]);
|
---|
| 139 | if (isOnSurface(check_me))
|
---|
| 140 | intersections.push_back(check_me);
|
---|
| 141 | }
|
---|
| 142 |
|
---|
| 143 | LineSegmentSet result(line);
|
---|
| 144 | if (intersections.size()==2)
|
---|
| 145 | result.insert(LineSegment(intersections[0], intersections[1]));
|
---|
| 146 | return result;
|
---|
[0eb8f4] | 147 | }
|
---|
| 148 |
|
---|
| 149 | std::string Cylinder_impl::toString() const
|
---|
| 150 | {
|
---|
| 151 | return "Cylinder()";
|
---|
| 152 | }
|
---|
| 153 |
|
---|
| 154 | enum ShapeType Cylinder_impl::getType() const
|
---|
| 155 | {
|
---|
| 156 | return CylinderType;
|
---|
| 157 | }
|
---|
| 158 |
|
---|
| 159 | std::vector<Vector> Cylinder_impl::getHomogeneousPointsOnSurface(const size_t N) const {
|
---|
[9e2737] | 160 | const double nz_float = sqrt(N/M_PI);
|
---|
| 161 | const int nu = round(N/nz_float);
|
---|
| 162 | const int nz = round(nz_float);
|
---|
[6f0507e] | 163 |
|
---|
| 164 | const double dphi = 2.0*M_PI/nu;
|
---|
| 165 | const double dz = 2.0/nz;
|
---|
| 166 |
|
---|
| 167 | std::vector<Vector> result;
|
---|
| 168 |
|
---|
| 169 | for(int useg=0; useg<nu; useg++)
|
---|
[9e2737] | 170 | for(int zseg=0; zseg<nz; zseg++)
|
---|
[6f0507e] | 171 | result.push_back(Vector(cos(useg*dphi), sin(useg*dphi), zseg*dz-1.0));
|
---|
| 172 |
|
---|
| 173 | return result;
|
---|
[0eb8f4] | 174 | }
|
---|
| 175 |
|
---|
| 176 | std::vector<Vector> Cylinder_impl::getHomogeneousPointsInVolume(const size_t N) const {
|
---|
[9e2737] | 177 | const double nz_float = pow(N/(2.0*M_PI), 1.0/3.0);
|
---|
| 178 | const int nu = round(nz_float*M_PI);
|
---|
| 179 | const int nr = round(nz_float*0.5);
|
---|
| 180 | const int nz = round(nz_float);
|
---|
| 181 |
|
---|
| 182 | const double dphi = 2.0*M_PI/nu;
|
---|
| 183 | const double dz = 2.0/nz;
|
---|
| 184 | const double dr = 1.0/nr;
|
---|
| 185 |
|
---|
| 186 | std::vector<Vector> result;
|
---|
| 187 |
|
---|
| 188 | for(int useg=0; useg<nu; useg++)
|
---|
| 189 | for(int zseg=0; zseg<nz; zseg++)
|
---|
| 190 | for(int rseg=0; rseg<nr; rseg++)
|
---|
| 191 | {
|
---|
[5d4179f] | 192 | const double r = dr+rseg*dr;
|
---|
[9e2737] | 193 | result.push_back(Vector(r*cos(useg*dphi), r*sin(useg*dphi), zseg*dz-1.0));
|
---|
| 194 | }
|
---|
| 195 |
|
---|
| 196 | return result;
|
---|
[0eb8f4] | 197 | }
|
---|
| 198 |
|
---|
| 199 | Shape Cylinder() {
|
---|
| 200 | Shape::impl_ptr impl = Shape::impl_ptr(new Cylinder_impl());
|
---|
| 201 | return Shape(impl);
|
---|
| 202 | }
|
---|
| 203 |
|
---|
| 204 | Shape Cylinder(const Vector ¢er, const double xrot, const double yrot,
|
---|
| 205 | const double height, const double radius)
|
---|
| 206 | {
|
---|
| 207 | RealSpaceMatrix rot;
|
---|
| 208 | rot.setRotation(xrot, yrot, 0.0);
|
---|
| 209 |
|
---|
| 210 | return translate(
|
---|
| 211 | transform(
|
---|
| 212 | stretch(
|
---|
| 213 | Cylinder(),
|
---|
| 214 | Vector(radius, radius, height*0.5)),
|
---|
| 215 | rot),
|
---|
| 216 | center);
|
---|
| 217 | }
|
---|
| 218 | // ----------------------------------------------------------------------------
|
---|
| 219 |
|
---|
[735940] | 220 | bool Sphere_impl::isInside(const Vector &point) const{
|
---|
| 221 | return point.NormSquared()<=1.;
|
---|
[e38447] | 222 | }
|
---|
| 223 |
|
---|
[735940] | 224 | bool Sphere_impl::isOnSurface(const Vector &point) const{
|
---|
| 225 | return fabs(point.NormSquared()-1.)<MYEPSILON;
|
---|
[5de9da] | 226 | }
|
---|
| 227 |
|
---|
[735940] | 228 | Vector Sphere_impl::getNormal(const Vector &point) const throw(NotOnSurfaceException){
|
---|
[5de9da] | 229 | if(!isOnSurface(point)){
|
---|
[b94634] | 230 | throw NotOnSurfaceException() << ShapeVector(&point);
|
---|
[5de9da] | 231 | }
|
---|
| 232 | return point;
|
---|
| 233 | }
|
---|
| 234 |
|
---|
[6acc2f3] | 235 | Vector Sphere_impl::getCenter() const
|
---|
| 236 | {
|
---|
| 237 | return Vector(0.,0.,0.);
|
---|
| 238 | }
|
---|
| 239 |
|
---|
| 240 | double Sphere_impl::getRadius() const
|
---|
| 241 | {
|
---|
| 242 | return 1.;
|
---|
| 243 | }
|
---|
| 244 |
|
---|
[c67c65] | 245 | double Sphere_impl::getVolume() const
|
---|
| 246 | {
|
---|
| 247 | return (4./3.)*M_PI; // 4/3 pi r^3
|
---|
| 248 | }
|
---|
| 249 |
|
---|
| 250 | double Sphere_impl::getSurfaceArea() const
|
---|
| 251 | {
|
---|
| 252 | return 2.*M_PI; // 2 pi r^2
|
---|
| 253 | }
|
---|
| 254 |
|
---|
[6acc2f3] | 255 |
|
---|
[735940] | 256 | LineSegmentSet Sphere_impl::getLineIntersections(const Line &line) const{
|
---|
[c6f395] | 257 | LineSegmentSet res(line);
|
---|
| 258 | std::vector<Vector> intersections = line.getSphereIntersections();
|
---|
| 259 | if(intersections.size()==2){
|
---|
| 260 | res.insert(LineSegment(intersections[0],intersections[1]));
|
---|
| 261 | }
|
---|
| 262 | return res;
|
---|
| 263 | }
|
---|
| 264 |
|
---|
[b92e4a] | 265 | std::string Sphere_impl::toString() const{
|
---|
[cfda65] | 266 | return "Sphere()";
|
---|
| 267 | }
|
---|
| 268 |
|
---|
[b92e4a] | 269 | enum ShapeType Sphere_impl::getType() const
|
---|
| 270 | {
|
---|
| 271 | return SphereType;
|
---|
| 272 | }
|
---|
| 273 |
|
---|
[c5186e] | 274 | /**
|
---|
| 275 | * algorithm taken from http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere
|
---|
| 276 | * \param N number of points on surface
|
---|
| 277 | */
|
---|
[f6ba43] | 278 | std::vector<Vector> Sphere_impl::getHomogeneousPointsOnSurface(const size_t N) const
|
---|
| 279 | {
|
---|
[c5186e] | 280 | std::vector<Vector> PointsOnSurface;
|
---|
[125841] | 281 | if (true) {
|
---|
| 282 | // Exactly N points but not symmetric.
|
---|
| 283 |
|
---|
| 284 | // This formula is derived by finding a curve on the sphere that spirals down from
|
---|
| 285 | // the north pole to the south pole keeping a constant distance between consecutive turns.
|
---|
| 286 | // The curve is then parametrized by arch length and evaluated in constant intervals.
|
---|
| 287 | double a = sqrt(N) * 2;
|
---|
| 288 | for (int i=0; i<N; i++){
|
---|
| 289 | double t0 = ((double)i + 0.5) / (double)N;
|
---|
| 290 | double t = (sqrt(t0) - sqrt(1.0 - t0) + 1.0) / 2.0 * M_PI;
|
---|
| 291 | Vector point;
|
---|
| 292 | point.Zero();
|
---|
| 293 | point[0] = sin(t) * sin(t * a);
|
---|
| 294 | point[1] = sin(t) * cos(t * a);
|
---|
| 295 | point[2] = cos(t);
|
---|
| 296 | PointsOnSurface.push_back(point);
|
---|
| 297 | }
|
---|
| 298 | ASSERT(PointsOnSurface.size() == N,
|
---|
| 299 | "Sphere_impl::getHomogeneousPointsOnSurface() did not create "
|
---|
| 300 | +::toString(N)+" but "+::toString(PointsOnSurface.size())+" points.");
|
---|
| 301 | } else {
|
---|
| 302 | // Symmetric but only approximately N points.
|
---|
| 303 | double a=4*M_PI/N;
|
---|
[faca99] | 304 | double d= sqrt(a);
|
---|
[125841] | 305 | int Mtheta=int(M_PI/d);
|
---|
| 306 | double dtheta=M_PI/Mtheta;
|
---|
[faca99] | 307 | double dphi=a/dtheta;
|
---|
| 308 | for (int m=0; m<Mtheta; m++)
|
---|
| 309 | {
|
---|
[125841] | 310 | double theta=M_PI*(m+0.5)/Mtheta;
|
---|
| 311 | int Mphi=int(2*M_PI*sin(theta)/dphi);
|
---|
| 312 | for (int n=0; n<Mphi;n++)
|
---|
| 313 | {
|
---|
| 314 | double phi= 2*M_PI*n/Mphi;
|
---|
| 315 | Vector point;
|
---|
| 316 | point.Zero();
|
---|
| 317 | point[0]=sin(theta)*cos(phi);
|
---|
| 318 | point[1]=sin(theta)*sin(phi);
|
---|
| 319 | point[2]=cos(theta);
|
---|
| 320 | PointsOnSurface.push_back(point);
|
---|
| 321 | }
|
---|
[faca99] | 322 | }
|
---|
[125841] | 323 | }
|
---|
[c5186e] | 324 | return PointsOnSurface;
|
---|
| 325 | }
|
---|
| 326 |
|
---|
[5a8d61] | 327 | std::vector<Vector> Sphere_impl::getHomogeneousPointsInVolume(const size_t N) const {
|
---|
| 328 | ASSERT(0,
|
---|
| 329 | "Sphere_impl::getHomogeneousPointsInVolume() - not implemented.");
|
---|
| 330 | return std::vector<Vector>();
|
---|
| 331 | }
|
---|
[c5186e] | 332 |
|
---|
[e38447] | 333 | Shape Sphere(){
|
---|
| 334 | Shape::impl_ptr impl = Shape::impl_ptr(new Sphere_impl());
|
---|
| 335 | return Shape(impl);
|
---|
| 336 | }
|
---|
| 337 |
|
---|
[f3526d] | 338 | Shape Sphere(const Vector ¢er,double radius){
|
---|
| 339 | return translate(resize(Sphere(),radius),center);
|
---|
| 340 | }
|
---|
| 341 |
|
---|
| 342 | Shape Ellipsoid(const Vector ¢er, const Vector &radius){
|
---|
| 343 | return translate(stretch(Sphere(),radius),center);
|
---|
| 344 | }
|
---|
| 345 |
|
---|
[735940] | 346 | bool Cuboid_impl::isInside(const Vector &point) const{
|
---|
[0d02fb] | 347 | return (point[0]>=0 && point[0]<=1) && (point[1]>=0 && point[1]<=1) && (point[2]>=0 && point[2]<=1);
|
---|
[5de9da] | 348 | }
|
---|
| 349 |
|
---|
[735940] | 350 | bool Cuboid_impl::isOnSurface(const Vector &point) const{
|
---|
[5de9da] | 351 | bool retVal = isInside(point);
|
---|
| 352 | // test all borders of the cuboid
|
---|
| 353 | // double fabs
|
---|
| 354 | retVal = retVal &&
|
---|
[6c438f] | 355 | (((fabs(point[0]-1.) < MYEPSILON) || (fabs(point[0]) < MYEPSILON)) ||
|
---|
| 356 | ((fabs(point[1]-1.) < MYEPSILON) || (fabs(point[1]) < MYEPSILON)) ||
|
---|
| 357 | ((fabs(point[2]-1.) < MYEPSILON) || (fabs(point[2]) < MYEPSILON)));
|
---|
[5de9da] | 358 | return retVal;
|
---|
| 359 | }
|
---|
| 360 |
|
---|
[735940] | 361 | Vector Cuboid_impl::getNormal(const Vector &point) const throw(NotOnSurfaceException){
|
---|
[5de9da] | 362 | if(!isOnSurface(point)){
|
---|
[b94634] | 363 | throw NotOnSurfaceException() << ShapeVector(&point);
|
---|
[5de9da] | 364 | }
|
---|
| 365 | Vector res;
|
---|
| 366 | // figure out on which sides the Vector lies (maximum 3, when it is in a corner)
|
---|
| 367 | for(int i=NDIM;i--;){
|
---|
[7de208] | 368 | if((fabs(point[i])<MYEPSILON) || (fabs(point[i]-1)<MYEPSILON)){
|
---|
[5de9da] | 369 | // add the scaled (-1/+1) Vector to the set of surface vectors
|
---|
[7de208] | 370 | res[i] = point[i] * 2.0 - 1.0;
|
---|
[5de9da] | 371 | }
|
---|
| 372 | }
|
---|
| 373 | ASSERT(res.NormSquared()>=1 && res.NormSquared()<=3,"To many or to few sides found for this Vector");
|
---|
| 374 |
|
---|
| 375 | res.Normalize();
|
---|
| 376 | return res;
|
---|
[e38447] | 377 | }
|
---|
| 378 |
|
---|
[6acc2f3] | 379 |
|
---|
| 380 | Vector Cuboid_impl::getCenter() const
|
---|
| 381 | {
|
---|
| 382 | return Vector(0.5,0.5,0.5);
|
---|
| 383 | }
|
---|
| 384 |
|
---|
| 385 | double Cuboid_impl::getRadius() const
|
---|
| 386 | {
|
---|
| 387 | return .5;
|
---|
| 388 | }
|
---|
| 389 |
|
---|
[c67c65] | 390 | double Cuboid_impl::getVolume() const
|
---|
| 391 | {
|
---|
| 392 | return 1.; // l^3
|
---|
| 393 | }
|
---|
| 394 |
|
---|
| 395 | double Cuboid_impl::getSurfaceArea() const
|
---|
| 396 | {
|
---|
| 397 | return 6.; // 6 * l^2
|
---|
| 398 | }
|
---|
| 399 |
|
---|
[735940] | 400 | LineSegmentSet Cuboid_impl::getLineIntersections(const Line &line) const{
|
---|
[c6f395] | 401 | LineSegmentSet res(line);
|
---|
| 402 | // get the intersection on each of the six faces
|
---|
[955b91] | 403 | std::vector<Vector> intersections;
|
---|
[c6f395] | 404 | intersections.resize(2);
|
---|
| 405 | int c=0;
|
---|
| 406 | int x[2]={-1,+1};
|
---|
| 407 | for(int i=NDIM;i--;){
|
---|
| 408 | for(int p=0;p<2;++p){
|
---|
| 409 | if(c==2) goto end; // I know this sucks, but breaking two loops is stupid
|
---|
| 410 | Vector base;
|
---|
| 411 | base[i]=x[p];
|
---|
| 412 | // base now points to the surface and is normal to it at the same time
|
---|
| 413 | Plane p(base,base);
|
---|
| 414 | Vector intersection = p.GetIntersection(line);
|
---|
| 415 | if(isInside(intersection)){
|
---|
| 416 | // if we have a point on the edge it might already be contained
|
---|
| 417 | if(c==1 && intersections[0]==intersection)
|
---|
| 418 | continue;
|
---|
| 419 | intersections[c++]=intersection;
|
---|
| 420 | }
|
---|
| 421 | }
|
---|
| 422 | }
|
---|
| 423 | end:
|
---|
| 424 | if(c==2){
|
---|
| 425 | res.insert(LineSegment(intersections[0],intersections[1]));
|
---|
| 426 | }
|
---|
| 427 | return res;
|
---|
| 428 | }
|
---|
| 429 |
|
---|
[b92e4a] | 430 | std::string Cuboid_impl::toString() const{
|
---|
[cfda65] | 431 | return "Cuboid()";
|
---|
| 432 | }
|
---|
| 433 |
|
---|
[b92e4a] | 434 | enum ShapeType Cuboid_impl::getType() const
|
---|
| 435 | {
|
---|
| 436 | return CuboidType;
|
---|
| 437 | }
|
---|
| 438 |
|
---|
[c5186e] | 439 | /**
|
---|
| 440 | * \param N number of points on surface
|
---|
| 441 | */
|
---|
[9c1c89] | 442 | std::vector<Vector> Cuboid_impl::getHomogeneousPointsOnSurface(const size_t N) const {
|
---|
[c5186e] | 443 | std::vector<Vector> PointsOnSurface;
|
---|
[bf8318] | 444 | // sides
|
---|
| 445 | int n = sqrt((N - 1) / 6) + 1;
|
---|
| 446 | for (int i=0; i<=n; i++){
|
---|
| 447 | double ii = (double)i / (double)n;
|
---|
| 448 | for (int k=0; k<n; k++){
|
---|
| 449 | double kk = (double)k / (double)n;
|
---|
| 450 | PointsOnSurface.push_back(Vector(ii, kk, 1));
|
---|
| 451 | PointsOnSurface.push_back(Vector(ii, 1, 1-kk));
|
---|
| 452 | PointsOnSurface.push_back(Vector(ii, 1-kk, 0));
|
---|
| 453 | PointsOnSurface.push_back(Vector(ii, 0, kk));
|
---|
| 454 | }
|
---|
| 455 | }
|
---|
| 456 | // top and bottom
|
---|
| 457 | for (int i=1; i<n; i++){
|
---|
| 458 | double ii = (double)i / (double)n;
|
---|
| 459 | for (int k=1; k<n; k++){
|
---|
| 460 | double kk = (double)k / (double)n;
|
---|
| 461 | PointsOnSurface.push_back(Vector(0, ii, kk));
|
---|
| 462 | PointsOnSurface.push_back(Vector(1, ii, kk));
|
---|
| 463 | }
|
---|
| 464 | }
|
---|
[c5186e] | 465 | return PointsOnSurface;
|
---|
| 466 | }
|
---|
| 467 |
|
---|
[5a8d61] | 468 | std::vector<Vector> Cuboid_impl::getHomogeneousPointsInVolume(const size_t N) const {
|
---|
| 469 | ASSERT(0,
|
---|
| 470 | "Cuboid_impl::getHomogeneousPointsInVolume() - not implemented.");
|
---|
| 471 | return std::vector<Vector>();
|
---|
| 472 | }
|
---|
| 473 |
|
---|
[e38447] | 474 | Shape Cuboid(){
|
---|
[5de9da] | 475 | Shape::impl_ptr impl = Shape::impl_ptr(new Cuboid_impl());
|
---|
[e38447] | 476 | return Shape(impl);
|
---|
| 477 | }
|
---|
[d76a7c] | 478 |
|
---|
| 479 | Shape Cuboid(const Vector &corner1, const Vector &corner2){
|
---|
| 480 | // make sure the two edges are upper left front and lower right back
|
---|
| 481 | Vector sortedC1;
|
---|
| 482 | Vector sortedC2;
|
---|
| 483 | for(int i=NDIM;i--;){
|
---|
[955b91] | 484 | sortedC1[i] = std::min(corner1[i],corner2[i]);
|
---|
| 485 | sortedC2[i] = std::max(corner1[i],corner2[i]);
|
---|
[d76a7c] | 486 | ASSERT(corner1[i]!=corner2[i],"Given points for cuboid edges did not define a valid space");
|
---|
| 487 | }
|
---|
| 488 | // get the middle point
|
---|
| 489 | Vector middle = (1./2.)*(sortedC1+sortedC2);
|
---|
| 490 | Vector factors = sortedC2-middle;
|
---|
| 491 | return translate(stretch(Cuboid(),factors),middle);
|
---|
| 492 | }
|
---|