Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Sbalzarini Lab
S
Software
P
Parallel Computing
OpenFPM
openfpm_pdata
Commits
ca16d911
Commit
ca16d911
authored
Apr 15, 2016
by
Pietro Incardona
Browse files
Fixing getPos
parent
5af36978
Changes
16
Hide whitespace changes
Inline
Side-by-side
CHANGELOG.md
View file @
ca16d911
...
...
@@ -4,7 +4,7 @@ All notable changes to this project will be documented in this file.
## [0.3.0] -
### Added
-
Nothing to report
-
Molacular Dynamic example
### Fixed
-
Nothing to report
...
...
@@ -13,6 +13,8 @@ All notable changes to this project will be documented in this file.
-
Eliminated global_v_cluster, init_global_v_cluster, delete_global_v_cluster,
substituted by
create_vcluster, openfpm_init, openfpm_delete
-
CartDecomposition parameter for the distributed structures is now optional
-
template getPos
<
0
>
(), substituted by getPos()
## [0.2.1] -
...
...
example/Numerics/PSE/0_Derivative_approx_1D/main.cpp
View file @
ca16d911
...
...
@@ -127,7 +127,7 @@ int main(int argc, char* argv[])
auto
key
=
it2
.
get
();
// set the position of the particles
vd
.
template
getPos
<
0
>
(
key
)[
0
]
=
(
key
.
getKey
()
+
base
)
*
spacing
;
vd
.
getPos
(
key
)[
0
]
=
(
key
.
getKey
()
+
base
)
*
spacing
;
//set the property of the particles
vd
.
template
getProp
<
0
>(
key
)
=
f_xex2
((
key
.
getKey
()
+
base
)
*
spacing
);
...
...
@@ -185,19 +185,19 @@ int main(int argc, char* argv[])
auto
key
=
it
.
get
();
// set the position of the particles
if
(
m_pad
.
isInsideNB
(
vd
.
template
getPos
<
0
>
(
key
))
==
true
)
if
(
m_pad
.
isInsideNB
(
vd
.
getPos
(
key
))
==
true
)
{
vd
.
add
();
vd
.
template
getLastPos
<
0
>
()[
0
]
=
-
vd
.
template
getPos
<
0
>
(
key
)[
0
];
vd
.
getLastPos
()[
0
]
=
-
vd
.
getPos
(
key
)[
0
];
vd
.
template
getLastProp
<
0
>()
=
-
vd
.
template
getProp
<
0
>(
key
);
}
// set the position of the particles
if
(
m_pad2
.
isInsideNB
(
vd
.
template
getPos
<
0
>
(
key
))
==
true
)
if
(
m_pad2
.
isInsideNB
(
vd
.
getPos
(
key
))
==
true
)
{
vd
.
add
();
vd
.
template
getLastPos
<
0
>
()[
0
]
=
2.0
*
box
.
getHigh
(
0
)
-
vd
.
template
getPos
<
0
>
(
key
)[
0
];
vd
.
template
getLastProp
<
0
>()
=
f_xex2
(
vd
.
template
getLastPos
<
0
>
()[
0
]);
vd
.
getLastPos
()[
0
]
=
2.0
*
box
.
getHigh
(
0
)
-
vd
.
getPos
(
key
)[
0
];
vd
.
template
getLastProp
<
0
>()
=
f_xex2
(
vd
.
getLastPos
()[
0
]);
}
++
it
;
...
...
@@ -237,7 +237,7 @@ int main(int argc, char* argv[])
vect_dist_key_dx
key
=
it_p
.
get
();
// Get the position of the particles
Point
<
1
,
double
>
p
=
vd
.
template
getPos
<
0
>
(
key
);
Point
<
1
,
double
>
p
=
vd
.
getPos
(
key
);
// We are not interested in calculating out the domain
// note added padding particle are considered domain particles
...
...
@@ -261,7 +261,7 @@ int main(int argc, char* argv[])
if
(
nnp
!=
key
.
getKey
())
{
// W(x-y)
double
ker
=
lker
.
value
(
p
,
vd
.
template
getPos
<
0
>
(
nnp
));
double
ker
=
lker
.
value
(
p
,
vd
.
getPos
(
nnp
));
// f(y)
double
prp_y
=
vd
.
template
getProp
<
0
>(
nnp
);
...
...
example/Numerics/PSE/1_Derivative_approx_1D_mp/main_float128.cpp
View file @
ca16d911
...
...
@@ -130,7 +130,7 @@ int main(int argc, char* argv[])
auto
key
=
it2
.
get
();
// set the position of the particles
vd
.
template
getPos
<
0
>
(
key
)[
0
]
=
(
key
.
getKey
()
+
base
)
*
spacing
;
vd
.
getPos
(
key
)[
0
]
=
(
key
.
getKey
()
+
base
)
*
spacing
;
//set the property of the particles
vd
.
template
getProp
<
0
>(
key
)
=
f_xex2
((
key
.
getKey
()
+
base
)
*
spacing
);
...
...
@@ -188,19 +188,19 @@ int main(int argc, char* argv[])
auto
key
=
it
.
get
();
// set the position of the particles
if
(
m_pad
.
isInsideNB
(
vd
.
template
getPos
<
0
>
(
key
))
==
true
)
if
(
m_pad
.
isInsideNB
(
vd
.
getPos
(
key
))
==
true
)
{
vd
.
add
();
vd
.
template
getLastPos
<
0
>
()[
0
]
=
-
vd
.
template
getPos
<
0
>
(
key
)[
0
];
vd
.
getLastPos
()[
0
]
=
-
vd
.
getPos
(
key
)[
0
];
vd
.
template
getLastProp
<
0
>()
=
-
vd
.
template
getProp
<
0
>(
key
);
}
// set the position of the particles
if
(
m_pad2
.
isInsideNB
(
vd
.
template
getPos
<
0
>
(
key
))
==
true
)
if
(
m_pad2
.
isInsideNB
(
vd
.
getPos
(
key
))
==
true
)
{
vd
.
add
();
vd
.
template
getLastPos
<
0
>
()[
0
]
=
2.0
*
box
.
getHigh
(
0
)
-
vd
.
template
getPos
<
0
>
(
key
)[
0
];
vd
.
template
getLastProp
<
0
>()
=
f_xex2
(
vd
.
template
getLastPos
<
0
>
()[
0
]);
vd
.
getLastPos
()[
0
]
=
2.0
*
box
.
getHigh
(
0
)
-
vd
.
getPos
(
key
)[
0
];
vd
.
template
getLastProp
<
0
>()
=
f_xex2
(
vd
.
getLastPos
()[
0
]);
}
++
it
;
...
...
@@ -241,7 +241,7 @@ int main(int argc, char* argv[])
vect_dist_key_dx
key
=
it_p
.
get
();
// Get the position of the particles
Point
<
1
,
float128
>
p
=
vd
.
template
getPos
<
0
>
(
key
);
Point
<
1
,
float128
>
p
=
vd
.
getPos
(
key
);
// We are not interested in calculating out the domain
// note added padding particle are considered domain particles
...
...
@@ -265,7 +265,7 @@ int main(int argc, char* argv[])
if
(
nnp
!=
key
.
getKey
())
{
// W(x-y)
float128
ker
=
lker
.
value
(
p
,
vd
.
template
getPos
<
0
>
(
nnp
));
float128
ker
=
lker
.
value
(
p
,
vd
.
getPos
(
nnp
));
// f(y)
float128
prp_y
=
vd
.
template
getProp
<
0
>(
nnp
);
...
...
example/Numerics/PSE/1_Diffusion_1D/main.cpp
View file @
ca16d911
...
...
@@ -74,7 +74,7 @@ template<typename CellL> double calcLap(Point<1,double> p, vect_dist_key_dx key,
if
(
nnp
!=
key
.
getKey
())
{
// W(x-y)
double
ker
=
lker
.
value
(
p
,
vd
.
template
getPos
<
0
>
(
nnp
));
double
ker
=
lker
.
value
(
p
,
vd
.
getPos
(
nnp
));
// f(y)
double
prp_y
=
vd
.
template
getProp
<
0
>(
nnp
);
...
...
@@ -180,7 +180,7 @@ int main(int argc, char* argv[])
auto
key
=
it2
.
get
();
// set the position of the particles
vd
.
template
getPos
<
0
>
(
key
)[
0
]
=
(
key
.
getKey
()
+
base
)
*
spacing
;
vd
.
getPos
(
key
)[
0
]
=
(
key
.
getKey
()
+
base
)
*
spacing
;
//set the property of the particles
vd
.
template
getProp
<
0
>(
key
)
=
f_xex2
((
key
.
getKey
()
+
base
)
*
spacing
);
...
...
@@ -238,19 +238,19 @@ int main(int argc, char* argv[])
auto
key
=
it
.
get
();
// set the position of the particles
if
(
m_pad
.
isInsideNB
(
vd
.
template
getPos
<
0
>
(
key
))
==
true
)
if
(
m_pad
.
isInsideNB
(
vd
.
getPos
(
key
))
==
true
)
{
vd
.
add
();
vd
.
template
getLastPos
<
0
>
()[
0
]
=
-
vd
.
template
getPos
<
0
>
(
key
)[
0
];
vd
.
getLastPos
()[
0
]
=
-
vd
.
getPos
(
key
)[
0
];
vd
.
template
getLastProp
<
0
>()
=
-
vd
.
template
getProp
<
0
>(
key
);
}
// set the position of the particles
if
(
m_pad2
.
isInsideNB
(
vd
.
template
getPos
<
0
>
(
key
))
==
true
)
if
(
m_pad2
.
isInsideNB
(
vd
.
getPos
(
key
))
==
true
)
{
vd
.
add
();
vd
.
template
getLastPos
<
0
>
()[
0
]
=
2.0
*
box
.
getHigh
(
0
)
-
vd
.
template
getPos
<
0
>
(
key
)[
0
];
vd
.
template
getLastProp
<
0
>()
=
f_xex2
(
vd
.
template
getLastPos
<
0
>
()[
0
]);
vd
.
getLastPos
()[
0
]
=
2.0
*
box
.
getHigh
(
0
)
-
vd
.
getPos
(
key
)[
0
];
vd
.
template
getLastProp
<
0
>()
=
f_xex2
(
vd
.
getLastPos
()[
0
]);
}
++
it
;
...
...
@@ -284,7 +284,7 @@ int main(int argc, char* argv[])
// key
vect_dist_key_dx
key
=
it_p
.
get
();
Point
<
1
,
double
>
p
=
vd
.
template
getPos
<
0
>
(
key
);
Point
<
1
,
double
>
p
=
vd
.
getPos
(
key
);
// We are not interested in calculating out the domain
// note added padding particle are considered domain particles
...
...
example/SE/0_classes/main.cpp
View file @
ca16d911
...
...
@@ -98,7 +98,7 @@ int main(int argc, char* argv[])
try
{
vect_dist_key_dx
vt
(
5048
);
auto
it
=
vd
.
getPos
<
0
>
(
vt
);
auto
it
=
vd
.
getPos
(
vt
);
}
catch
(
size_t
e
)
{
...
...
@@ -185,7 +185,7 @@ int main(int argc, char* argv[])
try
{
vect_dist_key_dx
vt
(
0
);
auto
it
=
vd1
->
getPos
<
0
>
(
vt
);
auto
it
=
vd1
->
getPos
(
vt
);
}
catch
(
size_t
e
)
{
...
...
example/Vector/0_simple/main.cpp
View file @
ca16d911
...
...
@@ -91,8 +91,8 @@ int main(int argc, char* argv[])
{
auto
key
=
it
.
get
();
vd
.
template
getPos
<
s
::
x
>
(
key
)[
0
]
=
ud
(
eg
);
vd
.
template
getPos
<
s
::
x
>
(
key
)[
1
]
=
ud
(
eg
);
vd
.
getPos
(
key
)[
0
]
=
ud
(
eg
);
vd
.
getPos
(
key
)[
1
]
=
ud
(
eg
);
++
it
;
}
...
...
@@ -121,7 +121,7 @@ int main(int argc, char* argv[])
auto
key
=
it
.
get
();
// The template parameter is unuseful and will probably disappear
if
(
ct
.
isLocal
(
vd
.
template
getPos
<
0
>
(
key
))
==
false
)
if
(
ct
.
isLocal
(
vd
.
getPos
(
key
))
==
false
)
std
::
cerr
<<
"Error particle is not local"
<<
"
\n
"
;
// set the all the properties to 0.0
...
...
example/Vector/1_celllist/main.cpp
View file @
ca16d911
...
...
@@ -103,9 +103,9 @@ int main(int argc, char* argv[])
auto
key
=
it
.
get
();
vd
.
template
getLastPos
<
0
>
()[
0
]
=
key
.
get
(
0
)
*
it
.
getSpacing
(
0
);
vd
.
template
getLastPos
<
0
>
()[
1
]
=
key
.
get
(
1
)
*
it
.
getSpacing
(
1
);
vd
.
template
getLastPos
<
0
>
()[
2
]
=
key
.
get
(
2
)
*
it
.
getSpacing
(
2
);
vd
.
getLastPos
()[
0
]
=
key
.
get
(
0
)
*
it
.
getSpacing
(
0
);
vd
.
getLastPos
()[
1
]
=
key
.
get
(
1
)
*
it
.
getSpacing
(
1
);
vd
.
getLastPos
()[
2
]
=
key
.
get
(
2
)
*
it
.
getSpacing
(
2
);
++
it
;
}
...
...
@@ -133,9 +133,9 @@ int main(int argc, char* argv[])
{
auto
p
=
it2
.
get
();
Point
<
3
,
float
>
xp
=
vd
.
getPos
<
0
>
(
p
);
Point
<
3
,
float
>
xp
=
vd
.
getPos
(
p
);
auto
Np
=
NN
.
getIterator
(
NN
.
getCell
(
vd
.
getPos
<
0
>
(
p
)));
auto
Np
=
NN
.
getIterator
(
NN
.
getCell
(
vd
.
getPos
(
p
)));
while
(
Np
.
isNext
())
{
...
...
@@ -143,7 +143,7 @@ int main(int argc, char* argv[])
// repulsive
Point
<
3
,
float
>
xq
=
vd
.
getPos
<
0
>
(
q
);
Point
<
3
,
float
>
xq
=
vd
.
getPos
(
q
);
Point
<
3
,
float
>
f
=
(
xp
-
xq
);
// we sum the distance of all the particles
...
...
example/Vector/1_verlet/main.cpp
View file @
ca16d911
...
...
@@ -102,9 +102,9 @@ int main(int argc, char* argv[])
auto
key
=
it
.
get
();
vd
.
template
getLastPos
<
0
>
()[
0
]
=
key
.
get
(
0
)
*
it
.
getSpacing
(
0
);
vd
.
template
getLastPos
<
0
>
()[
1
]
=
key
.
get
(
1
)
*
it
.
getSpacing
(
1
);
vd
.
template
getLastPos
<
0
>
()[
2
]
=
key
.
get
(
2
)
*
it
.
getSpacing
(
2
);
vd
.
getLastPos
()[
0
]
=
key
.
get
(
0
)
*
it
.
getSpacing
(
0
);
vd
.
getLastPos
()[
1
]
=
key
.
get
(
1
)
*
it
.
getSpacing
(
1
);
vd
.
getLastPos
()[
2
]
=
key
.
get
(
2
)
*
it
.
getSpacing
(
2
);
++
it
;
}
...
...
@@ -138,14 +138,14 @@ int main(int argc, char* argv[])
for
(
size_t
i
=
0
;
i
<
verlet
.
size
()
;
i
++
)
{
Point
<
3
,
float
>
p
=
vd
.
getPos
<
0
>
(
i
);
Point
<
3
,
float
>
p
=
vd
.
getPos
(
i
);
// for each neighborhood particle
for
(
size_t
j
=
0
;
j
<
verlet
.
get
(
i
).
size
()
;
j
++
)
{
auto
&
NN
=
verlet
.
get
(
i
);
Point
<
3
,
float
>
q
=
vd
.
getPos
<
0
>
(
NN
.
get
(
j
));
Point
<
3
,
float
>
q
=
vd
.
getPos
(
NN
.
get
(
j
));
// some non-sense calculation as usage demo
...
...
example/Vector/2_molecular_dynamic/Makefile
0 → 100644
View file @
ca16d911
include
../../example.mk
CC
=
mpic++
LDIR
=
OBJ
=
main.o
%.o
:
%.cpp
$(CC)
-O3
-g
-c
--std
=
c++11
-o
$@
$<
$(INCLUDE_PATH)
md_dyn
:
$(OBJ)
$(CC)
-o
$@
$^
$(CFLAGS)
$(LIBS_PATH)
$(LIBS)
all
:
cell
.PHONY
:
clean all
clean
:
rm
-f
*
.o
*
~ core md_dyn
example/Vector/2_molecular_dynamic/config.cfg
0 → 100644
View file @
ca16d911
[pack]
files = main.cpp Makefile
example/Vector/2_molecular_dynamic/gc_plot2_out.html
0 → 100644
View file @
ca16d911
<html>
<head>
<script
type=
"text/javascript"
src=
"https://www.gstatic.com/charts/loader.js"
></script>
<script
type=
"text/javascript"
>
google
.
charts
.
load
(
'
current
'
,
{
'
packages
'
:[
'
corechart
'
]});
google
.
charts
.
setOnLoadCallback
(
drawVisualization
);
function
drawVisualization
()
{
var
data0
=
new
google
.
visualization
.
DataTable
();
data0
.
addColumn
(
'
number
'
,
'
iteration
'
);
data0
.
addColumn
(
'
number
'
,
'
line0
'
);
data0
.
addRows
([
[
0
,
18010.3
],
[
100
,
18010.5
],
[
200
,
16751.9
],
[
300
,
23684.3
],
[
400
,
23433.6
],
[
500
,
22997.1
],
[
600
,
24500.7
],
[
700
,
23194.2
],
[
800
,
23649.7
],
[
900
,
23663.8
],
[
1000
,
23838.5
],
[
1100
,
24004.5
],
[
1200
,
23994.6
],
[
1300
,
23702.6
],
[
1400
,
23636.1
],
[
1500
,
23877.8
],
[
1600
,
24171.7
],
[
1700
,
23985
],
[
1800
,
23507
],
[
1900
,
23308.7
],
[
2000
,
23820.8
],
[
2100
,
23417.1
],
[
2200
,
23607.6
],
[
2300
,
23670.5
],
[
2400
,
23722
],
[
2500
,
23757.4
],
[
2600
,
23914.1
],
[
2700
,
23541.1
],
[
2800
,
23037.1
],
[
2900
,
23520
],
[
3000
,
24058.2
],
[
3100
,
23087.6
],
[
3200
,
22901
],
[
3300
,
23182.9
],
[
3400
,
23352.6
],
[
3500
,
22819.6
],
[
3600
,
22970.7
],
[
3700
,
22831.1
],
[
3800
,
22506.7
],
[
3900
,
22614.4
],
[
4000
,
22716.8
],
[
4100
,
22969.7
],
[
4200
,
22976.2
],
[
4300
,
22690.9
],
[
4400
,
22382
],
[
4500
,
22731.3
],
[
4600
,
22547.5
],
[
4700
,
22194
],
[
4800
,
22350.9
],
[
4900
,
22686.1
],
[
5000
,
22328.7
],
[
5100
,
22629.4
],
[
5200
,
22597.8
],
[
5300
,
22272.9
],
[
5400
,
22649.2
],
[
5500
,
22057.9
],
[
5600
,
22137.9
],
[
5700
,
22347.8
],
[
5800
,
22465.9
],
[
5900
,
22530.8
],
[
6000
,
22393.6
],
[
6100
,
22117.9
],
[
6200
,
22008.9
],
[
6300
,
21868.5
],
[
6400
,
21330.7
],
[
6500
,
21390.1
],
[
6600
,
21782.7
],
[
6700
,
21393.1
],
[
6800
,
21718.1
],
[
6900
,
21874.9
],
[
7000
,
21867.5
],
[
7100
,
21453.3
],
[
7200
,
21734.5
],
[
7300
,
21439.8
],
[
7400
,
21454.8
],
[
7500
,
21700.6
],
[
7600
,
21638.9
],
[
7700
,
20965.1
],
[
7800
,
20846.4
],
[
7900
,
20689.3
],
[
8000
,
20758.1
],
[
8100
,
21154.9
],
[
8200
,
20726.4
],
[
8300
,
20594.4
],
[
8400
,
20460.2
],
[
8500
,
20694
],
[
8600
,
20509.1
],
[
8700
,
20492.3
],
[
8800
,
20342.2
],
[
8900
,
20397.8
],
[
9000
,
20673.6
],
[
9100
,
20500.1
],
[
9200
,
20203.3
],
[
9300
,
20221.9
],
[
9400
,
20219
],
[
9500
,
20021.4
],
[
9600
,
20164
],
[
9700
,
20237.4
],
[
9800
,
20089.2
],
[
9900
,
20061.7
],
]);
var
options0
=
{
title
:
'
Energy with time
'
,
vAxis
:
{
title
:
'
Energy
'
},
hAxis
:
{
title
:
'
iteration
'
},
curveType
:
'
function
'
,
lineWidth
:
1
,
intervals
:
{
'
style
'
:
'
area
'
},
};
var
chart
=
new
google
.
visualization
.
ComboChart
(
document
.
getElementById
(
'
chart_div0
'
));
chart
.
draw
(
data0
,
options0
);
}
</script>
</head>
<body>
<div
id=
"chart_div0"
style=
"width: 900px; height: 500px;"
></div>
</body>
</html>
example/Vector/2_molecular_dynamic/main.cpp
0 → 100644
View file @
ca16d911
#include "Vector/vector_dist.hpp"
#include "Decomposition/CartDecomposition.hpp"
#include "data_type/aggregate.hpp"
#include "Plot/GoogleChart.hpp"
#include "Plot/util.hpp"
/*
* ### WIKI 1 ###
*
* ## Molecular Dynamic with Lennard-Jones potential
*
* This example show a simple Lennard-Jones molecular dynamic simulation in a stable regime
*
* ### WIKI END ###
*
*/
constexpr
int
velocity
=
0
;
constexpr
int
force
=
1
;
/* ### WIKI 11 ###
*
* The function to calculate the forces between particles. It require the vector of particles
* Cell list and scaling factor.
*
*/
void
calc_forces
(
vector_dist
<
3
,
double
,
aggregate
<
double
[
3
],
double
[
3
]
>
>
&
vd
,
CellList
<
3
,
double
,
FAST
,
shift
<
3
,
double
>
>
&
NN
,
double
L
)
{
// ### WIKI 12 ###
//
// Update the cell list from the actual particle configuration
//
//
vd
.
updateCellList
(
NN
);
// ### WIKI 13 ###
//
// Calculate the forces
//
auto
it2
=
vd
.
getDomainIterator
();
while
(
it2
.
isNext
())
{
auto
p
=
it2
.
get
();
Point
<
3
,
double
>
xp
=
vd
.
getPos
<
0
>
(
p
);
vd
.
template
getProp
<
force
>(
p
)[
0
]
=
0.0
;
vd
.
template
getProp
<
force
>(
p
)[
1
]
=
0.0
;
vd
.
template
getProp
<
force
>(
p
)[
2
]
=
0.0
;
// For each neighborhood particle
auto
Np
=
NN
.
template
getNNIterator
<
NO_CHECK
>(
NN
.
getCell
(
vd
.
getPos
<
0
>
(
p
)));
while
(
Np
.
isNext
())
{
// Neighborhood particle q
auto
q
=
Np
.
get
();
if
(
q
==
p
.
getKey
())
{
++
Np
;
continue
;};
// repulsive
Point
<
3
,
double
>
xq
=
vd
.
getPos
<
0
>
(
q
);
Point
<
3
,
double
>
r
=
xp
-
xq
;
// take the norm, normalize
float
rn
=
r
.
norm
();
r
/=
rn
;
rn
*=
L
;
// Calculate the force, using pow is slower
Point
<
3
,
double
>
f
=
24.0
*
(
2.0
/
(
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
)
-
1.0
/
(
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
))
*
r
;
// we sum the force produced by q on p
vd
.
template
getProp
<
force
>(
p
)[
0
]
+=
f
.
get
(
0
);
vd
.
template
getProp
<
force
>(
p
)[
1
]
+=
f
.
get
(
1
);
vd
.
template
getProp
<
force
>(
p
)[
2
]
+=
f
.
get
(
2
);
++
Np
;
}
++
it2
;
}
}
/* ### WIKI 14 ###
*
* The function to calculate the total energy. It require the same parameter as calculate forces
*
*/
double
calc_energy
(
vector_dist
<
3
,
double
,
aggregate
<
double
[
3
],
double
[
3
]
>
>
&
vd
,
CellList
<
3
,
double
,
FAST
,
shift
<
3
,
double
>
>
&
NN
,
double
L
)
{
// ### WIKI 15 ###
//
// Reset the counter for the energy and
// update the cell list from the actual particle configuration
//
//
double
E
=
0.0
;
vd
.
updateCellList
(
NN
);
// ### WIKI 16 ###
//
// Calculate the forces
//
auto
it2
=
vd
.
getDomainIterator
();
while
(
it2
.
isNext
())
{
auto
p
=
it2
.
get
();
Point
<
3
,
double
>
xp
=
vd
.
getPos
<
0
>
(
p
);
vd
.
template
getProp
<
force
>(
p
)[
0
]
=
0.0
;
vd
.
template
getProp
<
force
>(
p
)[
1
]
=
0.0
;
vd
.
template
getProp
<
force
>(
p
)[
2
]
=
0.0
;
// For each neighborhood particle
auto
Np
=
NN
.
template
getNNIterator
<
NO_CHECK
>(
NN
.
getCell
(
vd
.
getPos
<
0
>
(
p
)));
while
(
Np
.
isNext
())
{
// Neighborhood particle q
auto
q
=
Np
.
get
();
if
(
q
==
p
.
getKey
())
{
++
Np
;
continue
;};
Point
<
3
,
double
>
xq
=
vd
.
getPos
<
0
>
(
q
);
Point
<
3
,
double
>
r
=
xp
-
xq
;
// take the normalized direction
float
rn
=
r
.
norm
();
r
/=
rn
;
rn
*=
L
;
// potential energy (using pow is slower)
E
+=
4.0
*
(
1.0
/
(
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
*
rn
)
+
1.0
/
(
rn
*
rn
*
rn
*
rn
*
rn
*
rn
)
);
// Kinetic energy
E
+=
sqrt
(
vd
.
template
getProp
<
force
>(
p
)[
0
]
*
vd
.
template
getProp
<
force
>(
p
)[
0
]
+
vd
.
template
getProp
<
force
>(
p
)[
1
]
*
vd
.
template
getProp
<
force
>(
p
)[
1
]
+
vd
.
template
getProp
<
force
>(
p
)[
2
]
*
vd
.
template
getProp
<
force
>(
p
)[
2
]);
++
Np
;
}
++
it2
;
}