1 /* GridMG.c */
2 /**********************************************************************************************************
3 Copyright (c) 2002-2013 Abdul-Rahman Allouche. All rights reserved
4
5 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
6 documentation files (the Gabedit), to deal in the Software without restriction, including without limitation
7 the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
8 and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
9
10 The above copyright notice and this permission notice shall be included in all copies or substantial portions
11 of the Software.
12
13 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
14 TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
15 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
16 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
17 DEALINGS IN THE SOFTWARE.
18 ************************************************************************************************************/
19
20 #include "../../Config.h"
21 #include <math.h>
22 #include <stdio.h>
23 #include <string.h>
24 #include <stdlib.h>
25 #include <ctype.h>
26 #include <gtk/gtk.h>
27 #ifdef ENABLE_OMP
28 #include <omp.h>
29 #endif
30 #include "../Utils/Vector3d.h"
31 #include "../Utils/Transformation.h"
32 #include "../Utils/Constants.h"
33 #include "GridMG.h"
34
35 #define PRECISION 1e-10
36
37
38 /*********************************************************/
39 /* private methods for GridMG */
40 static void initAllGridMG(GridMG* g, gdouble);
41 static void initInteriorGridMG(GridMG* g, gdouble);
42 static void initBoundaryGridMG(GridMG* g, gdouble);
43
44 static void equalAllGridMG(GridMG* g, GridMG* src);
45 static void equalInteriorGridMG(GridMG* g, GridMG* src);
46 static void equalBoundaryGridMG(GridMG* g, GridMG* src);
47
48 static void plusEqualAllGridMG(GridMG* g, GridMG* src);
49 static void plusEqualInteriorGridMG(GridMG* g, GridMG* src);
50 static void plusEqualBoundaryGridMG(GridMG* g, GridMG* src);
51
52 static void moinsEqualAllGridMG(GridMG* g, GridMG* src);
53 static void moinsEqualInteriorGridMG(GridMG* g, GridMG* src);
54 static void moinsEqualBoundaryGridMG(GridMG* g, GridMG* src);
55
56 static void multEqualAllGridMG(GridMG* g, GridMG* src);
57 static void multEqualInteriorGridMG(GridMG* g, GridMG* src);
58 static void multEqualBoundaryGridMG(GridMG* g, GridMG* src);
59
60 static void multEqualAllRealGridMG(GridMG* g, gdouble a);
61 static void multEqualInteriorRealGridMG(GridMG* g, gdouble a);
62 static void multEqualBoundaryRealGridMG(GridMG* g, gdouble a);
63
64 static void divEqualAllRealGridMG(GridMG* g, gdouble a);
65 static void divEqualInteriorRealGridMG(GridMG* g, gdouble a);
66 static void divEqualBoundaryRealGridMG(GridMG* g, gdouble a);
67
68 static gdouble dotAllGridMG(GridMG* g, GridMG* src);
69 static gdouble dotInteriorGridMG(GridMG* g, GridMG* src);
70 static gdouble dotBoundaryGridMG(GridMG* g, GridMG* src);
71
72 static gdouble normAllGridMG(GridMG* g);
73 static gdouble normInteriorGridMG(GridMG* g);
74 static gdouble normBoundaryGridMG(GridMG* g);
75
76 static gdouble normDiffAllGridMG(GridMG* g, GridMG* src);
77 static gdouble normDiffInteriorGridMG(GridMG* g, GridMG* src);
78 static gdouble normDiffBoundaryGridMG(GridMG* g, GridMG* src);
79
80 static gdouble sommeAllGridMG(GridMG* g);
81 static gdouble sommeInteriorGridMG(GridMG* g);
82 static gdouble sommeBoundaryGridMG(GridMG* g);
83 static void tradesBoundaryPeriodicGridMG(GridMG* g);
84
85
86 static void printAllGridMG(GridMG* g);
87 static void printInteriorGridMG(GridMG* g);
88 static void printBoundaryGridMG(GridMG* g);
89 /*********************************************************/
destroyGridMG(GridMG * g)90 void destroyGridMG(GridMG* g)
91 {
92 if(g && g->values) g_free(g->values);
93 }
94 /*********************************************************/
getNewGridMG()95 GridMG* getNewGridMG()
96 {
97 GridMG* g=g_malloc(sizeof(GridMG));
98 g->domain.xSize = 0;
99 g->domain.ySize = 0;
100 g->domain.zSize = 0;
101 g->domain.size = 0;
102 g->operationType = GABEDIT_ALL;
103 g->values = NULL;
104 return g;
105 }
106 /*********************************************************/
getNewGridMGUsingDomain(DomainMG * domain)107 GridMG* getNewGridMGUsingDomain(DomainMG* domain)
108 {
109 glong i;
110 GridMG* g=g_malloc(sizeof(GridMG));
111 g->domain = *domain;
112 g->operationType = GABEDIT_ALL;
113 g->values = g_malloc(g->domain.size*sizeof(gdouble));
114 #ifdef ENABLE_OMP
115 #pragma omp parallel for private(i)
116 #endif
117 for(i = 0;i<g->domain.size;i++) g->values[i] = 0.0;
118 return g;
119 }
120 /*********************************************************/
getNewGridMGFromOldGrid(GridMG * src)121 GridMG* getNewGridMGFromOldGrid(GridMG* src)
122 {
123 GridMG* g=g_malloc(sizeof(GridMG));
124 g->domain = src->domain;
125 g->operationType = GABEDIT_ALL;
126 g->values = NULL;
127 equalAllGridMG(g,src);
128 return g;
129 }
130 /*********************************************************/
initAllGridMG(GridMG * g,gdouble value)131 void initAllGridMG(GridMG* g, gdouble value)
132 {
133 gint i;
134 #ifdef ENABLE_OMP
135 #pragma omp parallel for private(i)
136 #endif
137 for(i = 0;i<g->domain.size;i++) g->values[i] = value;
138 }
139 /*********************************************************/
initInteriorGridMG(GridMG * g,gdouble value)140 void initInteriorGridMG(GridMG* g, gdouble value)
141 {
142 DomainMG domain = g->domain;
143 if(domain.size<=0) return;
144
145 int ix;
146 int iy;
147 int iz;
148 int iXBegin = domain.iXBeginInterior;
149 int iXEnd = domain.iXEndInterior;
150 int iYBegin = domain.iYBeginInterior;
151 int iYEnd = domain.iYEndInterior;
152 int iZBegin = domain.iZBeginInterior;
153 int iZEnd = domain.iZEndInterior;
154
155 #ifdef ENABLE_OMP
156 #pragma omp parallel for private(ix,iy,iz)
157 #endif
158 for(ix = iXBegin;ix <=iXEnd;ix++)
159 for(iy = iYBegin;iy <=iYEnd;iy++)
160 for(iz = iZBegin;iz <=iZEnd;iz++)
161 setValGridMG(g,ix,iy,iz, value);
162 }
163 /*********************************************************/
initBoundaryGridMG(GridMG * g,gdouble value)164 void initBoundaryGridMG(GridMG* g, gdouble value)
165 {
166 int ix;
167 int iy;
168 int iz;
169 DomainMG domain = g->domain;
170
171 #ifdef ENABLE_OMP
172 #pragma omp parallel for private(ix,iy,iz)
173 #endif
174 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
175 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
176 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
177 setValGridMG(g,ix,iy,iz, value);
178 #ifdef ENABLE_OMP
179 #pragma omp parallel for private(ix,iy,iz)
180 #endif
181 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
182 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
183 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
184 setValGridMG(g,ix,iy,iz, value);
185
186 #ifdef ENABLE_OMP
187 #pragma omp parallel for private(ix,iy,iz)
188 #endif
189 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
190 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
191 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
192 setValGridMG(g,ix,iy,iz, value);
193
194 #ifdef ENABLE_OMP
195 #pragma omp parallel for private(ix,iy,iz)
196 #endif
197 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
198 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
199 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
200 setValGridMG(g,ix,iy,iz, value);
201 #ifdef ENABLE_OMP
202 #pragma omp parallel for private(ix,iy,iz)
203 #endif
204 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
205 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
206 {
207 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
208 setValGridMG(g,ix,iy,iz, value);
209 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
210 setValGridMG(g,ix,iy,iz, value);
211 }
212
213 }
214 /*********************************************************/
initGridMG(GridMG * g,gdouble sommeValue)215 void initGridMG(GridMG*g, gdouble sommeValue)
216 {
217 switch(g->operationType)
218 {
219 case GABEDIT_ALL: initAllGridMG(g, sommeValue);break;
220 case GABEDIT_INTERIOR: initInteriorGridMG(g, sommeValue);break;
221 case GABEDIT_BOUNDARY: initBoundaryGridMG(g, sommeValue);break;
222 }
223 }
224 /*********************************************************/
addGaussian(GridMG * g,gdouble Z,gdouble x1,gdouble y1,gdouble z1,gdouble sigma)225 void addGaussian(GridMG* g, gdouble Z, gdouble x1, gdouble y1, gdouble z1, gdouble sigma)
226 {
227 gdouble sigma2 = sigma*sigma;
228
229 DomainMG* domain = &g->domain;
230 GridMG* tmp = getNewGridMGUsingDomain(domain);
231
232 gdouble x0 = domain->x0;
233 gdouble y0 = domain->y0;
234 gdouble z0 = domain->z0;
235
236 gdouble xh = domain->xh;
237 gdouble yh = domain->yh;
238 gdouble zh = domain->zh;
239
240 gdouble r2x;
241 gdouble r2y;
242 gdouble r2z;
243
244 gdouble r20x=2*sigma2*xh*xh;
245 gdouble r20y=2*sigma2*yh*yh;
246 gdouble r20z=2*sigma2*zh*zh;
247
248 gdouble x, y , z;
249 int ix, iy, iz;
250 gdouble s = 0;
251 gdouble ex = 0;
252
253 for(ix=domain->iXBeginBoundaryLeft;ix<=domain->iXEndBoundaryRight;ix++)
254 {
255 x = x0 + ix*domain->xh;
256 r2x = (x-x1)*(x-x1);
257 for(iy = domain->iYBeginBoundaryLeft;iy <=domain->iYEndBoundaryRight;iy++)
258 {
259 y = y0 + iy*domain->yh;
260 r2y = (y-y1)*(y-y1);
261 for(iz = domain->iZBeginBoundaryLeft;iz <=domain->iZEndBoundaryRight;iz++)
262 {
263 z = z0 + iz*domain->zh;
264 r2z = (z-z1)*(z-z1);
265
266 ex = exp(-r2x/r20x)*exp(-r2y/r20y)*exp(-r2z/r20z);
267 setValGridMG(tmp, ix, iy, iz, ex);
268 s += ex;
269
270 }
271 }
272 }
273
274
275 s *= xh*yh*zh;
276
277 if(fabs(s)>PRECISION)
278 {
279 OperationTypeMG operation = getOperationGridMG(g);
280 setOperationGridMG(tmp, GABEDIT_INTERIOR);
281 setOperationGridMG(g, GABEDIT_INTERIOR);
282 multEqualRealGridMG(tmp,Z/s);
283 plusEqualGridMG( g, tmp);
284 setOperationGridMG(g, operation);
285 }
286 destroyGridMG(tmp);
287 }
288 /*********************************************************/
equalAllGridMG(GridMG * g,GridMG * src)289 static void equalAllGridMG(GridMG* g, GridMG* src)
290 {
291 if (g != src)
292 {
293 gint i;
294 destroyGridMG(g);
295 g->domain = src->domain;
296 g->operationType = src->operationType;
297 g->values = g_malloc(g->domain.size*sizeof(gdouble));
298 #ifdef ENABLE_OMP
299 #pragma omp parallel for private(i)
300 #endif
301 for(i = 0;i<g->domain.size;i++) g->values[i] = src->values[i];
302 }
303 }
304 /*********************************************************/
equalInteriorGridMG(GridMG * g,GridMG * src)305 static void equalInteriorGridMG(GridMG* g, GridMG* src)
306 {
307 DomainMG domain = src->domain;
308 if (g == src) return;
309 if(!ifEqualDomainMG(&g->domain,&src->domain))
310 {
311 destroyGridMG(g);
312 g->domain = src->domain;
313 g->values = g_malloc(domain.size*sizeof(gdouble));
314 }
315
316 g->operationType = src->operationType;
317
318 int ix;
319 int iy;
320 int iz;
321 int iXBegin = domain.iXBeginInterior;
322 int iXEnd = domain.iXEndInterior;
323 int iYBegin = domain.iYBeginInterior;
324 int iYEnd = domain.iYEndInterior;
325 int iZBegin = domain.iZBeginInterior;
326 int iZEnd = domain.iZEndInterior;
327
328 #ifdef ENABLE_OMP
329 #pragma omp parallel for private(ix,iy,iz)
330 #endif
331 for(ix = iXBegin ; ix <=iXEnd ; ix++)
332 for(iy = iYBegin ; iy <=iYEnd ; iy++)
333 {
334 for(iz = iZBegin ; iz <=iZEnd ; iz++)
335 setValGridMG(g,ix,iy,iz, getValGridMG(src,ix,iy,iz));
336 }
337 }
338 /*********************************************************/
equalBoundaryGridMG(GridMG * g,GridMG * src)339 static void equalBoundaryGridMG(GridMG* g, GridMG* src)
340 {
341 DomainMG domain = src->domain;
342 if (g == src) return;
343 if(!ifEqualDomainMG(&g->domain,&src->domain))return;
344
345 int ix;
346 int iy;
347 int iz;
348
349 #ifdef ENABLE_OMP
350 #pragma omp parallel for private(ix,iy,iz)
351 #endif
352 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
353 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
354 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
355 setValGridMG(g,ix,iy,iz, getValGridMG(src,ix,iy,iz));
356 #ifdef ENABLE_OMP
357 #pragma omp parallel for private(ix,iy,iz)
358 #endif
359 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
360 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
361 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
362 setValGridMG(g,ix,iy,iz, getValGridMG(src,ix,iy,iz));
363
364 #ifdef ENABLE_OMP
365 #pragma omp parallel for private(ix,iy,iz)
366 #endif
367 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
368 {
369 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
370 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
371 setValGridMG(g,ix,iy,iz, getValGridMG(src,ix,iy,iz));
372 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
373 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
374 setValGridMG(g,ix,iy,iz, getValGridMG(src,ix,iy,iz));
375 }
376 #ifdef ENABLE_OMP
377 #pragma omp parallel for private(ix,iy,iz)
378 #endif
379 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
380 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
381 {
382 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
383 setValGridMG(g,ix,iy,iz, getValGridMG(src,ix,iy,iz));
384 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
385 setValGridMG(g,ix,iy,iz, getValGridMG(src,ix,iy,iz));
386 }
387 }
388 /*********************************************************/
equalGridMG(GridMG * g,GridMG * src)389 void equalGridMG(GridMG*g, GridMG* src)
390 {
391 switch(g->operationType)
392 {
393 case GABEDIT_ALL: equalAllGridMG(g, src);break;
394 case GABEDIT_INTERIOR: equalInteriorGridMG(g, src);break;
395 case GABEDIT_BOUNDARY: equalBoundaryGridMG(g, src);break;
396 }
397 }
398 /*********************************************************/
copyGridMG(GridMG * g,GridMG * src)399 void copyGridMG(GridMG*g, GridMG* src)
400 {
401 equalGridMG(g,src);
402 }
403 /*********************************************************/
plusEqualAllGridMG(GridMG * g,GridMG * right)404 static void plusEqualAllGridMG(GridMG* g, GridMG* right)
405 {
406 gint i;
407 if(!ifEqualDomainMG(&g->domain, &right->domain)) return;
408 #ifdef ENABLE_OMP
409 #pragma omp parallel for private(i)
410 #endif
411 for(i = 0;i<g->domain.size;i++) g->values[i] += right->values[i];
412 }
413 /*********************************************************/
plusEqualInteriorGridMG(GridMG * g,GridMG * src)414 static void plusEqualInteriorGridMG(GridMG* g, GridMG* src)
415 {
416 DomainMG domain = g->domain;
417 int ix;
418 int iy;
419 int iz;
420 int iXBegin = domain.iXBeginInterior;
421 int iXEnd = domain.iXEndInterior;
422 int iYBegin = domain.iYBeginInterior;
423 int iYEnd = domain.iYEndInterior;
424 int iZBegin = domain.iZBeginInterior;
425 int iZEnd = domain.iZEndInterior;
426 if(!ifEqualDomainMG(&g->domain, &src->domain)) return;
427
428 #ifdef ENABLE_OMP
429 #pragma omp parallel for private(ix,iy,iz)
430 #endif
431 for(ix = iXBegin;ix <=iXEnd;ix++)
432 for(iy = iYBegin;iy <=iYEnd;iy++)
433 for(iz = iZBegin;iz <=iZEnd;iz++)
434 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)+getValGridMG(src,ix,iy,iz));
435 }
436 /*********************************************************/
plusEqualBoundaryGridMG(GridMG * g,GridMG * src)437 static void plusEqualBoundaryGridMG(GridMG* g, GridMG* src)
438 {
439 DomainMG domain = g->domain;
440 int ix;
441 int iy;
442 int iz;
443 if(!ifEqualDomainMG(&g->domain, &src->domain)) return;
444
445 #ifdef ENABLE_OMP
446 #pragma omp parallel for private(ix,iy,iz)
447 #endif
448 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
449 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
450 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
451 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)+getValGridMG(src,ix,iy,iz));
452 #ifdef ENABLE_OMP
453 #pragma omp parallel for private(ix,iy,iz)
454 #endif
455 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
456 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
457 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
458 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)+getValGridMG(src,ix,iy,iz));
459
460 #ifdef ENABLE_OMP
461 #pragma omp parallel for private(ix,iy,iz)
462 #endif
463 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
464 {
465 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
466 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
467 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)+getValGridMG(src,ix,iy,iz));
468 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
469 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
470 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)+getValGridMG(src,ix,iy,iz));
471 }
472 #ifdef ENABLE_OMP
473 #pragma omp parallel for private(ix,iy,iz)
474 #endif
475 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
476 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
477 {
478 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
479 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)+getValGridMG(src,ix,iy,iz));
480 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
481 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)+getValGridMG(src,ix,iy,iz));
482 }
483 }
484 /*********************************************************/
plusEqualGridMG(GridMG * g,GridMG * src)485 void plusEqualGridMG(GridMG*g, GridMG* src)
486 {
487 switch(g->operationType)
488 {
489 case GABEDIT_ALL: plusEqualAllGridMG(g, src);break;
490 case GABEDIT_INTERIOR: plusEqualInteriorGridMG(g, src);break;
491 case GABEDIT_BOUNDARY: plusEqualBoundaryGridMG(g, src);break;
492 }
493 }
494 /*********************************************************/
moinsEqualAllGridMG(GridMG * g,GridMG * src)495 static void moinsEqualAllGridMG(GridMG* g, GridMG* src)
496 {
497 gint i;
498 if(!ifEqualDomainMG(&g->domain, &src->domain)) return;
499 #ifdef ENABLE_OMP
500 #pragma omp parallel for private(i)
501 #endif
502 for(i = 0;i<g->domain.size;i++) g->values[i] -= src->values[i];
503 }
504 /*********************************************************/
moinsEqualInteriorGridMG(GridMG * g,GridMG * src)505 static void moinsEqualInteriorGridMG(GridMG* g, GridMG* src)
506 {
507 DomainMG domain = g->domain;
508 int ix;
509 int iy;
510 int iz;
511 int iXBegin = domain.iXBeginInterior;
512 int iXEnd = domain.iXEndInterior;
513 int iYBegin = domain.iYBeginInterior;
514 int iYEnd = domain.iYEndInterior;
515 int iZBegin = domain.iZBeginInterior;
516 int iZEnd = domain.iZEndInterior;
517 if(!ifEqualDomainMG(&g->domain, &src->domain)) return;
518
519 #ifdef ENABLE_OMP
520 #pragma omp parallel for private(ix,iy,iz)
521 #endif
522 for(ix = iXBegin;ix <=iXEnd;ix++)
523 for(iy = iYBegin;iy <=iYEnd;iy++)
524 for(iz = iZBegin;iz <=iZEnd;iz++)
525 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)-getValGridMG(src,ix,iy,iz));
526 }
527 /*********************************************************/
moinsEqualBoundaryGridMG(GridMG * g,GridMG * src)528 static void moinsEqualBoundaryGridMG(GridMG* g, GridMG* src)
529 {
530 DomainMG domain = g->domain;
531 int ix;
532 int iy;
533 int iz;
534 if(!ifEqualDomainMG(&g->domain, &src->domain)) return;
535
536 #ifdef ENABLE_OMP
537 #pragma omp parallel for private(ix,iy,iz)
538 #endif
539 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
540 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
541 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
542 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)-getValGridMG(src,ix,iy,iz));
543 #ifdef ENABLE_OMP
544 #pragma omp parallel for private(ix,iy,iz)
545 #endif
546 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
547 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
548 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
549 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)-getValGridMG(src,ix,iy,iz));
550
551 #ifdef ENABLE_OMP
552 #pragma omp parallel for private(ix,iy,iz)
553 #endif
554 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
555 {
556 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
557 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
558 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)-getValGridMG(src,ix,iy,iz));
559 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
560 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
561 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)-getValGridMG(src,ix,iy,iz));
562 }
563 #ifdef ENABLE_OMP
564 #pragma omp parallel for private(ix,iy,iz)
565 #endif
566 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
567 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
568 {
569 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
570 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)-getValGridMG(src,ix,iy,iz));
571 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
572 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)-getValGridMG(src,ix,iy,iz));
573 }
574 }
575 /*********************************************************/
moinsEqualGridMG(GridMG * g,GridMG * src)576 void moinsEqualGridMG(GridMG*g, GridMG* src)
577 {
578 switch(g->operationType)
579 {
580 case GABEDIT_ALL: moinsEqualAllGridMG(g, src);break;
581 case GABEDIT_INTERIOR: moinsEqualInteriorGridMG(g, src);break;
582 case GABEDIT_BOUNDARY: moinsEqualBoundaryGridMG(g, src);break;
583 }
584 }
585 /*********************************************************/
multEqualAllGridMG(GridMG * g,GridMG * src)586 static void multEqualAllGridMG(GridMG* g, GridMG* src)
587 {
588 gint i;
589 if(!ifEqualDomainMG(&g->domain, &src->domain)) return;
590 #ifdef ENABLE_OMP
591 #pragma omp parallel for private(i)
592 #endif
593 for(i= 0;i<g->domain.size;i++) g->values[i] *= src->values[i];
594 }
595 /*********************************************************/
multEqualInteriorGridMG(GridMG * g,GridMG * src)596 static void multEqualInteriorGridMG(GridMG* g, GridMG* src)
597 {
598 DomainMG domain = g->domain;
599 int ix;
600 int iy;
601 int iz;
602 int iXBegin = domain.iXBeginInterior;
603 int iXEnd = domain.iXEndInterior;
604 int iYBegin = domain.iYBeginInterior;
605 int iYEnd = domain.iYEndInterior;
606 int iZBegin = domain.iZBeginInterior;
607 int iZEnd = domain.iZEndInterior;
608 if(!ifEqualDomainMG(&g->domain, &src->domain)) return;
609
610 #ifdef ENABLE_OMP
611 #pragma omp parallel for private(ix,iy,iz)
612 #endif
613 for(ix = iXBegin;ix <=iXEnd;ix++)
614 for(iy = iYBegin;iy <=iYEnd;iy++)
615 for(iz = iZBegin;iz <=iZEnd;iz++)
616 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)*getValGridMG(src,ix,iy,iz));
617 }
618 /*********************************************************/
multEqualBoundaryGridMG(GridMG * g,GridMG * src)619 static void multEqualBoundaryGridMG(GridMG* g, GridMG* src)
620 {
621 DomainMG domain = g->domain;
622 int ix;
623 int iy;
624 int iz;
625 if(!ifEqualDomainMG(&g->domain, &src->domain)) return;
626
627 #ifdef ENABLE_OMP
628 #pragma omp parallel for private(ix,iy,iz)
629 #endif
630 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
631 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
632 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
633 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)*getValGridMG(src,ix,iy,iz));
634 #ifdef ENABLE_OMP
635 #pragma omp parallel for private(ix,iy,iz)
636 #endif
637 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
638 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
639 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
640 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)*getValGridMG(src,ix,iy,iz));
641
642 #ifdef ENABLE_OMP
643 #pragma omp parallel for private(ix,iy,iz)
644 #endif
645 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
646 {
647 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
648 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
649 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)*getValGridMG(src,ix,iy,iz));
650 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
651 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
652 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)*getValGridMG(src,ix,iy,iz));
653 }
654 #ifdef ENABLE_OMP
655 #pragma omp parallel for private(ix,iy,iz)
656 #endif
657 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
658 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
659 {
660 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
661 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)*getValGridMG(src,ix,iy,iz));
662 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
663 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)*getValGridMG(src,ix,iy,iz));
664 }
665 }
666 /*********************************************************/
multEqualGridMG(GridMG * g,GridMG * src)667 void multEqualGridMG(GridMG*g, GridMG* src)
668 {
669 switch(g->operationType)
670 {
671 case GABEDIT_ALL: multEqualAllGridMG(g, src);break;
672 case GABEDIT_INTERIOR: multEqualInteriorGridMG(g, src);break;
673 case GABEDIT_BOUNDARY: multEqualBoundaryGridMG(g, src);break;
674 }
675 }
676 /*********************************************************/
multEqualAllRealGridMG(GridMG * g,gdouble a)677 static void multEqualAllRealGridMG(GridMG* g, gdouble a)
678 {
679 gint i;
680 #ifdef ENABLE_OMP
681 #pragma omp parallel for private(i)
682 #endif
683 for(i = 0;i<g->domain.size;i++) g->values[i] *= a;
684 }
685 /*********************************************************/
multEqualInteriorRealGridMG(GridMG * g,gdouble a)686 static void multEqualInteriorRealGridMG(GridMG* g, gdouble a)
687 {
688 DomainMG domain = g->domain;
689 int ix;
690 int iy;
691 int iz;
692 int iXBegin = domain.iXBeginInterior;
693 int iXEnd = domain.iXEndInterior;
694 int iYBegin = domain.iYBeginInterior;
695 int iYEnd = domain.iYEndInterior;
696 int iZBegin = domain.iZBeginInterior;
697 int iZEnd = domain.iZEndInterior;
698
699 #ifdef ENABLE_OMP
700 #pragma omp parallel for private(ix,iy,iz)
701 #endif
702 for(ix = iXBegin;ix <=iXEnd;ix++)
703 for(iy = iYBegin;iy <=iYEnd;iy++)
704 for(iz = iZBegin;iz <=iZEnd;iz++)
705 multValGridMG(g,ix,iy,iz,a);
706 }
707 /*********************************************************/
multEqualBoundaryRealGridMG(GridMG * g,gdouble a)708 static void multEqualBoundaryRealGridMG(GridMG* g, gdouble a)
709 {
710 int ix;
711 int iy;
712 int iz;
713 DomainMG domain = g->domain;
714
715 #ifdef ENABLE_OMP
716 #pragma omp parallel for private(ix,iy,iz)
717 #endif
718 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
719 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
720 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
721 multValGridMG(g,ix,iy,iz,a);
722 #ifdef ENABLE_OMP
723 #pragma omp parallel for private(ix,iy,iz)
724 #endif
725 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
726 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
727 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
728 multValGridMG(g,ix,iy,iz,a);
729
730 #ifdef ENABLE_OMP
731 #pragma omp parallel for private(ix,iy,iz)
732 #endif
733 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
734 {
735 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
736 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
737 multValGridMG(g,ix,iy,iz,a);
738 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
739 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
740 multValGridMG(g,ix,iy,iz,a);
741 }
742 #ifdef ENABLE_OMP
743 #pragma omp parallel for private(ix,iy,iz)
744 #endif
745 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
746 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
747 {
748 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
749 multValGridMG(g,ix,iy,iz,a);
750 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
751 multValGridMG(g,ix,iy,iz,a);
752 }
753 }
754 /*********************************************************/
multEqualRealGridMG(GridMG * g,gdouble a)755 void multEqualRealGridMG(GridMG*g, gdouble a)
756 {
757 switch(g->operationType)
758 {
759 case GABEDIT_ALL: multEqualAllRealGridMG(g, a);break;
760 case GABEDIT_INTERIOR: multEqualInteriorRealGridMG(g, a);break;
761 case GABEDIT_BOUNDARY: multEqualBoundaryRealGridMG(g, a);break;
762 }
763 }
764 /*********************************************************/
divEqualAllRealGridMG(GridMG * g,gdouble a)765 static void divEqualAllRealGridMG(GridMG* g, gdouble a)
766 {
767 gint i;
768 for(i = 0;i<g->domain.size;i++) g->values[i] /= a;
769 }
770 /*********************************************************/
divEqualInteriorRealGridMG(GridMG * g,gdouble a)771 static void divEqualInteriorRealGridMG(GridMG* g, gdouble a)
772 {
773
774 DomainMG domain = g->domain;
775 int ix;
776 int iy;
777 int iz;
778 int iXBegin = domain.iXBeginInterior;
779 int iXEnd = domain.iXEndInterior;
780 int iYBegin = domain.iYBeginInterior;
781 int iYEnd = domain.iYEndInterior;
782 int iZBegin = domain.iZBeginInterior;
783 int iZEnd = domain.iZEndInterior;
784 if(a==0) return ;
785
786 #ifdef ENABLE_OMP
787 #pragma omp parallel for private(ix,iy,iz)
788 #endif
789 for(ix = iXBegin;ix <=iXEnd;ix++)
790 for(iy = iYBegin;iy <=iYEnd;iy++)
791 for(iz = iZBegin;iz <=iZEnd;iz++)
792 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)/a);
793 }
794 /*********************************************************/
divEqualBoundaryRealGridMG(GridMG * g,gdouble a)795 static void divEqualBoundaryRealGridMG(GridMG* g, gdouble a)
796 {
797 int ix;
798 int iy;
799 int iz;
800 DomainMG domain = g->domain;
801 if(a==0) return;
802
803 #ifdef ENABLE_OMP
804 #pragma omp parallel for private(ix,iy,iz)
805 #endif
806 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
807 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
808 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
809 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)/a);
810 #ifdef ENABLE_OMP
811 #pragma omp parallel for private(ix,iy,iz)
812 #endif
813 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
814 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
815 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
816 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)/a);
817
818 #ifdef ENABLE_OMP
819 #pragma omp parallel for private(ix,iy,iz)
820 #endif
821 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
822 {
823 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
824 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
825 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)/a);
826 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
827 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
828 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)/a);
829 }
830 #ifdef ENABLE_OMP
831 #pragma omp parallel for private(ix,iy,iz)
832 #endif
833 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
834 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
835 {
836 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
837 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)/a);
838 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
839 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)/a);
840 }
841 }
842 /*********************************************************/
divEqualRealGridMG(GridMG * g,gdouble a)843 void divEqualRealGridMG(GridMG*g, gdouble a)
844 {
845 switch(g->operationType)
846 {
847 case GABEDIT_ALL: divEqualAllRealGridMG(g, a);break;
848 case GABEDIT_INTERIOR: divEqualInteriorRealGridMG(g, a);break;
849 case GABEDIT_BOUNDARY: divEqualBoundaryRealGridMG(g, a);break;
850 }
851 }
852 /*********************************************************/
getDiagGridMG(GridMG * g)853 gdouble getDiagGridMG(GridMG* g)
854 {
855 return g->domain.diag;
856 }
857 /*********************************************************/
laplacianGridMG(GridMG * g,GridMG * src)858 gdouble laplacianGridMG(GridMG* g, GridMG* src)
859 {
860 int ix, iy, iz;
861 DomainMG domain = g->domain;
862
863
864 gdouble cc = domain.cc;
865 gdouble diag = domain.diag;
866 gdouble* fcx = domain.fLaplacinaX;
867 gdouble* fcy = domain.fLaplacinaY;
868 gdouble* fcz = domain.fLaplacinaZ;
869 int i;
870
871 int iXBegin = domain.iXBeginInterior;
872 int iXEnd = domain.iXEndInterior;
873 int iYBegin = domain.iYBeginInterior;
874 int iYEnd = domain.iYEndInterior;
875 int iZBegin = domain.iZBeginInterior;
876 int iZEnd = domain.iZEndInterior;
877 int nBoundary = domain.nBoundary;
878 gdouble v;
879
880 if(!ifEqualDomainMG(&g->domain,&src->domain))
881 {
882 destroyGridMG(g);
883 g->domain = src->domain;
884 g->operationType = src->operationType;
885 g->values = g_malloc(domain.size*sizeof(gdouble));
886 }
887 initBoundaryGridMG(g,0.0);
888
889 #ifdef ENABLE_OMP
890 #pragma omp parallel for private(ix,iy,iz,v)
891 #endif
892 for(ix = iXBegin;ix <= iXEnd;ix++)
893 for(iy = iYBegin;iy <= iYEnd;iy++)
894 for(iz = iZBegin;iz <= iZEnd;iz++)
895 {
896 {
897 v = cc * getValGridMG(src, ix,iy,iz);
898 for(i=1;i<=nBoundary;i++)
899 {
900 v += fcx[i] *(getValGridMG(src, ix-i,iy,iz)+getValGridMG(src, ix+i,iy,iz));
901 v += fcy[i] *(getValGridMG(src, ix,iy-i,iz)+getValGridMG(src, ix,iy+i,iz));
902 v += fcz[i] *(getValGridMG(src, ix,iy,iz-i)+getValGridMG(src, ix,iy,iz+i));
903 }
904 setValGridMG(g,ix,iy,iz, v);
905 }
906 }
907 return diag;
908 }
909 /*********************************************************/
plusLaplacianGridMG(GridMG * g,GridMG * src)910 gdouble plusLaplacianGridMG(GridMG* g, GridMG* src)
911 {
912 int ix, iy, iz;
913
914 DomainMG domain = g->domain;
915 gdouble cc = domain.cc;
916 gdouble diag = domain.diag;
917 gdouble* fcx = domain.fLaplacinaX;
918 gdouble* fcy = domain.fLaplacinaY;
919 gdouble* fcz = domain.fLaplacinaZ;
920 int i;
921
922 int iXBegin = domain.iXBeginInterior;
923 int iXEnd = domain.iXEndInterior;
924 int iYBegin = domain.iYBeginInterior;
925 int iYEnd = domain.iYEndInterior;
926 int iZBegin = domain.iZBeginInterior;
927 int iZEnd = domain.iZEndInterior;
928 int nBoundary = domain.nBoundary;
929 gdouble v;
930
931 if(!ifEqualDomainMG(&g->domain,&src->domain))
932 {
933 destroyGridMG(g);
934 g->domain = src->domain;
935 g->operationType = src->operationType;
936 g->values = g_malloc(domain.size*sizeof(gdouble));
937 initAllGridMG(g, 0.0);
938 }
939 else
940 initBoundaryGridMG(g, 0.0);
941
942 #ifdef ENABLE_OMP
943 #pragma omp parallel for private(ix,iy,iz,v)
944 #endif
945 for(ix = iXBegin;ix <= iXEnd;ix++)
946 for(iy = iYBegin;iy <= iYEnd;iy++)
947 for(iz = iZBegin;iz <= iZEnd;iz++)
948 {
949 v = cc * getValGridMG(src, ix,iy,iz);
950
951 for(i=1;i<=nBoundary;i++)
952 {
953 v += fcx[i] *(getValGridMG(src, ix-i,iy,iz)+getValGridMG(src, ix+i,iy,iz));
954 v += fcy[i] *(getValGridMG(src, ix,iy-i,iz)+getValGridMG(src, ix,iy+i,iz));
955 v += fcz[i] *(getValGridMG(src, ix,iy,iz-i)+getValGridMG(src, ix,iy,iz+i));
956 }
957 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)+v);
958 }
959
960 return diag;
961 }
962 /*********************************************************/
moinsLaplacianGridMG(GridMG * g,GridMG * src)963 gdouble moinsLaplacianGridMG(GridMG* g, GridMG* src)
964 {
965 int ix, iy, iz;
966
967 DomainMG domain = g->domain;
968 gdouble cc = domain.cc;
969 gdouble diag = domain.diag;
970 gdouble* fcx = domain.fLaplacinaX;
971 gdouble* fcy = domain.fLaplacinaY;
972 gdouble* fcz = domain.fLaplacinaZ;
973 int i;
974
975 int iXBegin = domain.iXBeginInterior;
976 int iXEnd = domain.iXEndInterior;
977 int iYBegin = domain.iYBeginInterior;
978 int iYEnd = domain.iYEndInterior;
979 int iZBegin = domain.iZBeginInterior;
980 int iZEnd = domain.iZEndInterior;
981 int nBoundary = domain.nBoundary;
982
983 if(!ifEqualDomainMG(&g->domain,&src->domain))
984 {
985 destroyGridMG(g);
986 g->domain = src->domain;
987 g->operationType = src->operationType;
988 g->values = g_malloc(domain.size*sizeof(gdouble));
989 initAllGridMG(g, 0.0);
990 }
991 else
992 initBoundaryGridMG(g, 0.0);
993
994 #ifdef ENABLE_OMP
995 #pragma omp parallel for private(ix,iy,iz)
996 #endif
997 for(ix = iXBegin;ix <= iXEnd;ix++)
998 for(iy = iYBegin;iy <= iYEnd;iy++)
999 for(iz = iZBegin;iz <= iZEnd;iz++)
1000 {
1001 gdouble v;
1002 gdouble vx, vy, vz;
1003 v = cc * getValGridMG(src, ix,iy,iz);
1004 vx = 0.0;
1005 vy = 0.0;
1006 vz = 0.0;
1007
1008 for(i=1;i<=nBoundary;i++)
1009 {
1010 vx += fcx[i] *(getValGridMG(src, ix-i,iy,iz)+getValGridMG(src, ix+i,iy,iz));
1011 vy += fcy[i] *(getValGridMG(src, ix,iy-i,iz)+getValGridMG(src, ix,iy+i,iz));
1012 vz += fcz[i] *(getValGridMG(src, ix,iy,iz-i)+getValGridMG(src, ix,iy,iz+i));
1013 }
1014 setValGridMG(g,ix,iy,iz, getValGridMG(g,ix,iy,iz)-(v + vx + vy + vz));
1015 }
1016
1017 return diag;
1018 }
1019 /*********************************************************/
averageGridMG(GridMG * g)1020 void averageGridMG(GridMG* g)
1021 {
1022 int ix, iy, iz;
1023 int x0, xp, xm, y0, yp, ym, z0, zp, zm;
1024 gdouble face, corner, edge;
1025 static gdouble scale = 1.0 / 64.0;
1026 DomainMG domain = g->domain;
1027 int iXBegin = domain.iXBeginInterior;
1028 int iXEnd = domain.iXEndInterior;
1029 int iYBegin = domain.iYBeginInterior;
1030 int iYEnd = domain.iYEndInterior;
1031 int iZBegin = domain.iZBeginInterior;
1032 int iZEnd = domain.iZEndInterior;
1033
1034 GridMG* src = getNewGridMGFromOldGrid(g);
1035
1036 #ifdef ENABLE_OMP
1037 #pragma omp parallel for private(ix,iy,iz,x0,xp,xm,y0,yp,ym,z0,zp,zm,face,corner,edge)
1038 #endif
1039 for(ix = iXBegin ; ix <= iXEnd ; ix++)
1040 {
1041 x0 = ix;
1042 xp = x0 + 1;
1043 xm = x0 - 1;
1044 for(iy = iYBegin ; iy <= iYEnd ; iy++)
1045 {
1046 y0 = iy;
1047 yp = y0 + 1;
1048 ym = y0 - 1;
1049 for(iz = iZBegin ; iz <= iZEnd ; iz++)
1050 {
1051 z0 = iz;
1052 zp = z0 + 1;
1053 zm = z0 - 1;
1054
1055 face =
1056 getValGridMG(src, xm , y0 , z0) +
1057 getValGridMG(src, xp , y0 , z0) +
1058 getValGridMG(src, x0 , ym , z0) +
1059 getValGridMG(src, x0 , yp , z0) +
1060 getValGridMG(src, x0 , y0 , zm) +
1061 getValGridMG(src, x0 , y0 , zp);
1062
1063 corner =
1064 getValGridMG(src, xm , ym , zm) +
1065 getValGridMG(src, xm , ym , zp) +
1066 getValGridMG(src, xm , yp , zm) +
1067 getValGridMG(src, xm , yp , zp) +
1068 getValGridMG(src, xp , ym , zm) +
1069 getValGridMG(src, xp , ym , zp) +
1070 getValGridMG(src, xp , yp , zm) +
1071 getValGridMG(src, xp , yp , zp);
1072
1073 edge =
1074 getValGridMG(src, xm , y0 , zm) +
1075 getValGridMG(src, xm , ym , z0) +
1076 getValGridMG(src, xm , yp , z0) +
1077 getValGridMG(src, xm , y0 , zp) +
1078 getValGridMG(src, x0 , ym , zm) +
1079 getValGridMG(src, x0 , yp , zm) +
1080 getValGridMG(src, x0 , ym , zp) +
1081 getValGridMG(src, x0 , yp , zp) +
1082 getValGridMG(src, xp , y0 , zm) +
1083 getValGridMG(src, xp , ym , z0) +
1084 getValGridMG(src, xp , yp , z0) +
1085 getValGridMG(src, xp , y0 , zp);
1086
1087
1088 setValGridMG(g,ix,iy,iz,
1089 scale * (
1090 8.0 * getValGridMG(src, x0 , y0 , z0) +
1091 4.0 * face +
1092 2.0 * edge +
1093 corner
1094 )
1095 );
1096 }
1097 }
1098 }
1099 destroyGridMG(src);
1100 }
1101 /*********************************************************/
resetLaplacianOrderGridMG(GridMG * g,LaplacianOrderMG order)1102 void resetLaplacianOrderGridMG(GridMG* g, LaplacianOrderMG order)
1103 {
1104 DomainMG domain = g->domain;
1105 DomainMG newDomain = getDomainMG(domain.xSize,domain.ySize,domain.xSize,
1106 domain.x0, domain.y0, domain.z0,
1107 domain.xLength, domain.yLength, domain.zLength,
1108 order);
1109
1110 int ix, iy, iz;
1111 int ixNew, iyNew, izNew;
1112
1113
1114 GridMG* newGrid = getNewGridMGUsingDomain(&newDomain);
1115
1116 int iXBegin = domain.iXBeginInterior;
1117 int iXEnd = domain.iXEndInterior;
1118 int iYBegin = domain.iYBeginInterior;
1119 int iYEnd = domain.iYEndInterior;
1120 int iZBegin = domain.iZBeginInterior;
1121 int iZEnd = domain.iZEndInterior;
1122
1123 int iXBeginNew = newDomain.iXBeginInterior;
1124 int iYBeginNew = newDomain.iYBeginInterior;
1125 int iZBeginNew = newDomain.iZBeginInterior;
1126
1127 #ifdef ENABLE_OMP
1128 #pragma omp parallel for private(ix,iy,iz,ixNew,iyNew,izNew)
1129 #endif
1130 for(ix = iXBegin ; ix <= iXEnd ; ix++)
1131 {
1132 ixNew = iXBeginNew+ix-iXBegin;
1133 iyNew = iYBeginNew ;
1134 for(iy = iYBegin ; iy <= iYEnd ; iy++, iyNew++)
1135 {
1136 izNew = iZBeginNew ;
1137 for(iz = iZBegin ; iz <= iZEnd ; iz++, izNew++)
1138 setValGridMG(newGrid, ixNew , iyNew , izNew, getValGridMG(g, ix , iy , iz));
1139 }
1140 }
1141
1142 equalAllGridMG(g,newGrid);
1143 destroyGridMG(newGrid);
1144 }
1145 /*********************************************************/
reAllocValuesTableGridMG(GridMG * g)1146 void reAllocValuesTableGridMG(GridMG* g)
1147 {
1148 glong i;
1149 if(g->domain.size <1) return;
1150 if(g->values) g_free(g->values);
1151 g->values = g_malloc(g->domain.size*sizeof(gdouble));
1152 #ifdef ENABLE_OMP
1153 #pragma omp parallel for private(i)
1154 #endif
1155 for(i = 0;i<g->domain.size;i++) g->values[i] = 0.0;
1156 }
1157 /*********************************************************/
levelUpGridMG(GridMG * g)1158 void levelUpGridMG(GridMG* g)
1159 {
1160 levelUpDomainMG(&g->domain);
1161 reAllocValuesTableGridMG(g);
1162 }
1163 /*********************************************************/
interpolationCubicSrcGridMG(GridMG * g,GridMG * src)1164 void interpolationCubicSrcGridMG(GridMG* g, GridMG* src)
1165 {
1166 gint ix, iy, iz;
1167 DomainMG domain;
1168
1169 g->domain = src->domain;
1170 levelUpGridMG(g);
1171
1172 domain = g->domain;
1173
1174 /*
1175 * transfer coarse grid pogints to
1176 * fine grid along with the
1177 * high side image pogint
1178 */
1179 /*printf("Je suis dans prolongation de Grid\n");*/
1180 /*printf("xSize =%d\n",xSize);*/
1181
1182 #ifdef ENABLE_OMP
1183 #pragma omp parallel for private(ix,iy,iz)
1184 #endif
1185 for(ix = src->domain.iXBeginInterior-1;ix <=src->domain.iXEndInterior+1;ix++)
1186 for(iy = src->domain.iYBeginInterior-1;iy <=src->domain.iYEndInterior+1;iy++)
1187 for(iz = src->domain.iZBeginInterior-1;iz <=src->domain.iZEndInterior+1;iz++)
1188 setValGridMG(g,2*ix,2*iy,2*iz, getValGridMG(src,ix,iy,iz));
1189
1190 /* ginterior center pogints */
1191 #ifdef ENABLE_OMP
1192 #pragma omp parallel for private(ix,iy,iz)
1193 #endif
1194 for(ix = domain.iXBeginInterior;ix <=domain.iXEndInterior;ix += 2)
1195 for(iy = domain.iYBeginInterior;iy <=domain.iYEndInterior;iy += 2)
1196 for(iz = domain.iZBeginInterior;iz <=domain.iZEndInterior;iz += 2)
1197 {
1198
1199 setValGridMG(g,ix,iy,iz,
1200 0.125 * getValGridMG(g,ix-1 , iy-1 , iz-1) +
1201 0.125 * getValGridMG(g, ix-1 , iy-1 , iz+1) +
1202 0.125 * getValGridMG(g, ix-1 , iy+1 , iz-1) +
1203 0.125 * getValGridMG(g, ix-1 , iy+1 , iz+1) +
1204 0.125 * getValGridMG(g, ix+1 , iy-1 , iz-1) +
1205 0.125 * getValGridMG(g, ix+1 , iy-1 , iz+1) +
1206 0.125 * getValGridMG(g, ix+1 , iy+1 , iz-1) +
1207 0.125 * getValGridMG(g, ix+1 , iy+1 , iz+1)
1208 );
1209 }
1210
1211 #ifdef ENABLE_OMP
1212 #pragma omp parallel for private(ix,iy,iz)
1213 #endif
1214 for(ix = domain.iXBeginInterior;ix <=domain.iXEndInterior;ix += 2)
1215 for(iy = domain.iYBeginInterior;iy <=domain.iYEndInterior;iy += 2)
1216 for(iz = domain.iZBeginInterior+1;iz <=domain.iZEndInterior;iz += 2)
1217 {
1218
1219 setValGridMG(g,ix,iy,iz,
1220 0.5 * getValGridMG(g, ix , iy , iz-1) +
1221 0.5 * getValGridMG(g, ix , iy , iz+1)
1222 );
1223 }
1224 #ifdef ENABLE_OMP
1225 #pragma omp parallel for private(ix,iy,iz)
1226 #endif
1227 for(ix = domain.iXBeginInterior;ix <=domain.iXEndInterior;ix += 2)
1228 for(iy = domain.iYBeginInterior+1;iy <=domain.iYEndInterior;iy += 2)
1229 for(iz = domain.iZBeginInterior;iz <=domain.iZEndInterior;iz += 2)
1230 {
1231 setValGridMG(g,ix,iy,iz,
1232 0.5 * getValGridMG(g, ix , (iy-1) , iz) +
1233 0.5 * getValGridMG(g, ix , (iy+1) , iz)
1234 );
1235 }
1236 #ifdef ENABLE_OMP
1237 #pragma omp parallel for private(ix,iy,iz)
1238 #endif
1239 for(ix = domain.iXBeginInterior+1;ix <=domain.iXEndInterior;ix += 2)
1240 for(iy = domain.iYBeginInterior;iy <=domain.iYEndInterior;iy += 2)
1241 for(iz = domain.iZBeginInterior;iz <=domain.iZEndInterior;iz += 2)
1242 {
1243 setValGridMG(g,ix,iy,iz,
1244 0.5 * getValGridMG(g, (ix-1) , iy , iz) +
1245 0.5 * getValGridMG(g, (ix+1) , iy , iz)
1246 );
1247 }
1248 #ifdef ENABLE_OMP
1249 #pragma omp parallel for private(ix,iy,iz)
1250 #endif
1251 for(ix = domain.iXBeginInterior;ix <=domain.iXEndInterior;ix += 2)
1252 for(iy = domain.iYBeginInterior+1;iy <=domain.iYEndInterior;iy += 2)
1253 for(iz = domain.iZBeginInterior+1;iz <=domain.iZEndInterior;iz += 2)
1254 {
1255 setValGridMG(g,ix,iy,iz,
1256 0.25 * getValGridMG(g, ix , (iy-1) , iz-1) +
1257 0.25 * getValGridMG(g, ix , (iy-1) , iz+1) +
1258 0.25 * getValGridMG(g, ix , (iy+1) , iz-1) +
1259 0.25 * getValGridMG(g, ix , (iy+1) , iz+1)
1260 );
1261 }
1262 #ifdef ENABLE_OMP
1263 #pragma omp parallel for private(ix,iy,iz)
1264 #endif
1265 for(ix = domain.iXBeginInterior+1;ix <=domain.iXEndInterior;ix += 2)
1266 for(iy = domain.iYBeginInterior;iy <=domain.iYEndInterior;iy += 2)
1267 for(iz = domain.iZBeginInterior+1;iz <=domain.iZEndInterior;iz += 2)
1268 {
1269 setValGridMG(g,ix,iy,iz,
1270 0.25 * getValGridMG(g, (ix-1) , iy , iz-1) +
1271 0.25 * getValGridMG(g, (ix-1) , iy , iz+1) +
1272 0.25 * getValGridMG(g, (ix+1) , iy , iz-1) +
1273 0.25 * getValGridMG(g, (ix+1) , iy , iz+1)
1274 );
1275 }
1276 #ifdef ENABLE_OMP
1277 #pragma omp parallel for private(ix,iy,iz)
1278 #endif
1279 for(ix = domain.iXBeginInterior+1;ix <=domain.iXEndInterior;ix += 2)
1280 for(iy = domain.iYBeginInterior+1;iy <=domain.iYEndInterior;iy += 2)
1281 for(iz = domain.iZBeginInterior;iz <=domain.iZEndInterior;iz += 2)
1282 {
1283 setValGridMG(g,ix,iy,iz,
1284 0.25 * getValGridMG(g, (ix-1) , (iy-1) , iz) +
1285 0.25 * getValGridMG(g, (ix+1) , (iy-1) , iz) +
1286 0.25 * getValGridMG(g, (ix-1) , (iy+1) , iz) +
1287 0.25 * getValGridMG(g, (ix+1) , (iy+1) , iz)
1288 );
1289 }
1290 }
1291 /*********************************************************/
interpolationCubicGridMG(GridMG * g)1292 void interpolationCubicGridMG(GridMG* g)
1293 {
1294 GridMG* newGrid = getNewGridMGFromOldGrid(g);
1295 interpolationCubicSrcGridMG(g, newGrid);
1296 destroyGridMG(newGrid);
1297 }
1298 /*********************************************************/
interpolationTriLinearSrcGridMG(GridMG * g,GridMG * src)1299 void interpolationTriLinearSrcGridMG(GridMG* g, GridMG* src)
1300 {
1301 gint ix, iy, iz;
1302 gdouble a1, a2, a3, a4;
1303
1304 gint iXBegin = src->domain.iXBeginInterior - 1;
1305 gint iYBegin = src->domain.iYBeginInterior - 1;
1306 gint iZBegin = src->domain.iZBeginInterior - 1;
1307
1308 gint iXEnd = src->domain.iXEndInterior + 1;
1309 gint iYEnd = src->domain.iYEndInterior + 1;
1310 gint iZEnd = src->domain.iZEndInterior + 1;
1311
1312 g->domain = src->domain;
1313 levelUpGridMG(g);
1314
1315 initAllGridMG(g,0.0);
1316
1317 addValGridMG(g, iXBegin,iYBegin,iZBegin, getValGridMG(src, iXBegin,iYBegin,iZBegin) );
1318 /* Interpolation of the first xy-plane where z = 0 */
1319 #ifdef ENABLE_OMP
1320 #pragma omp parallel for private(ix)
1321 #endif
1322 for(ix = iXBegin+1 ; ix <=iXEnd ; ix++)
1323 {
1324 addValGridMG(g, 2*ix-1, iYBegin, iZBegin, 0.5*( getValGridMG(src, ix-1,iYBegin, iZBegin) + getValGridMG(src, ix, iYBegin, iZBegin) ) );
1325 addValGridMG(g, 2*ix, iYBegin, iZBegin , getValGridMG(src, ix, iYBegin, iZBegin));
1326 }
1327 #ifdef ENABLE_OMP
1328 #pragma omp parallel for private(ix,iy,a1,a2)
1329 #endif
1330 for(iy = iYBegin+1 ; iy <=iYEnd ; iy++)
1331 {
1332 addValGridMG(g, iXBegin, 2*iy-1, iZBegin, 0.5*( getValGridMG(src, iXBegin, iy-1, iZBegin) + getValGridMG(src, iXBegin, iy, iZBegin) ) );
1333 addValGridMG(g, iXBegin, 2*iy, iZBegin, getValGridMG(src, iXBegin,iy, iZBegin));
1334 for(ix = iXBegin+1 ; ix <=iXEnd ; ix++)
1335 {
1336 a1 = 0.5*( getValGridMG(src, ix, iy-1, iZBegin) + getValGridMG(src, ix, iy, iZBegin) );
1337 a2 = 0.5*( getValGridMG(src, ix-1, iy-1, iZBegin) + getValGridMG(src, ix-1, iy, iZBegin) );
1338 addValGridMG(g, 2*ix-1, 2*iy-1, iZBegin, 0.5 * ( a1 + a2));
1339 addValGridMG(g, 2*ix, 2*iy-1, iZBegin, a1);
1340 addValGridMG(g, 2*ix-1, 2*iy, iZBegin,
1341 0.5*( getValGridMG(src, ix-1, iy, iZBegin) + getValGridMG(src, ix, iy, iZBegin)));
1342 addValGridMG(g, 2*ix, 2*iy, iZBegin, getValGridMG(src, ix, iy, iZBegin));
1343 }
1344 }
1345 /* Interpolation of other xy-plane where 0<z<xSize */
1346 #ifdef ENABLE_OMP
1347 #pragma omp parallel for private(ix,iy,iz,a1,a2)
1348 #endif
1349 for(iz = iZBegin+1 ; iz <= iZEnd ; iz++)
1350 {
1351 /* Interpolation on even planes */
1352 addValGridMG(g, iXBegin, iYBegin, 2*iz, getValGridMG(src, iXBegin, iYBegin, iz));
1353 for(ix = iXBegin+1 ; ix <=iXEnd ; ix++)
1354 {
1355 addValGridMG(g, 2*ix-1, iYBegin, 2*iz,
1356 0.5*( getValGridMG(src, ix-1, iYBegin, iz) + getValGridMG(src, ix, iYBegin, iz) ) );
1357 addValGridMG(g, 2*ix, iYBegin, 2*iz, getValGridMG(src, ix, iYBegin, iz));
1358 }
1359 for(iy = iYBegin+1 ; iy <=iYEnd ; iy++)
1360 {
1361 addValGridMG(g, iXBegin, 2*iy-1, 2*iz,
1362 0.5*( getValGridMG(src, iXBegin, iy-1, iz) + getValGridMG(src, iXBegin, iy, iz) ) );
1363 addValGridMG(g, iXBegin, 2*iy, 2*iz, getValGridMG(src, iXBegin, iy, iz));
1364 for(ix = iXBegin+1 ; ix <=iXEnd ; ix++)
1365 {
1366 a1 = 0.5*( getValGridMG(src, ix, iy-1, iz) + getValGridMG(src, ix, iy, iz) );
1367 a2 = 0.5*( getValGridMG(src, ix-1, iy-1, iz) + getValGridMG(src, ix-1, iy, iz) );
1368 addValGridMG(g, 2*ix-1, 2*iy-1, 2*iz, 0.5*( a1 + a2 ));
1369 addValGridMG(g, 2*ix, 2*iy-1, 2*iz, a1);
1370 addValGridMG(g, 2*ix-1, 2*iy, 2*iz, 0.5*( getValGridMG(src, ix-1,iy,iz) + getValGridMG(src, ix,iy,iz) ));
1371 addValGridMG(g, 2*ix, 2*iy, 2*iz, getValGridMG(src, ix, iy, iz));
1372 }
1373 }
1374 /* Interpolation on odd planes */
1375 addValGridMG(g, iXBegin, iYBegin, 2*iz-1,
1376 0.5*( getValGridMG(src, iXBegin, iYBegin, iz-1) + getValGridMG(src, iXBegin, iYBegin, iz) ));
1377 for(ix = iXBegin+1 ; ix <=iXEnd ; ix++)
1378 {
1379 addValGridMG(g, 2*ix-1, iYBegin, 2*iz-1, 0.25*(
1380 getValGridMG(src, ix-1, iYBegin, iz-1) + getValGridMG(src, ix, iYBegin, iz-1) +
1381 getValGridMG(src, ix-1, iYBegin, iz) + getValGridMG(src, ix, iYBegin, iz)
1382 ) );
1383 addValGridMG(g, 2*ix, iYBegin, 2*iz-1,
1384 0.5* (getValGridMG(src, ix, iYBegin, iz-1) + getValGridMG(src, ix, iYBegin, iz)));
1385 }
1386 for(iy = iYBegin+1 ; iy <=iYEnd ; iy++)
1387 {
1388 addValGridMG(g, iXBegin, 2*iy-1, 2*iz-1, 0.25*(
1389 getValGridMG(src, iXBegin, iy-1, iz) + getValGridMG(src, iXBegin, iy, iz)+
1390 getValGridMG(src, iXBegin, iy-1, iz-1) + getValGridMG(src, iXBegin, iy, iz-1)
1391 ) );
1392 addValGridMG(g, iXBegin, 2*iy, 2*iz-1,
1393 0.5* ( getValGridMG(src, iXBegin, iy, iz-1) + getValGridMG(src, iXBegin, iy, iz)));
1394
1395 for(ix = iXBegin+1 ; ix <=iXEnd ; ix++)
1396 {
1397 a1 = 0.5*( getValGridMG(src, ix-1, iy, iz) + getValGridMG(src, ix-1, iy, iz-1) );
1398 a2 = 0.5*( getValGridMG(src, ix, iy, iz) + getValGridMG(src, ix, iy, iz-1) );
1399 a3 = 0.5*( getValGridMG(src, ix, iy-1, iz) + getValGridMG(src, ix, iy-1, iz-1) );
1400 a4 = 0.5*( getValGridMG(src, ix-1, iy-1, iz) + getValGridMG(src, ix-1, iy-1, iz-1) );
1401 addValGridMG(g, 2*ix, 2*iy, 2*iz-1, a2);
1402 addValGridMG(g, 2*ix-1, 2*iy, 2*iz-1, 0.5*( a1 + a2 ));
1403 addValGridMG(g, 2*ix, 2*iy-1, 2*iz-1, 0.5*( a2 + a3 ));
1404 addValGridMG(g, 2*ix-1, 2*iy-1, 2*iz-1, 0.25*( a1 + a2 + a3 + a4 ));
1405 }
1406 }
1407 }
1408 }
1409 /*********************************************************/
interpolationTriLinearGridMG(GridMG * g)1410 void interpolationTriLinearGridMG(GridMG* g)
1411 {
1412 GridMG* newGrid = getNewGridMGFromOldGrid(g);
1413 interpolationTriLinearSrcGridMG(g, newGrid);
1414 destroyGridMG(newGrid);
1415 }
1416 /*********************************************************/
prolongationGridMG(GridMG * g)1417 void prolongationGridMG(GridMG* g)
1418 {
1419 interpolationTriLinearGridMG(g);
1420 }
1421 /*********************************************************/
levelDownGridMG(GridMG * g)1422 void levelDownGridMG(GridMG* g)
1423 {
1424 levelDownDomainMG(&g->domain);
1425 reAllocValuesTableGridMG(g);
1426 }
1427 /*********************************************************/
restrictionSrcGridMG(GridMG * g,GridMG * src)1428 void restrictionSrcGridMG(GridMG* g, GridMG* src)
1429 {
1430 gint ix, iy, iz;
1431 gint x0, xp, xm, y0, yp, ym, z0, zp, zm;
1432 gdouble face, corner, edge;
1433 static gdouble scale = 1.0 / 64.0;
1434 DomainMG domain;
1435
1436 /*printf("Begin restriction\n");*/
1437
1438 g->domain = src->domain;
1439 levelDownGridMG(g);
1440 domain = g->domain;
1441
1442 gint iXBegin = domain.iXBeginInterior;
1443 gint iXEnd = domain.iXEndInterior;
1444 gint iYBegin = domain.iYBeginInterior;
1445 gint iYEnd = domain.iYEndInterior;
1446 gint iZBegin = domain.iZBeginInterior;
1447 gint iZEnd = domain.iZEndInterior;
1448
1449 #ifdef ENABLE_OMP
1450 #pragma omp parallel for private(ix,iy,iz,x0,xp,xm,y0,yp,ym,z0,zp,zm,face,corner,edge)
1451 #endif
1452 for(ix = iXBegin ; ix <= iXEnd ; ix++)
1453 {
1454 x0 = 2 * ix;
1455 xp = x0 + 1;
1456 xm = x0 - 1;
1457 for(iy = iYBegin ; iy <= iYEnd ; iy++)
1458 {
1459 y0 = 2 * iy;
1460 yp = y0 + 1;
1461 ym = y0 - 1;
1462 for(iz = iZBegin ; iz <= iZEnd ; iz++)
1463 {
1464 z0 = 2 * iz;
1465 zp = z0 + 1;
1466 zm = z0 - 1;
1467
1468 face =
1469 getValGridMG(src, xm , y0 , z0) +
1470 getValGridMG(src, xp , y0 , z0) +
1471 getValGridMG(src, x0 , ym , z0) +
1472 getValGridMG(src, x0 , yp , z0) +
1473 getValGridMG(src, x0 , y0 , zm) +
1474 getValGridMG(src, x0 , y0 , zp);
1475
1476 corner =
1477 getValGridMG(src, xm , ym , zm) +
1478 getValGridMG(src, xm , ym , zp) +
1479 getValGridMG(src, xm , yp , zm) +
1480 getValGridMG(src, xm , yp , zp) +
1481 getValGridMG(src, xp , ym , zm) +
1482 getValGridMG(src, xp , ym , zp) +
1483 getValGridMG(src, xp , yp , zm) +
1484 getValGridMG(src, xp , yp , zp);
1485
1486 edge =
1487 getValGridMG(src, xm , y0 , zm) +
1488 getValGridMG(src, xm , ym , z0) +
1489 getValGridMG(src, xm , yp , z0) +
1490 getValGridMG(src, xm , y0 , zp) +
1491 getValGridMG(src, x0 , ym , zm) +
1492 getValGridMG(src, x0 , yp , zm) +
1493 getValGridMG(src, x0 , ym , zp) +
1494 getValGridMG(src, x0 , yp , zp) +
1495 getValGridMG(src, xp , y0 , zm) +
1496 getValGridMG(src, xp , ym , z0) +
1497 getValGridMG(src, xp , yp , z0) +
1498 getValGridMG(src, xp , y0 , zp);
1499
1500
1501 setValGridMG(g, ix , iy , iz,
1502 scale * (
1503 8.0 * getValGridMG(src, x0 , y0 , z0) +
1504 4.0 * face +
1505 2.0 * edge +
1506 corner
1507 ));
1508 }
1509 }
1510 }
1511 }
1512 /*********************************************************/
restrictionGridMG(GridMG * g)1513 void restrictionGridMG(GridMG* g)
1514 {
1515 GridMG* newGrid = getNewGridMGFromOldGrid(g);
1516 restrictionSrcGridMG(g, newGrid);
1517 destroyGridMG(newGrid);
1518 }
1519 /*********************************************************/
restrictionInjectionSrcGridMG(GridMG * g,GridMG * src)1520 void restrictionInjectionSrcGridMG(GridMG* g, GridMG* src)
1521 {
1522 gint ix, iy, iz;
1523 DomainMG domain;
1524 gint x0,y0,z0;
1525
1526 g->domain = src->domain;
1527 levelDownGridMG(g);
1528 domain = g->domain;
1529
1530 gint iXBegin = domain.iXBeginInterior - 1;
1531 gint iXEnd = domain.iXEndInterior + 1;
1532 gint iYBegin = domain.iYBeginInterior - 1;
1533 gint iYEnd = domain.iYEndInterior + 1;
1534 gint iZBegin = domain.iZBeginInterior - 1;
1535 gint iZEnd = domain.iZEndInterior + 1;
1536
1537 #ifdef ENABLE_OMP
1538 #pragma omp parallel for private(ix,iy,iz,x0,y0,z0)
1539 #endif
1540 for(ix = iXBegin ; ix <= iXEnd ; ix++)
1541 {
1542 x0 = 2 * ix;
1543 for(iy = iYBegin ; iy <= iYEnd ; iy++)
1544 {
1545 y0 = 2 * iy;
1546 for(iz = iZBegin ; iz <= iZEnd ; iz++)
1547 {
1548 z0 = 2 * iz;
1549 setValGridMG(g, ix , iy , iz, getValGridMG(src, x0 , y0 , z0));
1550 }
1551 }
1552 }
1553 }
1554 /*********************************************************/
restrictionInjectionGridMG(GridMG * g)1555 void restrictionInjectionGridMG(GridMG* g)
1556 {
1557 GridMG* newGrid = getNewGridMGFromOldGrid(g);
1558 restrictionInjectionSrcGridMG(g, newGrid);
1559 destroyGridMG(newGrid);
1560 }
1561 /*********************************************************/
getDomainGridMG(GridMG * g)1562 DomainMG getDomainGridMG(GridMG* g)
1563 {
1564 return g->domain;
1565 }
1566 /*********************************************************/
dotAllGridMG(GridMG * g,GridMG * src)1567 static gdouble dotAllGridMG(GridMG* g, GridMG* src)
1568 {
1569 gdouble p = 0.0;
1570 glong i;
1571 if(g->domain.size != src->domain.size)
1572 {
1573 printf(" Error in doAll\n ");
1574 return 0.0;
1575 }
1576 #ifdef ENABLE_OMP
1577 #pragma omp parallel for private(i) reduction(+:p)
1578 #endif
1579 for(i = 0 ; i < g->domain.size ; i++)
1580 p += g->values[i]*src->values[i];
1581
1582 p *= g->domain.cellVolume;
1583 return p;
1584 }
1585 /*********************************************************/
dotInteriorGridMG(GridMG * g,GridMG * src)1586 static gdouble dotInteriorGridMG(GridMG* g, GridMG* src)
1587 {
1588
1589 DomainMG domain = g->domain;
1590 gint ix;
1591 gint iy;
1592 gint iz;
1593 gint iXBegin = domain.iXBeginInterior;
1594 gint iXEnd = domain.iXEndInterior;
1595 gint iYBegin = domain.iYBeginInterior;
1596 gint iYEnd = domain.iYEndInterior;
1597 gint iZBegin = domain.iZBeginInterior;
1598 gint iZEnd = domain.iZEndInterior;
1599
1600 gdouble p = 0;
1601
1602 if(g->domain.size != src->domain.size)
1603 {
1604 printf(" Error in doInterior\n ");
1605 return 0.0;
1606 }
1607
1608 #ifdef ENABLE_OMP
1609 #pragma omp parallel for private(ix,iy,iz) reduction(+:p)
1610 #endif
1611 for(ix = iXBegin;ix <=iXEnd;ix++)
1612 for(iy = iYBegin;iy <=iYEnd;iy++)
1613 for(iz = iZBegin;iz <=iZEnd;iz++)
1614 {
1615 p += getValGridMG(g, ix,iy,iz)*getValGridMG(src, ix,iy,iz);
1616 }
1617
1618 p *= domain.cellVolume;
1619 return p;
1620 }
1621 /*********************************************************/
dotBoundaryGridMG(GridMG * g,GridMG * src)1622 static gdouble dotBoundaryGridMG(GridMG* g, GridMG* src)
1623 {
1624 gint ix;
1625 gint iy;
1626 gint iz;
1627 gdouble p = 0.0;
1628 DomainMG domain = g->domain;
1629 if(g->domain.size != src->domain.size)
1630 {
1631 printf(" Error in doBoundary\n ");
1632 return 0.0;
1633 }
1634
1635 #ifdef ENABLE_OMP
1636 #pragma omp parallel for private(ix,iy,iz) reduction(+:p)
1637 #endif
1638 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
1639 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
1640 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1641 {
1642 p += getValGridMG(g,ix,iy,iz)*getValGridMG(src, ix,iy,iz);
1643 }
1644 #ifdef ENABLE_OMP
1645 #pragma omp parallel for private(ix,iy,iz) reduction(+:p)
1646 #endif
1647 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
1648 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
1649 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1650 {
1651 p += getValGridMG(g,ix,iy,iz)*getValGridMG(src, ix,iy,iz);
1652 }
1653
1654 #ifdef ENABLE_OMP
1655 #pragma omp parallel for private(ix,iy,iz) reduction(+:p)
1656 #endif
1657 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
1658 {
1659 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
1660 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1661 {
1662 p += getValGridMG(g,ix,iy,iz)*getValGridMG(src, ix,iy,iz);
1663 }
1664 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
1665 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1666 {
1667 p += getValGridMG(g,ix,iy,iz)*getValGridMG(src, ix,iy,iz);
1668 }
1669 }
1670 #ifdef ENABLE_OMP
1671 #pragma omp parallel for private(ix,iy,iz) reduction(+:p)
1672 #endif
1673 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
1674 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
1675 {
1676 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
1677 {
1678 p += getValGridMG(g,ix,iy,iz)*getValGridMG(src, ix,iy,iz);
1679 }
1680 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
1681 {
1682 p += getValGridMG(g,ix,iy,iz)*getValGridMG(src, ix,iy,iz);
1683 }
1684 }
1685 p *= domain.cellVolume;
1686 return p;
1687
1688 }
1689 /*********************************************************/
dotGridMG(GridMG * g,GridMG * src)1690 gdouble dotGridMG(GridMG*g, GridMG* src)
1691 {
1692 switch(g->operationType)
1693 {
1694 case GABEDIT_ALL: return dotAllGridMG(g, src);break;
1695 case GABEDIT_INTERIOR: return dotInteriorGridMG(g, src);break;
1696 case GABEDIT_BOUNDARY: return dotBoundaryGridMG(g, src);break;
1697 }
1698 return 1.0;
1699 }
1700 /*********************************************************/
normAllGridMG(GridMG * g)1701 static gdouble normAllGridMG(GridMG* g)
1702 {
1703 glong i;
1704 gdouble n = 0;
1705
1706 #ifdef ENABLE_OMP
1707 #pragma omp parallel for private(i) reduction(+:n)
1708 #endif
1709 for(i = 0;i < g->domain.size ; i++)
1710 n += g->values[i]*g->values[i];
1711
1712 return sqrt(n*g->domain.cellVolume);
1713 }
1714 /*********************************************************/
normInteriorGridMG(GridMG * g)1715 static gdouble normInteriorGridMG(GridMG* g)
1716 {
1717
1718 DomainMG domain = g->domain;
1719 gdouble n = 0;
1720 gint ix;
1721 gint iy;
1722 gint iz;
1723 gint iXBegin = domain.iXBeginInterior;
1724 gint iXEnd = domain.iXEndInterior;
1725 gint iYBegin = domain.iYBeginInterior;
1726 gint iYEnd = domain.iYEndInterior;
1727 gint iZBegin = domain.iZBeginInterior;
1728 gint iZEnd = domain.iZEndInterior;
1729
1730 if(g->domain.size<=0) return 0;
1731
1732 #ifdef ENABLE_OMP
1733 #pragma omp parallel for private(ix,iy,iz) reduction(+:n)
1734 #endif
1735 for(ix = iXBegin;ix <=iXEnd;ix++)
1736 for(iy = iYBegin;iy <=iYEnd;iy++)
1737 for(iz = iZBegin;iz <=iZEnd;iz++)
1738 n += getValGridMG(g,ix,iy,iz)*getValGridMG(g,ix,iy,iz);
1739
1740 return sqrt(n*domain.cellVolume);
1741 }
1742 /*********************************************************/
normBoundaryGridMG(GridMG * g)1743 static gdouble normBoundaryGridMG(GridMG* g)
1744 {
1745 gdouble n = 0;
1746 DomainMG domain = g->domain;
1747 gint ix;
1748 gint iy;
1749 gint iz;
1750
1751 #ifdef ENABLE_OMP
1752 #pragma omp parallel for private(ix,iy,iz) reduction(+:n)
1753 #endif
1754 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
1755 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
1756 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1757 n += getValGridMG(g,ix,iy,iz)*getValGridMG(g,ix,iy,iz);
1758 #ifdef ENABLE_OMP
1759 #pragma omp parallel for private(ix,iy,iz) reduction(+:n)
1760 #endif
1761 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
1762 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
1763 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1764 n += getValGridMG(g,ix,iy,iz)*getValGridMG(g,ix,iy,iz);
1765
1766 #ifdef ENABLE_OMP
1767 #pragma omp parallel for private(ix,iy,iz) reduction(+:n)
1768 #endif
1769 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
1770 {
1771 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
1772 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1773 n += getValGridMG(g,ix,iy,iz)*getValGridMG(g,ix,iy,iz);
1774 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
1775 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1776 n += getValGridMG(g,ix,iy,iz)*getValGridMG(g,ix,iy,iz);
1777 }
1778 #ifdef ENABLE_OMP
1779 #pragma omp parallel for private(ix,iy,iz) reduction(+:n)
1780 #endif
1781 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
1782 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
1783 {
1784 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
1785 n += getValGridMG(g,ix,iy,iz)*getValGridMG(g,ix,iy,iz);
1786 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
1787 n += getValGridMG(g,ix,iy,iz)*getValGridMG(g,ix,iy,iz);
1788 }
1789 return sqrt(n*domain.cellVolume);
1790 }
1791 /*********************************************************/
normGridMG(GridMG * g)1792 gdouble normGridMG(GridMG* g)
1793 {
1794 gdouble n = 0;
1795 switch(g->operationType)
1796 {
1797 case GABEDIT_ALL: n = normAllGridMG(g);break;
1798 case GABEDIT_INTERIOR: n = normInteriorGridMG(g);break;
1799 case GABEDIT_BOUNDARY: n = normBoundaryGridMG(g);break;
1800 }
1801 return n;
1802 }
1803 /*********************************************************/
normDiffAllGridMG(GridMG * g,GridMG * src)1804 static gdouble normDiffAllGridMG(GridMG* g, GridMG* src)
1805 {
1806 glong i;
1807 gdouble n = 0;
1808 DomainMG domain = g->domain;
1809
1810 #ifdef ENABLE_OMP
1811 #pragma omp parallel for private(i) reduction(+:n)
1812 #endif
1813 for(i = 0;i < domain.size ; i++)
1814 n += (g->values[i]-src->values[i])*(g->values[i]-src->values[i]);
1815
1816 return sqrt(n/(domain.size));
1817 }
1818 /*********************************************************/
normDiffInteriorGridMG(GridMG * g,GridMG * src)1819 static gdouble normDiffInteriorGridMG(GridMG* g, GridMG* src)
1820 {
1821 DomainMG domain = g->domain;
1822
1823 gdouble n = 0;
1824 gint ix;
1825 gint iy;
1826 gint iz;
1827 gint iXBegin = domain.iXBeginInterior;
1828 gint iXEnd = domain.iXEndInterior;
1829 gint iYBegin = domain.iYBeginInterior;
1830 gint iYEnd = domain.iYEndInterior;
1831 gint iZBegin = domain.iZBeginInterior;
1832 gint iZEnd = domain.iZEndInterior;
1833 gdouble v;
1834 if(domain.size<=0)
1835 return 0;
1836
1837 #ifdef ENABLE_OMP
1838 #pragma omp parallel for private(ix,iy,iz,v) reduction(+:n)
1839 #endif
1840 for(ix = iXBegin;ix <=iXEnd;ix++)
1841 for(iy = iYBegin;iy <=iYEnd;iy++)
1842 for(iz = iZBegin;iz <=iZEnd;iz++)
1843 {
1844 v = getValGridMG(g,ix, iy, iz)-getValGridMG(src, ix, iy, iz );
1845 n += v*v ;
1846 }
1847
1848 return sqrt(n/((domain.xSize) * (domain.ySize)*(domain.zSize)));
1849 }
1850 /*********************************************************/
normDiffBoundaryGridMG(GridMG * g,GridMG * src)1851 static gdouble normDiffBoundaryGridMG(GridMG* g, GridMG* src)
1852 {
1853 gdouble n = 0;
1854
1855 gint ix;
1856 gint iy;
1857 gint iz;
1858 gdouble v;
1859 DomainMG domain = g->domain;
1860
1861 #ifdef ENABLE_OMP
1862 #pragma omp parallel for private(ix,iy,iz,v) reduction(+:n)
1863 #endif
1864 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
1865 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
1866 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1867 {
1868 v = getValGridMG(g,ix, iy, iz)-getValGridMG(src, ix, iy, iz );
1869 n += v*v ;
1870 }
1871 #ifdef ENABLE_OMP
1872 #pragma omp parallel for private(ix,iy,iz,v) reduction(+:n)
1873 #endif
1874 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
1875 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
1876 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1877 {
1878 v = getValGridMG(g,ix, iy, iz)-getValGridMG(src, ix, iy, iz );
1879 n += v*v ;
1880 }
1881
1882 #ifdef ENABLE_OMP
1883 #pragma omp parallel for private(ix,iy,iz,v) reduction(+:n)
1884 #endif
1885 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
1886 {
1887 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
1888 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1889 {
1890 v = getValGridMG(g,ix, iy, iz)-getValGridMG(src, ix, iy, iz );
1891 n += v*v ;
1892 }
1893 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
1894 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1895 {
1896 v = getValGridMG(g,ix, iy, iz)-getValGridMG(src, ix, iy, iz );
1897 n += v*v ;
1898 }
1899 }
1900 #ifdef ENABLE_OMP
1901 #pragma omp parallel for private(ix,iy,iz,v) reduction(+:n)
1902 #endif
1903 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
1904 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
1905 {
1906 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
1907 {
1908 v = getValGridMG(g,ix, iy, iz)-getValGridMG(src, ix, iy, iz );
1909 n += v*v ;
1910 }
1911 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
1912 {
1913 v = getValGridMG(g,ix, iy, iz)-getValGridMG(src, ix, iy, iz );
1914 n += v*v ;
1915 }
1916 }
1917 return sqrt(
1918 n/(
1919 (domain.ySize) * (domain.zSize )*2*domain.nBoundary
1920 + (domain.xSize) * (domain.zSize )*2*domain.nBoundary
1921 + (domain.ySize) * (domain.zSize )*2*domain.nBoundary
1922 )
1923 );
1924 }
1925 /*********************************************************/
normDiffGridMG(GridMG * g,GridMG * src)1926 gdouble normDiffGridMG(GridMG* g, GridMG* src)
1927 {
1928 gdouble n = 0;
1929 switch(g->operationType)
1930 {
1931 case GABEDIT_ALL: n = normDiffAllGridMG(g,src);break;
1932 case GABEDIT_INTERIOR: n = normDiffInteriorGridMG(g,src);break;
1933 case GABEDIT_BOUNDARY: n = normDiffBoundaryGridMG(g,src);break;
1934 }
1935 return n;
1936 }
1937 /*********************************************************/
sommeAllGridMG(GridMG * g)1938 static gdouble sommeAllGridMG(GridMG* g)
1939 {
1940 glong i;
1941 gdouble s = 0;
1942 DomainMG domain = g->domain;
1943
1944 #ifdef ENABLE_OMP
1945 #pragma omp parallel for private(i) reduction(+:s)
1946 #endif
1947 for(i = 0;i < domain.size ; i++)
1948 s += g->values[i];
1949
1950 s *= domain.cellVolume;
1951
1952 return s;
1953 }
1954 /*********************************************************/
sommeInteriorGridMG(GridMG * g)1955 static gdouble sommeInteriorGridMG(GridMG* g)
1956 {
1957
1958 gdouble s = 0;
1959 DomainMG domain = g->domain;
1960 gint ix;
1961 gint iy;
1962 gint iz;
1963 gint iXBegin = domain.iXBeginInterior;
1964 gint iXEnd = domain.iXEndInterior;
1965 gint iYBegin = domain.iYBeginInterior;
1966 gint iYEnd = domain.iYEndInterior;
1967 gint iZBegin = domain.iZBeginInterior;
1968 gint iZEnd = domain.iZEndInterior;
1969
1970 if(domain.size<=0) return 0;
1971
1972 #ifdef ENABLE_OMP
1973 #pragma omp parallel for private(ix,iy,iz) reduction(+:s)
1974 #endif
1975 for(ix = iXBegin;ix <=iXEnd;ix++)
1976 for(iy = iYBegin;iy <=iYEnd;iy++)
1977 for(iz = iZBegin;iz <=iZEnd;iz++)
1978 s += getValGridMG(g,ix,iy,iz);
1979
1980 s *= domain.cellVolume;
1981 return s;
1982 }
1983 /*********************************************************/
sommeBoundaryGridMG(GridMG * g)1984 static gdouble sommeBoundaryGridMG(GridMG* g)
1985 {
1986 gdouble s = 0;
1987
1988 gint ix;
1989 gint iy;
1990 gint iz;
1991 DomainMG domain = g->domain;
1992
1993 #ifdef ENABLE_OMP
1994 #pragma omp parallel for private(ix,iy,iz) reduction(+:s)
1995 #endif
1996 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
1997 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
1998 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
1999 s += getValGridMG(g,ix,iy,iz);
2000 #ifdef ENABLE_OMP
2001 #pragma omp parallel for private(ix,iy,iz) reduction(+:s)
2002 #endif
2003 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
2004 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
2005 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
2006 s += getValGridMG(g,ix,iy,iz);
2007 #ifdef ENABLE_OMP
2008 #pragma omp parallel for private(ix,iy,iz) reduction(+:s)
2009 #endif
2010 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
2011 {
2012 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
2013 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
2014 s += getValGridMG(g,ix,iy,iz);
2015 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
2016 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
2017 s += getValGridMG(g,ix,iy,iz);
2018 }
2019 #ifdef ENABLE_OMP
2020 #pragma omp parallel for private(ix,iy,iz) reduction(+:s)
2021 #endif
2022 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
2023 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
2024 {
2025 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
2026 s += getValGridMG(g,ix,iy,iz);
2027 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
2028 s += getValGridMG(g,ix,iy,iz);
2029 }
2030 s *= domain.cellVolume;
2031 return s;
2032 }
2033 /*********************************************************/
sommeGridMG(GridMG * g)2034 gdouble sommeGridMG(GridMG* g)
2035 {
2036 gdouble s = 0;
2037 switch(g->operationType)
2038 {
2039 case GABEDIT_ALL: s = sommeAllGridMG(g);break;
2040 case GABEDIT_INTERIOR: s = sommeInteriorGridMG(g);break;
2041 case GABEDIT_BOUNDARY: s = sommeBoundaryGridMG(g);break;
2042 }
2043 return s;
2044 }
2045 /*********************************************************/
normalizeGridMG(GridMG * g)2046 gdouble normalizeGridMG(GridMG* g)
2047 {
2048 gdouble sum2 = dotInteriorGridMG(g,g);
2049 sum2 = 1/sqrt(sum2);
2050 multEqualInteriorRealGridMG(g,sum2);
2051 return sum2;
2052 }
2053 /*********************************************************/
setOperationGridMG(GridMG * g,const OperationTypeMG operation)2054 void setOperationGridMG(GridMG* g, const OperationTypeMG operation)
2055 {
2056 g->operationType = operation;
2057 }
2058 /*********************************************************/
getOperationGridMG(GridMG * g)2059 OperationTypeMG getOperationGridMG(GridMG* g)
2060 {
2061 return g->operationType;
2062 }
2063 /*********************************************************/
tradesBoundaryPeriodicGridMG(GridMG * g)2064 void tradesBoundaryPeriodicGridMG(GridMG* g)
2065 {
2066 gint ix;
2067 gint iy;
2068 gint iz;
2069
2070 gint j;
2071 DomainMG domain = g->domain;
2072
2073 #ifdef ENABLE_OMP
2074 #pragma omp parallel for private(ix,iy,iz,j)
2075 #endif
2076 for(ix=domain.iXBeginBoundaryLeft ; ix <= domain.iXEndBoundaryLeft ; ix++)
2077 {
2078 for(iy = domain.iYBeginBoundaryLeft ; iy <= domain.iYEndBoundaryRight ; iy++)
2079 for(iz = domain.iZBeginBoundaryLeft ; iz <= domain.iZEndBoundaryRight ; iz++)
2080 {
2081 j= domain.iXEndInterior - domain.nBoundary+ix-domain.iXBeginBoundaryLeft;
2082 setValGridMG(g, ix, iy, iz, getValGridMG(g, j, iy, iz));
2083 }
2084 }
2085
2086 #ifdef ENABLE_OMP
2087 #pragma omp parallel for private(ix,iy,iz,j)
2088 #endif
2089 for(ix=domain.iXBeginBoundaryRight ; ix <= domain.iXEndBoundaryRight ; ix++)
2090 {
2091 for(iy = domain.iYBeginBoundaryLeft ; iy <= domain.iYEndBoundaryRight ; iy++)
2092 for(iz = domain.iZBeginBoundaryLeft ; iz <= domain.iZEndBoundaryRight ; iz++)
2093 {
2094 j = domain.iXBeginInterior+ix-domain.iXBeginBoundaryRight;
2095 setValGridMG(g, ix, iy, iz, getValGridMG(g, j, iy, iz));
2096 }
2097 }
2098
2099 #ifdef ENABLE_OMP
2100 #pragma omp parallel for private(ix,iy,iz,j)
2101 #endif
2102 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
2103 {
2104 j = domain.iYEndInterior - domain.nBoundary;
2105 for(iy = domain.iYBeginBoundaryLeft ; iy <=domain.iYEndBoundaryLeft ; iy++, j++)
2106 for(iz = domain.iZBeginBoundaryLeft ; iz <=domain.iZEndBoundaryRight ; iz++)
2107 setValGridMG(g, ix, iy, iz, getValGridMG(g, ix, j, iz));
2108
2109 j = domain.iYBeginInterior;
2110 for(iy = domain.iYBeginBoundaryRight ; iy <=domain.iYEndBoundaryRight; iy++, j++)
2111 for(iz = domain.iZBeginBoundaryLeft; iz <=domain.iZEndBoundaryRight; iz++)
2112 setValGridMG(g, ix, iy, iz, getValGridMG(g, ix, j, iz));
2113 }
2114 #ifdef ENABLE_OMP
2115 #pragma omp parallel for private(ix,iy,iz,j)
2116 #endif
2117 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
2118 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
2119 {
2120 j = domain.iZEndInterior - domain.nBoundary;
2121 for(iz = domain.iZBeginBoundaryLeft ; iz <=domain.iZEndBoundaryLeft ; iz++, j++)
2122 setValGridMG(g, ix, iy, iz, getValGridMG(g, ix, iy, j));
2123
2124 j = domain.iZBeginInterior;
2125 for(iz = domain.iZBeginBoundaryRight ; iz <=domain.iZEndBoundaryRight ; iz++, j++)
2126 setValGridMG(g, ix, iy, iz, getValGridMG(g, ix, iy, j));
2127 }
2128
2129 }
2130 /*********************************************************/
tradesBoundaryGridMG(GridMG * g,const Condition condition)2131 void tradesBoundaryGridMG(GridMG* g, const Condition condition)
2132 {
2133 switch(condition)
2134 {
2135 case GABEDIT_CONDITION_PERIODIC : tradesBoundaryPeriodicGridMG(g);break;
2136 case GABEDIT_CONDITION_CLUSTER : initBoundaryGridMG(g,0.0);break;
2137 case GABEDIT_CONDITION_MULTIPOL :
2138 case GABEDIT_CONDITION_EXTERNAL :
2139 case GABEDIT_CONDITION_EWALD :
2140 printf("Error(Grid Class), I can not set boundaris using MULTIPOL approximation or EXTERNAL\n");
2141 break;
2142 }
2143 }
2144 /**************************************************************************************/
printGridMG(GridMG * g,const gint ix,const gint iy,const gint iz)2145 void printGridMG(GridMG* g, const gint ix, const gint iy, const gint iz)
2146 {
2147 char t[BSIZE];
2148
2149 sprintf(t,"SLICE %4d %4d %4d",ix,iy,iz);
2150 printf("%20s %14.8f \n", t, getValGridMG(g, ix, iy, iz));
2151 }
2152 /*********************************************************/
printAllGridMG(GridMG * g)2153 void printAllGridMG(GridMG* g)
2154 {
2155 gint ix;
2156 gint iy;
2157 gint iz;
2158 char t1[BSIZE];
2159 DomainMG domain = g->domain;
2160
2161 for(ix=domain.iXBeginBoundaryLeft ; ix<=domain.iXEndBoundaryRight ; ix++)
2162 for(iy = domain.iYBeginBoundaryLeft ; iy <=domain.iYEndBoundaryRight ; iy++)
2163 for(iz = domain.iZBeginBoundaryLeft ; iz <=domain.iZEndBoundaryRight ; iz++)
2164 {
2165 sprintf(t1,"%d %d %d %14.8f\n",ix, iy, iz, getValGridMG(g,ix, iy, iz));
2166 printf("%s",t1);
2167 }
2168 }
2169 /*********************************************************/
printInteriorGridMG(GridMG * g)2170 void printInteriorGridMG(GridMG* g)
2171 {
2172 char t1[BSIZE];
2173
2174 DomainMG domain = g->domain;
2175 gint ix;
2176 gint iy;
2177 gint iz;
2178 gint iXBegin = domain.iXBeginInterior;
2179 gint iXEnd = domain.iXEndInterior;
2180 gint iYBegin = domain.iYBeginInterior;
2181 gint iYEnd = domain.iYEndInterior;
2182 gint iZBegin = domain.iZBeginInterior;
2183 gint iZEnd = domain.iZEndInterior;
2184
2185 for(ix = iXBegin;ix <=iXEnd;ix++)
2186 for(iy = iYBegin;iy <=iYEnd;iy++)
2187 for(iz = iZBegin;iz <=iZEnd;iz++)
2188 {
2189 sprintf(t1,"%d %d %d %14.8f\n",ix, iy, iz, getValGridMG(g,ix,iy,iz));
2190 /*
2191 sprintf(t1,"%d %d %d %f %f %f %14.8f\n",ix, iy, iz,
2192 domain.x0+ix*domain.xh,
2193 domain.y0+iy*domain.yh,
2194 domain.z0+iz*domain.zh,
2195 getValGridMG(g,ix,iy,iz));
2196 */
2197 printf("%s",t1);
2198 }
2199 }
2200 /*********************************************************/
printBoundaryGridMG(GridMG * g)2201 void printBoundaryGridMG(GridMG* g)
2202 {
2203 char t1[BSIZE];
2204 gint ix;
2205 gint iy;
2206 gint iz;
2207 DomainMG domain = g->domain;
2208
2209 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
2210 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
2211 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
2212 {
2213 sprintf(t1,"%d %d %d %14.8f\n",ix, iy, iz, getValGridMG(g,ix, iy, iz));
2214 printf("%s",t1);
2215 }
2216 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
2217 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
2218 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
2219 {
2220 sprintf(t1,"%d %d %d %14.8f\n",ix, iy, iz, getValGridMG(g,ix, iy, iz));
2221 printf("%s",t1);
2222 }
2223
2224 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
2225 {
2226 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
2227 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
2228 {
2229 sprintf(t1,"%d %d %d %14.8f\n",ix, iy, iz, getValGridMG(g,ix, iy, iz));
2230 printf("%s",t1);
2231 }
2232 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
2233 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
2234 {
2235 sprintf(t1,"%d %d %d %14.8f\n",ix, iy, iz, getValGridMG(g,ix, iy, iz));
2236 printf("%s",t1);
2237 }
2238 }
2239 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
2240 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
2241 {
2242 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
2243 {
2244 sprintf(t1,"%d %d %d %14.8f\n",ix, iy, iz, getValGridMG(g,ix, iy, iz));
2245 printf("%s",t1);
2246 }
2247 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
2248 {
2249 sprintf(t1,"%d %d %d %14.8f\n",ix, iy, iz, getValGridMG(g,ix, iy, iz));
2250 printf("%s",t1);
2251 }
2252 }
2253 }
2254 /*********************************************************/
printGridMGAll(GridMG * g)2255 void printGridMGAll(GridMG* g)
2256 {
2257 switch(g->operationType)
2258 {
2259 case GABEDIT_ALL: printAllGridMG(g);break;
2260 case GABEDIT_INTERIOR: printInteriorGridMG(g);break;
2261 case GABEDIT_BOUNDARY: printBoundaryGridMG(g);break;
2262 }
2263 }
2264 /*********************************************************/
printGridFileGridMG(GridMG * g,FILE * file)2265 void printGridFileGridMG(GridMG* g, FILE* file)
2266 {
2267 DomainMG domain = g->domain;
2268 gint ix;
2269 gint iy;
2270 gint iz;
2271 char t1[BSIZE];
2272
2273 fprintf(file,"GABEDIT\n");
2274 fprintf(file,"Poisson density\n");
2275 sprintf(t1,"%d %14.8f %14.8f %14.8f",-(8),domain.x0,domain.y0,domain.z0);
2276 fprintf(file,"%s\n",t1);
2277 sprintf(t1,"%d %14.8f %14.8f %14.8f",domain.xSize + 2,domain.xh,0.0,0.0);
2278 fprintf(file,"%s\n",t1);
2279 sprintf(t1,"%d %14.8f %14.8f %14.8f",domain.ySize + 2 ,0.0,domain.yh,0.0);
2280 fprintf(file,"%s\n",t1);
2281 sprintf(t1,"%d %14.8f %14.8f %14.8f",domain.zSize + 2 ,0.0,0.0,domain.zh);
2282 fprintf(file,"%s\n",t1);
2283
2284 gdouble x,y,z;
2285
2286 x = domain.x0 + domain.iXEndBoundaryLeft * domain.xh;
2287 y = domain.y0 + domain.iYEndBoundaryLeft * domain.yh;
2288 z = domain.z0 + domain.iZEndBoundaryLeft * domain.zh;
2289 sprintf(t1,"%d %14.8f %14.8f %14.8f %14.8f ", 4, 4.0, x, y, z);
2290 fprintf(file,"%s\n",t1);
2291
2292 x = domain.x0 + domain.iXEndBoundaryLeft * domain.xh;
2293 y = domain.y0 + domain.iYEndBoundaryLeft * domain.yh;
2294 z = domain.z0 + domain.iZBeginBoundaryRight * domain.zh;
2295 sprintf(t1,"%d %14.8f %14.8f %14.8f %14.8f ", 4, 4.0, x, y, z);
2296 fprintf(file,"%s\n",t1);
2297
2298 x = domain.x0 + domain.iXEndBoundaryLeft * domain.xh;
2299 y = domain.y0 + domain.iYBeginBoundaryRight * domain.yh;
2300 z = domain.z0 + domain.iZEndBoundaryLeft * domain.zh;
2301 sprintf(t1,"%d %14.8f %14.8f %14.8f %14.8f ", 4, 4.0, x, y, z);
2302 fprintf(file,"%s\n",t1);
2303
2304 x = domain.x0 + domain.iXBeginBoundaryRight * domain.xh;
2305 y = domain.y0 + domain.iYEndBoundaryLeft * domain.yh;
2306 z = domain.z0 + domain.iZEndBoundaryLeft * domain.zh;
2307 sprintf(t1,"%d %14.8f %14.8f %14.8f %14.8f ", 4, 4.0, x, y, z);
2308 fprintf(file,"%s\n",t1);
2309
2310 x = domain.x0 + domain.iXBeginBoundaryRight * domain.xh;
2311 y = domain.y0 + domain.iYBeginBoundaryRight * domain.yh;
2312 z = domain.z0 + domain.iZEndBoundaryLeft * domain.zh;
2313 sprintf(t1,"%d %14.8f %14.8f %14.8f %14.8f ", 4, 4.0, x, y, z);
2314 fprintf(file,"%s\n",t1);
2315
2316 x = domain.x0 + domain.iXBeginBoundaryRight * domain.xh;
2317 y = domain.y0 + domain.iYEndBoundaryLeft * domain.yh;
2318 z = domain.z0 + domain.iZBeginBoundaryRight * domain.zh;
2319 sprintf(t1,"%d %14.8f %14.8f %14.8f %14.8f ", 4, 4.0, x, y, z);
2320 fprintf(file,"%s\n",t1);
2321
2322 x = domain.x0 + domain.iXEndBoundaryLeft * domain.xh;
2323 y = domain.y0 + domain.iYBeginBoundaryRight * domain.yh;
2324 z = domain.z0 + domain.iZBeginBoundaryRight * domain.zh;
2325 sprintf(t1,"%d %14.8f %14.8f %14.8f %14.8f ", 4, 4.0, x, y, z);
2326 fprintf(file,"%s\n",t1);
2327
2328 x = domain.x0 + domain.iXBeginBoundaryRight * domain.xh;
2329 y = domain.y0 + domain.iYBeginBoundaryRight * domain.yh;
2330 z = domain.z0 + domain.iZBeginBoundaryRight * domain.zh;
2331 sprintf(t1,"%d %14.8f %14.8f %14.8f %14.8f ", 4, 4.0, x, y, z);
2332 fprintf(file,"%s\n",t1);
2333
2334 sprintf(t1,"%d %d",1,1);
2335 fprintf(file,"%s\n",t1);
2336
2337
2338 for(ix=domain.iXBeginInterior-1 ; ix<=domain.iXEndInterior + 1 ; ix++)
2339 for(iy = domain.iYBeginInterior - 1 ; iy <=domain.iYEndInterior + 1 ; iy++)
2340 {
2341 for(iz = domain.iZBeginInterior -1 ; iz <=domain.iZEndInterior + 1 ; iz++)
2342 {
2343 sprintf(t1,"%14.8f",getValGridMG(g,ix, iy, iz));
2344 fprintf(file,"%s ",t1);
2345 if((iz+1)%6==0)
2346 fprintf(file,"\n");
2347 }
2348 if((domain.iZEndBoundaryRight - domain.iZBeginBoundaryLeft +1)%6 != 0)
2349 fprintf(file,"\n");
2350 }
2351 }
2352 /*********************************************************/
printFileNameGridMG(GridMG * g,char * fileName)2353 void printFileNameGridMG(GridMG* g, char* fileName)
2354 {
2355 FILE* file;
2356
2357 file = fopen(fileName,"w");
2358 if(!file) return;
2359
2360 printGridFileGridMG(g,file);
2361 fclose(file);
2362 }
2363 /*********************************************************/
getMaxGridMG(GridMG * g)2364 gdouble getMaxGridMG(GridMG* g)
2365 {
2366
2367 DomainMG domain = g->domain;
2368 gint ix;
2369 gint iy;
2370 gint iz;
2371 gint ixMax = domain.iXBeginInterior;
2372 gint iyMax = domain.iYBeginInterior;
2373 gint izMax = domain.iZBeginInterior;
2374
2375 if(domain.size <1)
2376 {
2377 printf("ERROR Size =%d\n",domain.size);
2378 return -1;
2379 }
2380
2381
2382 #ifdef ENABLE_OMP
2383 #pragma omp parallel for private(ix,iy,iz)
2384 #endif
2385 for(ix=domain.iXBeginInterior ; ix<=domain.iXEndInterior ; ix++)
2386 for(iy = domain.iYBeginInterior ; iy <=domain.iYEndInterior ; iy++)
2387 for(iz = domain.iZBeginInterior ; iz <=domain.iZEndInterior ; iz++)
2388 if(getValGridMG(g,ixMax, iyMax, izMax) < getValGridMG(g,ix, iy, iz))
2389 {
2390 ixMax = ix;
2391 iyMax = iy;
2392 izMax = iz;
2393 }
2394 return getValGridMG(g,ixMax, iyMax, izMax);
2395 }
2396 /*********************************************************/
printMaxGridMG(GridMG * g)2397 void printMaxGridMG(GridMG* g)
2398 {
2399
2400 DomainMG domain = g->domain;
2401 gint ix;
2402 gint iy;
2403 gint iz;
2404 gint ixMax = domain.iXBeginInterior;
2405 gint iyMax = domain.iYBeginInterior;
2406 gint izMax = domain.iZBeginInterior;
2407
2408 if(domain.size <1)
2409 {
2410 printf("ERROR Size =%d\n",domain.size);
2411 return;
2412 }
2413
2414 for(ix=domain.iXBeginInterior ; ix<=domain.iXEndInterior ; ix++)
2415 for(iy = domain.iYBeginInterior ; iy <=domain.iYEndInterior ; iy++)
2416 for(iz = domain.iZBeginInterior ; iz <=domain.iZEndInterior ; iz++)
2417 if(getValGridMG(g,ixMax, iyMax, izMax) < getValGridMG(g,ix, iy, iz))
2418 {
2419 ixMax = ix;
2420 iyMax = iy;
2421 izMax = iz;
2422 }
2423 printf("MAX :");
2424 printGridMG(g, ixMax,iyMax,izMax);
2425 }
2426 /*********************************************************/
printMinGridMG(GridMG * g)2427 void printMinGridMG(GridMG* g)
2428 {
2429
2430 DomainMG domain = g->domain;
2431 gint ix;
2432 gint iy;
2433 gint iz;
2434 gint ixMin = domain.iXBeginInterior;
2435 gint iyMin = domain.iYBeginInterior;
2436 gint izMin = domain.iZBeginInterior;
2437
2438 if(domain.size <1)
2439 {
2440 printf("ERROR Size =%d\n",domain.size);
2441 return;
2442 }
2443
2444 for(ix=domain.iXBeginInterior ; ix<=domain.iXEndInterior ; ix++)
2445 for(iy = domain.iYBeginInterior ; iy <=domain.iYEndInterior ; iy++)
2446 for(iz = domain.iZBeginInterior ; iz <=domain.iZEndInterior ; iz++)
2447 if(getValGridMG(g,ixMin, iyMin, izMin) > getValGridMG(g,ix, iy, iz))
2448 {
2449 ixMin = ix;
2450 iyMin = iy;
2451 izMin = iz;
2452 }
2453 printf("MIN :");
2454 printGridMG(g, ixMin,iyMin,izMin);
2455 }
2456 /*********************************************************/
getValGridMG(GridMG * g,gint ix,gint iy,gint iz)2457 gdouble getValGridMG(GridMG* g, gint ix, gint iy, gint iz)
2458 {
2459 /*
2460 gint i =
2461 (ix+g->domain.nShift)*g->domain.incx
2462 + (iy+g->domain.nShift)*g->domain.incy
2463 + (iz+g->domain.nShift)*g->domain.incz;
2464 if(i>g->domain.size) printf("ERROR i>size , ix = %d iy = %d iz = %d size = %d i = %d\n",ix,iy,iz,g->domain.size,i);
2465 */
2466 return g->values[
2467 (ix+g->domain.nShift)*g->domain.incx
2468 + (iy+g->domain.nShift)*g->domain.incy
2469 + (iz+g->domain.nShift)*g->domain.incz
2470 ];
2471 }
2472 /*********************************************************/
setValGridMG(GridMG * g,gint ix,gint iy,gint iz,gdouble v)2473 void setValGridMG(GridMG* g, gint ix, gint iy, gint iz, gdouble v)
2474 {
2475 g->values[
2476 (ix+g->domain.nShift)*g->domain.incx
2477 + (iy+g->domain.nShift)*g->domain.incy
2478 + (iz+g->domain.nShift)*g->domain.incz
2479 ]=v;
2480 }
2481 /*********************************************************/
addValGridMG(GridMG * g,gint ix,gint iy,gint iz,gdouble v)2482 void addValGridMG(GridMG* g, gint ix, gint iy, gint iz, gdouble v)
2483 {
2484 g->values[
2485 (ix+g->domain.nShift)*g->domain.incx
2486 + (iy+g->domain.nShift)*g->domain.incy
2487 + (iz+g->domain.nShift)*g->domain.incz
2488 ]+=v;
2489 }
2490 /*********************************************************/
multValGridMG(GridMG * g,gint ix,gint iy,gint iz,gdouble v)2491 void multValGridMG(GridMG* g, gint ix, gint iy, gint iz, gdouble v)
2492 {
2493 g->values[
2494 (ix+g->domain.nShift)*g->domain.incx
2495 + (iy+g->domain.nShift)*g->domain.incy
2496 + (iz+g->domain.nShift)*g->domain.incz
2497 ]*=v;
2498 }
2499