SparseGrid.hpp 45.4 KB
Newer Older
incardon's avatar
incardon committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/*
 * SparseGrid.hpp
 *
 *  Created on: Oct 22, 2017
 *      Author: i-bird
 */

#ifndef OPENFPM_DATA_SRC_SPARSEGRID_SPARSEGRID_HPP_
#define OPENFPM_DATA_SRC_SPARSEGRID_SPARSEGRID_HPP_

#include "memory_ly/memory_array.hpp"
#include "memory_ly/memory_c.hpp"
#include "memory_ly/memory_conf.hpp"
#include "hash_map/hopscotch_map.h"
#include "hash_map/hopscotch_set.h"
#include "Vector/map_vector.hpp"
#include "util/variadic_to_vmpl.hpp"
#include "data_type/aggregate.hpp"
#include "SparseGridUtil.hpp"
#include "SparseGrid_iterator.hpp"
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
// We do not want parallel writer
#define NO_PARALLEL
#include "VTKWriter/VTKWriter.hpp"


template<typename Tsrc,typename Tdst>
class copy_prop_to_vector
{
	//! source
	Tsrc src;

	//! destination
	Tdst dst;

	size_t pos;

public:

	copy_prop_to_vector(Tsrc src, Tdst dst,size_t pos)
	:src(src),dst(dst),pos(pos)
	{}

	//! It call the copy function for each property
	template<typename T>
	inline void operator()(T& t) const
	{
		typedef typename std::remove_reference<decltype(dst.template get<T::value>())>::type copy_rtype;

		meta_copy<copy_rtype>::meta_copy_(src.template get<T::value>()[pos],dst.template get<T::value>());
	}

};
incardon's avatar
incardon committed
53

54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

template<unsigned int dim, typename Tsrc,typename Tdst>
class copy_sparse_to_sparse
{
	//! source
	const Tsrc & src;

	//! destination
	Tdst & dst;

	//! source position
	grid_key_dx<dim> & pos_src;

	//! destination position
	grid_key_dx<dim> & pos_dst;

public:

	copy_sparse_to_sparse(const Tsrc & src, Tdst & dst,grid_key_dx<dim> & pos_src, grid_key_dx<dim> & pos_dst)
	:src(src),dst(dst),pos_src(pos_src),pos_dst(pos_dst)
	{}

	//! It call the copy function for each property
	template<typename T>
	inline void operator()(T& t) const
	{
		typedef typename std::remove_reference<decltype(dst.template insert<T::value>(pos_dst))>::type copy_rtype;

		meta_copy<copy_rtype>::meta_copy_(src.template get<T::value>(pos_src),dst.template insert<T::value>(pos_dst));
	}

};

incardon's avatar
incardon committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
template< template<typename,typename> class op,unsigned int dim, typename Tsrc,typename Tdst, unsigned int ... prp>
class copy_sparse_to_sparse_op
{
	//! source
	const Tsrc & src;

	//! destination
	Tdst & dst;

	//! source position
	grid_key_dx<dim> & pos_src;

	//! destination position
	grid_key_dx<dim> & pos_dst;

	//! Convert the packed properties into an MPL vector
	typedef typename to_boost_vmpl<prp...>::type v_prp;

public:

	copy_sparse_to_sparse_op(const Tsrc & src, Tdst & dst,grid_key_dx<dim> & pos_src, grid_key_dx<dim> & pos_dst)
	:src(src),dst(dst),pos_src(pos_src),pos_dst(pos_dst)
	{}

	//! It call the copy function for each property
	template<typename T>
	inline void operator()(T& t) const
	{
		typedef typename boost::mpl::at<v_prp,boost::mpl::int_<T::value>>::type idx_type;
		typedef typename std::remove_reference<decltype(dst.template insert<idx_type::value>(pos_dst))>::type copy_rtype;

incardon's avatar
incardon committed
118
119
120
121
		if (dst.existPoint(pos_dst) == false)
		{meta_copy_op<replace_,copy_rtype>::meta_copy_op_(src.template get<idx_type::value>(pos_src),dst.template insert<idx_type::value>(pos_dst));}
		else
		{meta_copy_op<op,copy_rtype>::meta_copy_op_(src.template get<idx_type::value>(pos_src),dst.template insert<idx_type::value>(pos_dst));}
incardon's avatar
incardon committed
122
123
124
125
	}

};

incardon's avatar
incardon committed
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
/*! \brief this class is a functor for "for_each" algorithm
 *
 * This class is a functor for "for_each" algorithm. For each
 * element of the boost::vector the operator() is called.
 * Is mainly used to copy a boost::mpl::vector into runtime array
 *
 */

template<unsigned int dim, typename mpl_v>
struct copy_sz
{
	//! sz site_t
	size_t (& sz)[dim];


	/*! \brief constructor
	 *
	 * \param sz runtime sz to fill
	 *
	 */
	inline copy_sz(size_t (& sz)[dim])
	:sz(sz)
	{
	};

	//! It call the copy function for each property
	template<typename T>
	inline void operator()(T& t) const
	{
		sz[T::value] = boost::mpl::at<mpl_v,boost::mpl::int_<T::value>>::type::value;
	}
};

#define BIT_SHIFT_SIZE_T 6

template<unsigned int dim, typename T, typename S, typename chunking>
class sgrid_cpu<dim,T,S,typename memory_traits_lin<T>::type,chunking>
{
	//! cache pointer
165
	mutable size_t cache_pnt;
incardon's avatar
incardon committed
166
167
168
169
170

	//! background values
	T background;

	//! cache
171
	mutable long int cache[SGRID_CACHE];
incardon's avatar
incardon committed
172
173

	//! cached id
174
	mutable long int cached_id[SGRID_CACHE];
incardon's avatar
incardon committed
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196

	//! Map to convert from grid coordinates to chunk
	tsl::hopscotch_map<size_t, size_t> map;

	//! indicate which element in the chunk are really filled
	openfpm::vector<cheader<dim,chunking::size::value>> header;

	//Definition of the chunks
	typedef typename v_transform_two<Ft_chunk,boost::mpl::int_<chunking::size::value>,typename T::type>::type chunk_def;

	//! vector of chunks
	openfpm::vector<aggregate_bfv<chunk_def>> chunks;

	//! grid size information
	grid_sm<dim,void> g_sm;

	//! grid size information with shift
	grid_sm<dim,void> g_sm_shift;

	//! conversion position in the chunks
	grid_key_dx<dim> pos_chunk[chunking::size::value];

197
198
199
	//! size of the chunk
	size_t sz_cnk[dim];

200
201
	openfpm::vector<size_t> empty_v;

incardon's avatar
incardon committed
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
	/*! \brief Remove
	 *
	 *
	 */
	template<unsigned int n_ele>
	inline void remove_from_chunk(size_t sub_id,
			 	 	 	 	 	  size_t & nele,
								  size_t (& mask)[n_ele / (sizeof(size_t)*8) + (n_ele % (sizeof(size_t)*8) != 0) + 1])
	{
		// set the mask to null
		size_t mask_check = (size_t)1 << (sub_id & ((1 << BIT_SHIFT_SIZE_T) - 1));

		size_t swt = mask[sub_id >> BIT_SHIFT_SIZE_T] & mask_check;

		nele = (swt)?nele-1:nele;

		mask[sub_id >> BIT_SHIFT_SIZE_T] &= ~mask_check;
	}

	/*! \brief reconstruct the map
	 *
	 *
	 *
	 */
	inline void reconstruct_map()
	{
		// reconstruct map

		map.clear();
		for (size_t i = 0 ; i < header.size() ; i++)
		{
			grid_key_dx<dim> kh = header.get(i).pos;
			grid_key_dx<dim> kl;

			// shift the key
			key_shift<dim,chunking>::shift(kh,kl);

			long int lin_id = g_sm_shift.LinId(kh);

			map[lin_id] = i;
		}
	}

245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
	/*! \brief Eliminate empty chunks
	 *
	 * \warning Because this operation is time consuming it perform the operation once
	 *          we reach a critical size in the list of the empty chunks
	 *
	 */
	inline void remove_empty()
	{
		if (empty_v.size() >= FLUSH_REMOVE)
		{
			// eliminate double entry

			empty_v.sort();
			empty_v.unique();

			// Because chunks can be refilled the empty list can contain chunks that are
			// filled so before remove we have to check that they are really empty

			for (int i = empty_v.size() - 1 ; i >= 0  ; i--)
			{
				if (header.get(empty_v.get(i)).nele != 0)
				{empty_v.remove(i);}
			}

			header.remove(empty_v);
			chunks.remove(empty_v);

			// reconstruct map

incardon's avatar
incardon committed
274
			reconstruct_map();
275
276
277
278
279

			empty_v.clear();

			// cache must be cleared

incardon's avatar
incardon committed
280
			clear_cache();
281
282
283
		}
	}

incardon's avatar
incardon committed
284
285
286
287
288
289
	/*! \brief add on cache
	 *
	 * \param lin_id linearized id
	 * \param active_cnk active chunk
	 *
	 */
290
	inline void add_on_cache(size_t lin_id, size_t active_cnk) const
incardon's avatar
incardon committed
291
292
293
294
295
296
297
	{
		// Add on cache the chunk
		cache[cache_pnt] = lin_id;
		cached_id[cache_pnt] = active_cnk;
		cache_pnt++;
		cache_pnt = (cache_pnt >= SGRID_CACHE)?0:cache_pnt;
	}
incardon's avatar
incardon committed
298

incardon's avatar
incardon committed
299
	/*! \brief reset the cache
incardon's avatar
incardon committed
300
301
302
	 *
	 *
	 */
incardon's avatar
incardon committed
303
	inline void clear_cache()
incardon's avatar
incardon committed
304
	{
incardon's avatar
incardon committed
305
		cache_pnt = 0;
incardon's avatar
incardon committed
306
307
		for (size_t i = 0 ; i < SGRID_CACHE ; i++)
		{cache[i] = -1;}
incardon's avatar
incardon committed
308
	}
incardon's avatar
incardon committed
309

incardon's avatar
incardon committed
310
311
312
313
314
315
316
317
	/*! \brief set the grid shift from size
	 *
	 * \param sz size of the grid
	 * \param g_sm_shift chunks grid size
	 *
	 */
	void set_g_shift_from_size(const size_t (& sz)[dim], grid_sm<dim,void> & g_sm_shift)
	{
incardon's avatar
incardon committed
318
319
320
321
322
323
324
325
326
327
328
329
330
331
		grid_key_dx<dim> cs;
		grid_key_dx<dim> unused;

		for (size_t i = 0 ; i < dim ; i++)
		{cs.set_d(i,sz[i]);}

		key_shift<dim,chunking>::shift(cs,unused);

		size_t sz_i[dim];

		for (size_t i = 0 ; i < dim ; i++)
		{sz_i[i] = cs.get(i) + 1;}

		g_sm_shift.setDimensions(sz_i);
incardon's avatar
incardon committed
332
333
334
335
336
337
338
339
340
341
	}

	/*! \brief initialize
	 *
	 *
	 */
	void init()
	{
		for (size_t i = 0 ; i < SGRID_CACHE ; i++)
		{cache[i] = -1;}
incardon's avatar
incardon committed
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366

		// fill pos_g

		copy_sz<dim,typename chunking::type> cpsz(sz_cnk);
		boost::mpl::for_each_ref< boost::mpl::range_c<int,0,dim> >(cpsz);

		grid_sm<dim,void> gs(sz_cnk);

		grid_key_dx_iterator<dim> it(gs);
		size_t cnt = 0;

		while (it.isNext())
		{
			auto key = it.get();

			for (size_t i = 0 ; i < dim ; i++)
			{
				pos_chunk[cnt].set_d(i,key.get(i));
			}

			++cnt;
			++it;
		}
	}

incardon's avatar
incardon committed
367
368
369
370
371
372
373
	/*! \brief Before insert data you have to do this
	 *
	 * \param v1 grid key where you want to insert data
	 * \param active_cnk output which chunk is competent for this data
	 * \param sub_id element inside the chunk
	 *
	 */
incardon's avatar
incardon committed
374
	inline bool pre_insert(const grid_key_dx<dim> & v1, size_t & active_cnk, size_t & sub_id)
incardon's avatar
incardon committed
375
	{
incardon's avatar
incardon committed
376
		bool exist = true;
incardon's avatar
incardon committed
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
		active_cnk = 0;

		grid_key_dx<dim> kh = v1;
		grid_key_dx<dim> kl;

		// shift the key
		key_shift<dim,chunking>::shift(kh,kl);

		long int lin_id = g_sm_shift.LinId(kh);

		size_t id = 0;
		for (size_t k = 0 ; k < SGRID_CACHE; k++)
		{id += (cache[k] == lin_id)?k+1:0;}

		if (id == 0)
		{
			// we do not have it in cache we check if we have it in the map

			auto fnd = map.find(lin_id);
			if (fnd == map.end())
			{
				// we do not have it in the map create a chunk

				map[lin_id] = chunks.size();
				chunks.add();
				header.add();
				header.last().pos = kh;
				header.last().nele = 0;

				// set the mask to null
				auto & h = header.last().mask;

				for (size_t i = 0 ; i < chunking::size::value / (sizeof(size_t)*8) + (chunking::size::value % (sizeof(size_t)*8) != 0) + 1 ; i++)
				{h[i] = 0;}

				key_shift<dim,chunking>::cpos(header.last().pos);

				active_cnk = chunks.size() - 1;
			}
			else
			{
				// we have it in the map

				active_cnk = fnd->second;
			}

			// Add on cache the chunk
			cache[cache_pnt] = lin_id;
			cached_id[cache_pnt] = active_cnk;
			cache_pnt++;
			cache_pnt = (cache_pnt >= SGRID_CACHE)?0:cache_pnt;
		}
		else
		{
			active_cnk = cached_id[id-1];
			cache_pnt = id;
			cache_pnt = (cache_pnt == SGRID_CACHE)?0:cache_pnt;
		}

		sub_id = sublin<dim,typename chunking::shift_c>::lin(kl);

		// the chunk is in cache, solve

		// we notify that we added one element
		auto & h = header.get(active_cnk);

		// we set the mask

		size_t mask_check = (size_t)1 << (sub_id & ((1 << BIT_SHIFT_SIZE_T) - 1));

incardon's avatar
incardon committed
447
448
		exist = h.mask[sub_id >> BIT_SHIFT_SIZE_T] & mask_check;
		h.nele = (exist)?h.nele:h.nele + 1;
incardon's avatar
incardon committed
449
		h.mask[sub_id >> BIT_SHIFT_SIZE_T] |= mask_check;
incardon's avatar
incardon committed
450
451

		return exist;
incardon's avatar
incardon committed
452
453
	}

incardon's avatar
incardon committed
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
	inline void remove_point(const grid_key_dx<dim> & v1)
	{
		size_t active_cnk = 0;

		grid_key_dx<dim> kh = v1;
		grid_key_dx<dim> kl;

		// shift the key
		key_shift<dim,chunking>::shift(kh,kl);

		long int lin_id = g_sm_shift.LinId(kh);

		size_t id = 0;
		for (size_t k = 0 ; k < SGRID_CACHE; k++)
		{id += (cache[k] == lin_id)?k+1:0;}

		if (id == 0)
		{
			// we do not have it in cache we check if we have it in the map

			auto fnd = map.find(lin_id);
			if (fnd == map.end())
			{return;}
			else
			{active_cnk = fnd->second;}

			// Add on cache the chunk
			cache[cache_pnt] = lin_id;
			cached_id[cache_pnt] = active_cnk;
			cache_pnt++;
			cache_pnt = (cache_pnt >= SGRID_CACHE)?0:cache_pnt;
		}
		else
		{
			active_cnk = cached_id[id-1];
			cache_pnt = id;
			cache_pnt = (cache_pnt == SGRID_CACHE)?0:cache_pnt;
		}

		size_t sub_id = sublin<dim,typename chunking::shift_c>::lin(kl);

		// eliminate the element

		// set the mask to null
		size_t mask_check = (size_t)1 << (sub_id & ((1 << BIT_SHIFT_SIZE_T) - 1));

		auto & h = header.get(active_cnk);
		size_t swt = h.mask[sub_id >> BIT_SHIFT_SIZE_T] & mask_check;

		h.nele = (swt)?h.nele-1:h.nele;

		h.mask[sub_id >> BIT_SHIFT_SIZE_T] &= ~mask_check;

		if (h.nele == 0 && swt != 0)
		{
			// Add the chunks in the empty list
			empty_v.add(active_cnk);
		}
	}

incardon's avatar
incardon committed
514
515
516
517
518
public:

	//! it define that this data-structure is a grid
	typedef int yes_i_am_grid;

519
520
521
	//! base_key for the grid
	typedef grid_key_dx<dim> base_key;

incardon's avatar
incardon committed
522
523
524
	//! expose the dimansionality as a static const
	static constexpr unsigned int dims = dim;

525
526
527
528
529
530
531
	//! value that the grid store
	//! The object type the grid is storing
	typedef T value_type;

	//! sub-grid iterator type
	typedef grid_key_sparse_dx_iterator_sub<dim, chunking::size::value> sub_grid_iterator_type;

incardon's avatar
incardon committed
532
533
534
535
536
537
	/*! \brief Trivial constructor
	 *
	 *
	 */
	inline sgrid_cpu()
	:cache_pnt(0)
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
	{
		init();
	}

	/*! \brief create a sparse grid from another grid
	 *
	 * \param g the grid to copy
	 *
	 */
	inline sgrid_cpu(const sgrid_cpu & g) THROW
	{
		this->operator=(g);
	}

	/*! \brief create a sparse grid from another grid
	 *
	 * \param g the grid to copy
	 *
	 */
	inline sgrid_cpu(const sgrid_cpu && g) THROW
	{
		this->operator=(g);
	}
incardon's avatar
incardon committed
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578

	/*! \brief Constructor for sparse grid
	 *
	 * \param sz size in each direction of the sparse grid
	 *
	 */
	sgrid_cpu(const size_t (& sz)[dim])
	:cache_pnt(0),g_sm(sz)
	{
		// calculate the chunks grid

		set_g_shift_from_size(sz,g_sm_shift);

		// fill pos_g

		init();
	}

incardon's avatar
incardon committed
579
580
581
582
583
584
585
586
587
588
	/*! \brief Get the background value
	 *
	 * \return background value
	 *
	 */
	T & getBackgroundValue()
	{
		return background;
	}

incardon's avatar
incardon committed
589
590
591
592
593
	/*! \brief This is a multiresolution sparse grid so is a compressed format
	 *
	 * \return true
	 *
	 */
incardon's avatar
incardon committed
594
	static constexpr bool isCompressed()
incardon's avatar
incardon committed
595
596
597
598
	{
		return true;
	}

incardon's avatar
incardon committed
599
600
601
602
603
604
605
606
607
608
609
610
611
612
	/*! \brief Insert a full element (with all properties)
	 *
	 * \param v1 where you want to insert the element
	 *
	 */
	inline auto insert_o(const grid_key_dx<dim> & v1, size_t & ele_id) -> decltype(chunks.template get_o(0))
	{
		size_t active_cnk;

		pre_insert(v1,active_cnk,ele_id);

		return chunks.template get_o(active_cnk);
	}

incardon's avatar
incardon committed
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
	/*! \brief Get the reference of the selected element
	 *
	 * \param v1 grid_key that identify the element in the grid
	 *
	 * \return the reference of the element
	 *
	 */
	template <unsigned int p, typename r_type=decltype(chunks.template get<p>(0)[0])>
	inline r_type insert(const grid_key_dx<dim> & v1)
	{
		size_t active_cnk = 0;

		grid_key_dx<dim> kh = v1;
		grid_key_dx<dim> kl;

		// shift the key
		key_shift<dim,chunking>::shift(kh,kl);

		long int lin_id = g_sm_shift.LinId(kh);

		size_t id = 0;
		for (size_t k = 0 ; k < SGRID_CACHE; k++)
		{id += (cache[k] == lin_id)?k+1:0;}

		if (id == 0)
		{
			// we do not have it in cache we check if we have it in the map

			auto fnd = map.find(lin_id);
			if (fnd == map.end())
			{
				// we do not have it in the map create a chunk

				map[lin_id] = chunks.size();
				chunks.add();
				header.add();
				header.last().pos = kh;
				header.last().nele = 0;

				// set the mask to null
				auto & h = header.last().mask;

				for (size_t i = 0 ; i < chunking::size::value / (sizeof(size_t)*8) + (chunking::size::value % (sizeof(size_t)*8) != 0) + 1 ; i++)
				{h[i] = 0;}

				key_shift<dim,chunking>::cpos(header.last().pos);

				active_cnk = chunks.size() - 1;
			}
			else
			{
				// we have it in the map

				active_cnk = fnd->second;
			}

			// Add on cache the chunk
			cache[cache_pnt] = lin_id;
			cached_id[cache_pnt] = active_cnk;
			cache_pnt++;
			cache_pnt = (cache_pnt >= SGRID_CACHE)?0:cache_pnt;
		}
		else
		{
			active_cnk = cached_id[id-1];
incardon's avatar
incardon committed
678
679
			cache_pnt = id;
			cache_pnt = (cache_pnt == SGRID_CACHE)?0:cache_pnt;
incardon's avatar
incardon committed
680
681
682
683
684
685
686
687
688
689
690
		}

		size_t sub_id = sublin<dim,typename chunking::shift_c>::lin(kl);

		// the chunk is in cache, solve

		// we notify that we added one element
		auto & h = header.get(active_cnk);

		// we set the mask

691
692
693
694
		size_t mask_check = (size_t)1 << (sub_id & ((1 << BIT_SHIFT_SIZE_T) - 1));

		h.nele = (h.mask[sub_id >> BIT_SHIFT_SIZE_T] & mask_check)?h.nele:h.nele + 1;
		h.mask[sub_id >> BIT_SHIFT_SIZE_T] |= mask_check;
incardon's avatar
incardon committed
695
696
697
698

		return chunks.template get<p>(active_cnk)[sub_id];
	}

699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
	/*! \brief Get the reference of the selected element
	 *
	 * \param v1 grid_key that identify the element in the grid
	 *
	 * \return the reference of the element
	 *
	 */
	template <unsigned int p, typename r_type=decltype(chunks.template get<p>(0)[0])>
	inline r_type insert(const grid_key_sparse_lin_dx & v1)
	{
		size_t active_cnk = v1.getChunk();
		size_t sub_id = v1.getPos();

		// the chunk is in cache, solve

		// we notify that we added one element
		auto & h = header.get(active_cnk);

		// we set the mask

		size_t mask_check = (size_t)1 << (sub_id & ((1 << BIT_SHIFT_SIZE_T) - 1));

		h.nele = (h.mask[sub_id >> BIT_SHIFT_SIZE_T] & mask_check)?h.nele:h.nele + 1;
		h.mask[sub_id >> BIT_SHIFT_SIZE_T] |= mask_check;

		return chunks.template get<p>(active_cnk)[sub_id];
	}

incardon's avatar
incardon committed
727
728
729
730
731
732
733
	/*! \brief Get the reference of the selected element
	 *
	 * \param v1 grid_key that identify the element in the grid
	 *
	 * \return the reference of the element
	 *
	 */
734
735
	template <unsigned int p>
	inline auto get(const grid_key_dx<dim> & v1) const -> decltype(openfpm::as_const(chunks.template get<p>(0))[0])
incardon's avatar
incardon committed
736
	{
incardon's avatar
incardon committed
737
		size_t active_cnk;
incardon's avatar
incardon committed
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
		grid_key_dx<dim> kh = v1;
		grid_key_dx<dim> kl;

		// shift the key
		key_shift<dim,chunking>::shift(kh,kl);

		long int lin_id = g_sm_shift.LinId(kh);

		size_t id = 0;
		for (size_t k = 0 ; k < SGRID_CACHE; k++)
		{id += (cache[k] == lin_id)?k+1:0;}

		if (id == 0)
		{
			auto it = map.find(lin_id);

			if (it == map.end())
			{return background.template get<p>();}

757
			add_on_cache(lin_id,it->second);
incardon's avatar
incardon committed
758

incardon's avatar
incardon committed
759
760
761
762
763
			active_cnk = it->second;
		}
		else
		{
			active_cnk = cached_id[id-1];
incardon's avatar
incardon committed
764
765
766
767
		}

		size_t sub_id = sublin<dim,typename chunking::shift_c>::lin(kl);

incardon's avatar
incardon committed
768
		// we check the mask
incardon's avatar
incardon committed
769
		auto & h = header.get(active_cnk);
incardon's avatar
incardon committed
770

incardon's avatar
incardon committed
771
		// We check the mask
incardon's avatar
incardon committed
772
773
774
775
776
		size_t mask_check = (size_t)1 << (sub_id & ((1 << BIT_SHIFT_SIZE_T) - 1));

		if ((h.mask[sub_id >> BIT_SHIFT_SIZE_T] & mask_check) == 0)
		{return background.template get<p>();}

incardon's avatar
incardon committed
777
		return chunks.template get<p>(active_cnk)[sub_id];
incardon's avatar
incardon committed
778
779
780
781
782
783
784
785
786
	}

	/*! \brief Get the reference of the selected element
	 *
	 * \param v1 grid_key that identify the element in the grid
	 *
	 * \return the reference of the element
	 *
	 */
incardon's avatar
incardon committed
787
788
789
	template <unsigned int p>
	inline auto get(const grid_key_dx<dim> & v1) -> decltype(openfpm::as_const(chunks.template get<p>(0))[0])
	{
incardon's avatar
incardon committed
790
		size_t active_cnk;
incardon's avatar
incardon committed
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
		grid_key_dx<dim> kh = v1;
		grid_key_dx<dim> kl;

		// shift the key
		key_shift<dim,chunking>::shift(kh,kl);

		long int lin_id = g_sm_shift.LinId(kh);

		size_t id = 0;
		for (size_t k = 0 ; k < SGRID_CACHE; k++)
		{id += (cache[k] == lin_id)?k+1:0;}

		if (id == 0)
		{
			auto it = map.find(lin_id);

			if (it == map.end())
			{return background.template get<p>();}

810
			add_on_cache(lin_id,it->second);
incardon's avatar
incardon committed
811

incardon's avatar
incardon committed
812
813
814
815
816
			active_cnk = it->second;
		}
		else
		{
			active_cnk = cached_id[id-1];
incardon's avatar
incardon committed
817
818
819
820
		}

		size_t sub_id = sublin<dim,typename chunking::shift_c>::lin(kl);

incardon's avatar
incardon committed
821
		// we check the mask
incardon's avatar
incardon committed
822
		auto & h = header.get(active_cnk);
incardon's avatar
incardon committed
823
824
825
826
827
828
829

		// We check the mask
		size_t mask_check = (size_t)1 << (sub_id & ((1 << BIT_SHIFT_SIZE_T) - 1));

		if ((h.mask[sub_id >> BIT_SHIFT_SIZE_T] & mask_check) == 0)
		{return background.template get<p>();}

incardon's avatar
incardon committed
830
		return chunks.template get<p>(active_cnk)[sub_id];
incardon's avatar
incardon committed
831
832
	}

incardon's avatar
incardon committed
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
	/*! \brief Check if the point exist
	 *
	 * \param v1 grid_key that identify the element in the grid
	 *
	 * \return the true if the point exist
	 *
	 */
	inline bool existPoint(const grid_key_dx<dim> & v1) const
	{
		size_t active_cnk;
		grid_key_dx<dim> kh = v1;
		grid_key_dx<dim> kl;

		// shift the key
		key_shift<dim,chunking>::shift(kh,kl);

		long int lin_id = g_sm_shift.LinId(kh);

		size_t id = 0;
		for (size_t k = 0 ; k < SGRID_CACHE; k++)
		{id += (cache[k] == lin_id)?k+1:0;}

		if (id == 0)
		{
			auto it = map.find(lin_id);

			if (it == map.end())
			{return false;}

			add_on_cache(lin_id,it->second);

			active_cnk = it->second;
		}
		else
		{
			active_cnk = cached_id[id-1];
		}

		size_t sub_id = sublin<dim,typename chunking::shift_c>::lin(kl);

		// we check the mask
		auto & h = header.get(active_cnk);

		// We check the mask
		size_t mask_check = (size_t)1 << (sub_id & ((1 << BIT_SHIFT_SIZE_T) - 1));

		if ((h.mask[sub_id >> BIT_SHIFT_SIZE_T] & mask_check) == 0)
		{return false;}

		return true;
	}

incardon's avatar
incardon committed
885
886
887
888
889
890
891
892
893
	/*! \brief Get the reference of the selected element
	 *
	 * \param v1 grid_key that identify the element in the grid
	 *
	 * \return the reference of the element
	 *
	 */
	template <unsigned int p>
	inline auto get(const grid_key_sparse_lin_dx & v1) -> decltype(chunks.template get<p>(0)[0])
incardon's avatar
incardon committed
894
895
896
897
898
899
900
901
902
	{
		return chunks.template get<p>(v1.getChunk())[v1.getPos()];
	}

	/*! \brief Return a Domain iterator
	 *
	 * \return return the domain iterator
	 *
	 */
903
	grid_key_sparse_dx_iterator<dim,chunking::size::value>
incardon's avatar
incardon committed
904
	getIterator(size_t opt = 0) const
incardon's avatar
incardon committed
905
	{
incardon's avatar
incardon committed
906
		return grid_key_sparse_dx_iterator<dim,chunking::size::value>(&header,&pos_chunk);
incardon's avatar
incardon committed
907
	}
908

909
910
911
912
913
914
	/*! \brief Return an iterator over a sub-grid
	 *
	 * \return return an iterator over a sub-grid
	 *
	 */
	grid_key_sparse_dx_iterator_sub<dim,chunking::size::value>
incardon's avatar
incardon committed
915
	getIterator(const grid_key_dx<dim> & start, const grid_key_dx<dim> & stop, size_t opt = 0) const
916
917
918
919
	{
		return grid_key_sparse_dx_iterator_sub<dim,chunking::size::value>(header,pos_chunk,start,stop,sz_cnk);
	}

incardon's avatar
incardon committed
920
921
922
923
924
925
926
	/*! \brief Return the internal grid information
	 *
	 * Return the internal grid information
	 *
	 * \return the internal grid
	 *
	 */
927
	const grid_sm<dim,void> & getGrid() const
incardon's avatar
incardon committed
928
929
930
931
932
933
934
	{
#ifdef SE_CLASS2
		check_valid(this,8);
#endif
		return g_sm;
	}

935
936
937
938
939
940
941
	/*! \brief Remove the point
	 *
	 * \param v1 element to remove
	 *
	 */
	void remove(const grid_key_dx<dim> & v1)
	{
incardon's avatar
incardon committed
942
943
944
		remove_point(v1);
		remove_empty();
	}
945

incardon's avatar
incardon committed
946
947
948
949
950
951
952
953
954
	/*! \brief Remove the point but do not flush the remove
	 *
	 * \param v1 element to remove
	 *
	 */
	void remove_no_flush(const grid_key_dx<dim> & v1)
	{
		remove_point(v1);
	}
955

incardon's avatar
incardon committed
956
957
958
959
960
961
962
963
	/*! \brief Remove the point
	 *
	 * \param v1 element to remove
	 *
	 */
	void flush_remove()
	{
		remove_empty();
964
965
	}

incardon's avatar
incardon committed
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
	/*! \brief Resize the grid
	 *
	 * The old information is retained on the new grid if the new grid is bigger.
	 * When smaller the data are cropped
	 *
	 * \param sz reference to an array of dimension dim
	 *
	 */
	void resize(const size_t (& sz)[dim])
	{
		bool is_bigger = true;

		// we check if we are resizing bigger, because if is the case we do not have to do
		// much

		for (size_t i = 0 ; i < dim ; i++)
		{
			if (sz[i] < g_sm.size(i))
			{is_bigger = false;}
		}

		g_sm.setDimensions(sz);

		// set g_sm_shift

		set_g_shift_from_size(sz,g_sm_shift);

		clear_cache();

		if (is_bigger == true)
		{

			// if we resize bigger we do not have to do anything in the headers
			// and in the chunks we just have to update g_sm and reset the cache
			// and reconstruct the map. So we reconstruct the map and we just
			// finish

			reconstruct_map();

			return;
		}

		// create a box that is as big as the grid

		Box<dim,size_t> gs_box;

		for (size_t i = 0 ; i < dim ; i++)
		{
			gs_box.setLow(i,0);
			gs_box.setHigh(i,g_sm.size(i));
		}

		// we take a list of all chunks to remove
		openfpm::vector<size_t> rmh;

		// in this case we have to crop data, we go through all the headers

		for (size_t i = 0 ; i < header.size() ; i++)
		{
			Box<dim,size_t> cnk;

			for (size_t j = 0 ; j < dim ; j++)
			{
				cnk.setLow(j,header.get(i).pos.get(j));
				cnk.setHigh(j,sz_cnk[j] + header.get(i).pos.get(j));
			}

			// if the chunk is not fully contained in the new smaller sparse grid
			// we have to crop it
			if (!cnk.isContained(gs_box))
			{
				// We check if the chunks is fully out or only partially in
				// cheking the intersection between the new grid and the box
				// enclosing the chunk as it was before
				Box<dim,size_t> inte;

				if (gs_box.Intersect(cnk,inte))
				{
					// part of the chunk is in, part is out

					// shift P1 to the origin
					// this is the valid box everything out must me reset
					inte -= inte.getP1();

					size_t mask_nele;
					short unsigned int mask_it[chunking::size::value];

					auto & mask = header.get(i).mask;
					auto & n_ele = header.get(i).nele;

					// ok so the box is not fully contained so we must crop data

					fill_mask(mask_it,mask,mask_nele);

					// now we have the mask of all the filled elements

					for (size_t j = 0 ; j < mask_nele ; j++)
					{
						if (!inte.isInside(pos_chunk[mask_it[j]].toPoint()))
						{
							// if is not inside, the point must be deleted

							remove_from_chunk<chunking::size::value>(mask_it[j],n_ele,mask);
						}
					}
				}
				else
				{
					// the chunk is completely out and must be removed completely
					// we add it to the list of the chunks to remove

					rmh.add(i);
				}
			}
		}

		header.remove(rmh,0);
		chunks.remove(rmh,0);

		reconstruct_map();
	}
1087

1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
	/*! \brief Calculate the memory size required to pack n elements
	 *
	 * Calculate the total size required to store n-elements in a vector
	 *
	 * \param n number of elements
	 * \param e unused
	 *
	 * \return the size of the allocation number e
	 *
	 */
	template<int ... prp> static inline size_t packMem(size_t n, size_t e)
	{
		if (sizeof...(prp) == 0)
		{return n * sizeof(typename T::type);}

		typedef object<typename object_creator<typename T::type,prp...>::type> prp_object;

		return n * sizeof(prp_object);
	}

	/*! \brief Insert an allocation request
	 *
	 * \tparam prp set of properties to pack
	 *

	 * \param sub sub-grid iterator
	 * \param vector of requests
	 *
	 */
	template<int ... prp> inline
incardon's avatar
incardon committed
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
	void packRequest(size_t & req) const
	{
		grid_sm<dim,void> gs_cnk(sz_cnk);

		// For sure we have to pack the number of chunk we want to pack

		req += sizeof(size_t);
		req += dim*sizeof(size_t);

		// Here we have to calculate the number of points to pack

		for (size_t i = 0 ; i < header.size() ; i++)
		{
			auto & h = header.get(i);

			size_t mask_nele;
			short unsigned int mask_it[chunking::size::value];

			fill_mask(mask_it,h.mask,mask_nele);

			for (size_t j = 0 ; j < mask_nele ; j++)
			{
				// If all of the aggregate properties do not have a "pack()" member
				if (has_pack_agg<T,prp...>::result::value == false)
				{
					// here we count how many chunks must be sent

					size_t alloc_ele = this->packMem<prp...>(1,0);
					req += alloc_ele;
				}
				else
				{
					//Call a pack request
					call_aggregatePackRequestChunking<decltype(chunks.template get_o(i)),
																  S,prp ... >
																  ::call_packRequest(chunks.template get_o(i),mask_it[j],req);
				}
			}

			// There are point to send. So we have to save the mask chunk
			req += sizeof(header.get(i).mask);
			// the chunk position
			req += sizeof(header.get(i).pos);
			// and the number of element
			req += sizeof(header.get(i).nele);
		}
	}

	/*! \brief Insert an allocation request
	 *
	 * \tparam prp set of properties to pack
	 *

	 * \param sub sub-grid to pack
	 * \param vector of requests
	 *
	 */
	template<int ... prp> inline
1176
1177
1178
1179
1180
1181
1182
1183
	void packRequest(grid_key_sparse_dx_iterator_sub<dim,chunking::size::value> & sub_it,
					 size_t & req) const
	{
		grid_sm<dim,void> gs_cnk(sz_cnk);

		// For sure we have to pack the number of chunk we want to pack

		req += sizeof(size_t);
incardon's avatar
incardon committed
1184
		req += dim*sizeof(size_t);
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289

		// Here we have to calculate the number of points to pack

		Box<dim,size_t> section_to_pack;

		for (size_t i = 0; i < dim ; i++)
		{
			section_to_pack.setLow(i,sub_it.getStart().get(i));
			section_to_pack.setHigh(i,sub_it.getStop().get(i));
		}

		for (size_t i = 0 ; i < header.size() ; i++)
		{
			auto & h = header.get(i);

			Box<dim,size_t> bc;

			for (size_t j = 0 ; j < dim ; j++)
			{
				bc.setLow(j,header.get(i).pos.get(j));
				bc.setHigh(j,header.get(i).pos.get(j) + sz_cnk[j] - 1);
			}

			// now we intersect the chunk box with the box

			Box<dim,size_t> inte;
			bool stp = bc.Intersect(section_to_pack,inte);

			if (stp == true)
			{
				// If it is intersect ok we have to check if there are points to pack
				// we shift inte to be relative to the chunk origin

				inte -= header.get(i).pos.toPoint();

				// we iterate all the points

				size_t old_req = req;
				grid_key_dx_iterator_sub<dim> sit(gs_cnk,inte.getKP1(),inte.getKP2());

				while (sit.isNext())
				{
					auto key = sit.get();

					size_t sub_id = gs_cnk.LinId(key);

					size_t mask_check = (size_t)1 << (sub_id & ((1 << BIT_SHIFT_SIZE_T) - 1));

					if (h.mask[sub_id >> BIT_SHIFT_SIZE_T] & mask_check)
					{
						// If all of the aggregate properties do not have a "pack()" member
						if (has_pack_agg<T,prp...>::result::value == false)
						{
							// here we count how many chunks must be sent

							size_t alloc_ele = this->packMem<prp...>(1,0);
							req += alloc_ele;
						}
						//If at least one property has "pack()"
						else
						{
							//Call a pack request
							call_aggregatePackRequestChunking<decltype(chunks.template get_o(i)),
																  S,prp ... >
																  ::call_packRequest(chunks.template get_o(i),sub_id,req);
						}
					}

					++sit;
				}

				if (old_req != req)
				{
					// There are point to send. So we have to save the mask chunk
					req += sizeof(header.get(i).mask);
					// the chunk position
					req += sizeof(header.get(i).pos);
					// and the number of element
					req += sizeof(header.get(i).nele);
				}
			}
		}
	}

	/*! \brief Pack the object into the memory given an iterator
	 *
	 * \tparam prp properties to pack
	 *
	 * \param mem preallocated memory where to pack the objects
	 * \param sub_it sub grid iterator ( or the elements in the grid to pack )
	 * \param sts pack statistic
	 *
	 */
	template<int ... prp> void pack(ExtPreAlloc<S> & mem,
									grid_key_sparse_dx_iterator_sub<dims,chunking::size::value> & sub_it,
									Pack_stat & sts)
	{
		grid_sm<dim,void> gs_cnk(sz_cnk);

		// Here we allocate a size_t that indicate the number of chunk we are packing,
		// because we do not know a priory, we will fill it later

		mem.allocate(sizeof(size_t));
		size_t * number_of_chunks = (size_t *)mem.getPointer();

incardon's avatar
incardon committed
1290
1291
1292
1293
1294
		// Pack the size of the grid

		for (size_t i = 0 ; i < dim ; i++)
		{Packer<size_t,S>::pack(mem,getGrid().size(i),sts);}

1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
		// Here we have to calculate the number of points to pack

		Box<dim,size_t> section_to_pack;

		for (size_t i = 0; i < dim ; i++)
		{
			section_to_pack.setLow(i,sub_it.getStart().get(i));
			section_to_pack.setHigh(i,sub_it.getStop().get(i));
		}

		size_t n_packed_chunk = 0;

		for (size_t i = 0 ; i < header.size() ; i++)
		{
			auto & h = header.get(i);

			Box<dim,size_t> bc;

			for (size_t j = 0 ; j < dim ; j++)
			{
				bc.setLow(j,header.get(i).pos.get(j));
				bc.setHigh(j,header.get(i).pos.get(j) + sz_cnk[j] - 1);
			}

			// now we intersect the chunk box with the box

			Box<dim,size_t> inte;
			bool stp = bc.Intersect(section_to_pack,inte);

			if (stp == true)
			{
				// This flag indicate if something has been packed from this chunk
				bool has_packed = false;

				size_t mask_to_pack[chunking::size::value / (sizeof(size_t)*8) + (chunking::size::value % (sizeof(size_t)*8) != 0) + 1];
				memset(mask_to_pack,0,sizeof(mask_to_pack));
incardon's avatar
incardon committed
1331
				mem.allocate_nocheck(sizeof(header.get(i).mask) + sizeof(header.get(i).pos) + sizeof(header.get(i).nele));
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358

				// here we get the pointer of the memory in case we have to pack the header
				// and we also shift the memory pointer by an offset equal to the header
				// to pack
				unsigned char * ptr_start = (unsigned char *)mem.getPointer();

				// If it is intersect ok we have to check if there are points to pack
				// we shift inte intp the chunk origin

				inte -= header.get(i).pos.toPoint();

				// we iterate all the points

				grid_key_dx_iterator_sub<dim> sit(gs_cnk,inte.getKP1(),inte.getKP2());

				while (sit.isNext())
				{
					auto key = sit.get();

					size_t sub_id = gs_cnk.LinId(key);

					size_t mask_check = (size_t)1 << (sub_id & ((1 << BIT_SHIFT_SIZE_T) - 1));

					if (h.mask[sub_id >> BIT_SHIFT_SIZE_T] & mask_check)
					{
						Packer<decltype(chunks.template get_o(i)),
									S,
incardon's avatar
incardon committed
1359
									PACKER_ENCAP_OBJECTS_CHUNKING>::template pack<T,prp...>(mem,chunks.template get_o(i),sub_id,sts);
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404

						mask_to_pack[sub_id >> BIT_SHIFT_SIZE_T] |= mask_check;
						has_packed = true;

					}

					++sit;
				}

				if (has_packed == true)
				{
					unsigned char * ptr_final = (unsigned char *)mem.getPointer();
					unsigned char * ptr_final_for = (unsigned char *)mem.getPointerEnd();

					// Ok we packed something so we have to pack the header
					 size_t shift = ptr_final - ptr_start;

					 mem.shift_backward(shift);

					 // The position of the chunks

					 grid_key_dx<dim> pos = header.get(i).pos - sub_it.getStart();

					 Packer<decltype(header.get(i).mask),S>::pack(mem,mask_to_pack,sts);
					 Packer<decltype(header.get(i).pos),S>::pack(mem,pos,sts);
					 Packer<decltype(header.get(i).nele),S>::pack(mem,header.get(i).nele,sts);

					 size_t shift_for = ptr_final_for - (unsigned char *)mem.getPointer();

					 mem.shift_forward(shift_for);

					 n_packed_chunk++;
				}
				else
				{
					// This just reset the last allocation
					mem.shift_backward(0);
				}
			}
		}

		// Now we fill the number of packed chunks
		*number_of_chunks = n_packed_chunk;
	}

incardon's avatar
incardon committed
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
	/*! \brief Pack the object into the memory given an iterator
	 *
	 * \tparam prp properties to pack
	 *
	 * \param mem preallocated memory where to pack the objects
	 * \param sub_it sub grid iterator ( or the elements in the grid to pack )
	 * \param sts pack statistic
	 *
	 */
	template<int ... prp> void pack(ExtPreAlloc<S> & mem,
									Pack_stat & sts)
	{
		grid_sm<dim,void> gs_cnk(sz_cnk);

		// Here we allocate a size_t that indicate the number of chunk we are packing,
		// because we do not know a priory, we will fill it later

		Packer<size_t,S>::pack(mem,header.size(),sts);

		for (size_t i = 0 ; i < dim ; i++)
		{Packer<size_t,S>::pack(mem,getGrid().size(i),sts);}

		// Here we have to calculate the number of points to pack

		for (size_t i = 0 ; i < header.size() ; i++)
		{
			auto & h = header.get(i);

			Packer<decltype(header.get(i).mask),S>::pack(mem,h.mask,sts);
			Packer<decltype(header.get(i).pos),S>::pack(mem,h.pos,sts);
			Packer<decltype(header.get(i).nele),S>::pack(mem,h.nele,sts);

			// we iterate all the points

			size_t mask_nele;
			short unsigned int mask_it[chunking::size::value];

			fill_mask(mask_it,h.mask,mask_nele);

			for (size_t j = 0 ; j < mask_nele ; j++)
			{
				Packer<decltype(chunks.template get_o(i)),
								S,
incardon's avatar
incardon committed
1448
								PACKER_ENCAP_OBJECTS_CHUNKING>::template pack<T,prp...>(mem,chunks.template get_o(i),mask_it[j],sts);
incardon's avatar
incardon committed
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
			};
		}
	}

	/*! \brief It does materially nothing
	 *
	 */
	void setMemory()
	{}

1459

1460
1461
1462
1463
1464
1465
1466
	/*! \number of element inserted
	 *
	 * \warning this function is not as fast as the size in other structures
	 *
	 * \return the total number of elements inserted
	 *
	 */
1467
	size_t size() const
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
	{
		size_t tot = 0;

		for (size_t i = 0 ; i < header.size() ; i++)
		{
			tot += header.get(i).nele;
		}

		return tot;
	}

incardon's avatar
incardon committed
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
	/*! \brief Remove all the points in this region
	 *
	 * \param box_src box to kill the points
	 *
	 */
	void remove(Box<dim,size_t> & section_to_delete)
	{
		grid_sm<dim,void> gs_cnk(sz_cnk);

		for (size_t i = 0 ; i < header.size() ; i++)
		{
			auto & h = header.get(i);

			if (i == 9672)
			{
				int debug = 0;
				debug++;
			}

			Box<dim,size_t> bc;

			for (size_t j = 0 ; j < dim ; j++)
			{
				bc.setLow(j,header.get(i).pos.get(j));
				bc.setHigh(j,header.get(i).pos.get(j) + sz_cnk[j] - 1);
			}

			// now we intersect the chunk box with the box

			Box<dim,size_t> inte;
			bool stp = bc.Intersect(section_to_delete,inte);

			if (stp == true)
			{
				// If it is intersect ok we have to check if there are points to pack
				// we shift inte intp the chunk origin

				inte -= header.get(i).pos.toPoint();

				// we iterate all the points

				grid_key_dx_iterator_sub<dim> sit(gs_cnk,inte.getKP1(),inte.getKP2());

				while (sit.isNext())
				{
					auto key = sit.get();

					size_t sub_id = gs_cnk.LinId(key);
					size_t mask_check = (size_t)1 << (sub_id & ((1 << BIT_SHIFT_SIZE_T) - 1));

					size_t swt = header.get(i).mask[sub_id >> BIT_SHIFT_SIZE_T] & mask_check;

					h.nele = (swt)?h.nele-1:h.nele;
					h.mask[sub_id >> BIT_SHIFT_SIZE_T] &= ~mask_check;

					if (h.nele == 0 && swt != 0)
					{
						// Add the chunks in the empty list
						empty_v.add(i);
					}

					++sit;
				}
			}
		}

		remove_empty();
	}

1548
	void copy_to(const sgrid_cpu<dim,T,S,typename memory_traits_lin<T>::type,chunking> & grid_src,
incardon's avatar
incardon committed
1549
1550
1551
		         const Box<dim,size_t> & box_src,
			     const Box<dim,size_t> & box_dst)
	{
1552
		auto it = grid_src.getIterator(box_src.getKP1(),box_src.getKP2());
incardon's avatar
incardon committed
1553

1554
1555
1556
1557
		while (it.isNext())
		{
			auto key_src = it.get();
			grid_key_dx<dim> key_dst = key_src + box_dst.getKP1();
incardon's avatar
incardon committed
1558
			key_dst -= box_src.getKP1();
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568

			typedef typename std::remove_const<typename std::remove_reference<decltype(grid_src)>::type>::type gcopy;

			copy_sparse_to_sparse<dim,gcopy,gcopy> caps(grid_src,*this,key_src,key_dst);
			boost::mpl::for_each_ref< boost::mpl::range_c<int,0,T::max_prop> >(caps);

			++it;
		}
	}

incardon's avatar
incardon committed
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
	template<template <typename,typename> class op, unsigned int ... prp >
	void copy_to_op(const sgrid_cpu<dim,T,S,typename memory_traits_lin<T>::type,chunking> & grid_src,
		         const Box<dim,size_t> & box_src,
			     const Box<dim,size_t> & box_dst)
	{
		auto it = grid_src.getIterator(box_src.getKP1(),box_src.getKP2());

		while (it.isNext())
		{
			auto key_src = it.get();
			grid_key_dx<dim> key_dst = key_src + box_dst.getKP1();
			key_dst -= box_src.getKP1();

			typedef typename std::remove_const<typename std::remove_reference<decltype(grid_src)>::type>::type gcopy;

			copy_sparse_to_sparse_op<op,dim,gcopy,gcopy,prp ...> caps(grid_src,*this,key_src,key_dst);
			boost::mpl::for_each_ref< boost::mpl::range_c<int,0,sizeof...(prp)> >(caps);

			++it;
		}
	}
incardon's avatar
incardon committed
1590

1591
1592
1593
1594
1595
1596
1597
1598
1599
	/*! \brief unpack the sub-grid object
	 *
	 * \tparam prp properties to unpack
	 *
	 * \param mem preallocated memory from where to unpack the object
	 * \param sub sub-grid iterator
	 * \param obj object where to unpack
	 *
	 */
incardon's avatar
incardon committed
1600
1601
	template<unsigned int ... prp, typename S2>
	void unpack(ExtPreAlloc<S2> & mem,
1602
1603
1604
1605
				grid_key_sparse_dx_iterator_sub<dims,chunking::size::value> & sub_it,
				Unpack_stat & ps)
	{
		short unsigned int mask_it[chunking::size::value];
incardon's avatar
incardon committed
1606

1607
		// first we unpack the number of chunks
incardon's avatar
incardon committed
1608

1609
1610
		size_t n_chunks;

incardon's avatar
incardon committed
1611
		Unpacker<size_t,S2>::unpack(mem,n_chunks,ps);
1612

incardon's avatar
incardon committed
1613
1614
1615
1616
		size_t sz[dim];
		for (size_t i = 0 ; i < dim ; i++)
		{Unpacker<size_t,S2>::unpack(mem,sz[i],ps);}

incardon's avatar
incardon committed
1617
1618
1619
1620
1621
		openfpm::vector<cheader<dim,chunking::size::value>> header_tmp;
		openfpm::vector<aggregate_bfv<chunk_def>> chunks_tmp;

		header_tmp.resize(n_chunks);
		chunks_tmp.resize(n_chunks);
1622
1623

		for (size_t i = 0 ; i < n_chunks ; i++)
incardon's avatar
incardon committed
1624
		{
incardon's avatar
incardon committed
1625
			auto & h = header_tmp.get(i);
incardon's avatar
incardon committed
1626

incardon's avatar
incardon committed
1627
1628
1629
			Unpacker<decltype(header.get(i).mask),S2>::unpack(mem,h.mask,ps);
			Unpacker<decltype(header.get(i).pos),S2>::unpack(mem,h.pos,ps);
			Unpacker<decltype(header.get(i).nele),S2>::unpack(mem,h.nele,ps);
1630
1631
1632
1633
1634
1635

			// fill the mask_it

			fill_mask(mask_it,h.mask,h.nele);

			// now we unpack the information
incardon's avatar
incardon committed
1636
1637
			size_t active_cnk;
			size_t ele_id;
1638
1639

			for (size_t k = 0 ; k < h.nele ; k++)
incardon's avatar
incardon committed
1640
			{
incardon's avatar
incardon committed
1641
1642
1643
1644
1645
1646
1647
				// construct v1
				grid_key_dx<dim> v1;
				for (size_t i = 0 ; i < dim ; i++)
				{v1.set_d(i,h.pos.get(i) + pos_chunk[mask_it[k]].get(i) + sub_it.getStart().get(i));}

				pre_insert(v1,active_cnk,ele_id);

1648
				Unpacker<decltype(chunks.template get_o(mask_it[k])),
incardon's avatar
incardon committed
1649
							S2,
incardon's avatar
incardon committed
1650
							PACKER_ENCAP_OBJECTS_CHUNKING>::template unpack<T,prp...>(mem,chunks.template get_o(active_cnk),ele_id,ps);
1651

incardon's avatar
incardon committed
1652
			}
1653
		}
incardon's avatar
incardon committed
1654
	}
incardon's avatar
incardon committed
1655

incardon's avatar
incardon committed
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
	/*! \brief unpack the sub-grid object
	 *
	 * \tparam prp properties to unpack
	 *
	 * \param mem preallocated memory from where to unpack the object
	 * \param sub sub-grid iterator
	 * \param obj object where to unpack
	 *
	 */
	template<unsigned int ... prp, typename S2>
	void unpack(ExtPreAlloc<S2> & mem,
				Unpack_stat & ps)
	{
		this->clear();

		grid_key_dx<dim> start;
		grid_key_dx<dim> stop;

		// We preunpack some data
1675
		Unpack_stat ps_tmp = ps;
incardon's avatar
incardon committed
1676
1677
1678
1679
1680
1681
1682
1683

		size_t unused;
		Unpacker<size_t,S2>::unpack(mem,unused,ps_tmp);

		size_t sz[dim];
		for (size_t i = 0 ; i < dim ; i++)
		{Unpacker<size_t,S2>::unpack(mem,sz[i],ps_tmp);}

1684
		g_sm.setDimensions(sz);
incardon's avatar
incardon committed
1685
1686
1687
1688
1689
1690
1691
1692
		for (size_t i = 0 ; i < dim ; i++)
		{
			start.set_d(i,0);
			stop.set_d(i,getGrid().size(i)-1);
		}

		set_g_shift_from_size(sz,g_sm_shift);

1693
		auto sub_it = this->getIterator(start,stop);
incardon's avatar
incardon committed
1694

incardon's avatar
incardon committed
1695
		unpack<prp...>(mem,sub_it,ps);
incardon's avatar
incardon committed
1696
1697
	}

incardon's avatar
incardon committed
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
	/*! \brief unpack the sub-grid object applying an operation
	 *
	 * \tparam op operation
	 * \tparam prp properties to unpack
	 *
	 * \param mem preallocated memory from where to unpack the object
	 * \param sub sub-grid iterator
	 * \param obj object where to unpack
	 *
	 */
	template<template<typename,typename> class op, typename S2, unsigned int ... prp>
	void unpack_with_op(ExtPreAlloc<S2> & mem,
						grid_key_sparse_dx_iterator_sub<dim,chunking::size::value> & sub2,
						Unpack_stat & ps)
	{
		short unsigned int mask_it[chunking::size::value];

		// first we unpack the number of chunks

		size_t n_chunks;

		Unpacker<size_t,S2>::unpack(mem,n_chunks,ps);

incardon's avatar
incardon committed
1721
1722
1723
1724
		size_t sz[dim];
		for (size_t i = 0 ; i < dim ; i++)
		{Unpacker<size_t,S2>::unpack(mem,sz[i],ps);}

incardon's avatar
incardon committed
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
		openfpm::vector<cheader<dim,chunking::size::value>> header_tmp;
		openfpm::vector<aggregate_bfv<chunk_def>> chunks_tmp;

		header_tmp.resize(n_chunks);
		chunks_tmp.resize(n_chunks);

		for (size_t i = 0 ; i < n_chunks ; i++)
		{
			auto & h = header_tmp.get(i);

			Unpacker<decltype(header.get(i).mask),S2>::unpack(mem,h.mask,ps);
			Unpacker<decltype(header.get(i).pos),S2>::unpack(mem,h.pos,ps);
			Unpacker<decltype(header.get(i).nele),S2>::unpack(mem,h.nele,ps);

			// fill the mask_it

			fill_mask(mask_it,h.mask,h.nele);

			// now we unpack the information
			size_t active_cnk;
			size_t ele_id;

			for (size_t k = 0 ; k < h.nele ; k++)
			{
				// construct v1
				grid_key_dx<dim> v1;
				for (size_t i = 0 ; i < dim ; i++)
				{v1.set_d(i,h.pos.get(i) + pos_chunk[mask_it[k]].get(i) + sub2.getStart().get(i));}

incardon's avatar
incardon committed
1754
				bool exist = pre_insert(v1,active_cnk,ele_id);
incardon's avatar
incardon committed
1755

incardon's avatar
incardon committed
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
				if (exist == false)
				{
					Unpacker<decltype(chunks.template get_o(mask_it[k])),
								S2,
								PACKER_ENCAP_OBJECTS_CHUNKING>::template unpack_op<replace_,prp...>(mem,chunks.template get_o(active_cnk),ele_id,ps);
				}
				else
				{
					Unpacker<decltype(chunks.template get_o(mask_it[k])),
								S2,
								PACKER_ENCAP_OBJECTS_CHUNKING>::template unpack_op<op,prp...>(mem,chunks.template get_o(active_cnk),ele_id,ps);
				}
incardon's avatar
incardon committed
1768
1769
1770

			}
		}
1771
	}
incardon's avatar
incardon committed
1772

1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
	/*! \brief This is a meta-function return which type of sub iterator a grid produce
	 *
	 * \return the type of the sub-grid iterator
	 *
	 */
	template <typename stencil = no_stencil>
	static grid_key_sparse_dx_iterator_sub<dim, chunking::size::value>
	type_of_subiterator()
	{
		return  grid_key_sparse_dx_iterator_sub<dim, chunking::size::value>();
	}
incardon's avatar
incardon committed
1784

incardon's avatar
incardon committed
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
	/*! \brief This is a meta-function return which type of sub iterator a grid produce
	 *
	 * \return the type of the sub-grid iterator
	 *
	 */
	static grid_key_sparse_dx_iterator<dim, chunking::size::value>
	type_of_iterator()
	{
		return  grid_key_sparse_dx_iterator<dim, chunking::size::value>();
	}

1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
	/*! \brief Here we convert the linearized sparse key into the grid_key_dx
	 *
	 * \param key_out output key
	 * \param key_in input key
	 *
	 */
	void convert_key(grid_key_dx<dim> & key_out, const grid_key_sparse_lin_dx & key_in) const
	{
		auto & ph = header.get(key_in.getChunk()).pos;
		auto & pos_h = pos_chunk[key_in.getPos()];
incardon's avatar
incardon committed
1806

1807
1808
1809
1810
1811
		for (size_t i = 0 ; i < dim ; i++)
		{
			key_out.set_d(i,ph.get(i) + pos_h.get(i));
		}
	}
incardon's avatar
incardon committed
1812

1813
1814
1815
1816
1817
1818
1819
1820
1821
	/*! \brief copy an sparse grid
	 *
	 * \param tmp sparse grdi to copy
	 *
	 */
	sgrid_cpu & operator=(const sgrid_cpu & sg)
	{
		cache_pnt = sg.cache_pnt;
		meta_copy<T>::meta_copy_(sg.background,background);
incardon's avatar
incardon committed
1822

1823
1824
1825
1826
1827
		for (size_t i = 0 ; i < SGRID_CACHE ; i++)
		{
			cache[i] = sg.cache[i];
			cached_id[i] = sg.cached_id[i];
		}
incardon's avatar
incardon committed
1828

1829
1830
1831
1832
1833
1834
		//! Map to convert from grid coordinates to chunk
		map = sg.map;
		header = sg.header;
		chunks = sg.chunks;
		g_sm = sg.g_sm;
		g_sm_shift = sg.g_sm_shift;
incardon's avatar
incardon committed
1835

1836
1837
1838
1839
		for (size_t i = 0 ; i < chunking::size::value ; i++)
		{
			pos_chunk[i] = sg.pos_chunk[i];
		}
incardon's avatar
incardon committed
1840
1841


1842
1843
		for (size_t i = 0 ; i < dim ; i++)
		{sz_cnk[i] = sg.sz_cnk[i];}
incardon's avatar
incardon committed
1844

1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
		empty_v = sg.empty_v;

		return *this;
	}

	/*! \brief copy an sparse grid
	 *
	 * \param tmp sparse grdi to copy
	 *
	 */
	sgrid_cpu & operator=(sgrid_cpu && sg)
	{
		cache_pnt = sg.cache_pnt;
		meta_copy<T>::meta_copy_(sg.background,background);

		for (size_t i = 0 ; i < SGRID_CACHE ; i++)
		{
			cache[i] = sg.cache[i];
			cached_id[i] = sg.cached_id[i];
incardon's avatar
incardon committed
1864
		}
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884

		//! Map to convert from grid coordinates to chunk
		map.swap(sg.map);
		header.swap(sg.header);
		chunks.swap(sg.chunks);
		g_sm = sg.g_sm;
		g_sm_shift = sg.g_sm_shift;

		for (size_t i = 0 ; i < chunking::size::value ; i++)
		{
			pos_chunk[i] = sg.pos_chunk[i];
		}


		for (size_t i = 0 ; i < dim ; i++)
		{sz_cnk[i] = sg.sz_cnk[i];}

		empty_v = sg.empty_v;

		return *this;
incardon's avatar
incardon committed
1885
	}
1886

incardon's avatar
incardon committed
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
	/*! \brief Get the number of inserted points
	 *
	 * \return the number of inserted points
	 *
	 */
	size_t size_inserted()
	{
		return size();
	}

	/*! \brief delete all the points
	 *
	 *
	 */
	void clear()
	{
		header.clear();
		chunks.clear();

		clear_cache();
incardon's avatar
incardon committed
1907
		reconstruct_map();
incardon's avatar
incardon committed
1908
1909
	}

incardon's avatar
incardon committed
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
	/*! \brief This is an internal function to clear the cache
	 *
	 *
	 *
	 */
	void internal_clear_cache()
	{
		clear_cache();
	}

1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
	/*! \brief write the sparse grid into VTK
	 *
	 * \param out VTK output
	 *
	 */
	template<typename Tw = float> bool write(const std::string & output)
	{
		file_type ft = file_type::BINARY;

		openfpm::vector<Point<dim,Tw>> tmp_pos;
		openfpm::vector<T> tmp_prp;

		// copy position and properties

1934
		auto it = getIterator();