percolation/diamond.cpp at 0.2 · cmdawson/percolation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include "diamond.h"
#include "point3d.h"
#include <vector>
#include <map>


#include <iostream>

using namespace std;
using namespace boost;

Diamond::Diamond(int n) : Lattice(n)
{
    // cubic unit cell of the diamond 'lattice'
    std:: vector<point3d> ucell{
	point3d(0,0,0), point3d(0,2,2), point3d(2,0,2),
	point3d(2,2,0), point3d(3,3,3), point3d(3,1,1),
	point3d(1,3,1), point3d(1,1,3), point3d(4,0,0),
	point3d(0,4,0), point3d(0,0,4), point3d(4,4,0),
	point3d(0,4,4), point3d(4,0,4), point3d(4,4,4),
	point3d(4,2,2), point3d(2,4,2), point3d(2,2,4)
    };
    // valid edge vectors
    std::vector<point3d> evec{
	point3d(1,1,1), point3d(1,1,-1), point3d(1,-1,1), point3d(-1,1,1)
    };

    map<point3d,Vertex> vertex_map;
    int index = 0;

    // add vertices by translating unit cell
    for (int i=0;i<n;i++)
	for (int j=0;j<n;j++)
	    for (int k=0;k<n;k++)
	    {
		point3d tr(4*i,4*j,4*k);
		for (const auto& p : ucell)
		{
		    point3d newp = p+tr;
		    if (vertex_map.find(newp) == vertex_map.end())
		    {
		    // Jan 22: New addition, defining which vertices of the diamond
		    // are NEAR, FAR and NEITHER for isCrossable()
		    vertex_map[newp] = point3d::add_as_vertex(newp, index++, adjacency_list);

		    if(n==1)
		    {
			if(p.x==0){
			    //cout<<n<<" "<<p.x<<" "<<VertexT::NEAR<<endl;
			    adjacency_list[vertex_map[newp]].traverse= VertexT::NEAR;}
			else if(p.x==4){
			    adjacency_list[vertex_map[newp]].traverse= VertexT::FAR;}
			else {adjacency_list[vertex_map[newp]].traverse= VertexT::NEITHER;}
		    }
		else if (n>1)
		{
			if(i==0 && p.x==0){
				//cout<<n<<" "<<p.x<<" "<<VertexT::NEAR<<endl;
						adjacency_list[vertex_map[newp]].traverse= VertexT::NEAR;}
						else if(i==(n-1) && p.x==4)
							adjacency_list[vertex_map[newp]].traverse=VertexT::FAR;
						else adjacency_list[vertex_map[newp]].traverse=VertexT::NEITHER;
}
				}
			}
	    }

    // add edges by searching from each vertex
    for (const auto& p : vertex_map)
    {
	Vertex vp = vertex_map[p.first];
	for (const auto& e : evec)
	{
	    point3d nbr = p.first + e;
	    const auto& vn = vertex_map.find(nbr);
	    if (vn != vertex_map.end())
			add_edge(vp, vn->second, adjacency_list);
	}
    }
}


Diamond::~Diamond() {}

void Diamond::percolate(const double& p, const double& q)
{
    // TODO
}