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 | * Box.cpp
|
---|
10 | *
|
---|
11 | * Created on: Jun 30, 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 "Box.hpp"
|
---|
23 |
|
---|
24 | #include <cmath>
|
---|
25 | #include <iostream>
|
---|
26 | #include <cstdlib>
|
---|
27 |
|
---|
28 | #include "CodePatterns/Assert.hpp"
|
---|
29 | #include "CodePatterns/Log.hpp"
|
---|
30 | #include "CodePatterns/Observer/Channels.hpp"
|
---|
31 | #include "CodePatterns/Observer/Notification.hpp"
|
---|
32 | #include "CodePatterns/Verbose.hpp"
|
---|
33 | #include "Helpers/defs.hpp"
|
---|
34 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
35 | #include "LinearAlgebra/Vector.hpp"
|
---|
36 | #include "LinearAlgebra/Plane.hpp"
|
---|
37 | #include "Shapes/BaseShapes.hpp"
|
---|
38 | #include "Shapes/ShapeOps.hpp"
|
---|
39 |
|
---|
40 |
|
---|
41 | Box::Box() :
|
---|
42 | Observable("Box"),
|
---|
43 | M(new RealSpaceMatrix()),
|
---|
44 | Minv(new RealSpaceMatrix())
|
---|
45 | {
|
---|
46 | internal_list.reserve(pow(3,3));
|
---|
47 | coords.reserve(NDIM);
|
---|
48 | index.reserve(NDIM);
|
---|
49 |
|
---|
50 | // observable stuff
|
---|
51 | Channels *OurChannel = new Channels;
|
---|
52 | NotificationChannels.insert( std::make_pair(this, OurChannel) );
|
---|
53 | // add instance for each notification type
|
---|
54 | for (size_t type = 0; type < NotificationType_MAX; ++type)
|
---|
55 | OurChannel->addChannel(type);
|
---|
56 |
|
---|
57 | M->setIdentity();
|
---|
58 | Minv->setIdentity();
|
---|
59 | conditions.resize(3);
|
---|
60 | conditions[0] = conditions[1] = conditions[2] = Wrap;
|
---|
61 | }
|
---|
62 |
|
---|
63 | Box::Box(const Box& src) :
|
---|
64 | Observable("Box"),
|
---|
65 | conditions(src.conditions),
|
---|
66 | M(new RealSpaceMatrix(*src.M)),
|
---|
67 | Minv(new RealSpaceMatrix(*src.Minv))
|
---|
68 | {
|
---|
69 | internal_list.reserve(pow(3,3));
|
---|
70 | coords.reserve(NDIM);
|
---|
71 | index.reserve(NDIM);
|
---|
72 |
|
---|
73 | // observable stuff
|
---|
74 | Channels *OurChannel = new Channels;
|
---|
75 | NotificationChannels.insert( std::make_pair(this, OurChannel) );
|
---|
76 | // add instance for each notification type
|
---|
77 | for (size_t type = 0; type < NotificationType_MAX; ++type)
|
---|
78 | OurChannel->addChannel(type);
|
---|
79 | }
|
---|
80 |
|
---|
81 | Box::Box(RealSpaceMatrix _M) :
|
---|
82 | Observable("Box"),
|
---|
83 | M(new RealSpaceMatrix(_M)),
|
---|
84 | Minv(new RealSpaceMatrix())
|
---|
85 | {
|
---|
86 | internal_list.reserve(pow(3,3));
|
---|
87 | coords.reserve(NDIM);
|
---|
88 | index.reserve(NDIM);
|
---|
89 |
|
---|
90 | // observable stuff
|
---|
91 | Channels *OurChannel = new Channels;
|
---|
92 | NotificationChannels.insert( std::make_pair(this, OurChannel) );
|
---|
93 | // add instance for each notification type
|
---|
94 | for (size_t type = 0; type < NotificationType_MAX; ++type)
|
---|
95 | OurChannel->addChannel(type);
|
---|
96 |
|
---|
97 | ASSERT(M->determinant()!=0,"Matrix in Box construction was not invertible");
|
---|
98 | *Minv = M->invert();
|
---|
99 | conditions.resize(3);
|
---|
100 | conditions[0] = conditions[1] = conditions[2] = Wrap;
|
---|
101 | }
|
---|
102 |
|
---|
103 | Box::~Box()
|
---|
104 | {
|
---|
105 | // observable stuff
|
---|
106 | std::map<Observable *, Channels*>::iterator iter = NotificationChannels.find(this);
|
---|
107 | delete iter->second;
|
---|
108 | NotificationChannels.erase(iter);
|
---|
109 |
|
---|
110 | delete M;
|
---|
111 | delete Minv;
|
---|
112 | }
|
---|
113 |
|
---|
114 | const RealSpaceMatrix &Box::getM() const{
|
---|
115 | return *M;
|
---|
116 | }
|
---|
117 | const RealSpaceMatrix &Box::getMinv() const{
|
---|
118 | return *Minv;
|
---|
119 | }
|
---|
120 |
|
---|
121 | void Box::setM(RealSpaceMatrix _M){
|
---|
122 | ASSERT(_M.determinant()!=0,"Matrix in Box construction was not invertible");
|
---|
123 | OBSERVE;
|
---|
124 | NOTIFY(MatrixChanged);
|
---|
125 | *M =_M;
|
---|
126 | *Minv = M->invert();
|
---|
127 | }
|
---|
128 |
|
---|
129 | Vector Box::translateIn(const Vector &point) const{
|
---|
130 | return (*M) * point;
|
---|
131 | }
|
---|
132 |
|
---|
133 | Vector Box::translateOut(const Vector &point) const{
|
---|
134 | return (*Minv) * point;
|
---|
135 | }
|
---|
136 |
|
---|
137 | Vector Box::WrapPeriodically(const Vector &point) const{
|
---|
138 | Vector helper = translateOut(point);
|
---|
139 | for(int i=NDIM;i--;){
|
---|
140 |
|
---|
141 | switch(conditions[i]){
|
---|
142 | case Wrap:
|
---|
143 | helper.at(i)=fmod(helper.at(i),1);
|
---|
144 | helper.at(i)+=(helper.at(i)>=0)?0:1;
|
---|
145 | break;
|
---|
146 | case Bounce:
|
---|
147 | {
|
---|
148 | // there probably is a better way to handle this...
|
---|
149 | // all the fabs and fmod modf probably makes it very slow
|
---|
150 | double intpart,fracpart;
|
---|
151 | fracpart = modf(fabs(helper.at(i)),&intpart);
|
---|
152 | helper.at(i) = fabs(fracpart-fmod(intpart,2));
|
---|
153 | }
|
---|
154 | break;
|
---|
155 | case Ignore:
|
---|
156 | break;
|
---|
157 | default:
|
---|
158 | ASSERT(0,"No default case for this");
|
---|
159 | break;
|
---|
160 | }
|
---|
161 |
|
---|
162 | }
|
---|
163 | return translateIn(helper);
|
---|
164 | }
|
---|
165 |
|
---|
166 | bool Box::isInside(const Vector &point) const
|
---|
167 | {
|
---|
168 | bool result = true;
|
---|
169 | Vector tester = translateOut(point);
|
---|
170 |
|
---|
171 | for(int i=0;i<NDIM;i++)
|
---|
172 | result = result &&
|
---|
173 | ((conditions[i] == Ignore) ||
|
---|
174 | ((tester[i] >= -MYEPSILON) &&
|
---|
175 | ((tester[i] - 1.) < MYEPSILON)));
|
---|
176 |
|
---|
177 | return result;
|
---|
178 | }
|
---|
179 |
|
---|
180 |
|
---|
181 | VECTORSET(std::vector) Box::explode(const Vector &point,int n) const{
|
---|
182 | ASSERT(isInside(point),"Exploded point not inside Box");
|
---|
183 | internal_explode(point, n);
|
---|
184 | VECTORSET(std::vector) res(internal_list);
|
---|
185 | return res;
|
---|
186 | }
|
---|
187 |
|
---|
188 | void Box::internal_explode(const Vector &point,int n) const{
|
---|
189 | internal_list.clear();
|
---|
190 | size_t list_index = 0;
|
---|
191 |
|
---|
192 | Vector translater = translateOut(point);
|
---|
193 | Vector mask; // contains the ignored coordinates
|
---|
194 |
|
---|
195 | // count the number of coordinates we need to do
|
---|
196 | int dims = 0; // number of dimensions that are not ignored
|
---|
197 | coords.clear();
|
---|
198 | index.clear();
|
---|
199 | for(int i=0;i<NDIM;++i){
|
---|
200 | if(conditions[i]==Ignore){
|
---|
201 | mask[i]=translater[i];
|
---|
202 | continue;
|
---|
203 | }
|
---|
204 | coords.push_back(i);
|
---|
205 | index.push_back(-n);
|
---|
206 | dims++;
|
---|
207 | } // there are max vectors in total we need to create
|
---|
208 | internal_list.resize(pow(2*n+1, dims));
|
---|
209 |
|
---|
210 | if(!dims){
|
---|
211 | // all boundaries are ignored
|
---|
212 | internal_list[list_index++] = point;
|
---|
213 | return;
|
---|
214 | }
|
---|
215 |
|
---|
216 | bool done = false;
|
---|
217 | while(!done){
|
---|
218 | // create this vector
|
---|
219 | Vector helper;
|
---|
220 | for(int i=0;i<dims;++i){
|
---|
221 | switch(conditions[coords[i]]){
|
---|
222 | case Wrap:
|
---|
223 | helper[coords[i]] = index[i]+translater[coords[i]];
|
---|
224 | break;
|
---|
225 | case Bounce:
|
---|
226 | {
|
---|
227 | // Bouncing the coordinate x produces the series:
|
---|
228 | // 0 -> x
|
---|
229 | // 1 -> 2-x
|
---|
230 | // 2 -> 2+x
|
---|
231 | // 3 -> 4-x
|
---|
232 | // 4 -> 4+x
|
---|
233 | // the first number is the next bigger even number (n+n%2)
|
---|
234 | // the next number is the value with alternating sign (x-2*(n%2)*x)
|
---|
235 | // the negative numbers produce the same sequence reversed and shifted
|
---|
236 | int n = abs(index[i]) + ((index[i]<0)?-1:0);
|
---|
237 | int sign = (index[i]<0)?-1:+1;
|
---|
238 | int even = n%2;
|
---|
239 | helper[coords[i]]=n+even+translater[coords[i]]-2*even*translater[coords[i]];
|
---|
240 | helper[coords[i]]*=sign;
|
---|
241 | }
|
---|
242 | break;
|
---|
243 | case Ignore:
|
---|
244 | ASSERT(0,"Ignored coordinate handled in generation loop");
|
---|
245 | break;
|
---|
246 | default:
|
---|
247 | ASSERT(0,"No default case for this switch-case");
|
---|
248 | break;
|
---|
249 | }
|
---|
250 |
|
---|
251 | }
|
---|
252 | // add back all ignored coordinates (not handled in above loop)
|
---|
253 | helper+=mask;
|
---|
254 | ASSERT(list_index < internal_list.size(),
|
---|
255 | "Box::internal_explode() - we have estimated the number of vectors wrong: "
|
---|
256 | +toString(list_index) +" >= "+toString(internal_list.size())+".");
|
---|
257 | internal_list[list_index++] = translateIn(helper);
|
---|
258 | // set the new indexes
|
---|
259 | int pos=0;
|
---|
260 | ++index[pos];
|
---|
261 | while(index[pos]>n){
|
---|
262 | index[pos++]=-n;
|
---|
263 | if(pos>=dims) { // it's trying to increase one beyond array... all vectors generated
|
---|
264 | done = true;
|
---|
265 | break;
|
---|
266 | }
|
---|
267 | ++index[pos];
|
---|
268 | }
|
---|
269 | }
|
---|
270 | }
|
---|
271 |
|
---|
272 | VECTORSET(std::vector) Box::explode(const Vector &point) const{
|
---|
273 | ASSERT(isInside(point),"Exploded point not inside Box");
|
---|
274 | return explode(point,1);
|
---|
275 | }
|
---|
276 |
|
---|
277 | double Box::periodicDistanceSquared(const Vector &point1,const Vector &point2) const{
|
---|
278 | Vector helper1(!isInside(point1) ? WrapPeriodically(point1) : point1);
|
---|
279 | Vector helper2(!isInside(point2) ? WrapPeriodically(point2) : point2);
|
---|
280 | internal_explode(helper1,1);
|
---|
281 | double res = internal_list.minDistSquared(helper2);
|
---|
282 | return res;
|
---|
283 | }
|
---|
284 |
|
---|
285 | double Box::periodicDistance(const Vector &point1,const Vector &point2) const{
|
---|
286 | double res;
|
---|
287 | res = sqrt(periodicDistanceSquared(point1,point2));
|
---|
288 | return res;
|
---|
289 | }
|
---|
290 |
|
---|
291 | double Box::DistanceToBoundary(const Vector &point) const
|
---|
292 | {
|
---|
293 | std::map<double, Plane> DistanceSet;
|
---|
294 | std::vector<std::pair<Plane,Plane> > Boundaries = getBoundingPlanes();
|
---|
295 | for (int i=0;i<NDIM;++i) {
|
---|
296 | const double tempres1 = Boundaries[i].first.distance(point);
|
---|
297 | const double tempres2 = Boundaries[i].second.distance(point);
|
---|
298 | DistanceSet.insert( make_pair(tempres1, Boundaries[i].first) );
|
---|
299 | LOG(1, "Inserting distance " << tempres1 << " and " << tempres2 << ".");
|
---|
300 | DistanceSet.insert( make_pair(tempres2, Boundaries[i].second) );
|
---|
301 | }
|
---|
302 | ASSERT(!DistanceSet.empty(), "Box::DistanceToBoundary() - no distances in map!");
|
---|
303 | return (DistanceSet.begin())->first;
|
---|
304 | }
|
---|
305 |
|
---|
306 | Shape Box::getShape() const{
|
---|
307 | return transform(Cuboid(Vector(0,0,0),Vector(1,1,1)),(*M));
|
---|
308 | }
|
---|
309 |
|
---|
310 | const Box::Conditions_t Box::getConditions() const
|
---|
311 | {
|
---|
312 | return conditions;
|
---|
313 | }
|
---|
314 |
|
---|
315 | void Box::setCondition(int i,Box::BoundaryCondition_t condition)
|
---|
316 | {
|
---|
317 | OBSERVE;
|
---|
318 | NOTIFY(BoundaryConditionsChanged);
|
---|
319 | conditions[i]=condition;
|
---|
320 | }
|
---|
321 |
|
---|
322 | const std::vector<std::pair<Plane,Plane> > Box::getBoundingPlanes() const
|
---|
323 | {
|
---|
324 | std::vector<std::pair<Plane,Plane> > res;
|
---|
325 | for(int i=0;i<NDIM;++i){
|
---|
326 | Vector base1,base2,base3;
|
---|
327 | base2[(i+1)%NDIM] = 1.;
|
---|
328 | base3[(i+2)%NDIM] = 1.;
|
---|
329 | Plane p1(translateIn(base1),
|
---|
330 | translateIn(base2),
|
---|
331 | translateIn(base3));
|
---|
332 | Vector offset;
|
---|
333 | offset[i]=1;
|
---|
334 | Plane p2(translateIn(base1+offset),
|
---|
335 | translateIn(base2+offset),
|
---|
336 | translateIn(base3+offset));
|
---|
337 | res.push_back(make_pair(p1,p2));
|
---|
338 | }
|
---|
339 | ASSERT(res.size() == 3, "Box::getBoundingPlanes() - does not have three plane pairs!");
|
---|
340 | return res;
|
---|
341 | }
|
---|
342 |
|
---|
343 | void Box::setCuboid(const Vector &endpoint)
|
---|
344 | {
|
---|
345 | OBSERVE;
|
---|
346 | NOTIFY(MatrixChanged);
|
---|
347 | ASSERT(endpoint[0]>0 && endpoint[1]>0 && endpoint[2]>0,"Vector does not define a full cuboid");
|
---|
348 | M->setIdentity();
|
---|
349 | M->diagonal()=endpoint;
|
---|
350 | Vector &dinv = Minv->diagonal();
|
---|
351 | for(int i=NDIM;i--;)
|
---|
352 | dinv[i]=1/endpoint[i];
|
---|
353 | }
|
---|
354 |
|
---|
355 | Box &Box::operator=(const Box &src)
|
---|
356 | {
|
---|
357 | if(&src!=this){
|
---|
358 | OBSERVE;
|
---|
359 | // new matrix
|
---|
360 | NOTIFY(MatrixChanged);
|
---|
361 | delete M;
|
---|
362 | delete Minv;
|
---|
363 | M = new RealSpaceMatrix(*src.M);
|
---|
364 | Minv = new RealSpaceMatrix(*src.Minv);
|
---|
365 | // new boundary conditions
|
---|
366 | NOTIFY(BoundaryConditionsChanged);
|
---|
367 | conditions = src.conditions;
|
---|
368 | }
|
---|
369 | return *this;
|
---|
370 | }
|
---|
371 |
|
---|
372 | Box &Box::operator=(const RealSpaceMatrix &mat)
|
---|
373 | {
|
---|
374 | OBSERVE;
|
---|
375 | NOTIFY(MatrixChanged);
|
---|
376 | setM(mat);
|
---|
377 | return *this;
|
---|
378 | }
|
---|
379 |
|
---|
380 | std::ostream & operator << (std::ostream& ost, const Box &m)
|
---|
381 | {
|
---|
382 | ost << m.getM();
|
---|
383 | return ost;
|
---|
384 | }
|
---|