source: CIVL/examples/mpi-omp/AMG2013/struct_mv/struct_vector.c@ a389857

main test-branch
Last change on this file since a389857 was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 53.5 KB
Line 
1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13
14/******************************************************************************
15 *
16 * Member functions for hypre_StructVector class.
17 *
18 *****************************************************************************/
19
20#include "headers.h"
21
22/*--------------------------------------------------------------------------
23 * hypre_StructVectorCreate
24 *--------------------------------------------------------------------------*/
25
26hypre_StructVector *
27hypre_StructVectorCreate( MPI_Comm comm,
28 hypre_StructGrid *grid )
29{
30 hypre_StructVector *vector;
31
32 int i;
33
34 vector = hypre_CTAlloc(hypre_StructVector, 1);
35
36 hypre_StructVectorComm(vector) = comm;
37 hypre_StructGridRef(grid, &hypre_StructVectorGrid(vector));
38 hypre_StructVectorDataAlloced(vector) = 1;
39 hypre_StructVectorOffProcAdd(vector) = 0;
40 hypre_StructVectorRefCount(vector) = 1;
41
42 /* set defaults */
43 for (i = 0; i < 6; i++)
44 {
45 hypre_StructVectorNumGhost(vector)[i] = 1;
46 hypre_StructVectorAddNumGhost(vector)[i] = 1;
47 }
48
49 return vector;
50}
51
52/*--------------------------------------------------------------------------
53 * hypre_StructVectorRef
54 *--------------------------------------------------------------------------*/
55
56hypre_StructVector *
57hypre_StructVectorRef( hypre_StructVector *vector )
58{
59 hypre_StructVectorRefCount(vector) ++;
60
61 return vector;
62}
63
64/*--------------------------------------------------------------------------
65 * hypre_StructVectorDestroy
66 *--------------------------------------------------------------------------*/
67
68int
69hypre_StructVectorDestroy( hypre_StructVector *vector )
70{
71 int ierr = 0;
72
73 if (vector)
74 {
75 hypre_StructVectorRefCount(vector) --;
76 if (hypre_StructVectorRefCount(vector) == 0)
77 {
78 if (hypre_StructVectorDataAlloced(vector))
79 {
80 hypre_SharedTFree(hypre_StructVectorData(vector));
81 }
82 hypre_TFree(hypre_StructVectorDataIndices(vector));
83 hypre_BoxArrayDestroy(hypre_StructVectorDataSpace(vector));
84 hypre_StructGridDestroy(hypre_StructVectorGrid(vector));
85 hypre_TFree(vector);
86 }
87 }
88
89 return ierr;
90}
91
92/*--------------------------------------------------------------------------
93 * hypre_StructVectorInitializeShell
94 *--------------------------------------------------------------------------*/
95
96int
97hypre_StructVectorInitializeShell( hypre_StructVector *vector )
98{
99 int ierr = 0;
100
101 hypre_StructGrid *grid;
102
103 int *num_ghost;
104
105 hypre_BoxArray *data_space;
106 hypre_BoxArray *boxes;
107 hypre_Box *box;
108 hypre_Box *data_box;
109
110 int *data_indices;
111 int data_size;
112
113 int i, d;
114
115 /*-----------------------------------------------------------------------
116 * Set up data_space
117 *-----------------------------------------------------------------------*/
118
119 grid = hypre_StructVectorGrid(vector);
120
121 if (hypre_StructVectorDataSpace(vector) == NULL)
122 {
123 num_ghost = hypre_StructVectorNumGhost(vector);
124
125 boxes = hypre_StructGridBoxes(grid);
126 data_space = hypre_BoxArrayCreate(hypre_BoxArraySize(boxes));
127
128 hypre_ForBoxI(i, boxes)
129 {
130 box = hypre_BoxArrayBox(boxes, i);
131 data_box = hypre_BoxArrayBox(data_space, i);
132
133 hypre_CopyBox(box, data_box);
134 for (d = 0; d < 3; d++)
135 {
136 hypre_BoxIMinD(data_box, d) -= num_ghost[2*d];
137 hypre_BoxIMaxD(data_box, d) += num_ghost[2*d + 1];
138 }
139 }
140
141 hypre_StructVectorDataSpace(vector) = data_space;
142 }
143
144 /*-----------------------------------------------------------------------
145 * Set up data_indices array and data_size
146 *-----------------------------------------------------------------------*/
147
148 if (hypre_StructVectorDataIndices(vector) == NULL)
149 {
150 data_space = hypre_StructVectorDataSpace(vector);
151 data_indices = hypre_CTAlloc(int, hypre_BoxArraySize(data_space));
152
153 data_size = 0;
154 hypre_ForBoxI(i, data_space)
155 {
156 data_box = hypre_BoxArrayBox(data_space, i);
157
158 data_indices[i] = data_size;
159 data_size += hypre_BoxVolume(data_box);
160 }
161
162 hypre_StructVectorDataIndices(vector) = data_indices;
163 hypre_StructVectorDataSize(vector) = data_size;
164 }
165
166 /*-----------------------------------------------------------------------
167 * Set total number of nonzero coefficients
168 *-----------------------------------------------------------------------*/
169
170 hypre_StructVectorGlobalSize(vector) = hypre_StructGridGlobalSize(grid);
171
172 return ierr;
173}
174
175/*--------------------------------------------------------------------------
176 * hypre_StructVectorInitializeData
177 *--------------------------------------------------------------------------*/
178
179int
180hypre_StructVectorInitializeData( hypre_StructVector *vector,
181 double *data )
182{
183 int ierr = 0;
184
185 hypre_StructVectorData(vector) = data;
186 hypre_StructVectorDataAlloced(vector) = 0;
187
188 return ierr;
189}
190
191/*--------------------------------------------------------------------------
192 * hypre_StructVectorInitialize
193 *--------------------------------------------------------------------------*/
194
195int
196hypre_StructVectorInitialize( hypre_StructVector *vector )
197{
198 int ierr = 0;
199
200 double *data;
201
202 ierr = hypre_StructVectorInitializeShell(vector);
203
204 data = hypre_SharedCTAlloc(double, hypre_StructVectorDataSize(vector));
205 hypre_StructVectorInitializeData(vector, data);
206 hypre_StructVectorDataAlloced(vector) = 1;
207
208 return ierr;
209}
210
211/*--------------------------------------------------------------------------
212 * hypre_StructVectorSetValues
213 *--------------------------------------------------------------------------*/
214
215int
216hypre_StructVectorSetValues( hypre_StructVector *vector,
217 hypre_Index grid_index,
218 double values,
219 int add_to )
220{
221 int ierr = 0;
222
223 MPI_Comm comm= hypre_StructVectorComm(vector);
224 hypre_BoxArray *boxes;
225 hypre_Box *box;
226
227 double *vecp;
228
229 int i, found;
230 int true = 1;
231 int false= 0;
232
233 int nprocs;
234
235 MPI_Comm_size(comm, &nprocs);
236
237 boxes = hypre_StructGridBoxes(hypre_StructVectorGrid(vector));
238
239 found= false;
240 hypre_ForBoxI(i, boxes)
241 {
242 box = hypre_BoxArrayBox(boxes, i);
243
244 if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(box)) &&
245 (hypre_IndexX(grid_index) <= hypre_BoxIMaxX(box)) &&
246 (hypre_IndexY(grid_index) >= hypre_BoxIMinY(box)) &&
247 (hypre_IndexY(grid_index) <= hypre_BoxIMaxY(box)) &&
248 (hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(box)) &&
249 (hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(box)) )
250 {
251 vecp = hypre_StructVectorBoxDataValue(vector, i, grid_index);
252 if (add_to)
253 {
254 *vecp += values;
255 }
256 else
257 {
258 *vecp = values;
259 }
260 found= true;
261 }
262 }
263
264 /* to permit ADD values off myproc, but only glayers away from myproc's
265 grid, use the data_space boxes of vector instead of the grid boxes. */
266 if ((!found) && (add_to) && (nprocs > 1))
267 {
268 hypre_Box *orig_box;
269
270 int *add_num_ghost= hypre_StructVectorAddNumGhost(vector);
271 int j;
272
273 hypre_ForBoxI(i, boxes)
274 {
275 orig_box = hypre_BoxArrayBox(boxes, i);
276 box = hypre_BoxDuplicate(orig_box );
277 for (j= 0; j< 3; j++)
278 {
279 hypre_BoxIMin(box)[j]-= add_num_ghost[2*j];
280 hypre_BoxIMax(box)[j]+= add_num_ghost[2*j+1];
281 }
282
283 if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(box)) &&
284 (hypre_IndexX(grid_index) <= hypre_BoxIMaxX(box)) &&
285 (hypre_IndexY(grid_index) >= hypre_BoxIMinY(box)) &&
286 (hypre_IndexY(grid_index) <= hypre_BoxIMaxY(box)) &&
287 (hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(box)) &&
288 (hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(box)) )
289 {
290 vecp = hypre_StructVectorBoxDataValue(vector, i, grid_index);
291 if (add_to)
292 {
293 *vecp += values;
294 }
295 else
296 {
297 *vecp = values;
298 }
299 found= true;
300 }
301 hypre_BoxDestroy(box);
302
303 if (found) break;
304 }
305
306 /* set OffProcAdd for communication. Note that
307 we have an Add if only one point switches this on. */
308 if (found)
309 {
310 if (add_to)
311 {
312 hypre_StructVectorOffProcAdd(vector)= 1;
313 }
314 }
315 else
316 {
317 printf("not found- grid_index off the extended vector grid\n");
318 }
319 }
320
321 return ierr;
322}
323
324/*--------------------------------------------------------------------------
325 * hypre_StructVectorSetBoxValues
326 *--------------------------------------------------------------------------*/
327
328int
329hypre_StructVectorSetBoxValues( hypre_StructVector *vector,
330 hypre_Box *value_box,
331 double *values,
332 int add_to )
333{
334 int ierr = 0;
335
336 MPI_Comm comm = hypre_StructVectorComm(vector);
337 int *add_num_ghost= hypre_StructVectorAddNumGhost(vector);
338
339 hypre_BoxArray *grid_boxes;
340 hypre_Box *grid_box;
341 hypre_BoxArray *box_array1, *box_array2, *tmp_box_array;
342 hypre_BoxArray *value_boxarray;
343 hypre_BoxArrayArray *box_aarray;
344
345 hypre_Box *box, *tmp_box, *orig_box, *vbox;
346
347 hypre_BoxArray *data_space;
348 hypre_Box *data_box;
349 hypre_IndexRef data_start;
350 hypre_Index data_stride;
351 int datai;
352 double *datap;
353
354 hypre_Box *dval_box;
355 hypre_Index dval_start;
356 hypre_Index dval_stride;
357 int dvali;
358
359 hypre_Index loop_size;
360
361 int i, j, k, vol_vbox, vol_iboxes, vol_offproc;
362 int loopi, loopj, loopk;
363 int nprocs;
364
365 MPI_Comm_size(comm, &nprocs);
366
367 /*-----------------------------------------------------------------------
368 * Set up `box_array' by intersecting `box' with the grid boxes
369 *-----------------------------------------------------------------------*/
370
371 /* Find the intersecting boxes of the grid with value_box. Record the
372 volumes of the intersections for possible off_proc settings. */
373 vol_vbox = hypre_BoxVolume(value_box);
374 vol_iboxes= 0;
375
376 grid_boxes= hypre_StructGridBoxes(hypre_StructVectorGrid(vector));
377 box_array1= hypre_BoxArrayCreate(hypre_BoxArraySize(grid_boxes));
378 box = hypre_BoxCreate();
379 hypre_ForBoxI(i, grid_boxes)
380 {
381 grid_box = hypre_BoxArrayBox(grid_boxes, i);
382 hypre_IntersectBoxes(value_box, grid_box, box);
383 hypre_CopyBox(box, hypre_BoxArrayBox(box_array1, i));
384 vol_iboxes+= hypre_BoxVolume(box);
385 }
386
387 /* Check if possible off_proc setting */
388 vol_offproc= 0;
389 if ((vol_vbox > vol_iboxes) && (add_to) && (nprocs > 1))
390 {
391 box_aarray= hypre_BoxArrayArrayCreate(hypre_BoxArraySize(grid_boxes));
392
393 /* to prevent overlapping intersected boxes, we subtract the intersected
394 boxes from value_box. This requires a box_array structure. */
395 value_boxarray= hypre_BoxArrayCreate(0);
396 hypre_AppendBox(value_box, value_boxarray);
397
398 hypre_ForBoxI(i, grid_boxes)
399 {
400 tmp_box_array= hypre_BoxArrayCreate(0);
401
402 /* get ghostlayer boxes */
403 orig_box= hypre_BoxArrayBox(grid_boxes, i);
404 tmp_box = hypre_BoxDuplicate(orig_box );
405 for (j= 0; j< 3; j++)
406 {
407 hypre_BoxIMin(tmp_box)[j]-= add_num_ghost[2*j];
408 hypre_BoxIMax(tmp_box)[j]+= add_num_ghost[2*j+1];
409 }
410 hypre_SubtractBoxes(tmp_box, orig_box, tmp_box_array);
411 hypre_BoxDestroy(tmp_box);
412
413 box_array2= hypre_BoxArrayArrayBoxArray(box_aarray, i);
414 /* intersect the value_box and the ghostlayer boxes */
415 hypre_ForBoxI(j, tmp_box_array)
416 {
417 tmp_box= hypre_BoxArrayBox(tmp_box_array, j);
418 hypre_ForBoxI(k, value_boxarray)
419 {
420 vbox= hypre_BoxArrayBox(value_boxarray, k);
421 hypre_IntersectBoxes(vbox, tmp_box, box);
422 hypre_AppendBox(box, box_array2);
423
424 vol_offproc+= hypre_BoxVolume(box);
425 }
426 }
427
428 /* eliminate intersected boxes so that we do not get overlapping */
429 hypre_SubtractBoxArrays(value_boxarray, box_array2, tmp_box_array);
430 hypre_BoxArrayDestroy(tmp_box_array);
431
432 } /* hypre_ForBoxI(i, grid_boxes) */
433
434 /* if vol_offproc= 0, trying to set values away from ghostlayer */
435 if (!vol_offproc)
436 {
437 hypre_BoxArrayArrayDestroy(box_aarray);
438 }
439 else
440 {
441 /* set OffProcAdd for communication. Note that
442 we have an Add if only one point switches this on. */
443 hypre_StructVectorOffProcAdd(vector)= 1;
444 }
445 hypre_BoxArrayDestroy(value_boxarray);
446
447 }
448 hypre_BoxDestroy(box);
449
450 /*-----------------------------------------------------------------------
451 * Set the vector coefficients
452 *-----------------------------------------------------------------------*/
453
454 if (box_array1)
455 {
456 data_space = hypre_StructVectorDataSpace(vector);
457 hypre_SetIndex(data_stride, 1, 1, 1);
458
459 dval_box = hypre_BoxDuplicate(value_box);
460 hypre_SetIndex(dval_stride, 1, 1, 1);
461
462 hypre_ForBoxI(i, box_array1)
463 {
464 box = hypre_BoxArrayBox(box_array1, i);
465 data_box = hypre_BoxArrayBox(data_space, i);
466
467 /* if there was an intersection */
468 if (box)
469 {
470 data_start = hypre_BoxIMin(box);
471 hypre_CopyIndex(data_start, dval_start);
472
473 datap = hypre_StructVectorBoxData(vector, i);
474
475 hypre_BoxGetSize(box, loop_size);
476
477 if (add_to)
478 {
479 hypre_BoxLoop2Begin(loop_size,
480 data_box,data_start,data_stride,datai,
481 dval_box,dval_start,dval_stride,dvali);
482
483 hypre_BoxLoop2For(loopi, loopj, loopk, datai, dvali)
484 {
485 datap[datai] += values[dvali];
486 }
487 hypre_BoxLoop2End(datai, dvali);
488 }
489 else
490 {
491 hypre_BoxLoop2Begin(loop_size,
492 data_box,data_start,data_stride,datai,
493 dval_box,dval_start,dval_stride,dvali);
494
495 hypre_BoxLoop2For(loopi, loopj, loopk, datai, dvali)
496 {
497 datap[datai] = values[dvali];
498 }
499 hypre_BoxLoop2End(datai, dvali);
500 }
501 }
502 }
503
504 hypre_BoxDestroy(dval_box);
505 }
506 hypre_BoxArrayDestroy(box_array1);
507
508 if (vol_offproc) /* nonzero only when adding values. */
509 {
510 data_space = hypre_StructVectorDataSpace(vector);
511 hypre_SetIndex(data_stride, 1, 1, 1);
512
513 dval_box = hypre_BoxDuplicate(value_box);
514 hypre_SetIndex(dval_stride, 1, 1, 1);
515
516 hypre_ForBoxI(i, data_space)
517 {
518 data_box = hypre_BoxArrayBox(data_space, i);
519 box_array2= hypre_BoxArrayArrayBoxArray(box_aarray, i);
520
521 hypre_ForBoxI(j, box_array2)
522 {
523 box= hypre_BoxArrayBox(box_array2, j);
524
525 /* if there was an intersection */
526 if (box)
527 {
528 data_start = hypre_BoxIMin(box);
529 hypre_CopyIndex(data_start, dval_start);
530
531 datap = hypre_StructVectorBoxData(vector, i);
532
533 hypre_BoxGetSize(box, loop_size);
534
535 if (add_to) /* don't really need this conditional */
536 {
537 hypre_BoxLoop2Begin(loop_size,
538 data_box,data_start,data_stride,datai,
539 dval_box,dval_start,dval_stride,dvali);
540
541 hypre_BoxLoop2For(loopi, loopj, loopk, datai, dvali)
542 {
543 datap[datai] += values[dvali];
544 }
545 hypre_BoxLoop2End(datai, dvali);
546 }
547 } /* if (box) */
548 } /* hypre_ForBoxI(j, box_array2) */
549 } /* hypre_ForBoxI(i, data_space) */
550
551 hypre_BoxDestroy(dval_box);
552 hypre_BoxArrayArrayDestroy(box_aarray);
553 }
554
555 return ierr;
556}
557
558/*--------------------------------------------------------------------------
559 * hypre_StructVectorGetValues. OffProc values on the ghostlayer will be
560 * extracted out, and hence, the values_ptr must contain ghostlayers.
561 *--------------------------------------------------------------------------*/
562
563int
564hypre_StructVectorGetValues( hypre_StructVector *vector,
565 hypre_Index grid_index,
566 double *values_ptr )
567{
568 int ierr = 0;
569
570 int *add_num_ghost= hypre_StructVectorAddNumGhost(vector);
571 double values;
572
573 hypre_BoxArray *boxes;
574 hypre_Box *box, *orig_box;
575
576 double *vecp;
577
578 int i, j, found;
579 int true = 1;
580 int false= 0;
581
582 boxes = hypre_StructGridBoxes(hypre_StructVectorGrid(vector));
583
584 /* search first to see if it is in the box. If not then check
585 the ghostlayered boxes. */
586 found= false;
587 hypre_ForBoxI(i, boxes)
588 {
589 box = hypre_BoxArrayBox(boxes, i);
590
591 if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(box)) &&
592 (hypre_IndexX(grid_index) <= hypre_BoxIMaxX(box)) &&
593 (hypre_IndexY(grid_index) >= hypre_BoxIMinY(box)) &&
594 (hypre_IndexY(grid_index) <= hypre_BoxIMaxY(box)) &&
595 (hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(box)) &&
596 (hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(box)) )
597 {
598 vecp = hypre_StructVectorBoxDataValue(vector, i, grid_index);
599 values = *vecp;
600 found= true;
601 }
602 if (found) break;
603 }
604
605 /* now search if on the ghostlayer */
606 if (!found)
607 {
608 hypre_ForBoxI(i, boxes)
609 {
610 orig_box = hypre_BoxArrayBox(boxes, i);
611 box = hypre_BoxDuplicate(orig_box );
612 for (j= 0; j< 3; j++)
613 {
614 hypre_BoxIMin(box)[j]-= add_num_ghost[2*j];
615 hypre_BoxIMax(box)[j]+= add_num_ghost[2*j+1];
616 }
617
618 if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(box)) &&
619 (hypre_IndexX(grid_index) <= hypre_BoxIMaxX(box)) &&
620 (hypre_IndexY(grid_index) >= hypre_BoxIMinY(box)) &&
621 (hypre_IndexY(grid_index) <= hypre_BoxIMaxY(box)) &&
622 (hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(box)) &&
623 (hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(box)) )
624 {
625 vecp = hypre_StructVectorBoxDataValue(vector, i, grid_index);
626 values= *vecp;
627 found= true;
628 }
629 hypre_BoxDestroy(box);
630 if (found) break;
631 }
632 }
633
634 *values_ptr = values;
635
636 return ierr;
637}
638
639/*--------------------------------------------------------------------------
640 * hypre_StructVectorGetBoxValues. Loop over ghostlayer also.
641 *--------------------------------------------------------------------------*/
642
643int
644hypre_StructVectorGetBoxValues( hypre_StructVector *vector,
645 hypre_Box *value_box,
646 double *values )
647{
648 int ierr = 0;
649
650 int *add_num_ghost= hypre_StructVectorAddNumGhost(vector);
651
652 hypre_BoxArray *grid_boxes;
653 hypre_Box *grid_box;
654 hypre_BoxArray *box_array1, *box_array2, *tmp_box_array;
655 hypre_BoxArray *value_boxarray;
656 hypre_BoxArrayArray *box_aarray;
657
658 hypre_Box *box, *tmp_box, *orig_box, *vbox;
659
660 hypre_BoxArray *data_space;
661 hypre_Box *data_box;
662 hypre_IndexRef data_start;
663 hypre_Index data_stride;
664 int datai;
665 double *datap;
666
667 hypre_Box *dval_box;
668 hypre_Index dval_start;
669 hypre_Index dval_stride;
670 int dvali;
671
672 hypre_Index loop_size;
673
674 int i, j, k, vol_vbox, vol_iboxes, vol_offproc;
675 int loopi, loopj, loopk;
676
677 /*-----------------------------------------------------------------------
678 * Set up `box_array' by intersecting `box' with the grid boxes
679 *-----------------------------------------------------------------------*/
680 vol_vbox = hypre_BoxVolume(value_box);
681 vol_iboxes= 0;
682
683 grid_boxes = hypre_StructGridBoxes(hypre_StructVectorGrid(vector));
684 box_array1 = hypre_BoxArrayCreate(hypre_BoxArraySize(grid_boxes));
685 box = hypre_BoxCreate();
686 hypre_ForBoxI(i, grid_boxes)
687 {
688 grid_box = hypre_BoxArrayBox(grid_boxes, i);
689 hypre_IntersectBoxes(value_box, grid_box, box);
690 hypre_CopyBox(box, hypre_BoxArrayBox(box_array1, i));
691 vol_iboxes+= hypre_BoxVolume(box);
692 }
693
694 /* Check if possible off_proc setting */
695 vol_offproc= 0;
696 if (vol_vbox > vol_iboxes)
697 {
698 box_aarray= hypre_BoxArrayArrayCreate(hypre_BoxArraySize(grid_boxes));
699
700 /* to prevent overlapping intersected boxes, we subtract the intersected
701 boxes from value_box. This requires a box_array structure. */
702 value_boxarray= hypre_BoxArrayCreate(0);
703 hypre_AppendBox(value_box, value_boxarray);
704
705 hypre_ForBoxI(i, grid_boxes)
706 {
707 tmp_box_array= hypre_BoxArrayCreate(0);
708
709 /* get ghostlayer boxes */
710 orig_box= hypre_BoxArrayBox(grid_boxes, i);
711 tmp_box = hypre_BoxDuplicate(orig_box );
712 for (j= 0; j< 3; j++)
713 {
714 hypre_BoxIMin(tmp_box)[j]-= add_num_ghost[2*j];
715 hypre_BoxIMax(tmp_box)[j]+= add_num_ghost[2*j+1];
716 }
717 hypre_SubtractBoxes(tmp_box, orig_box, tmp_box_array);
718 hypre_BoxDestroy(tmp_box);
719
720 box_array2= hypre_BoxArrayArrayBoxArray(box_aarray, i);
721 /* intersect the value_box and the ghostlayer boxes */
722 hypre_ForBoxI(j, tmp_box_array)
723 {
724 tmp_box= hypre_BoxArrayBox(tmp_box_array, j);
725 hypre_ForBoxI(k, value_boxarray)
726 {
727 vbox= hypre_BoxArrayBox(value_boxarray, k);
728 hypre_IntersectBoxes(vbox, tmp_box, box);
729 hypre_AppendBox(box, box_array2);
730
731 vol_offproc+= hypre_BoxVolume(box);
732 }
733 }
734
735 /* eliminate intersected boxes so that we do not get overlapping */
736 hypre_SubtractBoxArrays(value_boxarray, box_array2, tmp_box_array);
737 hypre_BoxArrayDestroy(tmp_box_array);
738
739 } /* hypre_ForBoxI(i, grid_boxes) */
740 /* if vol_offproc= 0, trying to set values away from ghostlayer */
741
742 if (!vol_offproc)
743 {
744 hypre_BoxArrayArrayDestroy(box_aarray);
745 }
746 hypre_BoxArrayDestroy(value_boxarray);
747 }
748 hypre_BoxDestroy(box);
749
750 /*-----------------------------------------------------------------------
751 * Get the vector coefficients
752 *-----------------------------------------------------------------------*/
753 if (box_array1)
754 {
755 data_space = hypre_StructVectorDataSpace(vector);
756 hypre_SetIndex(data_stride, 1, 1, 1);
757
758 dval_box = hypre_BoxDuplicate(value_box);
759 hypre_SetIndex(dval_stride, 1, 1, 1);
760
761 hypre_ForBoxI(i, box_array1)
762 {
763 box = hypre_BoxArrayBox(box_array1, i);
764 data_box = hypre_BoxArrayBox(data_space, i);
765
766 /* if there was an intersection */
767 if (box)
768 {
769 data_start = hypre_BoxIMin(box);
770 hypre_CopyIndex(data_start, dval_start);
771
772 datap = hypre_StructVectorBoxData(vector, i);
773
774 hypre_BoxGetSize(box, loop_size);
775
776 hypre_BoxLoop2Begin(loop_size,
777 data_box, data_start, data_stride, datai,
778 dval_box, dval_start, dval_stride, dvali);
779
780 hypre_BoxLoop2For(loopi, loopj, loopk, datai, dvali)
781 {
782 values[dvali] = datap[datai];
783 }
784 hypre_BoxLoop2End(datai, dvali);
785 }
786 }
787
788 hypre_BoxDestroy(dval_box);
789 }
790 hypre_BoxArrayDestroy(box_array1);
791
792 if (vol_offproc)
793 {
794 data_space = hypre_StructVectorDataSpace(vector);
795 hypre_SetIndex(data_stride, 1, 1, 1);
796
797 dval_box = hypre_BoxDuplicate(value_box);
798 hypre_SetIndex(dval_stride, 1, 1, 1);
799
800 hypre_ForBoxI(i, data_space)
801 {
802 data_box = hypre_BoxArrayBox(data_space, i);
803 box_array2= hypre_BoxArrayArrayBoxArray(box_aarray, i);
804
805 hypre_ForBoxI(j, box_array2)
806 {
807 box= hypre_BoxArrayBox(box_array2, j);
808
809 /* if there was an intersection */
810 if (box)
811 {
812 data_start = hypre_BoxIMin(box);
813 hypre_CopyIndex(data_start, dval_start);
814
815 datap = hypre_StructVectorBoxData(vector, i);
816
817 hypre_BoxGetSize(box, loop_size);
818 hypre_BoxLoop2Begin(loop_size,
819 data_box, data_start, data_stride, datai,
820 dval_box, dval_start, dval_stride, dvali);
821
822 hypre_BoxLoop2For(loopi, loopj, loopk, datai, dvali)
823 {
824 values[dvali] = datap[datai];
825 }
826 hypre_BoxLoop2End(datai, dvali);
827 } /* if (box) */
828 } /* hypre_ForBoxI(j, box_array2) */
829 } /* hypre_ForBoxI(i, data_space) */
830
831 hypre_BoxDestroy(dval_box);
832 hypre_BoxArrayArrayDestroy(box_aarray);
833 }
834
835 return ierr;
836}
837
838/*--------------------------------------------------------------------------
839 * hypre_StructVectorSetNumGhost
840 *--------------------------------------------------------------------------*/
841
842int
843hypre_StructVectorSetNumGhost( hypre_StructVector *vector,
844 int *num_ghost )
845{
846 int ierr = 0;
847 int i;
848
849 for (i = 0; i < 6; i++)
850 hypre_StructVectorNumGhost(vector)[i] = num_ghost[i];
851
852 return ierr;
853}
854
855
856/*--------------------------------------------------------------------------
857 * hypre_StructVectorAssemble
858 * Before assembling the vector, all vector values added from
859 * off_procs are communicated to the correct proc. However, note that since
860 * the comm_pkg is created from an "inverted" comm_info derived from the
861 * vector, not all the communicated data is valid (i.e., we did not mark
862 * which values are actually set off_proc). Because the communicated values
863 * are added to existing values, the user is assumed to have set the
864 * values correctly on or off the proc.
865 *--------------------------------------------------------------------------*/
866
867int
868hypre_StructVectorAssemble( hypre_StructVector *vector )
869{
870 int ierr = 0;
871
872 int sum_OffProcAdd;
873 int OffProcAdd= hypre_StructVectorOffProcAdd(vector);
874
875 /* add_values may be off-proc. Communication needed, which is triggered
876 if one of the OffProcAdd is 1 */
877 sum_OffProcAdd= 0;
878 MPI_Allreduce(&OffProcAdd, &sum_OffProcAdd, 1, MPI_INT, MPI_SUM,
879 hypre_StructVectorComm(vector));
880
881 if (sum_OffProcAdd)
882 {
883 /* since the off_proc add_values are located on the ghostlayer, we
884 need "inverse" communication. */
885
886 hypre_CommInfo *comm_info;
887 hypre_CommInfo *inv_comm_info;
888 hypre_CommPkg *comm_pkg;
889 int *num_ghost = hypre_StructVectorAddNumGhost(vector);
890
891 hypre_BoxArrayArray *send_boxes;
892 hypre_BoxArrayArray *recv_boxes;
893 int **send_procs;
894 int **recv_procs;
895 int **send_rboxnums;
896 int **recv_rboxnums;
897 hypre_BoxArrayArray *send_rboxes;
898 hypre_BoxArray *box_array, *recv_array;
899
900 double *data;
901 hypre_CommHandle *comm_handle;
902
903 hypre_Box *data_box;
904 hypre_Box *box;
905
906 double *data_vec;
907 double *comm_data;
908 hypre_Index loop_size;
909 hypre_IndexRef start;
910 hypre_Index unit_stride;
911
912 int i, j, xi, loopi, loopj, loopk;
913
914 hypre_CreateCommInfoFromNumGhost(hypre_StructVectorGrid(vector),
915 num_ghost, &comm_info);
916
917 /* inverse CommInfo achieved by switching send & recv structures of
918 CommInfo */
919 send_boxes= hypre_BoxArrayArrayDuplicate(hypre_CommInfoRecvBoxes(comm_info));
920 recv_boxes= hypre_BoxArrayArrayDuplicate(hypre_CommInfoSendBoxes(comm_info));
921 send_rboxes= hypre_BoxArrayArrayDuplicate(send_boxes);
922
923 send_procs= hypre_CTAlloc(int *, hypre_BoxArrayArraySize(send_boxes));
924 recv_procs= hypre_CTAlloc(int *, hypre_BoxArrayArraySize(recv_boxes));
925 send_rboxnums= hypre_CTAlloc(int *, hypre_BoxArrayArraySize(send_boxes));
926 recv_rboxnums= hypre_CTAlloc(int *, hypre_BoxArrayArraySize(recv_boxes));
927
928 hypre_ForBoxArrayI(i, send_boxes)
929 {
930 box_array= hypre_BoxArrayArrayBoxArray(send_boxes, i);
931 send_procs[i]= hypre_CTAlloc(int, hypre_BoxArraySize(box_array));
932 memcpy(send_procs[i], hypre_CommInfoRecvProcesses(comm_info)[i],
933 hypre_BoxArraySize(box_array)*sizeof(int));
934
935 send_rboxnums[i]= hypre_CTAlloc(int, hypre_BoxArraySize(box_array));
936 memcpy(send_rboxnums[i], hypre_CommInfoRecvRBoxnums(comm_info)[i],
937 hypre_BoxArraySize(box_array)*sizeof(int));
938 }
939
940 hypre_ForBoxArrayI(i, recv_boxes)
941 {
942 box_array= hypre_BoxArrayArrayBoxArray(recv_boxes, i);
943 recv_procs[i]= hypre_CTAlloc(int, hypre_BoxArraySize(box_array));
944 memcpy(recv_procs[i], hypre_CommInfoSendProcesses(comm_info)[i],
945 hypre_BoxArraySize(box_array)*sizeof(int));
946
947 recv_rboxnums[i]= hypre_CTAlloc(int, hypre_BoxArraySize(box_array));
948 memcpy(recv_rboxnums[i], hypre_CommInfoSendRBoxnums(comm_info)[i],
949 hypre_BoxArraySize(box_array)*sizeof(int));
950 }
951
952 hypre_CommInfoCreate(send_boxes, recv_boxes, send_procs, recv_procs,
953 send_rboxnums, recv_rboxnums, send_rboxes,
954 &inv_comm_info);
955
956 hypre_CommPkgCreate(inv_comm_info,
957 hypre_StructVectorDataSpace(vector),
958 hypre_StructVectorDataSpace(vector),
959 1,
960 hypre_StructVectorComm(vector),
961 &comm_pkg);
962
963 /* communicate the add value entries */
964 data= hypre_CTAlloc(double, hypre_StructVectorDataSize(vector));
965 hypre_InitializeCommunication(comm_pkg,
966 hypre_StructVectorData(vector),
967 data,
968 &comm_handle);
969
970 hypre_FinalizeCommunication(comm_handle);
971
972 /* this proc will recved data in it's send_boxes of comm_info, or
973 equivalently, in the recv_boxes of inv_comm_info. Since inv_comm_info
974 has already been destroyed, we use the send_boxes of comm_info */
975 hypre_SetIndex(unit_stride, 1, 1, 1);
976 box_array= hypre_StructGridBoxes(hypre_StructVectorGrid(vector));
977 hypre_ForBoxI(i, box_array)
978 {
979 recv_array=
980 hypre_BoxArrayArrayBoxArray(hypre_CommInfoSendBoxes(comm_info), i);
981
982 data_box = hypre_BoxArrayBox(hypre_StructVectorDataSpace(vector), i);
983 data_vec = hypre_StructVectorBoxData(vector, i);
984 comm_data=(data + hypre_StructVectorDataIndices(vector)[i]);
985
986 hypre_ForBoxI(j, recv_array)
987 {
988 box = hypre_BoxArrayBox(recv_array, j);
989 start= hypre_BoxIMin(box);
990 hypre_BoxGetSize(box, loop_size);
991
992 /* note that every proc adds since we don't track which
993 ones should. */
994 hypre_BoxLoop1Begin(loop_size,
995 data_box, start, unit_stride, xi)
996
997 hypre_BoxLoop1For(loopi, loopj, loopk, xi)
998 {
999 data_vec[xi] += comm_data[xi];
1000 }
1001 hypre_BoxLoop1End(xi);
1002
1003 } /* hypre_ForBoxI(j, recv_array) */
1004 } /* hypre_ForBoxI(i, box_array) */
1005
1006 hypre_TFree(data);
1007 hypre_CommInfoDestroy(comm_info);
1008 hypre_CommPkgDestroy(comm_pkg);
1009
1010 }
1011
1012 return ierr;
1013}
1014
1015/*--------------------------------------------------------------------------
1016 * hypre_StructVectorCopy
1017 * copies data from x to y
1018 * y has its own data array, so this is a deep copy in that sense.
1019 * The grid and other size information are not copied - they are
1020 * assumed to have already been set up to be consistent.
1021 *--------------------------------------------------------------------------*/
1022
1023int
1024hypre_StructVectorCopy( hypre_StructVector *x,
1025 hypre_StructVector *y )
1026{
1027
1028
1029 int ierr = 0;
1030
1031 hypre_Box *x_data_box;
1032
1033 int vi;
1034 double *xp, *yp;
1035
1036 hypre_BoxArray *boxes;
1037 hypre_Box *box;
1038 hypre_Index loop_size;
1039 hypre_IndexRef start;
1040 hypre_Index unit_stride;
1041
1042 int i;
1043 int loopi, loopj, loopk;
1044
1045 /*-----------------------------------------------------------------------
1046 * Set the vector coefficients
1047 *-----------------------------------------------------------------------*/
1048
1049 hypre_SetIndex(unit_stride, 1, 1, 1);
1050
1051 boxes = hypre_StructGridBoxes( hypre_StructVectorGrid(x) );
1052 hypre_ForBoxI(i, boxes)
1053 {
1054 box = hypre_BoxArrayBox(boxes, i);
1055 start = hypre_BoxIMin(box);
1056
1057 x_data_box =
1058 hypre_BoxArrayBox(hypre_StructVectorDataSpace(x), i);
1059 xp = hypre_StructVectorBoxData(x, i);
1060 yp = hypre_StructVectorBoxData(y, i);
1061
1062 hypre_BoxGetSize(box, loop_size);
1063
1064 hypre_BoxLoop1Begin(loop_size,
1065 x_data_box, start, unit_stride, vi);
1066
1067 hypre_BoxLoop1For(loopi, loopj, loopk, vi)
1068 {
1069 yp[vi] = xp[vi];
1070 }
1071 hypre_BoxLoop1End(vi);
1072 }
1073
1074 return ierr;
1075}
1076
1077/*--------------------------------------------------------------------------
1078 * hypre_StructVectorSetConstantValues
1079 *--------------------------------------------------------------------------*/
1080
1081int
1082hypre_StructVectorSetConstantValues( hypre_StructVector *vector,
1083 double values )
1084{
1085 int ierr = 0;
1086
1087 hypre_Box *v_data_box;
1088
1089 int vi;
1090 double *vp;
1091
1092 hypre_BoxArray *boxes;
1093 hypre_Box *box;
1094 hypre_Index loop_size;
1095 hypre_IndexRef start;
1096 hypre_Index unit_stride;
1097
1098 int i;
1099 int loopi, loopj, loopk;
1100
1101 /*-----------------------------------------------------------------------
1102 * Set the vector coefficients
1103 *-----------------------------------------------------------------------*/
1104
1105 hypre_SetIndex(unit_stride, 1, 1, 1);
1106
1107 boxes = hypre_StructGridBoxes(hypre_StructVectorGrid(vector));
1108 hypre_ForBoxI(i, boxes)
1109 {
1110 box = hypre_BoxArrayBox(boxes, i);
1111 start = hypre_BoxIMin(box);
1112
1113 v_data_box =
1114 hypre_BoxArrayBox(hypre_StructVectorDataSpace(vector), i);
1115 vp = hypre_StructVectorBoxData(vector, i);
1116
1117 hypre_BoxGetSize(box, loop_size);
1118
1119 hypre_BoxLoop1Begin(loop_size,
1120 v_data_box, start, unit_stride, vi);
1121
1122 hypre_BoxLoop1For(loopi, loopj, loopk, vi)
1123 {
1124 vp[vi] = values;
1125 }
1126 hypre_BoxLoop1End(vi);
1127 }
1128
1129 return ierr;
1130}
1131
1132/*--------------------------------------------------------------------------
1133 * hypre_StructVectorSetFunctionValues
1134 *
1135 * Takes a function pointer of the form:
1136 *
1137 * double f(i,j,k)
1138 *--------------------------------------------------------------------------*/
1139
1140int
1141hypre_StructVectorSetFunctionValues( hypre_StructVector *vector,
1142 double (*fcn)(int, int, int) )
1143{
1144 int ierr = 0;
1145
1146 hypre_Box *v_data_box;
1147
1148 int vi;
1149 double *vp;
1150
1151 hypre_BoxArray *boxes;
1152 hypre_Box *box;
1153 hypre_Index loop_size;
1154 hypre_IndexRef start;
1155 hypre_Index unit_stride;
1156
1157 int b, i, j, k;
1158 int loopi, loopj, loopk;
1159
1160 /*-----------------------------------------------------------------------
1161 * Set the vector coefficients
1162 *-----------------------------------------------------------------------*/
1163
1164 hypre_SetIndex(unit_stride, 1, 1, 1);
1165
1166 boxes = hypre_StructGridBoxes(hypre_StructVectorGrid(vector));
1167 hypre_ForBoxI(b, boxes)
1168 {
1169 box = hypre_BoxArrayBox(boxes, b);
1170 start = hypre_BoxIMin(box);
1171
1172 v_data_box =
1173 hypre_BoxArrayBox(hypre_StructVectorDataSpace(vector), b);
1174 vp = hypre_StructVectorBoxData(vector, b);
1175
1176 hypre_BoxGetSize(box, loop_size);
1177
1178 hypre_BoxLoop1Begin(loop_size,
1179 v_data_box, start, unit_stride, vi);
1180 i = hypre_IndexX(start);
1181 j = hypre_IndexY(start);
1182 k = hypre_IndexZ(start);
1183
1184 hypre_BoxLoop1For(loopi, loopj, loopk, vi)
1185 {
1186 vp[vi] = fcn(i, j, k);
1187 i++;
1188 j++;
1189 k++;
1190 }
1191 hypre_BoxLoop1End(vi);
1192 }
1193
1194 return ierr;
1195}
1196
1197/*--------------------------------------------------------------------------
1198 * hypre_StructVectorClearGhostValues
1199 *--------------------------------------------------------------------------*/
1200
1201int
1202hypre_StructVectorClearGhostValues( hypre_StructVector *vector )
1203{
1204 int ierr = 0;
1205
1206 hypre_Box *v_data_box;
1207
1208 int vi;
1209 double *vp;
1210
1211 hypre_BoxArray *boxes;
1212 hypre_Box *box;
1213 hypre_BoxArray *diff_boxes;
1214 hypre_Box *diff_box;
1215 hypre_Index loop_size;
1216 hypre_IndexRef start;
1217 hypre_Index unit_stride;
1218
1219 int i, j;
1220 int loopi, loopj, loopk;
1221
1222 /*-----------------------------------------------------------------------
1223 * Set the vector coefficients
1224 *-----------------------------------------------------------------------*/
1225
1226 hypre_SetIndex(unit_stride, 1, 1, 1);
1227
1228 boxes = hypre_StructGridBoxes(hypre_StructVectorGrid(vector));
1229 diff_boxes = hypre_BoxArrayCreate(0);
1230 hypre_ForBoxI(i, boxes)
1231 {
1232 box = hypre_BoxArrayBox(boxes, i);
1233
1234 v_data_box =
1235 hypre_BoxArrayBox(hypre_StructVectorDataSpace(vector), i);
1236 vp = hypre_StructVectorBoxData(vector, i);
1237
1238 hypre_BoxArraySetSize(diff_boxes, 0);
1239 hypre_SubtractBoxes(v_data_box, box, diff_boxes);
1240 hypre_ForBoxI(j, diff_boxes)
1241 {
1242 diff_box = hypre_BoxArrayBox(diff_boxes, j);
1243 start = hypre_BoxIMin(diff_box);
1244
1245 hypre_BoxGetSize(diff_box, loop_size);
1246
1247 hypre_BoxLoop1Begin(loop_size,
1248 v_data_box, start, unit_stride, vi);
1249
1250 hypre_BoxLoop1For(loopi, loopj, loopk, vi)
1251 {
1252 vp[vi] = 0.0;
1253 }
1254 hypre_BoxLoop1End(vi);
1255 }
1256 }
1257 hypre_BoxArrayDestroy(diff_boxes);
1258
1259 return ierr;
1260}
1261
1262/*--------------------------------------------------------------------------
1263 * hypre_StructVectorClearBoundGhostValues
1264 * clears vector values on the physical boundaries
1265 *--------------------------------------------------------------------------*/
1266
1267int
1268hypre_StructVectorClearBoundGhostValues( hypre_StructVector *vector )
1269{
1270 int ierr = 0;
1271 int vi;
1272 double *vp;
1273 hypre_BoxArray *boxes;
1274 hypre_Box *box;
1275 hypre_Box *v_data_box;
1276 hypre_Index loop_size;
1277 hypre_IndexRef start;
1278 hypre_Index stride;
1279 hypre_Box *bbox;
1280 hypre_StructGrid *grid;
1281 hypre_BoxArray *boundary_boxes;
1282 hypre_BoxArray *array_of_box;
1283 hypre_BoxArray *work_boxarray;
1284
1285 int i, i2;
1286 int loopi, loopj, loopk;
1287
1288 /*-----------------------------------------------------------------------
1289 * Set the vector coefficients
1290 *-----------------------------------------------------------------------*/
1291
1292 grid = hypre_StructVectorGrid(vector);
1293 boxes = hypre_StructGridBoxes(grid);
1294 hypre_SetIndex(stride, 1, 1, 1);
1295
1296 hypre_ForBoxI(i, boxes)
1297 {
1298 box = hypre_BoxArrayBox(boxes, i);
1299 boundary_boxes = hypre_BoxArrayCreate( 0 );
1300 v_data_box =
1301 hypre_BoxArrayBox(hypre_StructVectorDataSpace(vector), i);
1302 ierr += hypre_BoxBoundaryG( v_data_box, grid, boundary_boxes );
1303 vp = hypre_StructVectorBoxData(vector, i);
1304
1305 /* box is a grid box, no ghost zones.
1306 v_data_box is vector data box, may or may not have ghost zones
1307 To get only ghost zones, subtract box from boundary_boxes. */
1308 work_boxarray = hypre_BoxArrayCreate( 0 );
1309 array_of_box = hypre_BoxArrayCreate( 1 );
1310 hypre_BoxArrayBoxes(array_of_box)[0] = *box;
1311 hypre_SubtractBoxArrays( boundary_boxes, array_of_box, work_boxarray );
1312
1313 hypre_ForBoxI(i2, boundary_boxes)
1314 {
1315 bbox = hypre_BoxArrayBox(boundary_boxes, i2);
1316 hypre_BoxGetSize(bbox, loop_size);
1317 start = hypre_BoxIMin(bbox);
1318 hypre_BoxLoop1Begin(loop_size,
1319 v_data_box, start, stride, vi);
1320
1321 hypre_BoxLoop1For(loopi, loopj, loopk, vi)
1322 {
1323 vp[vi] = 0.0;
1324 }
1325 hypre_BoxLoop1End(vi);
1326 }
1327 hypre_BoxArrayDestroy(boundary_boxes);
1328 hypre_BoxArrayDestroy(work_boxarray);
1329 hypre_BoxArrayDestroy(array_of_box);
1330 }
1331
1332 return ierr;
1333}
1334
1335/*--------------------------------------------------------------------------
1336 * hypre_StructVectorScaleValues
1337 *--------------------------------------------------------------------------*/
1338
1339int
1340hypre_StructVectorScaleValues( hypre_StructVector *vector, double factor )
1341{
1342 int ierr = 0;
1343
1344 int datai;
1345 double *data;
1346
1347 hypre_Index imin;
1348 hypre_Index imax;
1349 hypre_Box *box;
1350 hypre_Index loop_size;
1351
1352 int loopi, loopj, loopk;
1353
1354 /*-----------------------------------------------------------------------
1355 * Set the vector coefficients
1356 *-----------------------------------------------------------------------*/
1357
1358 box = hypre_BoxCreate();
1359 hypre_SetIndex(imin, 1, 1, 1);
1360 hypre_SetIndex(imax, hypre_StructVectorDataSize(vector), 1, 1);
1361 hypre_BoxSetExtents(box, imin, imax);
1362 data = hypre_StructVectorData(vector);
1363 hypre_BoxGetSize(box, loop_size);
1364
1365 hypre_BoxLoop1Begin(loop_size,
1366 box, imin, imin, datai);
1367
1368 hypre_BoxLoop1For(loopi, loopj, loopk, datai)
1369 {
1370 data[datai] *= factor;
1371 }
1372 hypre_BoxLoop1End(datai);
1373
1374 hypre_BoxDestroy(box);
1375
1376 return ierr;
1377}
1378
1379/*--------------------------------------------------------------------------
1380 * hypre_StructVectorClearAllValues
1381 *--------------------------------------------------------------------------*/
1382
1383int
1384hypre_StructVectorClearAllValues( hypre_StructVector *vector )
1385{
1386 int ierr = 0;
1387
1388 int datai;
1389 double *data;
1390
1391 hypre_Index imin;
1392 hypre_Index imax;
1393 hypre_Box *box;
1394 hypre_Index loop_size;
1395
1396 int loopi, loopj, loopk;
1397
1398 /*-----------------------------------------------------------------------
1399 * Set the vector coefficients
1400 *-----------------------------------------------------------------------*/
1401
1402 box = hypre_BoxCreate();
1403 hypre_SetIndex(imin, 1, 1, 1);
1404 hypre_SetIndex(imax, hypre_StructVectorDataSize(vector), 1, 1);
1405 hypre_BoxSetExtents(box, imin, imax);
1406 data = hypre_StructVectorData(vector);
1407 hypre_BoxGetSize(box, loop_size);
1408
1409 hypre_BoxLoop1Begin(loop_size,
1410 box, imin, imin, datai);
1411
1412 hypre_BoxLoop1For(loopi, loopj, loopk, datai)
1413 {
1414 data[datai] = 0.0;
1415 }
1416 hypre_BoxLoop1End(datai);
1417
1418 hypre_BoxDestroy(box);
1419
1420 return ierr;
1421}
1422
1423/*--------------------------------------------------------------------------
1424 * hypre_StructVectorGetMigrateCommPkg
1425 *--------------------------------------------------------------------------*/
1426
1427hypre_CommPkg *
1428hypre_StructVectorGetMigrateCommPkg( hypre_StructVector *from_vector,
1429 hypre_StructVector *to_vector )
1430{
1431 hypre_CommInfo *comm_info;
1432 hypre_CommPkg *comm_pkg;
1433
1434 /*------------------------------------------------------
1435 * Set up hypre_CommPkg
1436 *------------------------------------------------------*/
1437
1438 hypre_CreateCommInfoFromGrids(hypre_StructVectorGrid(from_vector),
1439 hypre_StructVectorGrid(to_vector),
1440 &comm_info);
1441 hypre_CommPkgCreate(comm_info,
1442 hypre_StructVectorDataSpace(from_vector),
1443 hypre_StructVectorDataSpace(to_vector), 1,
1444 hypre_StructVectorComm(from_vector), &comm_pkg);
1445 /* is this correct for periodic? */
1446
1447 return comm_pkg;
1448}
1449
1450/*--------------------------------------------------------------------------
1451 * hypre_StructVectorMigrate
1452 *--------------------------------------------------------------------------*/
1453
1454int
1455hypre_StructVectorMigrate( hypre_CommPkg *comm_pkg,
1456 hypre_StructVector *from_vector,
1457 hypre_StructVector *to_vector )
1458{
1459 hypre_CommHandle *comm_handle;
1460
1461 int ierr = 0;
1462
1463 /*-----------------------------------------------------------------------
1464 * Migrate the vector data
1465 *-----------------------------------------------------------------------*/
1466
1467 hypre_InitializeCommunication(comm_pkg,
1468 hypre_StructVectorData(from_vector),
1469 hypre_StructVectorData(to_vector),
1470 &comm_handle);
1471 hypre_FinalizeCommunication(comm_handle);
1472
1473 return ierr;
1474}
1475
1476/*--------------------------------------------------------------------------
1477 * hypre_StructVectorPrint
1478 *--------------------------------------------------------------------------*/
1479
1480int
1481hypre_StructVectorPrint( const char *filename,
1482 hypre_StructVector *vector,
1483 int all )
1484{
1485 int ierr = 0;
1486
1487 FILE *file;
1488 char new_filename[255];
1489
1490 hypre_StructGrid *grid;
1491 hypre_BoxArray *boxes;
1492
1493 hypre_BoxArray *data_space;
1494
1495 int myid;
1496
1497 /*----------------------------------------
1498 * Open file
1499 *----------------------------------------*/
1500
1501 MPI_Comm_rank(hypre_StructVectorComm(vector), &myid );
1502
1503 sprintf(new_filename, "%s.%05d", filename, myid);
1504
1505 if ((file = fopen(new_filename, "w")) == NULL)
1506 {
1507 printf("Error: can't open output file %s\n", new_filename);
1508 exit(1);
1509 }
1510
1511 /*----------------------------------------
1512 * Print header info
1513 *----------------------------------------*/
1514
1515 fprintf(file, "StructVector\n");
1516
1517 /* print grid info */
1518 fprintf(file, "\nGrid:\n");
1519 grid = hypre_StructVectorGrid(vector);
1520 hypre_StructGridPrint(file, grid);
1521
1522 /*----------------------------------------
1523 * Print data
1524 *----------------------------------------*/
1525
1526 data_space = hypre_StructVectorDataSpace(vector);
1527
1528 if (all)
1529 boxes = data_space;
1530 else
1531 boxes = hypre_StructGridBoxes(grid);
1532
1533 fprintf(file, "\nData:\n");
1534 hypre_PrintBoxArrayData(file, boxes, data_space, 1,
1535 hypre_StructVectorData(vector));
1536
1537 /*----------------------------------------
1538 * Close file
1539 *----------------------------------------*/
1540
1541 fflush(file);
1542 fclose(file);
1543
1544 return ierr;
1545}
1546
1547/*--------------------------------------------------------------------------
1548 * hypre_StructVectorRead
1549 *--------------------------------------------------------------------------*/
1550
1551hypre_StructVector *
1552hypre_StructVectorRead( MPI_Comm comm,
1553 const char *filename,
1554 int *num_ghost )
1555{
1556 FILE *file;
1557 char new_filename[255];
1558
1559 hypre_StructVector *vector;
1560
1561 hypre_StructGrid *grid;
1562 hypre_BoxArray *boxes;
1563
1564 hypre_BoxArray *data_space;
1565
1566 int myid;
1567
1568 /*----------------------------------------
1569 * Open file
1570 *----------------------------------------*/
1571
1572#ifdef HYPRE_USE_PTHREADS
1573#if MPI_Comm_rank == hypre_thread_MPI_Comm_rank
1574#undef MPI_Comm_rank
1575#endif
1576#endif
1577
1578 MPI_Comm_rank(comm, &myid );
1579
1580 sprintf(new_filename, "%s.%05d", filename, myid);
1581
1582 if ((file = fopen(new_filename, "r")) == NULL)
1583 {
1584 printf("Error: can't open output file %s\n", new_filename);
1585 exit(1);
1586 }
1587
1588 /*----------------------------------------
1589 * Read header info
1590 *----------------------------------------*/
1591
1592 fscanf(file, "StructVector\n");
1593
1594 /* read grid info */
1595 fscanf(file, "\nGrid:\n");
1596 hypre_StructGridRead(comm,file,&grid);
1597
1598 /*----------------------------------------
1599 * Initialize the vector
1600 *----------------------------------------*/
1601
1602 vector = hypre_StructVectorCreate(comm, grid);
1603 hypre_StructVectorSetNumGhost(vector, num_ghost);
1604 hypre_StructVectorInitialize(vector);
1605
1606 /*----------------------------------------
1607 * Read data
1608 *----------------------------------------*/
1609
1610 boxes = hypre_StructGridBoxes(grid);
1611 data_space = hypre_StructVectorDataSpace(vector);
1612
1613 fscanf(file, "\nData:\n");
1614 hypre_ReadBoxArrayData(file, boxes, data_space, 1,
1615 hypre_StructVectorData(vector));
1616
1617 /*----------------------------------------
1618 * Assemble the vector
1619 *----------------------------------------*/
1620
1621 hypre_StructVectorAssemble(vector);
1622
1623 /*----------------------------------------
1624 * Close file
1625 *----------------------------------------*/
1626
1627 fclose(file);
1628
1629 return vector;
1630}
1631
1632/*--------------------------------------------------------------------------
1633 * The following is used only as a debugging aid.
1634 *--------------------------------------------------------------------------*/
1635
1636int
1637hypre_StructVectorMaxValue( hypre_StructVector *vector,
1638 double *max_value, int *max_index,
1639 hypre_Index max_xyz_index )
1640/* Input: vector, and pointers to where to put returned data.
1641 Return value: error flag, 0 means ok.
1642 Finds the maximum value in a vector, puts it in max_value.
1643 The corresponding index is put in max_index.
1644 A hypre_Index corresponding to max_index is put in max_xyz_index.
1645 We assume that there is only one box to deal with. */
1646{
1647 int ierr = 0;
1648
1649 int datai;
1650 double *data;
1651
1652 hypre_Index imin;
1653 hypre_BoxArray *boxes;
1654 hypre_Box *box;
1655 hypre_Index loop_size;
1656 hypre_Index unit_stride;
1657
1658 int loopi, loopj, loopk, i;
1659 double maxvalue;
1660 int maxindex;
1661
1662 boxes = hypre_StructVectorDataSpace(vector);
1663 if ( hypre_BoxArraySize(boxes)!=1 ) {
1664 /* if more than one box, the return system max_xyz_index is too simple
1665 if needed, fix later */
1666 ierr = 1;
1667 return ierr;
1668 }
1669 hypre_SetIndex(unit_stride, 1, 1, 1);
1670 hypre_ForBoxI(i, boxes)
1671 {
1672 box = hypre_BoxArrayBox(boxes, i);
1673 /*v_data_box =
1674 hypre_BoxArrayBox(hypre_StructVectorDataSpace(vector), i);*/
1675 data = hypre_StructVectorBoxData(vector, i);
1676 hypre_BoxGetSize(box, loop_size);
1677 hypre_CopyIndex( hypre_BoxIMin(box), imin );
1678
1679 hypre_BoxLoop1Begin(loop_size,
1680 box, imin, unit_stride, datai);
1681
1682 maxindex = hypre_BoxIndexRank( box, imin );
1683 maxvalue = data[maxindex];
1684 hypre_CopyIndex( imin, max_xyz_index );
1685 hypre_BoxLoop1For(loopi, loopj, loopk, datai)
1686 {
1687 if ( data[datai] > maxvalue )
1688 {
1689 maxvalue = data[datai];
1690 maxindex = datai;
1691 hypre_SetIndex(max_xyz_index, loopi+hypre_IndexX(imin),
1692 loopj+hypre_IndexY(imin),
1693 loopk+hypre_IndexZ(imin) );
1694 }
1695 }
1696 hypre_BoxLoop1End(datai);
1697 }
1698
1699 *max_value = maxvalue;
1700 *max_index = maxindex;
1701
1702 return ierr;
1703}
1704
Note: See TracBrowser for help on using the repository browser.