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