1 | /*
|
---|
2 | * Box.cpp
|
---|
3 | *
|
---|
4 | * Created on: Jun 30, 2010
|
---|
5 | * Author: crueger
|
---|
6 | */
|
---|
7 |
|
---|
8 | //#include "Helpers/MemDebug.hpp"
|
---|
9 |
|
---|
10 | #include "Box.hpp"
|
---|
11 |
|
---|
12 | #include <cmath>
|
---|
13 |
|
---|
14 | #include "Matrix.hpp"
|
---|
15 | #include "vector.hpp"
|
---|
16 |
|
---|
17 | #include "Helpers/Assert.hpp"
|
---|
18 |
|
---|
19 | Box::Box()
|
---|
20 | {
|
---|
21 | M= new Matrix();
|
---|
22 | M->one();
|
---|
23 | Minv = new Matrix();
|
---|
24 | Minv->one();
|
---|
25 | }
|
---|
26 |
|
---|
27 | Box::Box(const Box& src){
|
---|
28 | M=new Matrix(*src.M);
|
---|
29 | Minv = new Matrix(*src.Minv);
|
---|
30 | }
|
---|
31 |
|
---|
32 | Box::~Box()
|
---|
33 | {
|
---|
34 | delete M;
|
---|
35 | delete Minv;
|
---|
36 | }
|
---|
37 |
|
---|
38 | const Matrix &Box::getM() const{
|
---|
39 | return *M;
|
---|
40 | }
|
---|
41 | const Matrix &Box::getMinv() const{
|
---|
42 | return *Minv;
|
---|
43 | }
|
---|
44 |
|
---|
45 | void Box::setM(Matrix _M){
|
---|
46 | ASSERT(_M.at(1,0)==_M.at(0,1),"Matrix used as cell_size was not symmetric");
|
---|
47 | ASSERT(_M.at(2,0)==_M.at(0,2),"Matrix used as cell_size was not symmetric");
|
---|
48 | ASSERT(_M.at(1,2)==_M.at(2,1),"Matrix used as cell_size was not symmetric");
|
---|
49 | *M =_M;
|
---|
50 | *Minv = M->invert();
|
---|
51 | }
|
---|
52 |
|
---|
53 | Vector Box::translateIn(const Vector &point) const{
|
---|
54 | return (*M) * point;
|
---|
55 | }
|
---|
56 |
|
---|
57 | Vector Box::translateOut(const Vector &point) const{
|
---|
58 | return (*Minv) * point;
|
---|
59 | }
|
---|
60 |
|
---|
61 | Vector Box::WrapPeriodically(const Vector &point) const{
|
---|
62 | Vector helper = translateOut(point);
|
---|
63 | for(int i=NDIM;i--;){
|
---|
64 | double intpart,fracpart;
|
---|
65 | fracpart = modf(helper.at(i),&intpart);
|
---|
66 | if(fracpart<0.)
|
---|
67 | fracpart+=1.;
|
---|
68 | helper.at(i)=fracpart;
|
---|
69 | }
|
---|
70 | return translateIn(helper);
|
---|
71 | }
|
---|
72 |
|
---|
73 | VECTORSET(std::list) Box::explode(const Vector &point) const{
|
---|
74 | VECTORSET(std::list) res;
|
---|
75 |
|
---|
76 | // translate the Vector into each of the 27 neighbourhoods
|
---|
77 |
|
---|
78 | // first create 27 translation Vectors
|
---|
79 | int N[NDIM];
|
---|
80 | for (N[0]=-1;N[0]<=1;N[0]++){
|
---|
81 | for (N[1]=-1;N[1]<=1;N[1]++){
|
---|
82 | for (N[2]=-1;N[2]<=1;N[2]++){
|
---|
83 | Vector translation = translateIn(Vector(N[0],N[1],N[2]));
|
---|
84 | res.push_back(translation);
|
---|
85 | }
|
---|
86 | }
|
---|
87 | }
|
---|
88 | // translate all the translation vector by the offset defined by point
|
---|
89 | res.translate(point);
|
---|
90 | return res;
|
---|
91 | }
|
---|
92 |
|
---|
93 | double Box::periodicDistanceSquared(const Vector &point1,const Vector &point2) const{
|
---|
94 | VECTORSET(std::list) expansion = explode(point1);
|
---|
95 | double res = expansion.minDistSquared(point2);
|
---|
96 | return res;
|
---|
97 | }
|
---|
98 |
|
---|
99 | double Box::periodicDistance(const Vector &point1,const Vector &point2) const{
|
---|
100 | double res;
|
---|
101 | res = sqrt(periodicDistanceSquared(point1,point2));
|
---|
102 | return res;
|
---|
103 | }
|
---|
104 |
|
---|
105 | Box &Box::operator=(const Box &src){
|
---|
106 | if(&src!=this){
|
---|
107 | delete M;
|
---|
108 | delete Minv;
|
---|
109 | M = new Matrix(*src.M);
|
---|
110 | Minv = new Matrix(*src.Minv);
|
---|
111 | }
|
---|
112 | return *this;
|
---|
113 | }
|
---|
114 |
|
---|
115 | Box &Box::operator=(const Matrix &mat){
|
---|
116 | setM(mat);
|
---|
117 | return *this;
|
---|
118 | }
|
---|