source: src/Shapes/BaseShapes.cpp@ 6f0507e

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 6f0507e was 6f0507e, checked in by Frederik Heber <heber@…>, 12 years ago

Implemented getHomogenousPointsOnSurface.

Incoming N is rounded to the next possible value to create evently
spaced points.

  • Property mode set to 100644
File size: 10.4 KB
Line 
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// ----------------------------------------------------------------------------
42bool 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
47bool 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
53Vector 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
64Vector Cylinder_impl::getCenter() const
65{
66 return Vector(0.0, 0.0, 0.0);
67}
68
69double Cylinder_impl::getRadius() const
70{
71 return 1.0;
72}
73
74double Cylinder_impl::getVolume() const
75{
76 return M_PI*2.0; // pi r^2 h
77}
78
79double Cylinder_impl::getSurfaceArea() const
80{
81 return 2.0*M_PI*2.0; // 2 pi r h
82}
83
84LineSegmentSet Cylinder_impl::getLineIntersections(const Line &line) const {
85 // TODO
86 ASSERT(0, "Cylinder_impl::getLineIntersections - not implemented.");
87}
88
89std::string Cylinder_impl::toString() const
90{
91 return "Cylinder()";
92}
93
94enum ShapeType Cylinder_impl::getType() const
95{
96 return CylinderType;
97}
98
99std::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
116std::vector<Vector> Cylinder_impl::getHomogeneousPointsInVolume(const size_t N) const {
117 // TODO
118 ASSERT(0, "Cylinder_impl::getHomogeneousPointsInVolume - not implemented.");
119}
120
121Shape Cylinder() {
122 Shape::impl_ptr impl = Shape::impl_ptr(new Cylinder_impl());
123 return Shape(impl);
124}
125
126Shape Cylinder(const Vector &center, 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
142bool Sphere_impl::isInside(const Vector &point) const{
143 return point.NormSquared()<=1.;
144}
145
146bool Sphere_impl::isOnSurface(const Vector &point) const{
147 return fabs(point.NormSquared()-1.)<MYEPSILON;
148}
149
150Vector Sphere_impl::getNormal(const Vector &point) const throw(NotOnSurfaceException){
151 if(!isOnSurface(point)){
152 throw NotOnSurfaceException() << ShapeVector(&point);
153 }
154 return point;
155}
156
157Vector Sphere_impl::getCenter() const
158{
159 return Vector(0.,0.,0.);
160}
161
162double Sphere_impl::getRadius() const
163{
164 return 1.;
165}
166
167double Sphere_impl::getVolume() const
168{
169 return (4./3.)*M_PI; // 4/3 pi r^3
170}
171
172double Sphere_impl::getSurfaceArea() const
173{
174 return 2.*M_PI; // 2 pi r^2
175}
176
177
178LineSegmentSet 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
187std::string Sphere_impl::toString() const{
188 return "Sphere()";
189}
190
191enum 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 */
200std::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
249std::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
255Shape Sphere(){
256 Shape::impl_ptr impl = Shape::impl_ptr(new Sphere_impl());
257 return Shape(impl);
258}
259
260Shape Sphere(const Vector &center,double radius){
261 return translate(resize(Sphere(),radius),center);
262}
263
264Shape Ellipsoid(const Vector &center, const Vector &radius){
265 return translate(stretch(Sphere(),radius),center);
266}
267
268bool 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
272bool 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
283Vector 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
302Vector Cuboid_impl::getCenter() const
303{
304 return Vector(0.5,0.5,0.5);
305}
306
307double Cuboid_impl::getRadius() const
308{
309 return .5;
310}
311
312double Cuboid_impl::getVolume() const
313{
314 return 1.; // l^3
315}
316
317double Cuboid_impl::getSurfaceArea() const
318{
319 return 6.; // 6 * l^2
320}
321
322LineSegmentSet 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
352std::string Cuboid_impl::toString() const{
353 return "Cuboid()";
354}
355
356enum ShapeType Cuboid_impl::getType() const
357{
358 return CuboidType;
359}
360
361/**
362 * \param N number of points on surface
363 */
364std::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
370std::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
376Shape Cuboid(){
377 Shape::impl_ptr impl = Shape::impl_ptr(new Cuboid_impl());
378 return Shape(impl);
379}
380
381Shape 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}
Note: See TracBrowser for help on using the repository browser.