1 | /*
|
---|
2 | * Project: MoleCuilder
|
---|
3 | * Description: creates and alters molecular systems
|
---|
4 | * Copyright (C) 2010 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 "Helpers/MemDebug.hpp"
|
---|
21 |
|
---|
22 | #include "Box.hpp"
|
---|
23 |
|
---|
24 | #include <cmath>
|
---|
25 | #include <iostream>
|
---|
26 | #include <cstdlib>
|
---|
27 |
|
---|
28 | #include "LinearAlgebra/Matrix.hpp"
|
---|
29 | #include "LinearAlgebra/Vector.hpp"
|
---|
30 | #include "LinearAlgebra/Plane.hpp"
|
---|
31 | #include "Shapes/BaseShapes.hpp"
|
---|
32 | #include "Shapes/ShapeOps.hpp"
|
---|
33 |
|
---|
34 | #include "Helpers/Assert.hpp"
|
---|
35 |
|
---|
36 | using namespace std;
|
---|
37 |
|
---|
38 | Box::Box()
|
---|
39 | {
|
---|
40 | M= new Matrix();
|
---|
41 | M->one();
|
---|
42 | Minv = new Matrix();
|
---|
43 | Minv->one();
|
---|
44 | conditions.resize(3);
|
---|
45 | conditions[0] = conditions[1] = conditions[2] = Wrap;
|
---|
46 | }
|
---|
47 |
|
---|
48 | Box::Box(const Box& src){
|
---|
49 | M=new Matrix(*src.M);
|
---|
50 | Minv = new Matrix(*src.Minv);
|
---|
51 | conditions = src.conditions;
|
---|
52 | }
|
---|
53 |
|
---|
54 | Box::~Box()
|
---|
55 | {
|
---|
56 | delete M;
|
---|
57 | delete Minv;
|
---|
58 | }
|
---|
59 |
|
---|
60 | const Matrix &Box::getM() const{
|
---|
61 | return *M;
|
---|
62 | }
|
---|
63 | const Matrix &Box::getMinv() const{
|
---|
64 | return *Minv;
|
---|
65 | }
|
---|
66 |
|
---|
67 | void Box::setM(Matrix _M){
|
---|
68 | ASSERT(_M.determinant()!=0,"Matrix in Box construction was not invertible");
|
---|
69 | *M =_M;
|
---|
70 | *Minv = M->invert();
|
---|
71 | }
|
---|
72 |
|
---|
73 | Vector Box::translateIn(const Vector &point) const{
|
---|
74 | return (*M) * point;
|
---|
75 | }
|
---|
76 |
|
---|
77 | Vector Box::translateOut(const Vector &point) const{
|
---|
78 | return (*Minv) * point;
|
---|
79 | }
|
---|
80 |
|
---|
81 | Vector Box::WrapPeriodically(const Vector &point) const{
|
---|
82 | Vector helper = translateOut(point);
|
---|
83 | for(int i=NDIM;i--;){
|
---|
84 |
|
---|
85 | switch(conditions[i]){
|
---|
86 | case Wrap:
|
---|
87 | helper.at(i)=fmod(helper.at(i),1);
|
---|
88 | helper.at(i)+=(helper.at(i)>=0)?0:1;
|
---|
89 | break;
|
---|
90 | case Bounce:
|
---|
91 | {
|
---|
92 | // there probably is a better way to handle this...
|
---|
93 | // all the fabs and fmod modf probably makes it very slow
|
---|
94 | double intpart,fracpart;
|
---|
95 | fracpart = modf(fabs(helper.at(i)),&intpart);
|
---|
96 | helper.at(i) = fabs(fracpart-fmod(intpart,2));
|
---|
97 | }
|
---|
98 | break;
|
---|
99 | case Ignore:
|
---|
100 | break;
|
---|
101 | default:
|
---|
102 | ASSERT(0,"No default case for this");
|
---|
103 | }
|
---|
104 |
|
---|
105 | }
|
---|
106 | return translateIn(helper);
|
---|
107 | }
|
---|
108 |
|
---|
109 | bool Box::isInside(const Vector &point) const
|
---|
110 | {
|
---|
111 | bool result = true;
|
---|
112 | Vector tester = translateOut(point);
|
---|
113 |
|
---|
114 | for(int i=0;i<NDIM;i++)
|
---|
115 | result = result &&
|
---|
116 | ((conditions[i] == Ignore) ||
|
---|
117 | ((tester[i] >= -MYEPSILON) &&
|
---|
118 | ((tester[i] - 1.) < MYEPSILON)));
|
---|
119 |
|
---|
120 | return result;
|
---|
121 | }
|
---|
122 |
|
---|
123 |
|
---|
124 | VECTORSET(std::list) Box::explode(const Vector &point,int n) const{
|
---|
125 | ASSERT(isInside(point),"Exploded point not inside Box");
|
---|
126 | VECTORSET(std::list) res;
|
---|
127 |
|
---|
128 | Vector translater = translateOut(point);
|
---|
129 | Vector mask; // contains the ignored coordinates
|
---|
130 |
|
---|
131 | // count the number of coordinates we need to do
|
---|
132 | int dims = 0; // number of dimensions that are not ignored
|
---|
133 | vector<int> coords;
|
---|
134 | vector<int> index;
|
---|
135 | for(int i=0;i<NDIM;++i){
|
---|
136 | if(conditions[i]==Ignore){
|
---|
137 | mask[i]=translater[i];
|
---|
138 | continue;
|
---|
139 | }
|
---|
140 | coords.push_back(i);
|
---|
141 | index.push_back(-n);
|
---|
142 | dims++;
|
---|
143 | } // there are max vectors in total we need to create
|
---|
144 |
|
---|
145 | if(!dims){
|
---|
146 | // all boundaries are ignored
|
---|
147 | res.push_back(point);
|
---|
148 | return res;
|
---|
149 | }
|
---|
150 |
|
---|
151 | bool done = false;
|
---|
152 | while(!done){
|
---|
153 | // create this vector
|
---|
154 | Vector helper;
|
---|
155 | for(int i=0;i<dims;++i){
|
---|
156 | switch(conditions[coords[i]]){
|
---|
157 | case Wrap:
|
---|
158 | helper[coords[i]] = index[i]+translater[coords[i]];
|
---|
159 | break;
|
---|
160 | case Bounce:
|
---|
161 | {
|
---|
162 | // Bouncing the coordinate x produces the series:
|
---|
163 | // 0 -> x
|
---|
164 | // 1 -> 2-x
|
---|
165 | // 2 -> 2+x
|
---|
166 | // 3 -> 4-x
|
---|
167 | // 4 -> 4+x
|
---|
168 | // the first number is the next bigger even number (n+n%2)
|
---|
169 | // the next number is the value with alternating sign (x-2*(n%2)*x)
|
---|
170 | // the negative numbers produce the same sequence reversed and shifted
|
---|
171 | int n = abs(index[i]) + ((index[i]<0)?-1:0);
|
---|
172 | int sign = (index[i]<0)?-1:+1;
|
---|
173 | int even = n%2;
|
---|
174 | helper[coords[i]]=n+even+translater[coords[i]]-2*even*translater[coords[i]];
|
---|
175 | helper[coords[i]]*=sign;
|
---|
176 | }
|
---|
177 | break;
|
---|
178 | case Ignore:
|
---|
179 | ASSERT(0,"Ignored coordinate handled in generation loop");
|
---|
180 | default:
|
---|
181 | ASSERT(0,"No default case for this switch-case");
|
---|
182 | }
|
---|
183 |
|
---|
184 | }
|
---|
185 | // add back all ignored coordinates (not handled in above loop)
|
---|
186 | helper+=mask;
|
---|
187 | res.push_back(translateIn(helper));
|
---|
188 | // set the new indexes
|
---|
189 | int pos=0;
|
---|
190 | ++index[pos];
|
---|
191 | while(index[pos]>n){
|
---|
192 | index[pos++]=-n;
|
---|
193 | if(pos>=dims) { // it's trying to increase one beyond array... all vectors generated
|
---|
194 | done = true;
|
---|
195 | break;
|
---|
196 | }
|
---|
197 | ++index[pos];
|
---|
198 | }
|
---|
199 | }
|
---|
200 | return res;
|
---|
201 | }
|
---|
202 |
|
---|
203 | VECTORSET(std::list) Box::explode(const Vector &point) const{
|
---|
204 | ASSERT(isInside(point),"Exploded point not inside Box");
|
---|
205 | return explode(point,1);
|
---|
206 | }
|
---|
207 |
|
---|
208 | double Box::periodicDistanceSquared(const Vector &point1,const Vector &point2) const{
|
---|
209 | Vector helper1 = WrapPeriodically(point1);
|
---|
210 | Vector helper2 = WrapPeriodically(point2);
|
---|
211 | VECTORSET(std::list) expansion = explode(helper1);
|
---|
212 | double res = expansion.minDistSquared(helper2);
|
---|
213 | return res;
|
---|
214 | }
|
---|
215 |
|
---|
216 | double Box::periodicDistance(const Vector &point1,const Vector &point2) const{
|
---|
217 | double res;
|
---|
218 | res = sqrt(periodicDistanceSquared(point1,point2));
|
---|
219 | return res;
|
---|
220 | }
|
---|
221 |
|
---|
222 | Shape Box::getShape() const{
|
---|
223 | return transform(Cuboid(Vector(0,0,0),Vector(1,1,1)),(*M));
|
---|
224 | }
|
---|
225 |
|
---|
226 | const Box::Conditions_t Box::getConditions(){
|
---|
227 | return conditions;
|
---|
228 | }
|
---|
229 |
|
---|
230 | void Box::setCondition(int i,Box::BoundaryCondition_t condition){
|
---|
231 | conditions[i]=condition;
|
---|
232 | }
|
---|
233 |
|
---|
234 | const vector<pair<Plane,Plane> > Box::getBoundingPlanes(){
|
---|
235 | vector<pair<Plane,Plane> > res;
|
---|
236 | for(int i=0;i<NDIM;++i){
|
---|
237 | Vector base1,base2,base3;
|
---|
238 | base2[(i+1)%NDIM] = 1.;
|
---|
239 | base3[(i+2)%NDIM] = 1.;
|
---|
240 | Plane p1(translateIn(base1),
|
---|
241 | translateIn(base2),
|
---|
242 | translateIn(base3));
|
---|
243 | Vector offset;
|
---|
244 | offset[i]=1;
|
---|
245 | Plane p2(translateIn(base1+offset),
|
---|
246 | translateIn(base2+offset),
|
---|
247 | translateIn(base3+offset));
|
---|
248 | res.push_back(make_pair(p1,p2));
|
---|
249 | }
|
---|
250 | return res;
|
---|
251 | }
|
---|
252 |
|
---|
253 | void Box::setCuboid(const Vector &endpoint){
|
---|
254 | ASSERT(endpoint[0]>0 && endpoint[1]>0 && endpoint[2]>0,"Vector does not define a full cuboid");
|
---|
255 | M->one();
|
---|
256 | M->diagonal()=endpoint;
|
---|
257 | Vector &dinv = Minv->diagonal();
|
---|
258 | for(int i=NDIM;i--;)
|
---|
259 | dinv[i]=1/endpoint[i];
|
---|
260 | }
|
---|
261 |
|
---|
262 | Box &Box::operator=(const Box &src){
|
---|
263 | if(&src!=this){
|
---|
264 | delete M;
|
---|
265 | delete Minv;
|
---|
266 | M = new Matrix(*src.M);
|
---|
267 | Minv = new Matrix(*src.Minv);
|
---|
268 | conditions = src.conditions;
|
---|
269 | }
|
---|
270 | return *this;
|
---|
271 | }
|
---|
272 |
|
---|
273 | Box &Box::operator=(const Matrix &mat){
|
---|
274 | setM(mat);
|
---|
275 | return *this;
|
---|
276 | }
|
---|