source: src/Shapes/BaseShapes.cpp@ 7de208

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since 7de208 was 7de208, checked in by Frederik Heber <heber@…>, 12 years ago

FIX: Cuboid_impl::getNormal used wrong coordinates

  • Property mode set to 100644
File size: 13.6 KB
RevLine 
[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// ----------------------------------------------------------------------------
57bool 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
62bool 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
68Vector 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
79Vector Cylinder_impl::getCenter() const
80{
81 return Vector(0.0, 0.0, 0.0);
82}
83
84double Cylinder_impl::getRadius() const
85{
86 return 1.0;
87}
88
89double Cylinder_impl::getVolume() const
90{
91 return M_PI*2.0; // pi r^2 h
92}
93
94double Cylinder_impl::getSurfaceArea() const
95{
96 return 2.0*M_PI*2.0; // 2 pi r h
97}
98
99LineSegmentSet 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
149std::string Cylinder_impl::toString() const
150{
151 return "Cylinder()";
152}
153
154enum ShapeType Cylinder_impl::getType() const
155{
156 return CylinderType;
157}
158
159std::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
176std::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
199Shape Cylinder() {
200 Shape::impl_ptr impl = Shape::impl_ptr(new Cylinder_impl());
201 return Shape(impl);
202}
203
204Shape Cylinder(const Vector &center, 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]220bool Sphere_impl::isInside(const Vector &point) const{
221 return point.NormSquared()<=1.;
[e38447]222}
223
[735940]224bool Sphere_impl::isOnSurface(const Vector &point) const{
225 return fabs(point.NormSquared()-1.)<MYEPSILON;
[5de9da]226}
227
[735940]228Vector 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]235Vector Sphere_impl::getCenter() const
236{
237 return Vector(0.,0.,0.);
238}
239
240double Sphere_impl::getRadius() const
241{
242 return 1.;
243}
244
[c67c65]245double Sphere_impl::getVolume() const
246{
247 return (4./3.)*M_PI; // 4/3 pi r^3
248}
249
250double Sphere_impl::getSurfaceArea() const
251{
252 return 2.*M_PI; // 2 pi r^2
253}
254
[6acc2f3]255
[735940]256LineSegmentSet 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]265std::string Sphere_impl::toString() const{
[cfda65]266 return "Sphere()";
267}
268
[b92e4a]269enum 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]278std::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]327std::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]333Shape Sphere(){
334 Shape::impl_ptr impl = Shape::impl_ptr(new Sphere_impl());
335 return Shape(impl);
336}
337
[f3526d]338Shape Sphere(const Vector &center,double radius){
339 return translate(resize(Sphere(),radius),center);
340}
341
342Shape Ellipsoid(const Vector &center, const Vector &radius){
343 return translate(stretch(Sphere(),radius),center);
344}
345
[735940]346bool 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]350bool 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]361Vector 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
380Vector Cuboid_impl::getCenter() const
381{
382 return Vector(0.5,0.5,0.5);
383}
384
385double Cuboid_impl::getRadius() const
386{
387 return .5;
388}
389
[c67c65]390double Cuboid_impl::getVolume() const
391{
392 return 1.; // l^3
393}
394
395double Cuboid_impl::getSurfaceArea() const
396{
397 return 6.; // 6 * l^2
398}
399
[735940]400LineSegmentSet 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]430std::string Cuboid_impl::toString() const{
[cfda65]431 return "Cuboid()";
432}
433
[b92e4a]434enum ShapeType Cuboid_impl::getType() const
435{
436 return CuboidType;
437}
438
[c5186e]439/**
440 * \param N number of points on surface
441 */
[9c1c89]442std::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]468std::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]474Shape Cuboid(){
[5de9da]475 Shape::impl_ptr impl = Shape::impl_ptr(new Cuboid_impl());
[e38447]476 return Shape(impl);
477}
[d76a7c]478
479Shape 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}
Note: See TracBrowser for help on using the repository browser.