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 1 : {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;
65 1950 : amrex::IndexType ixtype = a_field[lev]->ixType();
66 3900 : if (ixtype == amrex::IndexType::TheCellType()) bx = mfi.growntilebox();
67 0 : else if (ixtype == amrex::IndexType::TheNodeType()) bx = mfi.grownnodaltilebox();
68 :
69 1950 : int ncomp = a_field[lev]->nComp();
70 1950 : amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
71 1950 : amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
72 :
73 215272 : Set::Vector x;
74 215272 : AMREX_D_TERM(x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];,
75 : x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];,
76 : x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom[lev].CellSize()[2];);
77 :
78 215272 : amrex::Real min_distance = std::numeric_limits<amrex::Real>::infinity();
79 215272 : int min_grain_id = -1;
80 :
81 4520712 : for (int n = 0; n<number_of_grains; n++)
82 : {
83 4305440 : Set::Scalar d = (x - voronoi[n]).lpNorm<2>();
84 :
85 4305440 : if (geom[0].isPeriodic(0))
86 : {
87 4305440 : d = std::min(d,
88 4305440 : std::min( (x-voronoi[n] + size(0)*Set::Vector::Unit(0)).lpNorm<2>(),
89 8610880 : (x-voronoi[n] - size(0)*Set::Vector::Unit(0)).lpNorm<2>()));
90 : }
91 : #if AMREX_SPACEDIM>1
92 4305440 : if (geom[0].isPeriodic(1))
93 : {
94 4305440 : d = std::min(d,
95 4305440 : std::min( (x-voronoi[n] + size(0)*Set::Vector::Unit(1)).lpNorm<2>(),
96 8610880 : (x-voronoi[n] - size(0)*Set::Vector::Unit(1)).lpNorm<2>()));
97 : }
98 : #endif
99 : #if AMREX_SPACEDIM>2
100 0 : if (geom[0].isPeriodic(2))
101 : {
102 0 : d = std::min(d,
103 0 : std::min( (x-voronoi[n] + size(0)*Set::Vector::Unit(2)).lpNorm<2>(),
104 0 : (x-voronoi[n] - size(0)*Set::Vector::Unit(2)).lpNorm<2>()));
105 : }
106 : #endif
107 4305440 : if (d<min_distance)
108 : {
109 780384 : min_distance = d;
110 780384 : min_grain_id = n;
111 : }
112 : }
113 :
114 215272 : if (type == Type::Values) field(i,j,k) = alpha[min_grain_id];
115 430544 : else if (type == Type::Partition) field(i,j,k,min_grain_id % ncomp) = alpha[min_grain_id];
116 215272 : });
117 14 : }
118 14 : }
119 :
120 1 : static void Parse(Voronoi &value, IO::ParmParse &pp)
121 : {
122 1 : pp_query("number_of_grains",value.number_of_grains); // Number of grains
123 2 : if (pp.contains("alpha"))
124 : {
125 0 : pp_queryarr("alpha",value.alpha); // Value to take in the region [1.0]
126 0 : value.Define(value.number_of_grains,value.alpha);
127 : }
128 : else
129 : {
130 1 : value.Define(value.number_of_grains,1.0);
131 : }
132 1 : pp_query("seed",value.seed); // Random seed to use
133 1 : }
134 :
135 : private:
136 : int number_of_grains;
137 : int seed = 1;
138 : std::vector<Set::Scalar> alpha;
139 : std::vector<Set::Vector> voronoi;
140 : Type type;
141 : };
142 : }
143 : #endif
|
