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
      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           1 :     Voronoi (amrex::Vector<amrex::Geometry> &_geom,IO::ParmParse &pp, std::string name) : IC(_geom) 
      19           1 :     {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           1 :     void Define (int a_number_of_grains, Set::Scalar a_alpha)
      30             :     {
      31           1 :         Define(a_number_of_grains, std::vector<Set::Scalar>(a_number_of_grains,a_alpha),Type::Partition);
      32           1 :     }
      33             : 
      34           1 :     void Define (int a_number_of_grains,
      35             :                 std::vector<Set::Scalar> a_alpha,
      36             :                 Type a_type = Type::Values)
      37             :     {
      38           1 :         number_of_grains = a_number_of_grains;
      39           1 :         alpha = a_alpha;
      40           1 :         type = a_type;
      41             : 
      42           1 :         voronoi.resize(number_of_grains);
      43           1 :         srand(seed);
      44          21 :         for (int n = 0; n<number_of_grains; n++)
      45             :         {
      46          20 :             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           1 :     };
      51             :     
      52           7 :     void Add(const int &lev, Set::Field<Set::Scalar> &a_field, Set::Scalar)
      53             :     {
      54           7 :         Set::Vector size;
      55           7 :         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         982 :         for (amrex::MFIter mfi(*a_field[lev],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
      60             :         {
      61         975 :             amrex::Box bx = mfi.tilebox();
      62         975 :             bx.grow(a_field[lev]->nGrow());
      63         975 :             int ncomp = a_field[lev]->nComp();
      64         975 :             amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
      65         975 :             amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
      66             : 
      67      107636 :                 Set::Vector x;
      68      107636 :                 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      107636 :                 amrex::Real min_distance = std::numeric_limits<amrex::Real>::infinity();
      73      107636 :                 int min_grain_id = -1;
      74             : 
      75     2260360 :                 for (int n = 0; n<number_of_grains; n++)
      76             :                 {
      77     2152720 :                     Set::Scalar d = (x - voronoi[n]).lpNorm<2>();
      78             : 
      79     2152720 :                     if (geom[0].isPeriodic(0))
      80             :                         {
      81     2152720 :                             d = std::min(d,
      82     2152720 :                                         std::min( (x-voronoi[n] + size(0)*Set::Vector::Unit(0)).lpNorm<2>(),
      83     4305440 :                                                     (x-voronoi[n] - size(0)*Set::Vector::Unit(0)).lpNorm<2>()));
      84             :                         }
      85             : #if AMREX_SPACEDIM>1
      86     2152720 :                     if (geom[0].isPeriodic(1))
      87             :                         {
      88     2152720 :                             d = std::min(d,
      89     2152720 :                                         std::min( (x-voronoi[n] + size(0)*Set::Vector::Unit(1)).lpNorm<2>(),
      90     4305440 :                                                     (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     2152720 :                     if (d<min_distance)
     102             :                         {
     103      390192 :                             min_distance = d;
     104      390192 :                             min_grain_id = n;
     105             :                         }
     106             :                 }
     107             : 
     108      107636 :                 if (type == Type::Values) field(i,j,k) = alpha[min_grain_id];
     109      215272 :                 else if (type == Type::Partition) field(i,j,k,min_grain_id % ncomp) = alpha[min_grain_id];
     110      107636 :             });
     111             :         }
     112           7 :     }
     113             : 
     114           1 :     static void Parse(Voronoi &value, IO::ParmParse &pp)
     115             :     {
     116           1 :         pp_query("number_of_grains",value.number_of_grains); // Number of grains
     117           1 :         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           1 :             value.Define(value.number_of_grains,1.0);
     125             :         }
     126           1 :         pp_query("seed",value.seed); // Random seed to use
     127           1 :     }
     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