Line data Source code
1 : #ifndef IC_VORONOI_H_
2 : #define IC_VORONOI_H_
3 :
4 : #include "Set/Set.H"
5 : #include "IC/IC.H"
6 : #include "IO/ParmParse.H"
7 :
8 : namespace IC
9 : {
10 : class Voronoi : public IC<Set::Scalar>
11 : {
12 : public:
13 : static constexpr const char* name = "voronoi";
14 :
15 : enum Type {Partition, Values};
16 :
17 : Voronoi (amrex::Vector<amrex::Geometry> &_geom) : IC(_geom) {}
18 5 : Voronoi (amrex::Vector<amrex::Geometry> &_geom,IO::ParmParse &pp, std::string name) : IC(_geom)
19 20 : {pp_queryclass(name,*this);}
20 : Voronoi (amrex::Vector<amrex::Geometry> &_geom, int _number_of_grains, Set::Scalar a_alpha) : IC(_geom)
21 : {
22 : Define(_number_of_grains,a_alpha);
23 : }
24 : Voronoi (amrex::Vector<amrex::Geometry> &a_geom, int a_number_of_grains) : IC(a_geom)
25 : {
26 : Define(a_number_of_grains,1.0);
27 : }
28 :
29 5 : void Define (int a_number_of_grains, Set::Scalar a_alpha)
30 : {
31 5 : Define(a_number_of_grains, std::vector<Set::Scalar>(a_number_of_grains,a_alpha),Type::Partition);
32 5 : }
33 :
34 5 : void Define (int a_number_of_grains,
35 : std::vector<Set::Scalar> a_alpha,
36 : Type a_type = Type::Values)
37 : {
38 5 : number_of_grains = a_number_of_grains;
39 5 : alpha = a_alpha;
40 5 : type = a_type;
41 :
42 5 : voronoi.resize(number_of_grains);
43 5 : srand(seed);
44 245 : for (int n = 0; n<number_of_grains; n++)
45 : {
46 240 : AMREX_D_TERM(voronoi[n](0) = geom[0].ProbLo(0) + (geom[0].ProbHi(0)-geom[0].ProbLo(0))*Util::Random();,
47 : voronoi[n](1) = geom[0].ProbLo(1) + (geom[0].ProbHi(1)-geom[0].ProbLo(1))*Util::Random();,
48 : voronoi[n](2) = geom[0].ProbLo(2) + (geom[0].ProbHi(2)-geom[0].ProbLo(2))*Util::Random(););
49 : }
50 5 : };
51 :
52 56 : void Add(const int &lev, Set::Field<Set::Scalar> &a_field, Set::Scalar)
53 : {
54 56 : Set::Vector size;
55 56 : AMREX_D_TERM(size(0) = geom[0].ProbHi()[0] - geom[0].ProbLo()[0];,
56 : size(1) = geom[0].ProbHi()[1] - geom[0].ProbLo()[1];,
57 : size(2) = geom[0].ProbHi()[2] - geom[0].ProbLo()[2];)
58 :
59 10485 : for (amrex::MFIter mfi(*a_field[lev],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
60 : {
61 10430 : amrex::Box bx = mfi.tilebox();
62 10430 : bx.grow(a_field[lev]->nGrow());
63 10430 : int ncomp = a_field[lev]->nComp();
64 10430 : amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
65 10430 : amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
66 :
67 2068992 : Set::Vector x;
68 2068992 : AMREX_D_TERM(x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];,
69 : x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];,
70 : x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom[lev].CellSize()[2];);
71 :
72 2068992 : amrex::Real min_distance = std::numeric_limits<amrex::Real>::infinity();
73 2068992 : int min_grain_id = -1;
74 :
75 89290935 : for (int n = 0; n<number_of_grains; n++)
76 : {
77 87221944 : Set::Scalar d = (x - voronoi[n]).lpNorm<2>();
78 :
79 87221944 : if (geom[0].isPeriodic(0))
80 : {
81 18918344 : d = std::min(d,
82 18918344 : std::min( (x-voronoi[n] + size(0)*Set::Vector::Unit(0)).lpNorm<2>(),
83 37836688 : (x-voronoi[n] - size(0)*Set::Vector::Unit(0)).lpNorm<2>()));
84 : }
85 : #if AMREX_SPACEDIM>1
86 87221944 : if (geom[0].isPeriodic(1))
87 : {
88 18918343 : d = std::min(d,
89 18918344 : std::min( (x-voronoi[n] + size(0)*Set::Vector::Unit(1)).lpNorm<2>(),
90 37836687 : (x-voronoi[n] - size(0)*Set::Vector::Unit(1)).lpNorm<2>()));
91 : }
92 : #endif
93 : #if AMREX_SPACEDIM>2
94 0 : if (geom[0].isPeriodic(2))
95 : {
96 0 : d = std::min(d,
97 0 : std::min( (x-voronoi[n] + size(0)*Set::Vector::Unit(2)).lpNorm<2>(),
98 0 : (x-voronoi[n] - size(0)*Set::Vector::Unit(2)).lpNorm<2>()));
99 : }
100 : #endif
101 87221943 : if (d<min_distance)
102 : {
103 8874692 : min_distance = d;
104 8874692 : min_grain_id = n;
105 : }
106 : }
107 :
108 2068991 : if (type == Type::Values) field(i,j,k) = alpha[min_grain_id];
109 4137982 : else if (type == Type::Partition) field(i,j,k,min_grain_id % ncomp) = alpha[min_grain_id];
110 2068991 : });
111 55 : }
112 55 : }
113 :
114 5 : static void Parse(Voronoi &value, IO::ParmParse &pp)
115 : {
116 5 : pp_query("number_of_grains",value.number_of_grains); // Number of grains
117 10 : if (pp.contains("alpha"))
118 : {
119 0 : pp_queryarr("alpha",value.alpha); // Value to take in the region [1.0]
120 0 : value.Define(value.number_of_grains,value.alpha);
121 : }
122 : else
123 : {
124 5 : value.Define(value.number_of_grains,1.0);
125 : }
126 5 : pp_query("seed",value.seed); // Random seed to use
127 5 : }
128 :
129 : private:
130 : int number_of_grains;
131 : int seed = 1;
132 : std::vector<Set::Scalar> alpha;
133 : std::vector<Set::Vector> voronoi;
134 : Type type;
135 : };
136 : }
137 : #endif
|
