OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
DCPSE_op_Surface_tests.cpp
1//
2// Created by Abhinav Singh on 15.11.21.
3//
4
5#include "config.h"
6#ifdef HAVE_EIGEN
7#ifdef HAVE_PETSC
8
9
10#define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
11#define BOOST_MPL_LIMIT_VECTOR_SIZE 40
12
13#define BOOST_TEST_DYN_LINK
14
15#include "util/util_debug.hpp"
16#include <boost/test/unit_test.hpp>
17#include <iostream>
18#include "../DCPSE_surface_op.hpp"
19#include "../DCPSE_Solver.hpp"
20#include "Operators/Vector/vector_dist_operators.hpp"
21#include "Vector/vector_dist_subset.hpp"
22#include <iostream>
23#include "util/SphericalHarmonics.hpp"
24
25BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
26 BOOST_AUTO_TEST_CASE(dcpse_surface_simple) {
27 double boxP1{-1.5}, boxP2{1.5};
28 double boxSize{boxP2 - boxP1};
29 size_t n=256;
30 size_t sz[2] = {n,n};
31 double grid_spacing{boxSize/(sz[0]-1)};
32 double rCut{3.9 * grid_spacing};
33
34 Box<2,double> domain{{boxP1,boxP1},{boxP2,boxP2}};
35 size_t bc[2] = {NON_PERIODIC,NON_PERIODIC};
36 Ghost<2,double> ghost{rCut + grid_spacing/8.0};
37 auto &v_cl=create_vcluster();
38
39 vector_dist_ws<2, double, aggregate<double,double,double[2],double,double[2]>> Sparticles(0, domain,bc,ghost);
40 //Init_DCPSE(domain)
41 BOOST_TEST_MESSAGE("Init domain...");
42
43 // 1. Particles on a line
44 if (v_cl.rank() == 0) {
45 for (int i = 0; i < n; ++i) {
46 double xp = -1.5+i*grid_spacing;
47 Sparticles.add();
48 Sparticles.getLastPos()[0] = xp;
49 Sparticles.getLastPos()[1] = 0;
50 Sparticles.getLastProp<3>() = std::sin(xp);
51 Sparticles.getLastProp<2>()[0] = 0;
52 Sparticles.getLastProp<2>()[1] = 1.0;
53 Sparticles.getLastProp<1>() = -std::sin(xp);//sin(theta)*exp(-finalT/(radius*radius));
54 Sparticles.getLastSubset(0);
55 }
56 }
57 Sparticles.map();
58
59 BOOST_TEST_MESSAGE("Sync domain across processors...");
60
61 Sparticles.ghost_get<0,3>();
62 //Sparticles.write("Sparticles");
63 //Here template parameters are Normal property no.
64 SurfaceDerivative_xx<2> SDxx(Sparticles, 2, rCut,grid_spacing);
65 SurfaceDerivative_yy<2> SDyy(Sparticles, 2, rCut,grid_spacing);
66 //SurfaceDerivative_x<2> SDx(Sparticles, 4, rCut,grid_spacing);
67 //SurfaceDerivative_y<2> SDy(Sparticles, 4, rCut,grid_spacing);
68 auto INICONC = getV<3>(Sparticles);
69 auto CONC = getV<0>(Sparticles);
70 auto TEMP = getV<4>(Sparticles);
71 auto normal = getV<2>(Sparticles);
72 //auto ANASOL = getV<1>(domain);
73 CONC=SDxx(INICONC)+SDyy(INICONC);
74 //TEMP[0]=(-normal[0]*normal[0]+1.0) * SDx(INICONC) - normal[0]*normal[1] * SDy(INICONC);
75 //TEMP[1]=(-normal[1]*normal[1]+1.0) * SDy(INICONC) - normal[0]*normal[1] * SDx(INICONC);
76 //Sparticles.ghost_get<4>();
77 //CONC=SDxx(TEMP[0]) + SDyy(TEMP[1]);
78 auto it2 = Sparticles.getDomainIterator();
79 double worst = 0.0;
80 while (it2.isNext()) {
81 auto p = it2.get();
82 if (fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p)) > worst) {
83 worst = fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p));
84 }
85
86 ++it2;
87 }
88 Sparticles.deleteGhost();
89 //std::cout<<v_cl.rank()<<":WORST:"<<worst<<std::endl;
90 //Sparticles.write("Sparticles");
91 BOOST_REQUIRE(worst < 0.03);
92}
93 BOOST_AUTO_TEST_CASE(dcpse_surface_circle) {
94
95 double boxP1{-1.5}, boxP2{1.5};
96 double boxSize{boxP2 - boxP1};
97 size_t n=512;
98 auto &v_cl=create_vcluster();
99 //std::cout<<v_cl.rank()<<":Enter res: "<<std::endl;
100 //std::cin>>n;
101 size_t sz[2] = {n,n};
102 double grid_spacing{boxSize/(sz[0]-1)};
103 double rCut{5.1 * grid_spacing};
104
105 Box<2,double> domain{{boxP1,boxP1},{boxP2,boxP2}};
106 size_t bc[2] = {NON_PERIODIC,NON_PERIODIC};
107 Ghost<2,double> ghost{rCut + grid_spacing/8.0};
108 vector_dist_ws<2, double, aggregate<double,double,double[2],double,double[2],double>> Sparticles(0, domain,bc,ghost);
109 //Init_DCPSE(domain)
110 BOOST_TEST_MESSAGE("Init domain...");
111
112 // Surface prameters
113 const double radius{1.0};
114 std::array<double,2> center{0.0,0.0};
115 Point<2,double> coord;
116 const double pi{3.14159265358979323846};
117
118 // 1. Particles on surface
119 double theta{0.0};
120 double dtheta{2*pi/double(n)};
121 if (v_cl.rank() == 0) {
122 for (int i = 0; i < n; ++i) {
123 coord[0] = center[0] + radius * std::cos(theta);
124 coord[1] = center[1] + radius * std::sin(theta);
125 Sparticles.add();
126 Sparticles.getLastPos()[0] = coord[0];
127 Sparticles.getLastPos()[1] = coord[1];
128 Sparticles.getLastProp<3>() = std::sin(theta);
129 Sparticles.getLastProp<2>()[0] = std::cos(theta);
130 Sparticles.getLastProp<2>()[1] = std::sin(theta);
131 Sparticles.getLastProp<1>() = -std::sin(theta);;//sin(theta)*exp(-finalT/(radius*radius));
132 Sparticles.getLastSubset(0);
133 theta += dtheta;
134 }
135 }
136 Sparticles.map();
137
138 BOOST_TEST_MESSAGE("Sync domain across processors...");
139
140 Sparticles.ghost_get<0,3>();
141 //Sparticles.write("Sparticles");
142 //Here template parameters are Normal property no.
143 SurfaceDerivative_xx<2> SDxx(Sparticles, 2, rCut,grid_spacing);
144 SurfaceDerivative_yy<2> SDyy(Sparticles, 2, rCut,grid_spacing);
145 //SurfaceDerivative_xy<2> SDxy(Sparticles, 3, rCut,grid_spacing);
146 //SurfaceDerivative_x<2> SDx(Sparticles, 3, rCut,grid_spacing);
147 //SurfaceDerivative_y<2> SDy(Sparticles, 3, rCut,grid_spacing);
148 auto INICONC = getV<3>(Sparticles);
149 auto CONC = getV<0>(Sparticles);
150 auto TEMP = getV<4>(Sparticles);
151 auto normal = getV<2>(Sparticles);
152
153 CONC=SDxx(INICONC)+SDyy(INICONC);
154 //TEMP[0]=(-normal[0]*normal[0]+1.0) * SDx(INICONC) - normal[0]*normal[1] * SDy(INICONC);
155 //TEMP[1]=(-normal[1]*normal[1]+1.0) * SDy(INICONC) - normal[0]*normal[1] * SDx(INICONC);
156 //TEMP[0]=(-normal[0]*normal[0]+1.0);
157 //TEMP[1]=normal[0]*normal[1];
158 //Sparticles.ghost_get<2,4>();
159 //CONC=SDx(TEMP[0]) + SDy(TEMP[1]);
160 //CONC= (SDx(TEMP[0])*SDx(INICONC)+TEMP[0]*SDxx(INICONC)-(SDx(TEMP[1])*SDy(INICONC)+TEMP[1]*SDxy(INICONC)))+
161 // (SDy((-normal[1]*normal[1]+1.0))*SDy(INICONC)+(-normal[1]*normal[1]+1.0)*SDyy(INICONC)-(SDy(TEMP[1])*SDx(INICONC)+TEMP[1]*SDxy(INICONC)));
162 auto it2 = Sparticles.getDomainIterator();
163 double worst = 0.0;
164 while (it2.isNext()) {
165 auto p = it2.get();
166 Sparticles.getProp<5>(p) = fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p));
167 if (fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p)) > worst) {
168 worst = fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p));
169 }
170 ++it2;
171 }
172 Sparticles.deleteGhost();
173 //Sparticles.write("Sparticles");
174 //std::cout<<worst;
175 BOOST_REQUIRE(worst < 0.03);
176}
177 BOOST_AUTO_TEST_CASE(dcpse_surface_solver_circle) {
178 double boxP1{-1.5}, boxP2{1.5};
179 double boxSize{boxP2 - boxP1};
180 size_t n=512,k=2;
181 auto &v_cl=create_vcluster();
182 /*if(v_cl.rank()==0)
183 std::cout<<v_cl.rank()<<":Enter res: "<<std::endl;
184 std::cin>>n;
185 if(v_cl.rank()==0)
186 std::cout<<v_cl.rank()<<":Enter Freq: "<<std::endl;
187 std::cin>>k;*/
188 size_t sz[2] = {n,n};
189 double grid_spacing{boxSize/(sz[0]-1)};
190 double rCut{3.9 * grid_spacing};
191
192 Box<2,double> domain{{boxP1,boxP1},{boxP2,boxP2}};
193 size_t bc[2] = {NON_PERIODIC,NON_PERIODIC};
194 Ghost<2,double> ghost{rCut + grid_spacing/8.0};
195 vector_dist_ws<2, double, aggregate<double,double,double[2],double,double[2],double>> Sparticles(0, domain,bc,ghost);
196 //Init_DCPSE(domain)
197 BOOST_TEST_MESSAGE("Init domain...");
198
199 // Surface prameters
200 const double radius{1.0};
201 std::array<double,2> center{0.0,0.0};
202 Point<2,double> coord;
203 const double pi{3.14159265358979323846};
204
205 // 1. Particles on surface
206 double theta{0.0};
207 double dtheta{2*pi/double(n)};
208 if (v_cl.rank() == 0) {
209 for (int i = 0; i < n; ++i) {
210 coord[0] = center[0] + radius * std::cos(theta);
211 coord[1] = center[1] + radius * std::sin(theta);
212 Sparticles.add();
213 Sparticles.getLastPos()[0] = coord[0];
214 Sparticles.getLastPos()[1] = coord[1];
215 Sparticles.getLastProp<3>() = -openfpm::math::intpowlog(k,2)*std::sin(k*theta);
216 Sparticles.getLastProp<2>()[0] = std::cos(theta);
217 Sparticles.getLastProp<2>()[1] = std::sin(theta);
218 Sparticles.getLastProp<1>() = std::sin(k*theta);;//sin(theta)*exp(-finalT/(radius*radius));
219 Sparticles.getLastSubset(0);
220 if(coord[0]==1. && coord[1]==0.)
221 {Sparticles.getLastSubset(1);}
222 theta += dtheta;
223 }
224 }
225 Sparticles.map();
226
227 BOOST_TEST_MESSAGE("Sync domain across processors...");
228
229 Sparticles.ghost_get<0>();
230 //Sparticles.write("Sparticles");
231 vector_dist_subset<2, double, aggregate<double,double,double[2],double,double[2],double>> Sparticles_bulk(Sparticles,0);
232 vector_dist_subset<2, double, aggregate<double,double,double[2],double,double[2],double>> Sparticles_boundary(Sparticles,1);
233 auto & bulk=Sparticles_bulk.getIds();
234 auto & boundary=Sparticles_boundary.getIds();
235 //Here template parameters are Normal property no.
236 SurfaceDerivative_xx<2> SDxx(Sparticles, 2, rCut,grid_spacing);
237 SurfaceDerivative_yy<2> SDyy(Sparticles, 2, rCut,grid_spacing);
238 auto INICONC = getV<3>(Sparticles);
239 auto CONC = getV<0>(Sparticles);
240 auto TEMP = getV<4>(Sparticles);
241 auto normal = getV<2>(Sparticles);
242 auto ANASOL = getV<1>(Sparticles);
243 DCPSE_scheme<equations2d1,decltype(Sparticles)> Solver(Sparticles);
244 Solver.impose(SDxx(CONC)+SDyy(CONC), bulk, INICONC);
245 Solver.impose(CONC, boundary, ANASOL);
246 Solver.solve(CONC);
247 auto it2 = Sparticles.getDomainIterator();
248 double worst = 0.0;
249 while (it2.isNext()) {
250 auto p = it2.get();
251 Sparticles.getProp<5>(p) = fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p));
252 if (fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p)) > worst) {
253 worst = fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p));
254 }
255 ++it2;
256 }
257 Sparticles.deleteGhost();
258 //Sparticles.write("Sparticles");
259 //std::cout<<worst;
260 BOOST_REQUIRE(worst < 0.03);
261}
262BOOST_AUTO_TEST_CASE(dcpse_surface_sphere) {
263 auto & v_cl = create_vcluster();
264 timer tt;
265 tt.start();
266 size_t n=512;
267 size_t n_sp=n;
268 // Domain
269 double boxP1{-1.5}, boxP2{1.5};
270 double boxSize{boxP2 - boxP1};
271 size_t sz[3] = {n,n,n};
272 double grid_spacing{boxSize/(sz[0]-1)};
273 double grid_spacing_surf=grid_spacing*30;
274 double rCut{2.5 * grid_spacing_surf};
275
276 Box<3,double> domain{{boxP1,boxP1,boxP1},{boxP2,boxP2,boxP2}};
277 size_t bc[3] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
278 Ghost<3,double> ghost{rCut + grid_spacing/8.0};
279
280 constexpr int K = 1;
281 // particles
282 vector_dist_ws<3, double, aggregate<double,double,double[3],double,double[3],double>> Sparticles(0, domain,bc,ghost);
283 // 1. particles on the Spherical surface
284 double Golden_angle=M_PI * (3.0 - sqrt(5.0));
285 if (v_cl.rank() == 0) {
286 //std::vector<Vector3f> data;
287 //GenerateSphere(1,data);
288 for(int i=1;i<n_sp;i++)
289 {
290 double y = 1.0 - (i /double(n_sp - 1.0)) * 2.0;
291 double radius = sqrt(1 - y * y);
292 double Golden_theta = Golden_angle * i;
293 double x = cos(Golden_theta) * radius;
294 double z = sin(Golden_theta) * radius;
295 Sparticles.add();
296 Sparticles.getLastPos()[0] = x;
297 Sparticles.getLastPos()[1] = y;
298 Sparticles.getLastPos()[2] = z;
299 double rm=sqrt(x*x+y*y+z*z);
300 Sparticles.getLastProp<2>()[0] = x/rm;
301 Sparticles.getLastProp<2>()[1] = y/rm;
302 Sparticles.getLastProp<2>()[2] = z/rm;
303 Sparticles.getLastProp<4>()[0] = 1.0 ;
304 Sparticles.getLastProp<4>()[1] = std::atan2(sqrt(x*x+y*y),z);
305 Sparticles.getLastProp<4>()[2] = std::atan2(y,x);
306 if(i<=2*(K)+1)
307 {Sparticles.getLastSubset(1);}
308 else
309 {Sparticles.getLastSubset(0);}
310 }
311 //std::cout << "n: " << n << " - grid spacing: " << grid_spacing << " - rCut: " << rCut << "Surf Normal spacing" << grid_spacing<<std::endl;
312 }
313
314 Sparticles.map();
315 Sparticles.ghost_get<3>();
316
317 vector_dist_subset<3,double,aggregate<double,double,double[3],double,double[3],double>> Sparticles_bulk(Sparticles,0);
318 vector_dist_subset<3,double,aggregate<double,double,double[3],double,double[3],double>> Sparticles_boundary(Sparticles,1);
319 auto &bulkIds=Sparticles_bulk.getIds();
320 auto &bdrIds=Sparticles_boundary.getIds();
321 std::unordered_map<const lm,double,key_hash,key_equal> Alm;
322 //Setting max mode l_max
323 //Setting amplitudes to 1
324 for(int l=0;l<=K;l++){
325 for(int m=-l;m<=l;m++){
326 Alm[std::make_tuple(l,m)]=0;
327 }
328 }
329 Alm[std::make_tuple(1,0)]=1;
330 auto it2 = Sparticles.getDomainIterator();
331 while (it2.isNext()) {
332 auto p = it2.get();
333 Point<3, double> xP = Sparticles.getProp<4>(p);
334 /*double Sum=0;
335 for(int m=-spL;m<=spL;++m)
336 {
337 Sum+=openfpm::math::Y(spL,m,xP[1],xP[2]);
338 }*/
339 //Sparticles.getProp<ANADF>(p) = Sum;//openfpm::math::Y(K,K,xP[1],xP[2]);openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);;
340 Sparticles.getProp<3>(p)=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
341 Sparticles.getProp<1>(p)=-(K)*(K+1)*openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
342 ++it2;
343 }
344 auto f=getV<3>(Sparticles);
345 auto Df=getV<0>(Sparticles);
346
347 SurfaceDerivative_xx<2> Sdxx{Sparticles,2,rCut,grid_spacing_surf};
348 SurfaceDerivative_yy<2> Sdyy{Sparticles,2,rCut,grid_spacing_surf};
349 SurfaceDerivative_zz<2> Sdzz{Sparticles,2,rCut,grid_spacing_surf};
350 //Laplace_Beltrami<2> SLap{Sparticles,2,rCut,grid_spacing_surf};
351 //Sdyy.DrawKernel<5>(Sparticles,0);
352 //Sdzz.DrawKernel<5>(Sparticles,0);
353/* std::cout<<"SDXX:"<<std::endl;
354 Sdxx.checkMomenta(Sparticles);
355 std::cout<<"SDYY:"<<std::endl;
356 Sdyy.checkMomenta(Sparticles);
357 std::cout<<"SDZZ:"<<std::endl;
358 Sdzz.checkMomenta(Sparticles);*/
359
360 Sparticles.ghost_get<3>();
361 Df=(Sdxx(f)+Sdyy(f)+Sdzz(f));
362 //Df=SLap(f);
363 auto it3 = Sparticles.getDomainIterator();
364 double worst = 0.0;
365 while (it3.isNext()) {
366 auto p = it3.get();
367 //Sparticles.getProp<5>(p) = fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p));
368 if (fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p)) > worst) {
369 worst = fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p));
370 }
371 ++it3;
372 }
373 Sparticles.deleteGhost();
374 //Sparticles.write("Sparticles");
375 //std::cout<<worst;
376 BOOST_REQUIRE(worst < 0.03);
377}
378
379
380BOOST_AUTO_TEST_CASE(dcpse_surface_sphere_old) {
381 auto & v_cl = create_vcluster();
382 timer tt;
383 tt.start();
384 size_t n=512;
385 size_t n_sp=n;
386 // Domain
387 double boxP1{-1.5}, boxP2{1.5};
388 double boxSize{boxP2 - boxP1};
389 size_t sz[3] = {n,n,n};
390 double grid_spacing{boxSize/(sz[0]-1)};
391 double grid_spacing_surf=grid_spacing*30;
392 double rCut{2.5 * grid_spacing_surf};
393
394 Box<3,double> domain{{boxP1,boxP1,boxP1},{boxP2,boxP2,boxP2}};
395 size_t bc[3] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
396 Ghost<3,double> ghost{rCut + grid_spacing/8.0};
397
398 constexpr int K = 1;
399 // particles
400 vector_dist_ws<3, double, aggregate<double,double,double[3],double,double[3],double>> Sparticles(0, domain,bc,ghost);
401 // 1. particles on the Spherical surface
402 double Golden_angle=M_PI * (3.0 - sqrt(5.0));
403 if (v_cl.rank() == 0) {
404 //std::vector<Vector3f> data;
405 //GenerateSphere(1,data);
406 std::unordered_map<const lm,double,key_hash,key_equal> Alm;
407 //Setting max mode l_max
408 //Setting amplitudes to 1
409 for(int l=0;l<=K;l++){
410 for(int m=-l;m<=l;m++){
411 Alm[std::make_tuple(l,m)]=0;
412 }
413 }
414 Alm[std::make_tuple(1,0)]=1;
415 for(int i=1;i<n_sp;i++)
416 {
417 double y = 1.0 - (i /double(n_sp - 1.0)) * 2.0;
418 double radius = sqrt(1 - y * y);
419 double Golden_theta = Golden_angle * i;
420 double x = cos(Golden_theta) * radius;
421 double z = sin(Golden_theta) * radius;
422 Sparticles.add();
423 Sparticles.getLastPos()[0] = x;
424 Sparticles.getLastPos()[1] = y;
425 Sparticles.getLastPos()[2] = z;
426 double rm=sqrt(x*x+y*y+z*z);
427 Sparticles.getLastProp<2>()[0] = x/rm;
428 Sparticles.getLastProp<2>()[1] = y/rm;
429 Sparticles.getLastProp<2>()[2] = z/rm;
430 Sparticles.getLastProp<4>()[0] = 1.0 ;
431 Sparticles.getLastProp<4>()[1] = std::atan2(sqrt(x*x+y*y),z);
432 Sparticles.getLastProp<4>()[2] = std::atan2(y,x);
433 double m1=openfpm::math::sumY_Scalar<K>(1.0,std::atan2(sqrt(x*x+y*y),z),std::atan2(y,x),Alm);
434 double m2=-(K)*(K+1)*openfpm::math::sumY_Scalar<K>(1.0,std::atan2(sqrt(x*x+y*y),z),std::atan2(y,x),Alm);
435 Sparticles.getLastProp<3>()=m1;
436 Sparticles.getLastProp<1>()=m2;
437 Sparticles.getLastSubset(0);
438 for(int j=1;j<=2;++j){
439 Sparticles.add();
440 Sparticles.getLastPos()[0] = x+j*grid_spacing_surf*x/rm;
441 Sparticles.getLastPos()[1] = y+j*grid_spacing_surf*y/rm;
442 Sparticles.getLastPos()[2] = z+j*grid_spacing_surf*z/rm;
443 Sparticles.getLastProp<3>()=m1;
444 Sparticles.getLastSubset(1);
445 //Sparticles.getLastProp<1>(p)=m2;
446 Sparticles.add();
447 Sparticles.getLastPos()[0] = x-j*grid_spacing_surf*x/rm;
448 Sparticles.getLastPos()[1] = y-j*grid_spacing_surf*y/rm;
449 Sparticles.getLastPos()[2] = z-j*grid_spacing_surf*z/rm;
450 Sparticles.getLastProp<3>()=m1;
451 Sparticles.getLastSubset(1);
452 //Sparticles.getLastProp<1>(p)=m2;
453 }
454 }
455 //std::cout << "n: " << n << " - grid spacing: " << grid_spacing << " - rCut: " << rCut << "Surf Normal spacing" << grid_spacing<<std::endl;
456 }
457
458 Sparticles.map();
459 Sparticles.ghost_get<3>();
460 //Sparticles.write("SparticlesInit");
461
462 vector_dist_subset<3,double,aggregate<double,double,double[3],double,double[3],double>> Sparticles_bulk(Sparticles,0);
463 vector_dist_subset<3,double,aggregate<double,double,double[3],double,double[3],double>> Sparticles_boundary(Sparticles,1);
464 auto &bulkIds=Sparticles_bulk.getIds();
465 auto &bdrIds=Sparticles_boundary.getIds();
466 /*auto it2 = Sparticles.getDomainIterator();
467 while (it2.isNext()) {
468 auto p = it2.get();
469 Point<3, double> xP = Sparticles.getProp<4>(p);
470 *//*double Sum=0;
471 for(int m=-spL;m<=spL;++m)
472 {
473 Sum+=openfpm::math::Y(spL,m,xP[1],xP[2]);
474 }*//*
475 //Sparticles.getProp<ANADF>(p) = Sum;//openfpm::math::Y(K,K,xP[1],xP[2]);openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);;
476 Sparticles.getProp<3>(p)=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
477 Sparticles.getProp<1>(p)=-(K)*(K+1)*openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
478 ++it2;
479 }*/
480 auto f=getV<3>(Sparticles);
481 auto Df=getV<0>(Sparticles);
482
483 //SurfaceDerivative_xx<2> Sdxx{Sparticles,2,rCut,grid_spacing_surf};
484 //SurfaceDerivative_yy<2> Sdyy{Sparticles,2,rCut,grid_spacing_surf};
485 //SurfaceDerivative_zz<2> Sdzz{Sparticles,2,rCut,grid_spacing_surf};
486 Derivative_xx Sdxx{Sparticles,2,rCut};
487 //std::cout<<"Dxx Done"<<std::endl;
488 Derivative_yy Sdyy{Sparticles,2,rCut};
489 //std::cout<<"Dyy Done"<<std::endl;
490 Derivative_zz Sdzz{Sparticles,2,rCut};
491 //std::cout<<"Dzz Done"<<std::endl;
492
493 //Laplace_Beltrami<2> SLap{Sparticles,2,rCut,grid_spacing_surf};
494 //SLap.DrawKernel<5>(Sparticles,73);
495 //Sdxx.DrawKernel<5>(Sparticles,0);
496 //Sdyy.DrawKernel<5>(Sparticles,0);
497 //Sdzz.DrawKernel<5>(Sparticles,0);
498/* std::cout<<"SDXX:"<<std::endl;
499 Sdxx.checkMomenta(Sparticles);
500 std::cout<<"SDYY:"<<std::endl;
501 Sdyy.checkMomenta(Sparticles);
502 std::cout<<"SDZZ:"<<std::endl;
503 Sdzz.checkMomenta(Sparticles);*/
504
505 Sparticles.ghost_get<3>();
506 Df=(Sdxx(f)+Sdyy(f)+Sdzz(f));
507 //Df=SLap(f);
508 //auto it3 = Sparticles_bulk.getDomainIterator();
509 double worst = 0.0;
510 for (int j = 0; j < bulkIds.size(); j++) {
511 auto p = bulkIds.get<0>(j);
512 //Sparticles.getProp<5>(p) = fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p));
513 if (fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p)) > worst) {
514 worst = fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p));
515 }
516 }
517 Sparticles.deleteGhost();
518 //Sparticles.write("SparticlesNoo");
519 //std::cout<<"Worst: "<<worst<<std::endl;
520 BOOST_REQUIRE(worst < 0.03);
521}
522
523/*BOOST_AUTO_TEST_CASE(dcpse_surface_sphere_proj) {
524 auto & v_cl = create_vcluster();
525 timer tt;
526 tt.start();
527 size_t n=512;
528 size_t n_sp=n;
529 // Domain
530 double boxP1{-1.5}, boxP2{1.5};
531 double boxSize{boxP2 - boxP1};
532 size_t sz[3] = {n,n,n};
533 double grid_spacing{boxSize/(sz[0]-1)};
534 double grid_spacing_surf=grid_spacing*30;
535 double rCut{2.5 * grid_spacing_surf};
536
537 Box<3,double> domain{{boxP1,boxP1,boxP1},{boxP2,boxP2,boxP2}};
538 size_t bc[3] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
539 Ghost<3,double> ghost{rCut + grid_spacing/8.0};
540
541 constexpr int K = 1;
542 // particles
543 vector_dist_ws<3, double, aggregate<VectorS<3,double>,VectorS<3,double>,VectorS<3,double>,VectorS<3,double>,VectorS<3,double>,double>> Sparticles(0, domain,bc,ghost);
544 // 1. particles on the Spherical surface
545 double Golden_angle=M_PI * (3.0 - sqrt(5.0));
546 if (v_cl.rank() == 0) {
547 //std::vector<Vector3f> data;
548 //GenerateSphere(1,data);
549 for(int i=1;i<n_sp;i++)
550 {
551 double y = 1.0 - (i /double(n_sp - 1.0)) * 2.0;
552 double radius = sqrt(1 - y * y);
553 double Golden_theta = Golden_angle * i;
554 double x = cos(Golden_theta) * radius;
555 double z = sin(Golden_theta) * radius;
556 Sparticles.add();
557 Sparticles.getLastPos()[0] = x;
558 Sparticles.getLastPos()[1] = y;
559 Sparticles.getLastPos()[2] = z;
560 double rm=sqrt(x*x+y*y+z*z);
561 Sparticles.getLastProp<2>()[0] = x/rm;
562 Sparticles.getLastProp<2>()[1] = y/rm;
563 Sparticles.getLastProp<2>()[2] = z/rm;
564 Sparticles.getLastProp<4>()[0] = 1.0 ;
565 Sparticles.getLastProp<4>()[1] = std::atan2(sqrt(x*x+y*y),z);
566 Sparticles.getLastProp<4>()[2] = std::atan2(y,x);
567 if(i<=2*(K)+1)
568 {Sparticles.getLastSubset(1);}
569 else
570 {Sparticles.getLastSubset(0);}
571 }
572 //std::cout << "n: " << n << " - grid spacing: " << grid_spacing << " - rCut: " << rCut << "Surf Normal spacing" << grid_spacing<<std::endl;
573 }
574
575 Sparticles.map();
576 Sparticles.ghost_get<3>();
577
578 vector_dist_subset<3,double,aggregate<VectorS<3,double>,VectorS<3,double>,VectorS<3,double>,VectorS<3,double>,VectorS<3,double>,double>> Sparticles_bulk(Sparticles,0);
579 vector_dist_subset<3,double,aggregate<VectorS<3,double>,VectorS<3,double>,VectorS<3,double>,VectorS<3,double>,VectorS<3,double>,double>> Sparticles_boundary(Sparticles,1);
580 auto &bulkIds=Sparticles_bulk.getIds();
581 auto &bdrIds=Sparticles_boundary.getIds();
582 std::unordered_map<const lm,double,key_hash,key_equal> Alm;
583 //Setting max mode l_max
584 //Setting amplitudes to 1
585 for(int l=0;l<=K;l++){
586 for(int m=-l;m<=l;m++){
587 Alm[std::make_tuple(l,m)]=0;
588 }
589 }
590 Alm[std::make_tuple(1,0)]=1;
591 auto it2 = Sparticles.getDomainIterator();
592 while (it2.isNext()) {
593 auto p = it2.get();
594 Point<3, double> xP = Sparticles.getProp<4>(p);
595 *//*double Sum=0;
596 for(int m=-spL;m<=spL;++m)
597 {
598 Sum+=openfpm::math::Y(spL,m,xP[1],xP[2]);
599 }*//*
600 //Sparticles.getProp<ANADF>(p) = Sum;//openfpm::math::Y(K,K,xP[1],xP[2]);openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);;
601 Sparticles.getProp<3>(p)[0]=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
602 Sparticles.getProp<3>(p)[1]=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
603 Sparticles.getProp<3>(p)[2]=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
604 Sparticles.getProp<1>(p)[0]=-(K)*(K+1)*openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
605 Sparticles.getProp<1>(p)[1]=-(K)*(K+1)*openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
606 Sparticles.getProp<1>(p)[2]=-(K)*(K+1)*openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
607 ++it2;
608 }
609 auto f=getV<3>(Sparticles);
610 auto Df=getV<0>(Sparticles);
611
612 SurfaceProjectedGradient<2> SGP{Sparticles,2,rCut,grid_spacing_surf};
613
614 Sparticles.ghost_get<3>();
615 Df=SGP(f);
616 //Df=SLap(f);
617 auto it3 = Sparticles.getDomainIterator();
618 double worst = 0.0;
619 while (it3.isNext()) {
620 auto p = it3.get();
621 //Sparticles.getProp<5>(p) = fabs(Sparticles.getProp<1>(p) - Sparticles.getProp<0>(p));
622 if (fabs(Sparticles.getProp<1>(p)[0] - Sparticles.getProp<0>(p)[0]) > worst) {
623 worst = fabs(Sparticles.getProp<1>(p)[0] - Sparticles.getProp<0>(p)[0]);
624 }
625 ++it3;
626 }
627 Sparticles.deleteGhost();
628 //Sparticles.write("Sparticles");
629 //std::cout<<worst;
630 BOOST_REQUIRE(worst < 0.03);
631}*/
632
633
634 BOOST_AUTO_TEST_CASE(dcpse_surface_adaptive) {
635// int rank;
636// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
637 const size_t sz[3] = {81,81,1};
638 Box<3, double> box({0, 0,-5}, {0.5, 0.5,5});
639 size_t bc[3] = {NON_PERIODIC, NON_PERIODIC,NON_PERIODIC};
640 double spacing = box.getHigh(0) / (sz[0] - 1);
641 Ghost<3, double> ghost(spacing * 3.1);
642 double rCut = 3.1 * spacing;
643 BOOST_TEST_MESSAGE("Init vector_dist...");
644
645 vector_dist<3, double, aggregate<double,double,double,double,double,double,double[3]>> domain(0, box, bc, ghost);
646
647
648 //Init_DCPSE(domain)
649 BOOST_TEST_MESSAGE("Init domain...");
650
651 auto it = domain.getGridIterator(sz);
652 while (it.isNext()) {
653 domain.add();
654
655 auto key = it.get();
656 double x = key.get(0) * it.getSpacing(0);
657 domain.getLastPos()[0] = x;
658 double y = key.get(1) * it.getSpacing(1);
659 domain.getLastPos()[1] = y;
660
661 domain.getLastPos()[2] = 0;
662
663 ++it;
664 }
665
666 // Add multi res patch 1
667
668 {
669 const size_t sz2[3] = {40,40,1};
670 Box<3,double> bx({0.25 + it.getSpacing(0)/4.0,0.25 + it.getSpacing(0)/4.0,-0.5},{sz2[0]*it.getSpacing(0)/2.0 + 0.25 + it.getSpacing(0)/4.0, sz2[1]*it.getSpacing(0)/2.0 + 0.25 + it.getSpacing(0)/4.0,0.5});
671 openfpm::vector<size_t> rem;
672
673 auto it = domain.getDomainIterator();
674
675 while (it.isNext())
676 {
677 auto k = it.get();
678
679 Point<3,double> xp = domain.getPos(k);
680
681 if (bx.isInside(xp) == true)
682 {
683 rem.add(k.getKey());
684 }
685
686 ++it;
687 }
688
689 domain.remove(rem);
690
691 auto it2 = domain.getGridIterator(sz2);
692 while (it2.isNext()) {
693 domain.add();
694
695 auto key = it2.get();
696 double x = key.get(0) * spacing/2.0 + 0.25 + spacing/4.0;
697 domain.getLastPos()[0] = x;
698 double y = key.get(1) * spacing/2.0 + 0.25 + spacing/4.0;
699 domain.getLastPos()[1] = y;
700 domain.getLastPos()[2] = 0;
701
702 ++it2;
703 }
704 }
705
706 // Add multi res patch 2
707
708 {
709 const size_t sz2[3] = {40,40,1};
710 Box<3,double> bx({0.25 + 21.0*spacing/8.0,0.25 + 21.0*spacing/8.0,-5},{sz2[0]*spacing/4.0 + 0.25 + 21.0*spacing/8.0, sz2[1]*spacing/4.0 + 0.25 + 21*spacing/8.0,5});
711 openfpm::vector<size_t> rem;
712
713 auto it = domain.getDomainIterator();
714
715 while (it.isNext())
716 {
717 auto k = it.get();
718
719 Point<3,double> xp = domain.getPos(k);
720
721 if (bx.isInside(xp) == true)
722 {
723 rem.add(k.getKey());
724 }
725
726 ++it;
727 }
728
729 domain.remove(rem);
730
731 auto it2 = domain.getGridIterator(sz2);
732 while (it2.isNext()) {
733 domain.add();
734 auto key = it2.get();
735 double x = key.get(0) * spacing/4.0 + 0.25 + 21*spacing/8.0;
736 domain.getLastPos()[0] = x;
737 double y = key.get(1) * spacing/4.0 + 0.25 + 21*spacing/8.0;
738 domain.getLastPos()[1] = y;
739 domain.getLastPos()[2] = 0;
740
741 ++it2;
742 }
743 }
744
746
747 BOOST_TEST_MESSAGE("Sync domain across processors...");
748
749 domain.map();
750 domain.ghost_get<0>();
751
752 openfpm::vector<aggregate<int>> bulk;
753 openfpm::vector<aggregate<int>> up_p;
754 openfpm::vector<aggregate<int>> dw_p;
755 openfpm::vector<aggregate<int>> l_p;
756 openfpm::vector<aggregate<int>> r_p;
757 openfpm::vector<aggregate<int>> ref_p;
758
759 auto v = getV<0>(domain);
760 auto RHS=getV<1>(domain);
761 auto sol = getV<2>(domain);
762 auto anasol = getV<3>(domain);
763 auto err = getV<4>(domain);
764 auto DCPSE_sol=getV<5>(domain);
765
766 // Here fill me
767
768 Box<3, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0,-5},
769 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0,5});
770
771 Box<3, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0,-5},
772 {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0,5});
773
774 Box<3, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0,-5},
775 {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0,5});
776
777 Box<3, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0,-5},
778 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0,5});
779
780 openfpm::vector<Box<3, double>> boxes;
781 boxes.add(up);
782 boxes.add(down);
783 boxes.add(left);
784 boxes.add(right);
785
786 // Create a writer and write
787 VTKWriter<openfpm::vector<Box<3, double>>, VECTOR_BOX> vtk_box;
788 vtk_box.add(boxes);
789 //vtk_box.write("vtk_box.vtk");
790
791 auto it2 = domain.getDomainIterator();
792
793 while (it2.isNext()) {
794 auto p = it2.get();
795 Point<3, double> xp = domain.getPos(p);
796 if (up.isInside(xp) == true) {
797 up_p.add();
798 up_p.last().get<0>() = p.getKey();
799 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
800 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
801 } else if (down.isInside(xp) == true) {
802 dw_p.add();
803 dw_p.last().get<0>() = p.getKey();
804 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
805 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
806
807 } else if (left.isInside(xp) == true) {
808 l_p.add();
809 l_p.last().get<0>() = p.getKey();
810 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
811 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
812
813 } else if (right.isInside(xp) == true) {
814 r_p.add();
815 r_p.last().get<0>() = p.getKey();
816 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
817 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
818
819 } else {
820 bulk.add();
821 bulk.last().get<0>() = p.getKey();
822 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
823 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
824 }
825 domain.getProp<6>(p)[0] = 0;
826 domain.getProp<6>(p)[1] = 0;
827 domain.getProp<6>(p)[2] = 1;
828
829 ++it2;
830 }
831
832 domain.ghost_get<1,2,3>();
833 SurfaceDerivative_xx<6> Dxx(domain, 2, rCut,3.9,support_options::ADAPTIVE);
834
835/* v=0;
836 auto itNNN=domain.getDomainIterator();
837 while(itNNN.isNext()){
838 auto p=itNNN.get().getKey();
839 Dxx.DrawKernel<0,decltype(domain)>(domain,p);
840 domain.write_frame("Kernel",p);
841 v=0;
842 ++itNNN;
843 }
844
845*/
846 //Dxx.DrawKernel<5,decltype(domain)>(domain,6161);
847 //domain.write_frame("Kernel",6161);
848
849 SurfaceDerivative_yy<6> Dyy(domain, 2, rCut,3.9,support_options::ADAPTIVE);
850 SurfaceDerivative_zz<6> Dzz(domain, 2, rCut,3.9,support_options::ADAPTIVE);
851
852 Dxx.save(domain,"Sdxx_test");
853 Dyy.save(domain,"Sdyy_test");
854 Dzz.save(domain,"Sdzz_test");
855
856 domain.ghost_get<2>();
857 sol=Dxx(anasol)+Dyy(anasol)+Dzz(anasol);
858 domain.ghost_get<5>();
859
860 double worst1 = 0.0;
861
862 for(int j=0;j<bulk.size();j++)
863 { auto p=bulk.get<0>(j);
864 if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) >= worst1) {
865 worst1 = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
866 }
867 domain.getProp<4>(p) = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
868
869 }
870 //std::cout << "Maximum Analytic Error: " << worst1 << std::endl;
871 //domain.ghost_get<4>();
872 //domain.write("Robin_anasol");
873 BOOST_REQUIRE(worst1 < 0.03);
874
875 }
876
877
878 BOOST_AUTO_TEST_CASE(dcpse_surface_adaptive_load) {
879// int rank;
880// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
881 const size_t sz[3] = {81,81,1};
882 Box<3, double> box({0, 0,-5}, {0.5, 0.5,5});
883 size_t bc[3] = {NON_PERIODIC, NON_PERIODIC,NON_PERIODIC};
884 double spacing = box.getHigh(0) / (sz[0] - 1);
885 Ghost<3, double> ghost(spacing * 3.1);
886 double rCut = 3.1 * spacing;
887 BOOST_TEST_MESSAGE("Init vector_dist...");
888
889 vector_dist<3, double, aggregate<double,double,double,double,double,double,double[3]>> domain(0, box, bc, ghost);
890
891
892 //Init_DCPSE(domain)
893 BOOST_TEST_MESSAGE("Init domain...");
894
895 auto it = domain.getGridIterator(sz);
896 while (it.isNext()) {
897 domain.add();
898
899 auto key = it.get();
900 double x = key.get(0) * it.getSpacing(0);
901 domain.getLastPos()[0] = x;
902 double y = key.get(1) * it.getSpacing(1);
903 domain.getLastPos()[1] = y;
904
905 domain.getLastPos()[2] = 0;
906
907 ++it;
908 }
909
910 // Add multi res patch 1
911
912 {
913 const size_t sz2[3] = {40,40,1};
914 Box<3,double> bx({0.25 + it.getSpacing(0)/4.0,0.25 + it.getSpacing(0)/4.0,-0.5},{sz2[0]*it.getSpacing(0)/2.0 + 0.25 + it.getSpacing(0)/4.0, sz2[1]*it.getSpacing(0)/2.0 + 0.25 + it.getSpacing(0)/4.0,0.5});
915 openfpm::vector<size_t> rem;
916
917 auto it = domain.getDomainIterator();
918
919 while (it.isNext())
920 {
921 auto k = it.get();
922
923 Point<3,double> xp = domain.getPos(k);
924
925 if (bx.isInside(xp) == true)
926 {
927 rem.add(k.getKey());
928 }
929
930 ++it;
931 }
932
933 domain.remove(rem);
934
935 auto it2 = domain.getGridIterator(sz2);
936 while (it2.isNext()) {
937 domain.add();
938
939 auto key = it2.get();
940 double x = key.get(0) * spacing/2.0 + 0.25 + spacing/4.0;
941 domain.getLastPos()[0] = x;
942 double y = key.get(1) * spacing/2.0 + 0.25 + spacing/4.0;
943 domain.getLastPos()[1] = y;
944 domain.getLastPos()[2] = 0;
945
946 ++it2;
947 }
948 }
949
950 // Add multi res patch 2
951
952 {
953 const size_t sz2[3] = {40,40,1};
954 Box<3,double> bx({0.25 + 21.0*spacing/8.0,0.25 + 21.0*spacing/8.0,-5},{sz2[0]*spacing/4.0 + 0.25 + 21.0*spacing/8.0, sz2[1]*spacing/4.0 + 0.25 + 21*spacing/8.0,5});
955 openfpm::vector<size_t> rem;
956
957 auto it = domain.getDomainIterator();
958
959 while (it.isNext())
960 {
961 auto k = it.get();
962
963 Point<3,double> xp = domain.getPos(k);
964
965 if (bx.isInside(xp) == true)
966 {
967 rem.add(k.getKey());
968 }
969
970 ++it;
971 }
972
973 domain.remove(rem);
974
975 auto it2 = domain.getGridIterator(sz2);
976 while (it2.isNext()) {
977 domain.add();
978
979 auto key = it2.get();
980 double x = key.get(0) * spacing/4.0 + 0.25 + 21*spacing/8.0;
981 domain.getLastPos()[0] = x;
982 double y = key.get(1) * spacing/4.0 + 0.25 + 21*spacing/8.0;
983 domain.getLastPos()[1] = y;
984 domain.getLastPos()[2] = 0;
985
986 ++it2;
987 }
988 }
989
991
992 BOOST_TEST_MESSAGE("Sync domain across processors...");
993
994 domain.map();
995 domain.ghost_get<0>();
996
997 openfpm::vector<aggregate<int>> bulk;
998 openfpm::vector<aggregate<int>> up_p;
999 openfpm::vector<aggregate<int>> dw_p;
1000 openfpm::vector<aggregate<int>> l_p;
1001 openfpm::vector<aggregate<int>> r_p;
1002 openfpm::vector<aggregate<int>> ref_p;
1003
1004 auto v = getV<0>(domain);
1005 auto RHS=getV<1>(domain);
1006 auto sol = getV<2>(domain);
1007 auto anasol = getV<3>(domain);
1008 auto err = getV<4>(domain);
1009 auto DCPSE_sol=getV<5>(domain);
1010
1011 // Here fill me
1012
1013 Box<3, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0,-5},
1014 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0,5});
1015
1016 Box<3, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0,-5},
1017 {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0,5});
1018
1019 Box<3, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0,-5},
1020 {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0,5});
1021
1022 Box<3, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0,-5},
1023 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0,5});
1024
1025 openfpm::vector<Box<3, double>> boxes;
1026 boxes.add(up);
1027 boxes.add(down);
1028 boxes.add(left);
1029 boxes.add(right);
1030
1031 // Create a writer and write
1032 VTKWriter<openfpm::vector<Box<3, double>>, VECTOR_BOX> vtk_box;
1033 vtk_box.add(boxes);
1034 //vtk_box.write("vtk_box.vtk");
1035
1036
1037 auto it2 = domain.getDomainIterator();
1038
1039 while (it2.isNext()) {
1040 auto p = it2.get();
1041 Point<3, double> xp = domain.getPos(p);
1042 if (up.isInside(xp) == true) {
1043 up_p.add();
1044 up_p.last().get<0>() = p.getKey();
1045 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1046 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1047 } else if (down.isInside(xp) == true) {
1048 dw_p.add();
1049 dw_p.last().get<0>() = p.getKey();
1050 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1051 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1052
1053 } else if (left.isInside(xp) == true) {
1054 l_p.add();
1055 l_p.last().get<0>() = p.getKey();
1056 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1057 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1058
1059 } else if (right.isInside(xp) == true) {
1060 r_p.add();
1061 r_p.last().get<0>() = p.getKey();
1062 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1063 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1064
1065 } else {
1066 bulk.add();
1067 bulk.last().get<0>() = p.getKey();
1068 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1069 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1070 }
1071 domain.getProp<6>(p)[0] = 0;
1072 domain.getProp<6>(p)[1] = 0;
1073 domain.getProp<6>(p)[2] = 1;
1074
1075 ++it2;
1076 }
1077
1078 domain.ghost_get<1,2,3>();
1079 SurfaceDerivative_xx<6> Dxx(domain, 2, rCut,3.9,support_options::LOAD);
1080 SurfaceDerivative_yy<6> Dyy(domain, 2, rCut,3.9,support_options::LOAD);
1081 SurfaceDerivative_zz<6> Dzz(domain, 2, rCut,3.9,support_options::LOAD);
1082 Dxx.load(domain,"Sdxx_test");
1083 Dyy.load(domain,"Sdyy_test");
1084 Dzz.load(domain,"Sdzz_test");
1085
1086 domain.ghost_get<2>();
1087 sol=Dxx(anasol)+Dyy(anasol)+Dzz(anasol);
1088 domain.ghost_get<5>();
1089
1090 double worst1 = 0.0;
1091
1092 for(int j=0;j<bulk.size();j++)
1093 { auto p=bulk.get<0>(j);
1094 if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) >= worst1) {
1095 worst1 = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
1096 }
1097 domain.getProp<4>(p) = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
1098
1099 }
1100 //std::cout << "Maximum Analytic Error: " << worst1 << std::endl;
1101
1102 //domain.ghost_get<4>();
1103 //domain.write("Robin_anasol");
1104 BOOST_REQUIRE(worst1 < 0.03);
1105
1106 }
1107
1108
1109BOOST_AUTO_TEST_SUITE_END()
1110
1111#endif
1112#endif
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
Point
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28
Point::get
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition Point.hpp:172
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition map_vector.hpp:203
openfpm::vector::size
size_t size()
Stub size.
Definition map_vector.hpp:211
timer
Class for cpu time benchmarking.
Definition timer.hpp:28
timer::start
void start()
Start the timer.
Definition timer.hpp:90
vector_dist_subset
Definition vector_dist_subset.hpp:84
vector_dist_ws
Definition vector_dist_subset.hpp:17
vector_dist
Distributed vector.
Definition vector_dist.hpp:305
equations2d1
Specify the general characteristic of system to solve.
Definition EqnsStruct.hpp:14