container.hh

Go to the documentation of this file.

// Copyright Epic Games, Inc. All Rights Reserved.
// Modified version of Voro++'s source file
 
// Voro++, a 3D cell-based Voronoi library
//
// Author   : Chris H. Rycroft (LBL / UC Berkeley)
// Email    : chr@alum.mit.edu
// Date     : August 30th 2011
 
#ifndef VOROPP_CONTAINER_HH
#define VOROPP_CONTAINER_HH
 
#include <cstdio>
#include <vector>
 
#include "config.hh"
#include "common.hh"
#include "v_base.hh"
#include "cell.hh"
#include "c_loops.hh"
#include "v_compute.hh"
#include "rad_option.hh"
 
namespace voro {
 
class wall {
    public:
        virtual ~wall() {}
        virtual bool point_inside(double x,double y,double z) const = 0;
        virtual bool cut_cell(voronoicell &c,double x,double y,double z) const = 0;
        virtual bool cut_cell(voronoicell_neighbor &c,double x,double y,double z) const = 0;
};
 
class wall_list {
    public:
        wall **walls;
        wall **wep;
        wall_list();
        wall_list(const wall_list& Other) = delete;
        wall_list& operator=(const wall_list& Other) = delete;
        wall_list(wall_list&& Other) noexcept : walls(Other.walls), wep(Other.wep), wel(Other.wel), current_wall_size(Other.current_wall_size)
        {
            Other.walls = nullptr;
        }
        wall_list& operator=(wall_list&& Other)
        {
            if (this != &Other)
            {
                walls = Other.walls;
                wep = Other.wep;
                wel = Other.wel;
                current_wall_size = Other.current_wall_size;
                Other.walls = nullptr;
            }
            return *this;
        }
        ~wall_list();
        inline void add_wall(wall *w) {
            if(wep==wel) increase_wall_memory();
            *(wep++)=w;
        }
        inline void add_wall(wall &w) {add_wall(&w);}
        void add_wall(wall_list &wl);
        inline bool point_inside_walls(double x,double y,double z) const {
            for(wall **wp=walls;wp<wep;wp++) if(!((*wp)->point_inside(x,y,z))) return false;
            return true;
        }
        template<class c_class>
        bool apply_walls(c_class &c,double x,double y,double z) const {
            for(wall **wp=walls;wp<wep;wp++) if(!((*wp)->cut_cell(c,x,y,z))) return false;
            return true;
        }
        void deallocate();
    protected:
        void increase_wall_memory();
        wall **wel;
        int current_wall_size;
};
 
class container_base : public voro_base, public wall_list {
    public:
        const double ax;
        const double bx;
        const double ay;
        const double by;
        const double az;
        const double bz;
        const bool xperiodic;
        const bool yperiodic;
        const bool zperiodic;
        int **id;
        double **p;
        int *co;
        int *mem;
        const int ps;
        container_base(double ax_,double bx_,double ay_,double by_,double az_,double bz_,
                int nx_,int ny_,int nz_,bool xperiodic_,bool yperiodic_,bool zperiodic_,
                int init_mem,int ps_);
        container_base(container_base& Other) = delete;
        container_base& operator=(const container_base& Other) = delete;
        container_base(container_base&& Other) noexcept : voro_base(MoveTemp(Other)), wall_list(MoveTemp(Other)),
            ax(Other.ax), bx(Other.bx), ay(Other.ay), by(Other.by), az(Other.az), bz(Other.bz),
            xperiodic(Other.xperiodic), yperiodic(Other.yperiodic), zperiodic(Other.zperiodic),
            id(Other.id), p(Other.p), co(Other.co), mem(Other.mem), ps(Other.ps)
        {
            Other.id = nullptr;
            Other.p = nullptr;
            Other.co = nullptr;
            Other.mem = nullptr;
        }
        // Note: Move assignment not implemented due to const members on this class (and it has not been needed yet)
        container_base& operator=(container_base&& Other) = delete;
 
        ~container_base();
        bool point_inside(double x,double y,double z) const;
        void region_count() const;
        template<class v_cell>
        inline bool initialize_voronoicell(v_cell &c,int ijk,int q,int ci,int cj,int ck,
                int &i,int &j,int &k,double &x,double &y,double &z,int &disp) const {
            double x1,x2,y1,y2,z1,z2,*pp=p[ijk]+ps*q;
            x=*(pp++);y=*(pp++);z=*pp;
            if(xperiodic) {x1=-(x2=0.5*(bx-ax));i=nx;} else {x1=ax-x;x2=bx-x;i=ci;}
            if(yperiodic) {y1=-(y2=0.5*(by-ay));j=ny;} else {y1=ay-y;y2=by-y;j=cj;}
            if(zperiodic) {z1=-(z2=0.5*(bz-az));k=nz;} else {z1=az-z;z2=bz-z;k=ck;}
            c.init(x1,x2,y1,y2,z1,z2);
            if(!apply_walls(c,x,y,z)) return false;
            disp=ijk-i-nx*(j+ny*k);
            return true;
        }
        inline void initialize_search(int ci,int cj,int ck,int ijk,int &i,int &j,int &k,int &disp) const {
            i=xperiodic?nx:ci;
            j=yperiodic?ny:cj;
            k=zperiodic?nz:ck;
            disp=ijk-i-nx*(j+ny*k);
        }
        inline void frac_pos(double x,double y,double z,double ci,double cj,double ck,
                double &fx,double &fy,double &fz) const {
            fx=x-ax-boxx*ci;
            fy=y-ay-boxy*cj;
            fz=z-az-boxz*ck;
        }
        inline int region_index(int ci,int cj,int ck,int ei,int ej,int ek,double &qx,double &qy,double &qz,int &disp) const {
            if(xperiodic) {if(ci+ei<nx) {ei+=nx;qx=-(bx-ax);} else if(ci+ei>=(nx<<1)) {ei-=nx;qx=bx-ax;} else qx=0;}
            if(yperiodic) {if(cj+ej<ny) {ej+=ny;qy=-(by-ay);} else if(cj+ej>=(ny<<1)) {ej-=ny;qy=by-ay;} else qy=0;}
            if(zperiodic) {if(ck+ek<nz) {ek+=nz;qz=-(bz-az);} else if(ck+ek>=(nz<<1)) {ek-=nz;qz=bz-az;} else qz=0;}
            return disp+ei+nx*(ej+ny*ek);
        }
        inline int total_particles() const {
            int tp=*co;
            for(int *cop=co+1;cop<co+nxyz;cop++) tp+=*cop;
            return tp;
        }
    protected:
        void add_particle_memory(int i);
        bool put_locate_block(int &ijk,double &x,double &y,double &z);
        inline bool put_remap(int &ijk,double &x,double &y,double &z) const;
        inline bool remap(int &ai,int &aj,int &ak,int &ci,int &cj,int &ck,double &x,double &y,double &z,int &ijk) const;
};
 
class container : public container_base, public radius_mono {
    public:
        container(double ax_,double bx_,double ay_,double by_,double az_,double bz_,
                int nx_,int ny_,int nz_,bool xperiodic_,bool yperiodic_,bool zperiodic_,int init_mem);
        voro_compute<container> make_compute() const;
 
        void clear();
        void put(int n,double x,double y,double z);
        void put(particle_order &vo,int n,double x,double y,double z);
        void compute_all_cells(voro_compute<container>& vc);
        double sum_cell_volumes(voro_compute<container>& vc);
        bool find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid,voro_compute<container>& vc) const;
        
        template<class v_cell,class c_loop>
        inline bool compute_cell(v_cell &c,c_loop &vl,voro_compute<container>& vc) const {
            return vc.compute_cell(c,vl.ijk,vl.q,vl.i,vl.j,vl.k);
        }
        template<class v_cell>
        inline bool compute_cell(v_cell &c,int ijk,int q,voro_compute<container>& vc) const {
            int k=ijk/nxy,ijkt=ijk-nxy*k,j=ijkt/nx,i=ijkt-j*nx;
            return vc.compute_cell(c,ijk,q,i,j,k);
        }
        template<class v_cell>
        inline bool compute_ghost_cell(v_cell &c,double x,double y,double z,voro_compute<container>& vc) {
            int ijk;
            if(put_locate_block(ijk,x,y,z)) {
                double *pp=p[ijk]+3*co[ijk]++;
                *(pp++)=x;*(pp++)=y;*pp=z;
                bool q=compute_cell(c,ijk,co[ijk]-1,vc);
                co[ijk]--;
                return q;
            }
            return false;
        }
    private:
        friend class voro_compute<container>;
};
 
class container_poly : public container_base, public radius_poly {
    public:
        container_poly(double ax_,double bx_,double ay_,double by_,double az_,double bz_,
                int nx_,int ny_,int nz_,bool xperiodic_,bool yperiodic_,bool zperiodic_,int init_mem);
        voro_compute<container_poly> make_compute() const;
 
        void clear();
        void put(int n,double x,double y,double z,double r);
        void put(particle_order &vo,int n,double x,double y,double z,double r);
        void compute_all_cells(voro_compute<container>& vc);
        double sum_cell_volumes(voro_compute<container>& vc);
        template<class v_cell,class c_loop>
        inline bool compute_cell(v_cell &c,c_loop &vl,voro_compute<container>& vc) const {
            return vc.compute_cell(c,vl.ijk,vl.q,vl.i,vl.j,vl.k);
        }
        template<class v_cell>
        inline bool compute_cell(v_cell &c,int ijk,int q,voro_compute<container>& vc) const {
            int k=ijk/nxy,ijkt=ijk-nxy*k,j=ijkt/nx,i=ijkt-j*nx;
            return vc.compute_cell(c,ijk,q,i,j,k);
        }
        template<class v_cell>
        inline bool compute_ghost_cell(v_cell &c,double x,double y,double z,double r,voro_compute<container>& vc) {
            int ijk;
            if(put_locate_block(ijk,x,y,z)) {
                double *pp=p[ijk]+4*co[ijk]++,tm=max_radius;
                *(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
                if(r>max_radius) max_radius=r;
                bool q=compute_cell(c,ijk,co[ijk]-1);
                co[ijk]--;max_radius=tm;
                return q;
            }
            return false;
        }
        bool find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid,voro_compute<container>& vc) const;
    private:
        friend class voro_compute<container_poly>;
};
 
}
 
#endif