diff --git a/CHANGELOG.md b/CHANGELOG.md index 2d2795d5911b05abb948be856d683ac3eba49826..df303bd7b4a003b75227da6c1a1beacba0f24dbb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,8 +20,10 @@ All notable changes to this project will be documented in this file. - In case of miss compilation ignore system wide installation - Bug in VTK writer binary in case of vectors - Bug in VTK writer binary: long int are not supported removing output +- Bug in the constructor with stencil bigger than one - +### Changed +- CellList types has changed ## [0.8.0] 28 February 2016 diff --git a/example/Grid/3_gray_scott/main.cpp b/example/Grid/3_gray_scott/main.cpp index 932a15a344d334fac6856c49bcfad017d2caf543..3bb0eff8397731c57fc222a968199a1fa1faa3c8 100644 --- a/example/Grid/3_gray_scott/main.cpp +++ b/example/Grid/3_gray_scott/main.cpp @@ -10,10 +10,10 @@ * This example show the usage of periodic grid with ghost part given in grid units to solve * the following system of equations * - * \f$\frac{\partial u}{\partial t} = D_u \nabla u -uv^2 + F(1-u)\f$ + * \f$\frac{\partial u}{\partial t} = D_u \nabla u - uv^2 + F(1-u)\f$ * * - * \f$\frac{\partial v}{\partial t} = D_v \nabla v -uv^2 - (F + k)v\f$ + * \f$\frac{\partial v}{\partial t} = D_v \nabla v + uv^2 - (F + k)v\f$ * * ## Constants and functions ## * @@ -201,7 +201,7 @@ int main(int argc, char* argv[]) grid_dist_id<2, double, aggregate> Old(sz,domain,g,bc); // New grid with the decomposition of the old grid - grid_dist_id<2, double, aggregate> New(Old.getDecomposition(),sz,domain,g); + grid_dist_id<2, double, aggregate> New(Old.getDecomposition(),sz,g); // spacing of the grid on x and y diff --git a/example/Grid/3_gray_scott_3d/main.cpp b/example/Grid/3_gray_scott_3d/main.cpp index 73ea16464f9bc85cde87a13baddd44b21bdaaff7..91ec575f70c4bba39f08c08f14cf92a8e416d1be 100644 --- a/example/Grid/3_gray_scott_3d/main.cpp +++ b/example/Grid/3_gray_scott_3d/main.cpp @@ -46,8 +46,16 @@ void init(grid_dist_id<3,double,aggregate > & Old, grid_dist_id<3 ++it; } - grid_key_dx<3> start({(long int)std::floor(Old.size(0)*1.55f/domain.getHigh(0)),(long int)std::floor(Old.size(1)*1.55f/domain.getHigh(1)),(long int)std::floor(Old.size(1)*1.55f/domain.getHigh(2))}); - grid_key_dx<3> stop ({(long int)std::ceil (Old.size(0)*1.85f/domain.getHigh(0)),(long int)std::ceil (Old.size(1)*1.85f/domain.getHigh(1)),(long int)std::floor(Old.size(1)*1.85f/domain.getHigh(1))}); + long int x_start = Old.size(0)*1.55f/domain.getHigh(0); + long int y_start = Old.size(1)*1.55f/domain.getHigh(1); + long int z_start = Old.size(1)*1.55f/domain.getHigh(2); + + long int x_stop = Old.size(0)*1.85f/domain.getHigh(0); + long int y_stop = Old.size(1)*1.85f/domain.getHigh(1); + long int z_stop = Old.size(1)*1.85f/domain.getHigh(2); + + grid_key_dx<3> start({x_start,y_start,z_start}); + grid_key_dx<3> stop ({x_stop,y_stop,z_stop}); auto it_init = Old.getSubDomainIterator(start,stop); while (it_init.isNext()) @@ -72,7 +80,7 @@ int main(int argc, char* argv[]) Box<3,double> domain({0.0,0.0},{2.5,2.5,2.5}); // grid size - size_t sz[3] = {128,128,128}; + size_t sz[3] = {128,128,128}; // Define periodicity of the grid periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC}; @@ -90,11 +98,11 @@ int main(int argc, char* argv[]) double dv = 1*1e-5; // Number of timesteps - size_t timeSteps = 17000; + size_t timeSteps = 5000; // K and F (Physical constant in the equation) - double K = 0.065; - double F = 0.034; + double K = 0.014; + double F = 0.053; //! \cond [init lib] \endcond @@ -186,7 +194,7 @@ int main(int argc, char* argv[]) // visualization if (i % 60 == 0) { - Old.write("output",count,VTK_WRITER | FORMAT_BINARY); + Old.write_frame("output",count,VTK_WRITER | FORMAT_BINARY); count++; } } diff --git a/example/Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp b/example/Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp index 27b9cbf7d9bd7ce849dc54d80652234d84cfe2c3..6f5ae88b62bfcad1e73a46f07f028dae422b56e8 100644 --- a/example/Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp +++ b/example/Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp @@ -404,7 +404,7 @@ int main(int argc, char* argv[]) * * Once we have the solution we copy it on the grid * - * \snippet Numerics/Stoke_flow/0_2D_incompressible/main_eigen.cpp copy write + * \snippet Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp copy write * */ @@ -424,7 +424,7 @@ int main(int argc, char* argv[]) * * At the very end of the program we have always to de-initialize the library * - * \snippet Numerics/Stoke_flow/0_2D_incompressible/main_eigen.cpp fin lib + * \snippet Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp fin lib * */ @@ -440,7 +440,7 @@ int main(int argc, char* argv[]) * * # Full code # {#num_sk_inc_2D_ps_code} * - * \include Numerics/Stoke_flow/0_2D_incompressible/main_eigen.cpp + * \include Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp * */ } diff --git a/example/Vector/1_celllist/main.cpp b/example/Vector/1_celllist/main.cpp index 28c4a6fe1cbc4912dad3a63aee98f31030f0df6a..c244c46899c53a18f29c83fe0c6f4a80578c0196 100644 --- a/example/Vector/1_celllist/main.cpp +++ b/example/Vector/1_celllist/main.cpp @@ -247,7 +247,7 @@ int main(int argc, char* argv[]) Point<3,float> xp = vd.getPos(p); // Get an iterator of all the particles neighborhood of p - auto Np = NN.getIterator(NN.getCell(vd.getPos(p))); + auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p))); // For each particle near p while (Np.isNext()) diff --git a/example/Vector/3_molecular_dynamic/main.cpp b/example/Vector/3_molecular_dynamic/main.cpp index 8d70b8c9a1dca9d8ed1f7ca0efe7840893041d20..7a2e86db5b8ecd9fd73e16410abae4d768782f6d 100644 --- a/example/Vector/3_molecular_dynamic/main.cpp +++ b/example/Vector/3_molecular_dynamic/main.cpp @@ -6,8 +6,6 @@ */ #include "Vector/vector_dist.hpp" -#include "Decomposition/CartDecomposition.hpp" -#include "data_type/aggregate.hpp" #include "Plot/GoogleChart.hpp" #include "Plot/util.hpp" #include "timer.hpp" @@ -19,7 +17,10 @@ * * # Molecular Dynamic with Lennard-Jones potential # {#e3_md} * - * This example show a simple Lennard-Jones molecular dynamic simulation in a stable regime + * This example show a simple Lennard-Jones molecular dynamic simulation in a stable regime. + * Particle feel each other by the potential. + * + * \f$V(x_p,x_q) = 4( (\frac{\sigma}{r})^{12} - (\frac{\sigma}{r})^6 ) \f$ * * ## Constants ## * @@ -55,7 +56,7 @@ constexpr int force = 1; //! \cond [calc forces] \endcond -void calc_forces(vector_dist<3,double, aggregate > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut2) +template void calc_forces(vector_dist<3,double, aggregate > & vd, CellList & NN, double sigma12, double sigma6, double r_cut2) { //! \cond [calc forces] \endcond @@ -81,12 +82,12 @@ void calc_forces(vector_dist<3,double, aggregate > & vd, Ce /*! * - * \page Vector_3_md Vector 3 molecular dynamic with cell-list + * \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list * * Get an iterator over the particles and get its position. For each particle p iterate in its neighborhood q * and calculate the force based on the Lennard-Jhones potential given by * - * \f$F(x_p,x_q) = 24(\frac{2}{r^{13}} - \frac{1}{r^{7}}) r \f$ + * \f$F(x_p,x_q) = 24( \frac{2 \sigma^{12}}{r^{13}} - \frac{\sigma^6}{r^{7}}) \hat{r} \f$ * * \see \ref e0_s_assign_pos * @@ -108,7 +109,7 @@ void calc_forces(vector_dist<3,double, aggregate > & vd, Ce // Get the position xp of the particle Point<3,double> xp = vd.getPos(p); - // Reset the forice counter + // Reset the force counter vd.template getProp(p)[0] = 0.0; vd.template getProp(p)[1] = 0.0; vd.template getProp(p)[2] = 0.0; @@ -175,7 +176,7 @@ void calc_forces(vector_dist<3,double, aggregate > & vd, Ce //! \cond [calc energy] \endcond -double calc_energy(vector_dist<3,double, aggregate > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut2) +template double calc_energy(vector_dist<3,double, aggregate > & vd, CellList & NN, double sigma12, double sigma6, double r_cut2) { double rc = r_cut2; double shift = 2.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) ); @@ -208,7 +209,7 @@ double calc_energy(vector_dist<3,double, aggregate > & vd, * First we get an iterator over the particles and get its position. For each particle p iterate in its neighborhood q * and calculate the energy based on the Lennard-Jhones potential given by * - * \f$V(x_p,x_q) = 4(\frac{1}{r^{12}} - \frac{1}{r^{6}}) r \f$ + * \f$V(x_p,x_q) = 4(\frac{1}{r^{12}} - \frac{1}{r^{6}}) \f$ * * \see \ref e0_s_assign_pos * @@ -284,43 +285,24 @@ int main(int argc, char* argv[]) * * ## Initialization ## {#e3_md_init} * - * After we defined the two main function calc forces and calc energy, we define - * important parameters of the simulation, time step integration, + * After we defined the two main function calc forces and calc energy, we Initialize + * the library, we create a Box that define our domain, boundary conditions and ghost. + * Than we define important parameters of the simulation, time step integration, * size of the box, and cut-off radius of the interaction. We also define 2 vectors * x and y (they are like std::vector) used for statistic * + * \see \ref e0_s_init + * * \snippet Vector/3_molecular_dynamic/main.cpp constants run * */ //! \cond [constants run] \endcond - double dt = 0.0005; + openfpm_init(&argc,&argv); + double sigma = 0.1; double r_cut = 3.0*sigma; - double sigma12 = pow(sigma,12); - double sigma6 = pow(sigma,6); - - openfpm::vector x; - openfpm::vector> y; - - //! \cond [constants run] \endcond - - /*! - * \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list - * - * Here we Initialize the library, we create a Box that define our domain, boundary conditions and ghost - * - * \see \ref e0_s_init - * - * \snippet Vector/3_molecular_dynamic/main.cpp init - * - */ - - //! \cond [init] \endcond - - openfpm_init(&argc,&argv); - Vcluster & v_cl = create_vcluster(); // we will use it do place particles on a 10x10x10 Grid like size_t sz[3] = {10,10,10}; @@ -334,7 +316,14 @@ int main(int argc, char* argv[]) // ghost, big enough to contain the interaction radius Ghost<3,float> ghost(r_cut); - //! \cond [init] \endcond + double dt = 0.0005; + double sigma12 = pow(sigma,12); + double sigma6 = pow(sigma,6); + + openfpm::vector x; + openfpm::vector> y; + + //! \cond [constants run] \endcond /*! * \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list @@ -370,18 +359,24 @@ int main(int argc, char* argv[]) //! \cond [vect grid] \endcond + // We create the grid iterator auto it = vd.getGridIterator(sz); while (it.isNext()) { + // Create a new particle vd.add(); + // key contain (i,j,k) index of the grid auto key = it.get(); + // The index of the grid can be accessed with key.get(0) == i, key.get(1) == j ... + // We use getLastPos to set the position of the last particle added vd.getLastPos()[0] = key.get(0) * it.getSpacing(0); vd.getLastPos()[1] = key.get(1) * it.getSpacing(1); vd.getLastPos()[2] = key.get(2) * it.getSpacing(2); + // We use getLastProp to set the property value of the last particle we added vd.template getLastProp()[0] = 0.0; vd.template getLastProp()[1] = 0.0; vd.template getLastProp()[2] = 0.0; @@ -404,8 +399,8 @@ int main(int argc, char* argv[]) * * The verlet integration stepping look like this * - * \f[ \vec{v}(t_{n+1/2}) = \vec{v}_p(t_n) + \frac{1}{2} \delta t \vec{a}(t_n) \f] - * \f[ \vec{x}(t_{n}) = \vec{x}_p(t_n) + \delta t \vec{v}(t_n+1/2) \f] + * \f[ \vec{v}(t_{n}+1/2) = \vec{v}_p(t_n) + \frac{1}{2} \delta t \vec{a}(t_n) \f] + * \f[ \vec{x}(t_{n}+1) = \vec{x}_p(t_n) + \delta t \vec{v}(t_n+1/2) \f] * * calculate the forces from \f$\vec{a} (t_{n}) \f$ finally * @@ -465,7 +460,7 @@ int main(int argc, char* argv[]) ++it3; } - // Because we mooved the particles in space we have to map them and re-sync the ghost + // Because we moved the particles in space we have to map them and re-sync the ghost vd.map(); vd.template ghost_get<>(); @@ -488,7 +483,7 @@ int main(int argc, char* argv[]) ++it4; } - // After every iteration collect some statistic about the confoguration + // After every iteration collect some statistic about the configuration if (i % 100 == 0) { // We write the particle position for visualization (Without ghost) @@ -554,6 +549,15 @@ int main(int argc, char* argv[]) // width of the line options.lineWidth = 1.0; + // Resolution in x + options.width = 1280; + + // Resolution in y + options.heigh = 720; + + // Add zoom capability + options.more = GC_ZOOM; + // Object that draw the X Y graph GoogleChart cg; diff --git a/example/Vector/4_multiphase_celllist_verlet/main.cpp b/example/Vector/4_multiphase_celllist_verlet/main.cpp index de4093d5ddf09a92a91ad307cfb085fb2575ab8c..cfa51c3139d674de2c9820ab4708c868fd06460f 100644 --- a/example/Vector/4_multiphase_celllist_verlet/main.cpp +++ b/example/Vector/4_multiphase_celllist_verlet/main.cpp @@ -149,6 +149,8 @@ int main(int argc, char* argv[]) //! \cond [create multi-phase verlet] \endcond + { + // Get the cell list of the phase1 auto CL_phase1 = phases.get(1).getCellList(r_cut); @@ -204,6 +206,8 @@ int main(int argc, char* argv[]) ++it; } + } + //! \cond [count part from phase0 to 1] \endcond /*! @@ -311,13 +315,14 @@ int main(int argc, char* argv[]) * */ + { //! \cond [compute sym multi-phase two phase] \endcond // Get the cell list of the phase1 - CL_phase1 = phases.get(1).getCellListSym(r_cut); + auto CL_phase1 = phases.get(1).getCellListSym(r_cut); // This function create a Verlet-list between phases 0 and 1 - NN_ver01 = createVerletSym(phases.get(0),phases.get(1),CL_phase1,r_cut); + auto NN_ver01 = createVerletSym(phases.get(0),phases.get(1),CL_phase1,r_cut); // Get an iterator over the real and ghost particles it = phases.get(0).getDomainAndGhostIterator(); @@ -371,6 +376,8 @@ int main(int argc, char* argv[]) phases.get(0).ghost_put(); phases.get(1).ghost_put(); + } + //! \cond [compute sym multi-phase two phase] \endcond /*! diff --git a/example/Vector/4_reorder/energy_force.hpp b/example/Vector/4_reorder/energy_force.hpp index 5836b0fc8587a5f5e8bb0a2bb4f813af22531449..2108b811aaaf8b1c4363b9c89b2020dffa261bc2 100644 --- a/example/Vector/4_reorder/energy_force.hpp +++ b/example/Vector/4_reorder/energy_force.hpp @@ -12,7 +12,7 @@ constexpr int velocity = 0; constexpr int force = 1; -void calc_forces(vector_dist<3,double, aggregate > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6) +template void calc_forces(vector_dist<3,double, aggregate > & vd, CellList & NN, double sigma12, double sigma6) { // update the Cell-list vd.updateCellList(NN); @@ -73,7 +73,7 @@ void calc_forces(vector_dist<3,double, aggregate > & vd, Ce } -double calc_energy(vector_dist<3,double, aggregate > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6) +template double calc_energy(vector_dist<3,double, aggregate > & vd, CellList & NN, double sigma12, double sigma6) { double E = 0.0; diff --git a/example/Vector/4_reorder/main_comp_ord.cpp b/example/Vector/4_reorder/main_comp_ord.cpp index 9a9a1abf48e2597408482e2c1c1abb56cb17b851..77f915a464809300eb473282c42439c404a8f2a7 100644 --- a/example/Vector/4_reorder/main_comp_ord.cpp +++ b/example/Vector/4_reorder/main_comp_ord.cpp @@ -45,7 +45,7 @@ constexpr int force = 1; //! \cond [calc forces] \endcond -void calc_forces(vector_dist<3,double, aggregate > & vd, CellList_hilb<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6) +template void calc_forces(vector_dist<3,double, aggregate > & vd, CellList & NN, double sigma12, double sigma6) { //! \cond [calc forces] \endcond @@ -140,7 +140,7 @@ void calc_forces(vector_dist<3,double, aggregate > & vd, Ce //! \cond [calc energy all] \endcond -double calc_energy(vector_dist<3,double, aggregate > & vd, CellList_hilb<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6) +template double calc_energy(vector_dist<3,double, aggregate > & vd, CellList & NN, double sigma12, double sigma6) { double E = 0.0; diff --git a/example/Vector/5_molecular_dynamic_sym_crs/main.cpp b/example/Vector/5_molecular_dynamic_sym_crs/main.cpp index fe6bec6ba6215df5bab4fa418542dc06edd351a5..a8eb56155e0f78d6edc83825d752fe9a06de809c 100644 --- a/example/Vector/5_molecular_dynamic_sym_crs/main.cpp +++ b/example/Vector/5_molecular_dynamic_sym_crs/main.cpp @@ -148,7 +148,7 @@ void calc_forces(vector_dist<3,double, aggregate > & vd, Ve //! \cond [real and ghost] \endcond // Get an iterator over particles - auto it2 = vd.getParticleIteratorCRS(NN.getInternalCellList()); + auto it2 = vd.getParticleIteratorCRS(NN); //! \cond [real and ghost] \endcond @@ -250,7 +250,7 @@ double calc_energy(vector_dist<3,double, aggregate > & vd, double shift = 4.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) ); // Get an iterator over particles - auto it2 = vd.getParticleIteratorCRS(NN.getInternalCellList()); + auto it2 = vd.getParticleIteratorCRS(NN); // For each particle ... while (it2.isNext()) diff --git a/example/Vector/7_SPH_dlb_opt/main.cpp b/example/Vector/7_SPH_dlb_opt/main.cpp index cbb408053bd2042f3b78c5bd5b8bbee314305ae7..6c54fbd9129fca5466047a03ae560fed4f606036 100644 --- a/example/Vector/7_SPH_dlb_opt/main.cpp +++ b/example/Vector/7_SPH_dlb_opt/main.cpp @@ -912,7 +912,7 @@ int main(int argc, char* argv[]) cbar = coeff_sound * sqrt(gravity * h_swl); // for each particle inside the fluid box ... -/* while (fluid_it.isNext()) + while (fluid_it.isNext()) { // ... add a particle ... vd.add(); @@ -946,8 +946,7 @@ int main(int argc, char* argv[]) // next fluid particle ++fluid_it; - }*/ - + } // Recipient Box<3,double> recipient1({0.0,0.0,0.0},{1.6+dp/2.0,0.67+dp/2.0,0.4+dp/2.0}); @@ -1010,10 +1009,6 @@ int main(int argc, char* argv[]) vd.map(); - vd.write("Recipient"); - - return 0; - // Now that we fill the vector with particles ModelCustom md; diff --git a/openfpm_data b/openfpm_data index 3a732c98d6404fd81529c25af4fbf33f8ca3de9c..2c31320305af04b0621d2a10818354d064ec1bbf 160000 --- a/openfpm_data +++ b/openfpm_data @@ -1 +1 @@ -Subproject commit 3a732c98d6404fd81529c25af4fbf33f8ca3de9c +Subproject commit 2c31320305af04b0621d2a10818354d064ec1bbf diff --git a/openfpm_devices b/openfpm_devices index add2acfe9623e6045da1e41adbab6dc57ef50e5f..9817279ee3ac36b435d2178e994ab314a01aab30 160000 --- a/openfpm_devices +++ b/openfpm_devices @@ -1 +1 @@ -Subproject commit add2acfe9623e6045da1e41adbab6dc57ef50e5f +Subproject commit 9817279ee3ac36b435d2178e994ab314a01aab30 diff --git a/openfpm_io b/openfpm_io index a6bd0b797704f6faba662a10f5b101c218257490..9b0500dd1347d9f3d2898f2a9206220eea7907d5 160000 --- a/openfpm_io +++ b/openfpm_io @@ -1 +1 @@ -Subproject commit a6bd0b797704f6faba662a10f5b101c218257490 +Subproject commit 9b0500dd1347d9f3d2898f2a9206220eea7907d5 diff --git a/openfpm_vcluster b/openfpm_vcluster index 1856e665949ecca8dbcf68e08984aace853aad8d..7030e7f76e5d80a297e4cc9b2d4f3ad98beb8d0f 160000 --- a/openfpm_vcluster +++ b/openfpm_vcluster @@ -1 +1 @@ -Subproject commit 1856e665949ecca8dbcf68e08984aace853aad8d +Subproject commit 7030e7f76e5d80a297e4cc9b2d4f3ad98beb8d0f diff --git a/src/Decomposition/CartDecomposition.hpp b/src/Decomposition/CartDecomposition.hpp index c097931ddc9a8fef0ea0800acc071655064eeb7e..728292d04deab6e3621eededd970c05a65490575 100755 --- a/src/Decomposition/CartDecomposition.hpp +++ b/src/Decomposition/CartDecomposition.hpp @@ -135,7 +135,7 @@ protected: //! Structure that store the cartesian grid information grid_sm gr_dist; - //! Structure that decompose your structure into cell without creating them + //! Structure that decompose the space into cells without creating them //! useful to convert positions to CellId or sub-domain id in this case CellDecomposer_sm> cd; @@ -1724,6 +1724,16 @@ public: return dist.get_ndec(); } + /*! \brief Get the cell decomposer of the decomposition + * + * \return the cell decomposer + * + */ + const CellDecomposer_sm> & getCellDecomposer() + { + return cd; + } + //! friend classes friend extended_type; diff --git a/src/Decomposition/common.hpp b/src/Decomposition/common.hpp index 36575f125be68775f75b61b73cd5fcaed61f749d..f9cdfeb61108b4d7c52edd10bbe14544b7c6399a 100755 --- a/src/Decomposition/common.hpp +++ b/src/Decomposition/common.hpp @@ -114,7 +114,7 @@ struct lBox_dom { //! Intersection between the local sub-domain enlarged by the ghost and the contiguous processor //! sub-domains (External ghost) - openfpm::vector_std< Box_sub > ebx; + openfpm::vector_std< Box_sub_k > ebx; //! Intersection between the contiguous processor sub-domain enlarged by the ghost with the //! local sub-domain (Internal ghost) diff --git a/src/Decomposition/ie_ghost.hpp b/src/Decomposition/ie_ghost.hpp index cf74851475b13662fd810a65ae82c7625b20dbcb..a7b14d91e18f6956b156cf55bf92523d4a38754b 100755 --- a/src/Decomposition/ie_ghost.hpp +++ b/src/Decomposition/ie_ghost.hpp @@ -991,7 +991,7 @@ public: return true; } - /*! \brief Check if the ie_loc_ghosts contain the same information + /*! \brief Check if the ie_ghosts contain the same information * * \param ig Element to check * @@ -1065,52 +1065,6 @@ public: */ bool is_equal_ng(ie_ghost & ig) { - Box bt; - - if (getNEGhostBox() != ig.getNEGhostBox()) - return false; - - if (getNIGhostBox() != ig.getNIGhostBox()) - return false; - - for (size_t i = 0 ; i < proc_int_box.size() ; i++) - { - if (getProcessorNIGhost(i) != ig.getProcessorNIGhost(i)) - return false; - for (size_t j = 0 ; j < getProcessorNIGhost(i) ; j++) - { - if (getProcessorIGhostBox(i,j).Intersect(ig.getProcessorIGhostBox(i,j),bt) == false) - return false; - if (getProcessorIGhostId(i,j) != ig.getProcessorIGhostId(i,j)) - return false; - if (getProcessorIGhostSub(i,j) != ig.getProcessorIGhostSub(i,j)) - return false; - } - if (getIGhostBox(i).Intersect(ig.getIGhostBox(i),bt) == false) - return false; - if (getIGhostBoxProcessor(i) != ig.getIGhostBoxProcessor(i)) - return false; - } - - for (size_t i = 0 ; i < proc_int_box.size() ; i++) - { - if (getProcessorNEGhost(i) != ig.getProcessorNEGhost(i)) - return false; - for (size_t j = 0 ; j < getProcessorNEGhost(i) ; j++) - { - if (getProcessorEGhostBox(i,j).Intersect(ig.getProcessorEGhostBox(i,j),bt) == false) - return false; - if (getProcessorEGhostId(i,j) != ig.getProcessorEGhostId(i,j)) - return false; - if (getProcessorEGhostSub(i,j) != ig.getProcessorEGhostSub(i,j)) - return false; - } - if (getEGhostBox(i).Intersect(ig.getEGhostBox(i),bt) == false) - return false; - if (getEGhostBoxProcessor(i) != ig.getEGhostBoxProcessor(i)) - return false; - } - return true; } diff --git a/src/Decomposition/ie_loc_ghost.hpp b/src/Decomposition/ie_loc_ghost.hpp index 1d9bd57064d4d28d73a158be71511bceac7c46fe..727b1dca1d063db2be6ff4b4b88d8e0405320493 100755 --- a/src/Decomposition/ie_loc_ghost.hpp +++ b/src/Decomposition/ie_loc_ghost.hpp @@ -70,7 +70,7 @@ class ie_loc_ghost if (intersect == true) { - Box_sub b; + Box_sub_k b; b.sub = rj; b.bx = bi; b.cmb = sub_domains_prc.get(j).cmb; @@ -84,6 +84,7 @@ class ie_loc_ghost if (loc_ghost_box.get(rj).ibx.get(k).sub == i && loc_ghost_box.get(rj).ibx.get(k).cmb == sub_domains_prc.get(j).cmb.operator-()) { loc_ghost_box.get(rj).ibx.get(k).k = loc_ghost_box.get(i).ebx.size()-1; + loc_ghost_box.get(i).ebx.last().k = k; break; } } @@ -346,9 +347,11 @@ public: return loc_ghost_box.get(id).ibx.size(); } - /*! \brief For the sub-domain i intersected with the sub-domain j enlarged, the associated - * external ghost box is located in getLocalEGhostBox(j,k) with - * getLocalIGhostSub(j,k) == i, this function return k + /*! \brief For the sub-domain i intersected with a surrounding sub-domain enlarged. Produce a internal ghost box from + * the prospecive of i and an associated external ghost box from the prospective of j. + * In order to retrieve the information about the external ghost box we have to use getLocalEGhostBox(x,k). + * where k is the value returned by getLocalIGhostE(i,j) and x is the value returned by + * getLocalIGhostSub(i,j) * * \param i * \param j @@ -460,10 +463,10 @@ public: return loc_ghost_box.get(i).ebx.get(j).cmb; } - /*! \brief Considering that sub-domain has N internal local ghost box identified + /*! \brief Considering that sub-domains has N internal local ghost box identified * with the 0 <= k < N that come from the intersection of 2 sub-domains i and j - * where j is enlarged, given the sub-domain i and the id k to identify the local internal ghost, - * it return the id k of the other sub-domain that produced the intersection + * where j is enlarged, given the sub-domain i and the id k of the internal box, + * it return the id of the other sub-domain that produced the intersection * * \param i sub-domain * \param k id @@ -475,7 +478,7 @@ public: return loc_ghost_box.get(i).ibx.get(k).sub; } - /*! \brief Considering that sub-domain has N external local ghost box identified + /*! \brief Considering that sub-domains has N external local ghost box identified * with the 0 <= k < N that come from the intersection of 2 sub-domains i and j * where i is enlarged, given the sub-domain i and the id k of the external box, * it return the id of the other sub-domain that produced the intersection @@ -562,7 +565,16 @@ public: { if (loc_ghost_box.get(i).ibx.get(j).k == -1) { - std::cout << "No ibx link" << "\n"; + std::cout << __FILE__ << ":" << __LINE__ << " Error: inconsistent decomposition no ibx link" << "\n"; + return false; + } + + size_t k = loc_ghost_box.get(i).ibx.get(j).k; + size_t sub = loc_ghost_box.get(i).ibx.get(j).sub; + + if (loc_ghost_box.get(sub).ebx.get(k).k != (long int)j) + { + std::cout << __FILE__ << ":" << __LINE__ << " Error: inconsistent link between an external ghost box and an internal ghost box" << "\n"; return false; } } @@ -627,39 +639,6 @@ public: */ bool is_equal_ng(ie_loc_ghost & ilg) { - Box bt; - - if (ilg.loc_ghost_box.size() != loc_ghost_box.size()) - return false; - - // Explore all the subdomains - for (size_t i = 0 ; i < loc_ghost_box.size() ; i++) - { - if (getLocalNIGhost(i) != ilg.getLocalNIGhost(i)) - return false; - - if (getLocalNEGhost(i) != ilg.getLocalNEGhost(i)) - return false; - - for (size_t j = 0 ; j < getLocalNIGhost(i) ; j++) - { - if (getLocalIGhostE(i,j) != ilg.getLocalIGhostE(i,j)) - return false; - if (getLocalIGhostBox(i,j).Intersect(ilg.getLocalIGhostBox(i,j),bt) == false) - return false; - if (getLocalIGhostSub(i,j) != ilg.getLocalIGhostSub(i,j)) - return false; - } - for (size_t j = 0 ; j < getLocalNEGhost(i) ; j++) - { - if (getLocalEGhostBox(i,j).Intersect(ilg.getLocalEGhostBox(i,j),bt) == false) - return false; - if (getLocalEGhostSub(i,j) != ilg.getLocalEGhostSub(i,j)) - return false; - } - - } - return true; } diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp index 13239393c4695668bec48ed741ebb9f94517fbbd..88cb2af8ff1774fff7e415512dfc47862df9377a 100644 --- a/src/Grid/grid_dist_id.hpp +++ b/src/Grid/grid_dist_id.hpp @@ -376,73 +376,14 @@ class grid_dist_id : public grid_dist_id_comm void ghost_get_local() - { - //! For all the sub-domains - for (size_t i = 0 ; i < loc_ig_box.size() ; i++) - { - //! For all the internal ghost boxes of each sub-domain - for (size_t j = 0 ; j < loc_ig_box.get(i).bid.size() ; j++) - { - Box bx_src = loc_ig_box.get(i).bid.get(j).box; - // convert into local - bx_src -= gdb_ext.get(i).origin; - - // sub domain connected with external box - size_t sub_id_dst = loc_ig_box.get(i).bid.get(j).sub; - - // local internal ghost box connected - size_t k = loc_ig_box.get(i).bid.get(j).k; - - Box bx_dst = loc_eg_box.get(sub_id_dst).bid.get(k).box; - - // convert into local - bx_dst -= gdb_ext.get(sub_id_dst).origin; - - // create 2 sub grid iterator - - if (bx_dst.isValid() == false) - continue; - - grid_key_dx_iterator_sub sub_src(loc_grid.get(i).getGrid(),bx_src.getKP1(),bx_src.getKP2()); - grid_key_dx_iterator_sub sub_dst(loc_grid.get(sub_id_dst).getGrid(),bx_dst.getKP1(),bx_dst.getKP2()); - -#ifdef SE_CLASS1 - - if (loc_eg_box.get(sub_id_dst).bid.get(k).sub != i) - std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " source and destination are not correctly linked" << "\n"; - - if (sub_src.getVolume() != sub_dst.getVolume()) - std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " source and destination does not match in size" << "\n"; - -#endif - - const auto & gs = loc_grid.get(i); - auto & gd = loc_grid.get(sub_id_dst); - - while (sub_src.isNext()) - { - // Option 1 - gd.set(sub_dst.get(),gs,sub_src.get()); - - ++sub_src; - ++sub_dst; - } - } - } - } - /*! \brief Check the grid has a valid size * * \param g_sz size of the grid @@ -625,8 +566,8 @@ protected: // set the ghost for (size_t i = 0 ; i < dim ; i++) { - gc.setLow(i,-sp.getHigh(i)); - gc.setHigh(i,sp.getHigh(i)); + gc.setLow(i,gd.getLow(i)*(sp.getHigh(i))); + gc.setHigh(i,gd.getHigh(i)*(sp.getHigh(i))); } return gc; @@ -1018,6 +959,8 @@ public: return total; } + + /*! \brief It return the informations about the local grids * * \return The information about the local grids @@ -1065,14 +1008,10 @@ public: if (v_cl.getProcessUnitID() == 0) { for (size_t i = 1; i < v_cl.getProcessingUnits(); i++) - { v_cl.send(i,0,gdb_ext_global); - } } else - { v_cl.recv(0,0,gdb_ext_global); - } v_cl.execute(); } @@ -1208,6 +1147,33 @@ public: return false; } + /*! \brief Get the reference of the selected element + * + * \tparam p property to get (is an integer) + * \param v1 grid_key that identify the element in the grid + * + * \return the selected element + * + */ + template inline auto getProp(const grid_dist_key_dx & v1) const -> decltype(this->template get

(v1)) + { + return this->template get

(v1); + } + + /*! \brief Get the reference of the selected element + * + * \tparam p property to get (is an integer) + * \param v1 grid_key that identify the element in the grid + * + * \return the selected element + * + */ + template inline auto getProp(const grid_dist_key_dx & v1) -> decltype(this->template get

(v1)) + { + return this->template get

(v1); + } + + /*! \brief Get the reference of the selected element * * \tparam p property to get (is an integer) @@ -1240,12 +1206,6 @@ public: return loc_grid.get(v1.getSub()).template get