[petsc-dev] DMDA_*PERIODIC and DMDA_XYZGHOSTED
Ethan Coon
ecoon at lanl.gov
Thu Mar 10 14:30:36 CST 2011
Attached is a patchq for enabling DMDA_*GHOSTED to be used, in which
case ghost nodes are included even at domain boundaries of nonperiodic
dimensions. The DMDAPeriodicType enum was redone, but includes
DMDA_XYPERIODIC/etc to make it compatible with the old version of the
enum. To use, simply compose the DMDAPeriodicType with bitwise or:
DMDA_XPERIODIC | DMDA_YGHOSTED, for example.
This just slightly contributes to the ugliness of the 2d/3d SetUp
methods, but not much. The only place I've introduced negative indices
at this point is into the DMLocalToGlobalMapping, where there must be
global numbers for the local ghost cells which don't exist in the global
vec. Allowing negative indices to the VecScatters would simplify a lot
of this code, but I don't really have time to take on another
fundamental piece of PETSc at this point...
In the process, I've cleaned up all the old stuff where the index sets
were generated in true "dof-strided" indices (instead of block indices),
and removed all the code to patch in that change to ISCreateBlock(). At
the end, all the x-component pieces of the DMDA are multiplied by the #
of dofs, as it seems much of PETSc depends upon this (though it wasn't
clear why except if for historical reasons).
All the dm/examples/tests pass, and I redid the scripts in
dm/examples/tests/scripts so that they both worked and test all
combinations of periodicities/ghosted and stencils. Some regression
tests of my own code work as well. I checked the Interp operators for
MG, and they look fine (and seem to pass existing regression tests), but
that could use some checking. Let me know if I'm missing another set of
tests that I should be running... not sure what the standards are for
contributions like this.
Thanks,
Ethan
On Wed, 2011-03-02 at 18:33 -0600, Barry Smith wrote:
> On Mar 2, 2011, at 6:11 PM, Ethan Coon wrote:
>
> > On Wed, 2011-03-02 at 16:26 -0600, Barry Smith wrote:
> >> On Mar 2, 2011, at 1:47 PM, Ethan Coon wrote:
> >>
> >>>
> >>>>
> >>>> Ethan,
> >>>>
> >>>> MatSetValues() and VecSetValues() handle negative indices as "ignore these entries". Currently VecScatterCreate() does not handle negative indices as "ignore these entries" (at least it is not documented and I did not write it), likely it will either crash or generate an error.
> >>>>
> >>>> It sounds like you are proposing that if the from or to entry in a particular "slot" is negative you would like VecScatterCreate() to just ignore that slot? This seems like an ok proposal if you are willing to update VecScatterCreate() to handle it and add to VecScatterCreate() manual page this feature. If this truly simplifies all the horrible if () code in the DA construction to handle corner stuff then it would be worth doing.
> >>>>
> >>>
> >>> Hmm, I was proposing that, but because I thought that was the case
> >>> already.
> >>>
> >>> It would clean up the DASetUp code, but not as much as I thought
> >>> initially. Currently VecSetValuesLocal(), when used with the L2G
> >>> mapping from a STAR_STENCIL DMDA, will happily add/insert values from
> >>> the ghost cells in the corner to the global vector (why you would add
> >>> values to a ghost node on which you don't want to get information back
> >>> from I don't know, but someone made an effort to implement it...).
> >>
> >> I think that is a "bug" or "unintended feature" I hope nobody worked hard to get it to work. I think it is just that way because that is the way it worked out. Probably DMGetMatrix_DA() should call MatSetOption(mat,MAT_NO_NEW_NONZEROS,PETSC_TRUE); to prevent people of accidently using those slots (if they really want to for some perverse reason they could call MatSetOption() themselves to reset it.
> >>
> >> Barry
> >>
> >>> To
> >>> keep that feature, the L2GMapping must be different from the IS used for
> >>> the DMDAGlobalToLocal scatter, and all the ugly if crap has to stay.
> >>
> >> Even in the star case the L2G has to contain slots for those stencil points (though you could fill those slots with negative entries I guess) to get the VecSetValuesLocal() to work. Is that what you propose, putting negative numbers there?? Seems possibly ok to me. The one danger is that if they user intends to set that corner value with VecSetValuesLocal() or MatSetValuesLocal() it will be silently discarded (of course they should never try to set it but they may).
> >>
> >>
> >>
> >> Barry
> >>
> >
> > Right, got it, there has to be local number for the gxs,gys,gzs entry
> > even in the star case. The question is if it has to get mapped
> > someplace.
> >
> > Using VecSetValuesLocal() on an entire local domain with ADD_VALUES, the
> > current implementation will sum the entries from ghost nodes and
> > interior nodes, while in the INSERT_VALUES, it would lead to a race
> > condition, where the processor who owns the value may not win (and it
> > does so silently). (Note this is true for any stencil and any ghost
> > node, not just corners in the star stencil.)
>
> Using INSERT_VALUES with VecSetValues() also as well as the MatSetValues versions always has the condition that it is undefined who wins when several processes try to put in the same location. This is just a fact of (PETSc) life. I don't think DM or SetValuesLocal() really makes it any worse.
>
>
> >
> > If instead we put -1 in all ghost indices of the l2g mapping, then
> > VecSetValuesLocal() with INSERT_VALUES functions as expected (the value
> > in the global vec is the value given by the process owning the node),
> > while ADD_VALUES will only add in the value given by the process owning
> > the node (and do so silently). This behavior for ADD_VALUES is exactly
> > what is described in the "notes" section of DMLocalToGlobalBegin's man
> > page. Basically it just doesn't allow you to VecSetLocalValues() on
> > ghost nodes.
> >
> > Compare this to the result of the operation: DMDAVecGetArray(da, local),
> > assignment to the array, restore the array, then call
> > DMDALocalToGlobal(). With DMDALocalToGlobal() and INSERT_VALUES, it
> > does the l2g scatter, and so there is no race condition and the global
> > vec gets the value in the array of the process owning that entry. With
> > ADD_VALUES, it does the reverse of the g2l scatter, so all values get
> > added in.
> >
> > I guess I'm more concerned about the undocumented race condition than
> > not adding values from ghost cells, which could easily be explained in a
> > man page. So yes, I'm proposing to put -1 as the global index of all
> > ghost nodes in the local to global mapping.
> >
> > Note that I have to put something in the ghosted (nonperiodic) local
> > number spots as well -- they have no corresponding global number, so it
> > has to be -1. At least this is then consistent that all ghost nodes in
> > a VecSetValuesLocal() get ignored. And if you really want to do this,
> > it's likely you're using the (safer) DMLocalToGlobal() way anyway.
> >
>
> Go ahead and add support for VecScatterCreate() to handle negative indices and simplify (greatly) the DACreates if you are up for it.
>
> Thanks
>
> Barry
>
>
> > Ethan
> >
> >
> >>> I can just graft the Ghosted case on to that code (making it only
> >>> slightly more ugly). It will still depend upon the VecSetValues()
> >>> accepting and ignoring negative global indices, but that's already the
> >>> case.
> >>>
> >>> Ethan
> >>>
> >>>> Performance is not an issue since you would just discard those slots in the VecScatterCreate() phase and they would never appear in the actual scatter operations.
> >>>>
> >>>>
> >>>> Barry
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>> Ethan
> >>>>>
> >>>>>
> >>>>>
> >>>>> On Tue, 2010-12-07 at 15:41 -0700, Ethan Coon wrote:
> >>>>>> On Tue, 2010-12-07 at 13:45 -0600, Barry Smith wrote:
> >>>>>>> DMDA_XYZGHOSTED does not exist for 2d and 3d it was added, I'm guessing, as an experiment and was never in the initial design of DMDA. To fully support it one needs to go back tot he design of DMDA and see how to have it properly done and not just bolt it on. Some people like to use these types of ghost nodes so I agree it is a useful thing to have but who is going to properly add it?
> >>>>>>>
> >>>>>>
> >>>>>> At some point in the not-too-distant future I'll get frustrated enough
> >>>>>> to look into this, but I don't have the time at the moment. At first
> >>>>>> glance it looks like:
> >>>>>>
> >>>>>> - Ensure DMDA{X,Y,Z}Periodic() macros are used everywhere instead of
> >>>>>> direct comparisons to dd->wrap (they aren't used everywhere currently).
> >>>>>>
> >>>>>> - Define macros DMDA{X,Y,Z}Ghosted() to (in some places) replace
> >>>>>> DMDA{X,Y,Z}Periodic() and then choosing the appropriate macro in the
> >>>>>> right places.
> >>>>>>
> >>>>>> - This probably doesn't merit a change in the DMDACreate* API (it would
> >>>>>> affect a very large amount of user code). The most obvious alternative
> >>>>>> to an API change would be a larger, somewhat convoluted enum for the
> >>>>>> PeriodicType (DMDA_XPERIODIC_YGHOSTED, DMDA_XYGHOSTED, etc) which could
> >>>>>> at least be made backward compatible.
> >>>>>>
> >>>>>> At least all of the functionality should be there already (since it's
> >>>>>> needed in the periodic case)... it's just higher level code that would
> >>>>>> need to change.
> >>>>>>
> >>>>>> Ethan
> >>>>>>
> >>>>>>>
> >>>>>>> On Dec 7, 2010, at 1:30 PM, Jed Brown wrote:
> >>>>>>>
> >>>>>>>> On Tue, Dec 7, 2010 at 20:21, Ethan Coon <ecoon at lanl.gov> wrote:
> >>>>>>>> 'd like a DA where there are ghost cells on every boundary, and some of
> >>>>>>>> those ghost cells (but not all) are filled in with periodic values.
> >>>>>>>>
> >>>>>>>> It would be useful to people doing explicit stuff if there was a way to get ghost nodes in the local vector without implying periodic communication (and weird coordinate management).
> >>>>>>>>
> >>>>>>>> A related issue for purely explicit is to have a way to VecAXPY without needing to copy to and from a global vector. (TSSSP has low-memory schemes, paying for an extra vector or two is actually significant in that context, and (less significant) I'm certain I can cook up a realistic benchmark where the memcpy costs more than 10%.) I think I know how to implement this sharing transparently (more-or-less using VecNest) so we could make it non-default but be able to activate it at runtime.
> >>>>>>>
> >>>>>>> Why can you not use VecAXPY() on the local Vecs?
> >>>>>>>
> >>>>>>> Barry
> >>>>>>>
> >>>>>>>
> >>>>>>>
> >>>>>>>>
> >>>>>>>> Jed
> >>>>>>>
> >>>>>>
> >>>>>
> >>>>> --
> >>>>> ------------------------------------
> >>>>> Ethan Coon
> >>>>> Post-Doctoral Researcher
> >>>>> Applied Mathematics - T-5
> >>>>> Los Alamos National Laboratory
> >>>>> 505-665-8289
> >>>>>
> >>>>> http://www.ldeo.columbia.edu/~ecoon/
> >>>>> ------------------------------------
> >>>>>
> >>>>
> >>>
> >>> --
> >>> ------------------------------------
> >>> Ethan Coon
> >>> Post-Doctoral Researcher
> >>> Applied Mathematics - T-5
> >>> Los Alamos National Laboratory
> >>> 505-665-8289
> >>>
> >>> http://www.ldeo.columbia.edu/~ecoon/
> >>> ------------------------------------
> >>>
> >>
> >
> > --
> > ------------------------------------
> > Ethan Coon
> > Post-Doctoral Researcher
> > Applied Mathematics - T-5
> > Los Alamos National Laboratory
> > 505-665-8289
> >
> > http://www.ldeo.columbia.edu/~ecoon/
> > ------------------------------------
> >
>
--
------------------------------------
Ethan Coon
Post-Doctoral Researcher
Applied Mathematics - T-5
Los Alamos National Laboratory
505-665-8289
http://www.ldeo.columbia.edu/~ecoon/
------------------------------------
-------------- next part --------------
addition of DMDA_*GHOSTED as an option for DMDAPeriodicType. This enables the inclusion of ghost cells at (nonperiodic) domain boundaries in all dimensions.
diff -r 1bde78661028 include/finclude/petscdm.h
--- a/include/finclude/petscdm.h Thu Mar 10 12:10:19 2011 -0600
+++ b/include/finclude/petscdm.h Thu Mar 10 13:13:17 2011 -0700
@@ -16,6 +16,7 @@
! Types of periodicity
!
PetscEnum DMDA_NONPERIODIC
+ PetscEnum DMDA_NONGHOSTED
PetscEnum DMDA_XPERIODIC
PetscEnum DMDA_YPERIODIC
PetscEnum DMDA_XYPERIODIC
@@ -23,14 +24,24 @@
PetscEnum DMDA_XZPERIODIC
PetscEnum DMDA_YZPERIODIC
PetscEnum DMDA_ZPERIODIC
+ PetscEnum DMDA_XGHOSTED
+ PetscEnum DMDA_YGHOSTED
+ PetscEnum DMDA_ZGHOSTED
PetscEnum DMDA_XYZGHOSTED
- parameter (DMDA_NONPERIODIC = 0,DMDA_XPERIODIC = 1)
- parameter (DMDA_YPERIODIC = 2)
- parameter (DMDA_XYPERIODIC = 3,DMDA_XYZPERIODIC = 4)
- parameter (DMDA_XZPERIODIC = 5,DMDA_YZPERIODIC = 6)
- parameter (DMDA_ZPERIODIC = 7)
- parameter (DMDA_XYZGHOSTED = 8)
+ parameter (DMDA_NONPERIODIC = Z'0')
+ parameter (DMDA_NONGHOSTED = Z'0')
+ parameter (DMDA_XPERIODIC = Z'3')
+ parameter (DMDA_YPERIODIC = Z'C')
+ parameter (DMDA_XYPERIODIC = Z'F')
+ parameter (DMDA_XYZPERIODIC = Z'3F')
+ parameter (DMDA_XZPERIODIC = Z'33')
+ parameter (DMDA_YZPERIODIC = Z'3C')
+ parameter (DMDA_ZPERIODIC = Z'30')
+ parameter (DMDA_XGHOSTED = Z'1')
+ parameter (DMDA_YGHOSTED = Z'4')
+ parameter (DMDA_ZGHOSTED = Z'10')
+ parameter (DMDA_XYZGHOSTED = Z'15')
!
! DMDAInterpolationType
diff -r 1bde78661028 include/petscdm.h
--- a/include/petscdm.h Thu Mar 10 12:10:19 2011 -0600
+++ b/include/petscdm.h Thu Mar 10 13:13:17 2011 -0700
@@ -8,7 +8,6 @@
PETSC_EXTERN_CXX_BEGIN
extern PetscErrorCode DMInitializePackage(const char[]);
-
/*S
DM - Abstract PETSc object that manages an abstract grid object
@@ -55,17 +54,33 @@
M*/
/*E
- DMDAPeriodicType - Is the domain periodic in one or more directions
+ DMDAPeriodicType - Is the domain periodic or ghosted in one or more directions
Level: beginner
- DMDA_XYZGHOSTED means that ghost points are put around all the physical boundaries
- in the local representation of the Vec (i.e. DMDACreate/GetLocalVector().
+ Each dimension may be non-periodic, ghosted (meaning ghost nodes are added outside
+ of the boundary, but not filled by DMDAGlobalToLocal()), or periodic. Dimensions
+ may be composed using the bitwise or operator, i.e.:
+
+ DMDA_XYPERIODIC = DMDA_XPERIODIC | DMDA_YPERIODIC
.seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDACreate()
E*/
-typedef enum { DMDA_NONPERIODIC,DMDA_XPERIODIC,DMDA_YPERIODIC,DMDA_XYPERIODIC,
- DMDA_XYZPERIODIC,DMDA_XZPERIODIC,DMDA_YZPERIODIC,DMDA_ZPERIODIC,DMDA_XYZGHOSTED} DMDAPeriodicType;
+typedef enum {
+ DMDA_NONGHOSTED = 0x0,
+ DMDA_NONPERIODIC = 0x0, /* kept for backwards compatiblity, though it's not precise */
+ DMDA_XGHOSTED = 0x1,
+ DMDA_XPERIODIC = 0x3,
+ DMDA_YGHOSTED = 0x4,
+ DMDA_YPERIODIC = 0xc,
+ DMDA_ZGHOSTED = 0x10,
+ DMDA_ZPERIODIC = 0x30,
+ DMDA_XYPERIODIC = 0xf,
+ DMDA_XZPERIODIC = 0x33,
+ DMDA_YZPERIODIC = 0x3c,
+ DMDA_XYZPERIODIC = 0x3f,
+ DMDA_XYZGHOSTED = 0x15} DMDAPeriodicType;
+
extern const char *DMDAPeriodicTypes[];
/*E
@@ -94,9 +109,12 @@
extern PetscErrorCode DMDASetElementType(DM,DMDAElementType);
extern PetscErrorCode DMDAGetElementType(DM,DMDAElementType*);
-#define DMDAXPeriodic(pt) ((pt)==DMDA_XPERIODIC||(pt)==DMDA_XYPERIODIC||(pt)==DMDA_XZPERIODIC||(pt)==DMDA_XYZPERIODIC)
-#define DMDAYPeriodic(pt) ((pt)==DMDA_YPERIODIC||(pt)==DMDA_XYPERIODIC||(pt)==DMDA_YZPERIODIC||(pt)==DMDA_XYZPERIODIC)
-#define DMDAZPeriodic(pt) ((pt)==DMDA_ZPERIODIC||(pt)==DMDA_XZPERIODIC||(pt)==DMDA_YZPERIODIC||(pt)==DMDA_XYZPERIODIC)
+#define DMDAXPeriodic(pt) ((pt) & 0x2) /* (DMDA_XPERIODIC ^ DMDA_XGHOSTED)) */
+#define DMDAYPeriodic(pt) ((pt) & 0x8) /* (DMDA_YPERIODIC ^ DMDA_YGHOSTED)) */
+#define DMDAZPeriodic(pt) ((pt) & 0x20) /* (DMDA_ZPERIODIC ^ DMDA_ZGHOSTED)) */
+#define DMDAXGhosted(pt) ((pt) & DMDA_XGHOSTED)
+#define DMDAYGhosted(pt) ((pt) & DMDA_YGHOSTED)
+#define DMDAZGhosted(pt) ((pt) & DMDA_ZGHOSTED)
typedef enum { DMDA_X,DMDA_Y,DMDA_Z } DMDADirection;
@@ -137,7 +155,8 @@
extern PetscErrorCode DMDAGetNeighbors(DM,const PetscMPIInt**);
extern PetscErrorCode DMDAGetAO(DM,AO*);
-extern PetscErrorCode DMDASetCoordinates(DM,Vec);
+extern PetscErrorCode DMDASetCoordinates(DM,Vec);
+extern PetscErrorCode DMDASetGhostedCoordinates(DM,Vec);
extern PetscErrorCode DMDAGetCoordinates(DM,Vec *);
extern PetscErrorCode DMDAGetGhostedCoordinates(DM,Vec *);
extern PetscErrorCode DMDAGetCoordinateDA(DM,DM *);
diff -r 1bde78661028 src/characteristic/impls/da/slda.c
--- a/src/characteristic/impls/da/slda.c Thu Mar 10 12:10:19 2011 -0600
+++ b/src/characteristic/impls/da/slda.c Thu Mar 10 13:13:17 2011 -0700
@@ -120,11 +120,11 @@
if ( periodic_type == DMDA_NONPERIODIC ) {
ierr = 0;
} else {
- if (periodic_type==DMDA_XPERIODIC || periodic_type==DMDA_XYPERIODIC) {
+ if (DMDAXPeriodic(periodic_type)) {
while (*x >= ( PetscScalar ) gx ) { *x -= ( PetscScalar ) gx; }
while (*x < 0.0 ) { *x += ( PetscScalar ) gx; }
}
- if (periodic_type==DMDA_YPERIODIC || periodic_type==DMDA_XYPERIODIC) {
+ if (DMDAYPeriodic(periodic_type)) {
while (*y >= ( PetscScalar ) gy ) { *y -= ( PetscScalar ) gy; }
while (*y < 0.0 ) { *y += ( PetscScalar ) gy; }
}
diff -r 1bde78661028 src/characteristic/interface/characteristic.c
--- a/src/characteristic/interface/characteristic.c Thu Mar 10 12:10:19 2011 -0600
+++ b/src/characteristic/interface/characteristic.c Thu Mar 10 13:13:17 2011 -0700
@@ -800,10 +800,10 @@
ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
ierr = DMDAGetInfo(da, 0, 0, 0, 0, &PI,&PJ, 0, 0, 0, &periodic_type, 0);
- if (periodic_type==DMDA_XPERIODIC || periodic_type==DMDA_XYPERIODIC) {
+ if (DMDAXPeriodic(periodic_type)) {
IPeriodic = PETSC_TRUE;
}
- if (periodic_type==DMDA_YPERIODIC || periodic_type==DMDA_XYPERIODIC) {
+ if (DMDAYPeriodic(periodic_type)) {
JPeriodic = PETSC_TRUE;
}
diff -r 1bde78661028 src/dm/examples/tests/ex4.c
--- a/src/dm/examples/tests/ex4.c Thu Mar 10 12:10:19 2011 -0600
+++ b/src/dm/examples/tests/ex4.c Thu Mar 10 13:13:17 2011 -0700
@@ -26,20 +26,25 @@
ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",300,0,400,400,&viewer);CHKERRQ(ierr);
/* Readoptions */
- ierr = PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);CHKERRQ(ierr);
- ierr = PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);CHKERRQ(ierr);
+ ierr = PetscOptionsGetInt(PETSC_NULL,"-NX",&M,PETSC_NULL);CHKERRQ(ierr);
+ ierr = PetscOptionsGetInt(PETSC_NULL,"-NY",&N,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-s",&s,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-w",&w,PETSC_NULL);CHKERRQ(ierr);
+ wrap = 0x0;
flg = PETSC_FALSE;
- ierr = PetscOptionsGetBool(PETSC_NULL,"-xwrap",&flg,PETSC_NULL);CHKERRQ(ierr); if (flg) wrap = DMDA_XPERIODIC;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-xperiodic",&flg,PETSC_NULL);CHKERRQ(ierr); if (flg) wrap = wrap | DMDA_XPERIODIC;
flg = PETSC_FALSE;
- ierr = PetscOptionsGetBool(PETSC_NULL,"-ywrap",&flg,PETSC_NULL);CHKERRQ(ierr); if (flg) wrap = DMDA_YPERIODIC;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-yperiodic",&flg,PETSC_NULL);CHKERRQ(ierr); if (flg) wrap = wrap | DMDA_YPERIODIC;
flg = PETSC_FALSE;
- ierr = PetscOptionsGetBool(PETSC_NULL,"-xywrap",&flg,PETSC_NULL);CHKERRQ(ierr); if (flg) wrap = DMDA_XYPERIODIC;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-xghosted",&flg,PETSC_NULL);CHKERRQ(ierr); if (flg) wrap = wrap | DMDA_XGHOSTED;
flg = PETSC_FALSE;
- ierr = PetscOptionsGetBool(PETSC_NULL,"-star",&flg,PETSC_NULL);CHKERRQ(ierr); if (flg) st = DMDA_STENCIL_STAR;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-yghosted",&flg,PETSC_NULL);CHKERRQ(ierr); if (flg) wrap = wrap | DMDA_YGHOSTED;
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-star",&flg,PETSC_NULL);CHKERRQ(ierr); if (flg) st = DMDA_STENCIL_STAR;
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-box",&flg,PETSC_NULL);CHKERRQ(ierr); if (flg) st = DMDA_STENCIL_BOX;
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(PETSC_NULL,"-testorder",&testorder,PETSC_NULL);CHKERRQ(ierr);
/*
@@ -142,7 +147,7 @@
}
}
ierr = PetscFree(iglobal);CHKERRQ(ierr);
- }
+ }
/* Free memory */
ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
diff -r 1bde78661028 src/dm/examples/tests/ex6.c
--- a/src/dm/examples/tests/ex6.c Thu Mar 10 12:10:19 2011 -0600
+++ b/src/dm/examples/tests/ex6.c Thu Mar 10 13:13:17 2011 -0700
@@ -1,4 +1,3 @@
-
static char help[] = "Tests various 3-dimensional DMDA routines.\n\n";
#include "petscdm.h"
@@ -18,7 +17,7 @@
PetscViewer viewer;
Vec local,global;
PetscScalar value;
- DMDAPeriodicType wrap = DMDA_XYPERIODIC;
+ DMDAPeriodicType wrap = DMDA_NONPERIODIC;
DMDAStencilType stencil_type = DMDA_STENCIL_BOX;
AO ao;
PetscBool flg = PETSC_FALSE;
@@ -26,18 +25,50 @@
ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);
ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",300,0,400,300,&viewer);CHKERRQ(ierr);
- /* Read options */
- ierr = PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);CHKERRQ(ierr);
- ierr = PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);CHKERRQ(ierr);
- ierr = PetscOptionsGetInt(PETSC_NULL,"-P",&P,PETSC_NULL);CHKERRQ(ierr);
+ /* Read options */
+ ierr = PetscOptionsGetInt(PETSC_NULL,"-NX",&M,PETSC_NULL);CHKERRQ(ierr);
+ ierr = PetscOptionsGetInt(PETSC_NULL,"-NY",&N,PETSC_NULL);CHKERRQ(ierr);
+ ierr = PetscOptionsGetInt(PETSC_NULL,"-NZ",&P,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-p",&p,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-s",&s,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-w",&w,PETSC_NULL);CHKERRQ(ierr);
+ flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(PETSC_NULL,"-star",&flg,PETSC_NULL);CHKERRQ(ierr);
if (flg) stencil_type = DMDA_STENCIL_STAR;
- ierr = PetscOptionsGetBool(PETSC_NULL,"-test_order",&test_order,PETSC_NULL);CHKERRQ(ierr);
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-box",&flg,PETSC_NULL);CHKERRQ(ierr);
+ if (flg) stencil_type = DMDA_STENCIL_BOX;
+
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-xperiodic",&flg,PETSC_NULL);CHKERRQ(ierr);
+ if (flg) wrap = wrap | DMDA_XPERIODIC;
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-xghosted",&flg,PETSC_NULL);CHKERRQ(ierr);
+ if (flg) wrap = wrap | DMDA_XPERIODIC;
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-xnonghosted",&flg,PETSC_NULL);CHKERRQ(ierr);
+
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-yperiodic",&flg,PETSC_NULL);CHKERRQ(ierr);
+ if (flg) wrap = wrap | DMDA_YPERIODIC;
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-yghosted",&flg,PETSC_NULL);CHKERRQ(ierr);
+ if (flg) wrap = wrap | DMDA_YPERIODIC;
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-ynonghosted",&flg,PETSC_NULL);CHKERRQ(ierr);
+
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-zperiodic",&flg,PETSC_NULL);CHKERRQ(ierr);
+ if (flg) wrap = wrap | DMDA_ZPERIODIC;
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-zghosted",&flg,PETSC_NULL);CHKERRQ(ierr);
+ if (flg) wrap = wrap | DMDA_ZPERIODIC;
+ flg = PETSC_FALSE;
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-znonghosted",&flg,PETSC_NULL);CHKERRQ(ierr);
+
+ ierr = PetscOptionsGetBool(PETSC_NULL,"-testorder",&test_order,PETSC_NULL);CHKERRQ(ierr);
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(PETSC_NULL,"-distribute",&flg,PETSC_NULL);CHKERRQ(ierr);
diff -r 1bde78661028 src/dm/examples/tests/scripts/ex4.script
--- a/src/dm/examples/tests/scripts/ex4.script Thu Mar 10 12:10:19 2011 -0600
+++ b/src/dm/examples/tests/scripts/ex4.script Thu Mar 10 13:13:17 2011 -0700
@@ -4,25 +4,29 @@
#
# w = degress of freedom per node
# s = stencil width
-# M,N = global dimension in each direction of array
+# NX,NY = global dimension in each direction of array
# m,n = number of processors in each dimension
# np = number of processors
#
foreach m (1 2 3)
foreach n (1 2 3)
- foreach M (16 17)
- foreach N (14 15)
+ foreach NX (16 17)
+ foreach NY (14 15)
foreach w (1 2 3)
foreach s (1 2)
- set np = `expr $m \* $n`
- echo 'mpirun -np' $np' ex4 -testorder -nox -w' $w '-s' $s '-M' $M '-N' $N '-m' $m '-n' $n
- mpirun -np $np ex4 -testorder -nox -w $w -s $s -M $M -N $N -m $m -n $n
+ foreach stencil ('-star' '-box')
+ foreach xwrap ('-xperiodic' '-xghosted' '-xnonghosted')
+ foreach ywrap ('-yperiodic' '-yghosted' '-ynonghosted')
+ set np = `expr $m \* $n`
+ echo 'mpiexec -n' $np' ./ex4 -testorder -nox -w' $w '-s' $s '-NX' $NX '-NY' $NY '-m' $m '-n' $n $xwrap $ywrap $stencil
+ mpiexec -n $np ./ex4 -testorder -nox -w $w -s $s -NX $NX -NY $NY -m $m -n $n $xwrap $ywrap $stencil
+ end
+ end
+ end
end
end
end
end
end
end
-
-
diff -r 1bde78661028 src/dm/examples/tests/scripts/ex6.script
--- a/src/dm/examples/tests/scripts/ex6.script Thu Mar 10 12:10:19 2011 -0600
+++ b/src/dm/examples/tests/scripts/ex6.script Thu Mar 10 13:13:17 2011 -0700
@@ -12,14 +12,22 @@
foreach m (1 2 3)
foreach n (1 2 3)
foreach p (1 2 3)
- foreach M (16 17)
- foreach N (14 15)
- foreach P (13 12)
+ foreach NX (16 17)
+ foreach NY (14 15)
+ foreach NZ (13 12)
foreach w (1 2 3)
foreach s (1 2)
- set np = `expr $m \* $n \* $p`
- echo 'mpirun -np' $np' ex6 -testorder -nox -w' $w '-s' $s '-M' $M '-N' $N '-P' $P '-m' $m '-n' $n '-p' $p
- mpirun -np $np ex6 -testorder -nox -w $w -s $s -M $M -N $N -P $P -m $m -n $n -p $p
+ foreach stencil ('-box' '-star')
+ foreach xwrap ('-xperiodic' '-xghosted' '-xnonghosted')
+ foreach ywrap ('-yperiodic' '-yghosted' '-ynonghosted')
+ foreach zwrap ('-zperiodic' '-zghosted' '-znonghosted')
+ set np = `expr $m \* $n \* $p`
+ echo 'mpiexec -n' $np' ./ex6 -testorder -nox -w' $w '-s' $s '-NX' $NX '-NY' $NY '-NZ' $NZ '-m' $m '-n' $n '-p' $p $stencil $xwrap $ywrap $zwrap
+ mpiexec -n $np ./ex6 -testorder -nox -w $w -s $s -NX $NX -NY $NY -NZ $NZ -m $m -n $n -p $p $stencil $xwrap $ywrap $zwrap
+ end
+ end
+ end
+ end
end
end
end
diff -r 1bde78661028 src/dm/impls/da/da1.c
--- a/src/dm/impls/da/da1.c Thu Mar 10 12:10:19 2011 -0600
+++ b/src/dm/impls/da/da1.c Thu Mar 10 13:13:17 2011 -0700
@@ -7,7 +7,8 @@
#include "private/daimpl.h" /*I "petscdm.h" I*/
const char *DMDAPeriodicTypes[] = {"NONPERIODIC","XPERIODIC","YPERIODIC","XYPERIODIC",
- "XYZPERIODIC","XZPERIODIC","YZPERIODIC","ZPERIODIC","XYZGHOSTED","DMDAPeriodicType","DMDA_",0};
+ "XYZPERIODIC","XZPERIODIC","YZPERIODIC","ZPERIODIC",
+ "XYZGHOSTED","DMDAPeriodicType","DMDA_",0};
#undef __FUNCT__
#define __FUNCT__ "DMView_DA_1d"
@@ -143,7 +144,7 @@
IS to, from;
PetscBool flg1 = PETSC_FALSE, flg2 = PETSC_FALSE;
PetscMPIInt rank, size;
- PetscInt i,*idx,nn,left,xs,xe,x,Xs,Xe,start,end,m;
+ PetscInt i,*idx,nn,left,xs,xe,x,Xs,Xe,start,end,m,IXs,IXe;
PetscErrorCode ierr;
PetscFunctionBegin;
@@ -206,12 +207,20 @@
xe = xs + x;
/* determine ghost region */
- if (wrap == DMDA_XPERIODIC || wrap == DMDA_XYZGHOSTED) {
- Xs = xs - sDist;
+ if (DMDAXGhosted(wrap)) {
+ Xs = xs - sDist;
Xe = xe + sDist;
} else {
- if ((xs-sDist) >= 0) Xs = xs-sDist; else Xs = 0;
- if ((xe+sDist) <= M*dof) Xe = xe+sDist; else Xe = M*dof;
+ if ((xs-sDist) >= 0) Xs = xs-sDist; else Xs = 0;
+ if ((xe+sDist) <= M*dof) Xe = xe+sDist; else Xe = M*dof;
+ }
+
+ if (DMDAXPeriodic(wrap)) {
+ IXs = xs - sDist;
+ IXe = xe + sDist;
+ } else {
+ if ((xs-sDist) >= 0) IXs = xs-sDist; else IXs = 0;
+ if ((xe+sDist) <= M*dof) IXe = xe+sDist; else IXe = M*dof;
}
/* allocate the base parallel and sequential vectors */
@@ -221,7 +230,7 @@
dd->nlocal = (Xe-Xs);
ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dd->nlocal,0,&local);CHKERRQ(ierr);
ierr = VecSetBlockSize(local,dof);CHKERRQ(ierr);
-
+
/* Create Local to Global Vector Scatter Context */
/* local to global inserts non-ghost point region into global */
ierr = VecGetOwnershipRange(global,&start,&end);CHKERRQ(ierr);
@@ -234,57 +243,37 @@
/* Create Global to Local Vector Scatter Context */
/* global to local must retrieve ghost points */
- if (wrap == DMDA_XYZGHOSTED) {
- if (size == 1) {
- ierr = ISCreateStride(comm,(xe-xs),sDist,1,&to);CHKERRQ(ierr);
- } else if (!rank) {
- ierr = ISCreateStride(comm,(Xe-xs),sDist,1,&to);CHKERRQ(ierr);
- } else if (rank == size-1) {
- ierr = ISCreateStride(comm,(xe-Xs),0,1,&to);CHKERRQ(ierr);
- } else {
- ierr = ISCreateStride(comm,(Xe-Xs),0,1,&to);CHKERRQ(ierr);
- }
- } else {
- ierr = ISCreateStride(comm,(Xe-Xs),0,1,&to);CHKERRQ(ierr);
- }
-
+ ierr = ISCreateStride(comm,(IXe-IXs),IXs-Xs,1,&to);CHKERRQ(ierr);
+
ierr = PetscMalloc((x+2*sDist)*sizeof(PetscInt),&idx);CHKERRQ(ierr);
ierr = PetscLogObjectMemory(da,(x+2*sDist)*sizeof(PetscInt));CHKERRQ(ierr);
- nn = 0;
- if (wrap == DMDA_XPERIODIC) { /* Handle all cases with wrap first */
+ for (i=0; i<IXs-Xs; i++) {idx[i] = -1; } /* prepend with -1s if needed for ghosted case*/
+ nn = IXs-Xs;
+ if DMDAXPeriodic(wrap) { /* Handle all cases with wrap first */
for (i=0; i<sDist; i++) { /* Left ghost points */
if ((xs-sDist+i)>=0) { idx[nn++] = xs-sDist+i;}
else { idx[nn++] = M*dof+(xs-sDist+i);}
}
for (i=0; i<x; i++) { idx [nn++] = xs + i;} /* Non-ghost points */
-
+
for (i=0; i<sDist; i++) { /* Right ghost points */
if ((xe+i)<M*dof) { idx [nn++] = xe+i; }
else { idx [nn++] = (xe+i) - M*dof;}
}
- } else if (wrap == DMDA_XYZGHOSTED) {
-
- if (sDist <= xs) {for (i=0; i<sDist; i++) {idx[nn++] = xs - sDist + i;}}
-
- for (i=0; i<x; i++) { idx [nn++] = xs + i;}
-
- if ((xe+sDist)<=M*dof) {for (i=0; i<sDist; i++) {idx[nn++]=xe+i;}}
-
} else { /* Now do all cases with no wrapping */
-
if (sDist <= xs) {for (i=0; i<sDist; i++) {idx[nn++] = xs - sDist + i;}}
else {for (i=0; i<xs; i++) {idx[nn++] = i;}}
for (i=0; i<x; i++) { idx [nn++] = xs + i;}
-
+
if ((xe+sDist)<=M*dof) {for (i=0; i<sDist; i++) {idx[nn++]=xe+i;}}
else {for (i=xe; i<(M*dof); i++) {idx[nn++]=i;}}
}
- ierr = ISCreateGeneral(comm,nn,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
+ ierr = ISCreateGeneral(comm,nn-IXs+Xs,&idx[IXs-Xs],PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr);
ierr = PetscLogObjectParent(da,to);CHKERRQ(ierr);
ierr = PetscLogObjectParent(da,from);CHKERRQ(ierr);
@@ -306,32 +295,8 @@
Set the local to global ordering in the global vector, this allows use
of VecSetValuesLocal().
*/
- if (wrap == DMDA_XYZGHOSTED) {
- PetscInt *tmpidx;
- if (size == 1) {
- ierr = PetscMalloc((nn+2*sDist)*sizeof(PetscInt),&tmpidx);CHKERRQ(ierr);
- for (i=0; i<sDist; i++) tmpidx[i] = -1;
- ierr = PetscMemcpy(tmpidx+sDist,idx,nn*sizeof(PetscInt));CHKERRQ(ierr);
- for (i=nn+sDist; i<nn+2*sDist; i++) tmpidx[i] = -1;
- ierr = PetscFree(idx);CHKERRQ(ierr);
- idx = tmpidx;
- nn += 2*sDist;
- } else if (!rank) { /* must preprend -1 marker for ghost location that have no global value */
- ierr = PetscMalloc((nn+sDist)*sizeof(PetscInt),&tmpidx);CHKERRQ(ierr);
- for (i=0; i<sDist; i++) tmpidx[i] = -1;
- ierr = PetscMemcpy(tmpidx+sDist,idx,nn*sizeof(PetscInt));CHKERRQ(ierr);
- ierr = PetscFree(idx);CHKERRQ(ierr);
- idx = tmpidx;
- nn += sDist;
- } else if (rank == size-1) { /* must postpend -1 marker for ghost location that have no global value */
- ierr = PetscMalloc((nn+sDist)*sizeof(PetscInt),&tmpidx);CHKERRQ(ierr);
- ierr = PetscMemcpy(tmpidx,idx,nn*sizeof(PetscInt));CHKERRQ(ierr);
- for (i=nn; i<nn+sDist; i++) tmpidx[i] = -1;
- ierr = PetscFree(idx);CHKERRQ(ierr);
- idx = tmpidx;
- nn += sDist;
- }
- }
+ for (i=0; i<Xe-IXe; i++) {idx[nn++] = -1; } /* pad with -1s if needed for ghosted case*/
+
ierr = ISLocalToGlobalMappingCreate(comm,nn,idx,PETSC_OWN_POINTER,&da->ltogmap);CHKERRQ(ierr);
ierr = ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);CHKERRQ(ierr);
ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr);
@@ -353,8 +318,8 @@
Input Parameters:
+ comm - MPI communicator
-. wrap - type of periodicity should the array have, if any. Use
- either DMDA_NONPERIODIC or DMDA_XPERIODIC
+. wrap - type of ghost cells at the boundary the array should have, if any. Use
+ DMDA_NONGHOSTED, DMDA_XGHOSTED, or DMDA_XPERIODIC.
. M - global dimension of the array (use -M to indicate that it may be set to a different value
from the command line with -da_grid_x <M>)
. dof - number of degrees of freedom per node
diff -r 1bde78661028 src/dm/impls/da/da2.c
--- a/src/dm/impls/da/da2.c Thu Mar 10 12:10:19 2011 -0600
+++ b/src/dm/impls/da/da2.c Thu Mar 10 13:13:17 2011 -0700
@@ -1242,14 +1242,15 @@
PetscInt *ly = dd->ly;
MPI_Comm comm;
PetscMPIInt rank,size;
- PetscInt xs,xe,ys,ye,x,y,Xs,Xe,Ys,Ye,start,end;
- PetscInt up,down,left,i,n0,n1,n2,n3,n5,n6,n7,n8,*idx,nn;
- PetscInt xbase,*bases,*ldims,j,x_t,y_t,s_t,base,count;
+ PetscInt xs,xe,ys,ye,x,y,Xs,Xe,Ys,Ye,start,end,IXs,IXe,IYs,IYe;
+ PetscInt up,down,left,right,i,n0,n1,n2,n3,n5,n6,n7,n8,*idx,nn,*idx_cpy;
+ const PetscInt *idx_full;
+ PetscInt xbase,*bases,*ldims,j,x_t,y_t,s_t,base,count,count_dbg;
PetscInt s_x,s_y; /* s proportionalized to w */
PetscInt sn0 = 0,sn2 = 0,sn6 = 0,sn8 = 0;
Vec local,global;
VecScatter ltog,gtol;
- IS to,from;
+ IS to,from,ltogis;
PetscErrorCode ierr;
PetscFunctionBegin;
@@ -1349,32 +1350,21 @@
xe = xs + x;
ye = ys + y;
- /* determine ghost region */
+ /* determine ghost region (Xs) and region scattered into (IXs) */
/* Assume No Periodicity */
- if (xs-s > 0) Xs = xs - s; else Xs = 0;
- if (ys-s > 0) Ys = ys - s; else Ys = 0;
- if (xe+s <= M) Xe = xe + s; else Xe = M;
- if (ye+s <= N) Ye = ye + s; else Ye = N;
+ if (xs-s > 0) { Xs = xs - s; IXs = xs - s; } else { Xs = 0; IXs = 0; }
+ if (xe+s <= M) { Xe = xe + s; IXe = xe + s; } else { Xe = M; IXe = M; }
+ if (ys-s > 0) { Ys = ys - s; IYs = ys - s; } else { Ys = 0; IYs = 0; }
+ if (ye+s <= N) { Ye = ye + s; IYe = ye + s; } else { Ye = N; IYe = N; }
- /* X Periodic */
- if (DMDAXPeriodic(wrap)){
- Xs = xs - s;
- Xe = xe + s;
- }
-
- /* Y Periodic */
- if (DMDAYPeriodic(wrap)){
- Ys = ys - s;
- Ye = ye + s;
- }
+ /* fix for periodicity/ghosted */
+ if (DMDAXGhosted(wrap)) { Xs = xs - s; Xe = xe + s; }
+ if (DMDAXPeriodic(wrap)) { IXs = xs - s; IXe = xe + s; }
+ if (DMDAYGhosted(wrap)) { Ys = ys - s; Ye = ye + s; }
+ if (DMDAYPeriodic(wrap)) { IYs = ys - s; IYe = ye + s; }
/* Resize all X parameters to reflect w */
- x *= dof;
- xs *= dof;
- xe *= dof;
- Xs *= dof;
- Xe *= dof;
- s_x = s*dof;
+ s_x = s;
s_y = s;
/* determine starting point of each processor */
@@ -1388,73 +1378,101 @@
for (i=1; i<=size; i++) {
bases[i] += bases[i-1];
}
+ base = bases[rank]*dof;
/* allocate the base parallel and sequential vectors */
- dd->Nlocal = x*y;
+ dd->Nlocal = x*y*dof;
ierr = VecCreateMPIWithArray(comm,dd->Nlocal,PETSC_DECIDE,0,&global);CHKERRQ(ierr);
ierr = VecSetBlockSize(global,dof);CHKERRQ(ierr);
- dd->nlocal = (Xe-Xs)*(Ye-Ys);
+ dd->nlocal = (Xe-Xs)*(Ye-Ys)*dof;
ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dd->nlocal,0,&local);CHKERRQ(ierr);
ierr = VecSetBlockSize(local,dof);CHKERRQ(ierr);
-
/* generate appropriate vector scatters */
/* local to global inserts non-ghost point region into global */
ierr = VecGetOwnershipRange(global,&start,&end);CHKERRQ(ierr);
- ierr = ISCreateStride(comm,x*y,start,1,&to);CHKERRQ(ierr);
+ ierr = ISCreateStride(comm,x*y*dof,start,1,&to);CHKERRQ(ierr);
- left = xs - Xs; down = ys - Ys; up = down + y;
- ierr = PetscMalloc(x*(up - down)*sizeof(PetscInt),&idx);CHKERRQ(ierr);
+ count_dbg = x*y;
+ ierr = PetscMalloc(x*y*sizeof(PetscInt),&idx);CHKERRQ(ierr);
+ left = xs - Xs; right = left + x;
+ down = ys - Ys; up = down + y;
count = 0;
for (i=down; i<up; i++) {
- for (j=0; j<x/dof; j++) {
- idx[count++] = (left + i*(Xe-Xs) + j*dof)/dof;
+ for (j=left; j<right; j++) {
+ idx[count++] = i*(Xe-Xs) + j;
}
}
+ if (count != count_dbg) {
+ SETERRQ(((PetscObject)da)->comm, PETSC_ERR_SUP,"count != count_dbg");
+ PetscFunctionReturn(1);
+ }
+
ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&from);CHKERRQ(ierr);
-
ierr = VecScatterCreate(local,from,global,to,<og);CHKERRQ(ierr);
ierr = PetscLogObjectParent(dd,ltog);CHKERRQ(ierr);
ierr = ISDestroy(from);CHKERRQ(ierr);
ierr = ISDestroy(to);CHKERRQ(ierr);
- /* global to local must include ghost points */
+ /* global to local must include ghost points within the domain,
+ but not ghost points outside the domain that aren't periodic */
if (stencil_type == DMDA_STENCIL_BOX) {
- ierr = ISCreateStride(comm,(Xe-Xs)*(Ye-Ys),0,1,&to);CHKERRQ(ierr);
+ count_dbg = (IXe-IXs)*(IYe-IYs);
+ ierr = PetscMalloc(count_dbg*sizeof(PetscInt),&idx);CHKERRQ(ierr);
+
+ left = IXs - Xs; right = left + (IXe-IXs);
+ down = IYs - Ys; up = down + (IYe-IYs);
+ count = 0;
+ for (i=down; i<up; i++) {
+ for (j=left; j<right; j++) {
+ idx[count++] = j + i*(Xe-Xs);
+ }
+ }
+ if (count != count_dbg) {
+ SETERRQ(((PetscObject)da)->comm, PETSC_ERR_SUP,"count != count_dbg");
+ PetscFunctionReturn(1);
+ }
+ ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr);
+
} else {
/* must drop into cross shape region */
/* ---------|
| top |
- |--- ---|
+ |--- ---| up
| middle |
| |
- ---- ----
+ ---- ---- down
| bottom |
-----------
- Xs xs xe Xe */
+ Xs xs xe Xe */
+ count_dbg = (ys-IYs)*x + y*(IXe-IXs) + (IYe-ye)*x;
+ ierr = PetscMalloc(count_dbg*sizeof(PetscInt),&idx);CHKERRQ(ierr);
+
+ left = xs - Xs; right = left + x;
+ down = ys - Ys; up = down + y;
+ count = 0;
/* bottom */
- left = xs - Xs; down = ys - Ys; up = down + y;
- count = down*(xe-xs) + (up-down)*(Xe-Xs) + (Ye-Ys-up)*(xe-xs);
- ierr = PetscMalloc(count*sizeof(PetscInt)/dof,&idx);CHKERRQ(ierr);
- count = 0;
- for (i=0; i<down; i++) {
- for (j=0; j<xe-xs; j += dof) {
- idx[count++] = left + i*(Xe-Xs) + j;
+ for (i=(IYs-Ys); i<down; i++) {
+ for (j=left; j<right; j++) {
+ idx[count++] = j + i*(Xe-Xs);
}
}
/* middle */
for (i=down; i<up; i++) {
- for (j=0; j<Xe-Xs; j += dof) {
- idx[count++] = i*(Xe-Xs) + j;
+ for (j=(IXs-Xs); j<(IXe-Xs); j++) {
+ idx[count++] = j + i*(Xe-Xs);
}
}
/* top */
- for (i=up; i<Ye-Ys; i++) {
- for (j=0; j<xe-xs; j += dof) {
- idx[count++] = (left + i*(Xe-Xs) + j);
+ for (i=up; i<up+IYe-ye; i++) {
+ for (j=left; j<right; j++) {
+ idx[count++] = j + i*(Xe-Xs);
}
}
- for (i=0; i<count; i++) idx[i] = idx[i]/dof;
+ if (count != count_dbg) {
+ SETERRQ(((PetscObject)da)->comm, PETSC_ERR_SUP,"count != count_dbg");
+ PetscFunctionReturn(1);
+ }
ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr);
}
@@ -1465,9 +1483,9 @@
*/
/* Assume the Non-Periodic Case */
- n1 = rank - m;
+ n1 = rank - m;
if (rank % m) {
- n0 = n1 - 1;
+ n0 = n1 - 1;
} else {
n0 = -1;
}
@@ -1486,29 +1504,13 @@
}
n7 = rank + m; if (n7 >= m*n) n7 = -1;
-
+ if (DMDAXPeriodic(wrap) && DMDAYPeriodic(wrap)) {
/* Modify for Periodic Cases */
- if (wrap == DMDA_YPERIODIC) { /* Handle Top and Bottom Sides */
- if (n1 < 0) n1 = rank + m * (n-1);
- if (n7 < 0) n7 = rank - m * (n-1);
- if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
- if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
- if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
- if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
- } else if (wrap == DMDA_XPERIODIC) { /* Handle Left and Right Sides */
- if (n3 < 0) n3 = rank + (m-1);
- if (n5 < 0) n5 = rank - (m-1);
- if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
- if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
- if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
- if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
- } else if (wrap == DMDA_XYPERIODIC) {
-
/* Handle all four corners */
if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1;
if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m;
- if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1;
+ if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1;
/* Handle Top and Bottom Sides */
if (n1 < 0) n1 = rank + m * (n-1);
@@ -1525,7 +1527,22 @@
if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
+ } else if (DMDAYPeriodic(wrap)) { /* Handle Top and Bottom Sides */
+ if (n1 < 0) n1 = rank + m * (n-1);
+ if (n7 < 0) n7 = rank - m * (n-1);
+ if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
+ if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
+ if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
+ if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
+ } else if (DMDAXPeriodic(wrap)) { /* Handle Left and Right Sides */
+ if (n3 < 0) n3 = rank + (m-1);
+ if (n5 < 0) n5 = rank - (m-1);
+ if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
+ if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
+ if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
+ if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
}
+
ierr = PetscMalloc(9*sizeof(PetscInt),&dd->neighbors);CHKERRQ(ierr);
dd->neighbors[0] = n0;
dd->neighbors[1] = n1;
@@ -1543,14 +1560,14 @@
n0 = n2 = n6 = n8 = -1;
}
- ierr = PetscMalloc((x+2*s_x)*(y+2*s_y)*sizeof(PetscInt),&idx);CHKERRQ(ierr);
- ierr = PetscLogObjectMemory(da,(x+2*s_x)*(y+2*s_y)*sizeof(PetscInt));CHKERRQ(ierr);
+ ierr = PetscMalloc((Xe-Xs)*(Ye-Ys)*sizeof(PetscInt),&idx);CHKERRQ(ierr);
+ ierr = PetscLogObjectMemory(da,(Xe-Xs)*(Ye-Ys)*sizeof(PetscInt));CHKERRQ(ierr);
+
nn = 0;
-
xbase = bases[rank];
for (i=1; i<=s_y; i++) {
if (n0 >= 0) { /* left below */
- x_t = lx[n0 % m]*dof;
+ x_t = lx[n0 % m];
y_t = ly[(n0/m)];
s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
@@ -1562,7 +1579,7 @@
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
}
if (n2 >= 0) { /* right below */
- x_t = lx[n2 % m]*dof;
+ x_t = lx[n2 % m];
y_t = ly[(n2/m)];
s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
@@ -1571,7 +1588,7 @@
for (i=0; i<y; i++) {
if (n3 >= 0) { /* directly left */
- x_t = lx[n3 % m]*dof;
+ x_t = lx[n3 % m];
/* y_t = y; */
s_t = bases[n3] + (i+1)*x_t - s_x;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
@@ -1580,7 +1597,7 @@
for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */
if (n5 >= 0) { /* directly right */
- x_t = lx[n5 % m]*dof;
+ x_t = lx[n5 % m];
/* y_t = y; */
s_t = bases[n5] + (i)*x_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
@@ -1589,7 +1606,7 @@
for (i=1; i<=s_y; i++) {
if (n6 >= 0) { /* left above */
- x_t = lx[n6 % m]*dof;
+ x_t = lx[n6 % m];
/* y_t = ly[(n6/m)]; */
s_t = bases[n6] + (i)*x_t - s_x;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
@@ -1601,119 +1618,145 @@
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
}
if (n8 >= 0) { /* right above */
- x_t = lx[n8 % m]*dof;
+ x_t = lx[n8 % m];
/* y_t = ly[(n8/m)]; */
s_t = bases[n8] + (i-1)*x_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
}
}
- base = bases[rank];
- {
- PetscInt nnn = nn/dof,*iidx;
- ierr = PetscMalloc(nnn*sizeof(PetscInt),&iidx);CHKERRQ(ierr);
- for (i=0; i<nnn; i++) {
- iidx[i] = idx[dof*i]/dof;
- }
- ierr = ISCreateBlock(comm,dof,nnn,iidx,PETSC_OWN_POINTER,&from);CHKERRQ(ierr);
- }
+ ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr);
ierr = PetscLogObjectParent(da,gtol);CHKERRQ(ierr);
ierr = ISDestroy(to);CHKERRQ(ierr);
ierr = ISDestroy(from);CHKERRQ(ierr);
if (stencil_type == DMDA_STENCIL_STAR) {
+ n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8;
+ }
+
+ if ((stencil_type == DMDA_STENCIL_STAR) ||
+ (!DMDAXPeriodic(wrap) && DMDAXGhosted(wrap)) ||
+ (!DMDAYPeriodic(wrap) && DMDAYGhosted(wrap))) {
/*
Recompute the local to global mappings, this time keeping the
- information about the cross corner processor numbers.
+ information about the cross corner processor numbers and any ghosted
+ but not periodic indices.
*/
- n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8;
nn = 0;
xbase = bases[rank];
for (i=1; i<=s_y; i++) {
if (n0 >= 0) { /* left below */
- x_t = lx[n0 % m]*dof;
+ x_t = lx[n0 % m];
y_t = ly[(n0/m)];
s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xs-Xs > 0 && ys-Ys > 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
if (n1 >= 0) { /* directly below */
x_t = x;
y_t = ly[(n1/m)];
s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
+ } else if (ys-Ys > 0) {
+ for (j=0; j<x; j++) { idx[nn++] = -1;}
}
if (n2 >= 0) { /* right below */
- x_t = lx[n2 % m]*dof;
+ x_t = lx[n2 % m];
y_t = ly[(n2/m)];
s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xe-xe> 0 && ys-Ys > 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
for (i=0; i<y; i++) {
if (n3 >= 0) { /* directly left */
- x_t = lx[n3 % m]*dof;
+ x_t = lx[n3 % m];
/* y_t = y; */
s_t = bases[n3] + (i+1)*x_t - s_x;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xs-Xs > 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */
if (n5 >= 0) { /* directly right */
- x_t = lx[n5 % m]*dof;
+ x_t = lx[n5 % m];
/* y_t = y; */
s_t = bases[n5] + (i)*x_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xe-xe > 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
for (i=1; i<=s_y; i++) {
if (n6 >= 0) { /* left above */
- x_t = lx[n6 % m]*dof;
+ x_t = lx[n6 % m];
/* y_t = ly[(n6/m)]; */
s_t = bases[n6] + (i)*x_t - s_x;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xs-Xs > 0 && Ye-ye > 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
if (n7 >= 0) { /* directly above */
x_t = x;
/* y_t = ly[(n7/m)]; */
s_t = bases[n7] + (i-1)*x_t;
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
+ } else if (Ye-ye > 0) {
+ for (j=0; j<x; j++) { idx[nn++] = -1;}
}
if (n8 >= 0) { /* right above */
- x_t = lx[n8 % m]*dof;
+ x_t = lx[n8 % m];
/* y_t = ly[(n8/m)]; */
s_t = bases[n8] + (i-1)*x_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xe-xe > 0 && Ye-ye > 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
}
- ierr = PetscFree2(bases,ldims);CHKERRQ(ierr);
+ /*
+ Set the local to global ordering in the global vector, this allows use
+ of VecSetValuesLocal().
+ */
+ if (nn != (Xe-Xs)*(Ye-Ys)) {
+ SETERRQ(((PetscObject)da)->comm, PETSC_ERR_SUP,"nn != count_dbg");
+ PetscFunctionReturn(1);
+ }
+ ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_OWN_POINTER,<ogis);CHKERRQ(ierr);
+ ierr = PetscMalloc(nn*dof*sizeof(PetscInt),&idx_cpy);CHKERRQ(ierr);
+ /* ierr = PetscLogObjectMemory(da,nn*dof*sizeof(PetscInt));CHKERRQ(ierr);*/
+ ierr = ISGetIndices(ltogis, &idx_full);
+ ierr = PetscMemcpy(idx_cpy,idx_full,nn*dof*sizeof(PetscInt));CHKERRQ(ierr);
+ CHKMEMQ;
+ ierr = ISRestoreIndices(ltogis, &idx_full);
+ ierr = ISLocalToGlobalMappingCreateIS(ltogis,&da->ltogmap);CHKERRQ(ierr);
+ ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr);
+ ierr = ISDestroy(ltogis);CHKERRQ(ierr);
+ ierr = ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);CHKERRQ(ierr);
+ ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr);
+ ierr = PetscFree2(bases,ldims);CHKERRQ(ierr);
dd->m = m; dd->n = n;
- dd->xs = xs; dd->xe = xe; dd->ys = ys; dd->ye = ye; dd->zs = 0; dd->ze = 1;
- dd->Xs = Xs; dd->Xe = Xe; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = 0; dd->Ze = 1;
+ /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
+ dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = 0; dd->ze = 1;
+ dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = 0; dd->Ze = 1;
ierr = VecDestroy(local);CHKERRQ(ierr);
ierr = VecDestroy(global);CHKERRQ(ierr);
dd->gtol = gtol;
dd->ltog = ltog;
- dd->idx = idx;
- dd->Nl = nn;
+ dd->idx = idx_cpy;
+ dd->Nl = nn*dof;
dd->base = base;
da->ops->view = DMView_DA_2d;
-
- /*
- Set the local to global ordering in the global vector, this allows use
- of VecSetValuesLocal().
- */
- ierr = ISLocalToGlobalMappingCreate(comm,nn,idx,PETSC_OWN_POINTER,&da->ltogmap);CHKERRQ(ierr);
- ierr = ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);CHKERRQ(ierr);
- ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr);
-
dd->ltol = PETSC_NULL;
dd->ao = PETSC_NULL;
diff -r 1bde78661028 src/dm/impls/da/da3.c
--- a/src/dm/impls/da/da3.c Thu Mar 10 12:10:19 2011 -0600
+++ b/src/dm/impls/da/da3.c Thu Mar 10 13:13:17 2011 -0700
@@ -180,17 +180,19 @@
PetscInt *lz = dd->lz;
MPI_Comm comm;
PetscMPIInt rank,size;
- PetscInt xs = 0,xe,ys = 0,ye,zs = 0,ze,x = 0,y = 0,z = 0,Xs,Xe,Ys,Ye,Zs,Ze,start,end,pm;
- PetscInt left,up,down,bottom,top,i,j,k,*idx,nn;
+ PetscInt xs = 0,xe,ys = 0,ye,zs = 0,ze,x = 0,y = 0,z = 0;
+ PetscInt Xs,Xe,Ys,Ye,Zs,Ze,IXs,IXe,IYs,IYe,IZs,IZe,start,end,pm;
+ PetscInt left,right,up,down,bottom,top,i,j,k,*idx,*idx_cpy,nn;
+ const PetscInt *idx_full;
PetscInt n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n14;
PetscInt n15,n16,n17,n18,n19,n20,n21,n22,n23,n24,n25,n26;
- PetscInt *bases,*ldims,x_t,y_t,z_t,s_t,base,count,s_x,s_y,s_z;
+ PetscInt *bases,*ldims,base,x_t,y_t,z_t,s_t,count,count_dbg,s_x,s_y,s_z;
PetscInt sn0 = 0,sn1 = 0,sn2 = 0,sn3 = 0,sn5 = 0,sn6 = 0,sn7 = 0;
PetscInt sn8 = 0,sn9 = 0,sn11 = 0,sn15 = 0,sn24 = 0,sn25 = 0,sn26 = 0;
PetscInt sn17 = 0,sn18 = 0,sn19 = 0,sn20 = 0,sn21 = 0,sn23 = 0;
Vec local,global;
VecScatter ltog,gtol;
- IS to,from;
+ IS to,from,ltogis;
PetscErrorCode ierr;
PetscFunctionBegin;
@@ -329,38 +331,23 @@
/* determine ghost region */
/* Assume No Periodicity */
- if (xs-s > 0) Xs = xs - s; else Xs = 0;
- if (ys-s > 0) Ys = ys - s; else Ys = 0;
- if (zs-s > 0) Zs = zs - s; else Zs = 0;
- if (xe+s <= M) Xe = xe + s; else Xe = M;
- if (ye+s <= N) Ye = ye + s; else Ye = N;
- if (ze+s <= P) Ze = ze + s; else Ze = P;
+ if (xs-s > 0) { Xs = xs - s; IXs = xs - s; } else { Xs = 0; IXs = 0; }
+ if (xe+s <= M) { Xe = xe + s; IXe = xe + s; } else { Xe = M; IXe = M; }
+ if (ys-s > 0) { Ys = ys - s; IYs = ys - s; } else { Ys = 0; IYs = 0; }
+ if (ye+s <= N) { Ye = ye + s; IYe = ye + s; } else { Ye = N; IYe = N; }
+ if (zs-s > 0) { Zs = zs - s; IZs = zs - s; } else { Zs = 0; IZs = 0; }
+ if (ze+s <= P) { Ze = ze + s; IZe = ze + s; } else { Ze = P; IZe = P; }
- /* X Periodic */
- if (DMDAXPeriodic(wrap)){
- Xs = xs - s;
- Xe = xe + s;
- }
-
- /* Y Periodic */
- if (DMDAYPeriodic(wrap)){
- Ys = ys - s;
- Ye = ye + s;
- }
-
- /* Z Periodic */
- if (DMDAZPeriodic(wrap)){
- Zs = zs - s;
- Ze = ze + s;
- }
+ /* fix for periodicity/ghosted */
+ if (DMDAXGhosted(wrap)) { Xs = xs - s; Xe = xe + s; }
+ if (DMDAXPeriodic(wrap)) { IXs = xs - s; IXe = xe + s; }
+ if (DMDAYGhosted(wrap)) { Ys = ys - s; Ye = ye + s; }
+ if (DMDAYPeriodic(wrap)) { IYs = ys - s; IYe = ye + s; }
+ if (DMDAZGhosted(wrap)) { Zs = zs - s; Ze = ze + s; }
+ if (DMDAZPeriodic(wrap)) { IZs = zs - s; IZe = ze + s; }
/* Resize all X parameters to reflect w */
- x *= dof;
- xs *= dof;
- xe *= dof;
- Xs *= dof;
- Xe *= dof;
- s_x = s*dof;
+ s_x = s;
s_y = s;
s_z = s;
@@ -375,80 +362,108 @@
for (i=1; i<=size; i++) {
bases[i] += bases[i-1];
}
+ base = bases[rank]*dof;
/* allocate the base parallel and sequential vectors */
- dd->Nlocal = x*y*z;
+ dd->Nlocal = x*y*z*dof;
ierr = VecCreateMPIWithArray(comm,dd->Nlocal,PETSC_DECIDE,0,&global);CHKERRQ(ierr);
ierr = VecSetBlockSize(global,dof);CHKERRQ(ierr);
- dd->nlocal = (Xe-Xs)*(Ye-Ys)*(Ze-Zs);
- ierr = VecCreateSeqWithArray(MPI_COMM_SELF,dd->nlocal,0,&local);CHKERRQ(ierr);
+ dd->nlocal = (Xe-Xs)*(Ye-Ys)*(Ze-Zs)*dof;
+ ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dd->nlocal,0,&local);CHKERRQ(ierr);
ierr = VecSetBlockSize(local,dof);CHKERRQ(ierr);
/* generate appropriate vector scatters */
/* local to global inserts non-ghost point region into global */
ierr = VecGetOwnershipRange(global,&start,&end);CHKERRQ(ierr);
- ierr = ISCreateStride(comm,x*y*z,start,1,&to);CHKERRQ(ierr);
+ ierr = ISCreateStride(comm,x*y*z*dof,start,1,&to);CHKERRQ(ierr);
- left = xs - Xs;
+ count_dbg = x*y*z;
+ ierr = PetscMalloc(x*y*z*sizeof(PetscInt),&idx);CHKERRQ(ierr);
+ left = xs - Xs; right = left + x;
bottom = ys - Ys; top = bottom + y;
down = zs - Zs; up = down + z;
- count = x*(top-bottom)*(up-down);
- ierr = PetscMalloc(count*sizeof(PetscInt)/dof,&idx);CHKERRQ(ierr);
count = 0;
for (i=down; i<up; i++) {
for (j=bottom; j<top; j++) {
- for (k=0; k<x; k += dof) {
- idx[count++] = ((left+j*(Xe-Xs))+i*(Xe-Xs)*(Ye-Ys) + k)/dof;
+ for (k=left; k<right; k++) {
+ idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
}
}
}
+ if (count != count_dbg) {
+ SETERRQ(((PetscObject)da)->comm, PETSC_ERR_SUP,"count != count_dbg");
+ PetscFunctionReturn(1);
+ }
ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&from);CHKERRQ(ierr);
-
ierr = VecScatterCreate(local,from,global,to,<og);CHKERRQ(ierr);
ierr = PetscLogObjectParent(da,ltog);CHKERRQ(ierr);
ierr = ISDestroy(from);CHKERRQ(ierr);
ierr = ISDestroy(to);CHKERRQ(ierr);
- /* global to local must include ghost points */
+ /* global to local must include ghost points within the domain,
+ but not ghost points outside the domain that aren't periodic */
if (stencil_type == DMDA_STENCIL_BOX) {
- ierr = ISCreateStride(comm,(Xe-Xs)*(Ye-Ys)*(Ze-Zs),0,1,&to);CHKERRQ(ierr);
+ count_dbg = (IXe-IXs)*(IYe-IYs)*(IZe-IZs);
+ ierr = PetscMalloc(count_dbg*sizeof(PetscInt),&idx);CHKERRQ(ierr);
+
+ left = IXs - Xs; right = left + (IXe-IXs);
+ bottom = IYs - Ys; top = bottom + (IYe-IYs);
+ down = IZs - Zs; up = down + (IZe-IZs);
+ count = 0;
+ for (i=down; i<up; i++) {
+ for (j=bottom; j<top; j++) {
+ for (k=left; k<right; k++) {
+ idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
+ }
+ }
+ }
+ if (count != count_dbg) {
+ SETERRQ(((PetscObject)da)->comm, PETSC_ERR_SUP,"count != count_dbg");
+ PetscFunctionReturn(1);
+ }
+ ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr);
+
} else {
/* This is way ugly! We need to list the funny cross type region */
- /* the bottom chunck */
- left = xs - Xs;
+ count_dbg = ((ys-IYs) + (IYe-ye))*x*z + ((xs-IXs) + (IXe-xe))*y*z + ((zs-IZs) + (IZe-ze))*x*y + x*y*z;
+ ierr = PetscMalloc(count_dbg*sizeof(PetscInt),&idx);CHKERRQ(ierr);
+
+ left = xs - Xs; right = left + x;
bottom = ys - Ys; top = bottom + y;
down = zs - Zs; up = down + z;
- count = down*(top-bottom)*x + (up-down)*(bottom*x + (top-bottom)*(Xe-Xs) + (Ye-Ys-top)*x) + (Ze-Zs-up)*(top-bottom)*x;
- ierr = PetscMalloc(count*sizeof(PetscInt)/dof,&idx);CHKERRQ(ierr);
count = 0;
- for (i=0; i<down; i++) {
+ /* the bottom chunck */
+ for (i=(IZs-Zs); i<down; i++) {
for (j=bottom; j<top; j++) {
- for (k=0; k<x; k += dof) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
+ for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
}
}
/* the middle piece */
for (i=down; i<up; i++) {
/* front */
- for (j=0; j<bottom; j++) {
- for (k=0; k<x; k += dof) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
+ for (j=(IYs-Ys); j<bottom; j++) {
+ for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
}
/* middle */
for (j=bottom; j<top; j++) {
- for (k=0; k<Xe-Xs; k += dof) idx[count++] = j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
+ for (k=IXs-Xs; k<IXe-Xs; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
}
/* back */
- for (j=top; j<Ye-Ys; j++) {
- for (k=0; k<x; k += dof) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
+ for (j=top; j<top+IYe-ye; j++) {
+ for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
}
}
/* the top piece */
- for (i=up; i<Ze-Zs; i++) {
+ for (i=up; i<up+IZe-ze; i++) {
for (j=bottom; j<top; j++) {
- for (k=0; k<x; k += dof) idx[count++] = left+j*(Xe-Xs)+i*(Xe-Xs)*(Ye-Ys)+k;
+ for (k=left; k<right; k++) idx[count++] = (i*(Ye-Ys) + j)*(Xe-Xs) + k;
}
}
- for (i=0; i<count; i++) idx[i] = idx[i]/dof;
+ if (count != count_dbg) {
+ SETERRQ(((PetscObject)da)->comm, PETSC_ERR_SUP,"count != count_dbg");
+ PetscFunctionReturn(1);
+ }
ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr);
}
@@ -464,11 +479,9 @@
n3 n4 n5
n0 n1 n2
*/
-
+
/* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */
-
/* Assume Nodes are Internal to the Cube */
-
n0 = rank - m*n - m - 1;
n1 = rank - m*n - m;
n2 = rank - m*n - m + 1;
@@ -512,7 +525,7 @@
n24 = rank + 2*m -1 + (m*n);
}
- if (xe == M*dof) { /* First assume not corner or edge */
+ if (xe == M) { /* First assume not corner or edge */
n2 = rank -2*m +1 - (m*n);
n5 = rank - m +1 - (m*n);
n8 = rank +1 - (m*n);
@@ -596,25 +609,25 @@
n24 = rank - m*(n-1) + m-1 + m*n;
}
- if ((xe==M*dof) && (zs==0)) { /* Assume an edge, not corner */
+ if ((xe==M) && (zs==0)) { /* Assume an edge, not corner */
n2 = size - (m*n-rank) - (m-1) - m;
n5 = size - (m*n-rank) - (m-1);
n8 = size - (m*n-rank) - (m-1) + m;
}
- if ((xe==M*dof) && (ze==P)) { /* Assume an edge, not corner */
+ if ((xe==M) && (ze==P)) { /* Assume an edge, not corner */
n20 = m*n - (size - rank) - (m-1) - m;
n23 = m*n - (size - rank) - (m-1);
n26 = m*n - (size - rank) - (m-1) + m;
}
- if ((xe==M*dof) && (ys==0)) { /* Assume an edge, not corner */
+ if ((xe==M) && (ys==0)) { /* Assume an edge, not corner */
n2 = rank + m*(n-1) - (m-1) - m*n;
n11 = rank + m*(n-1) - (m-1);
n20 = rank + m*(n-1) - (m-1) + m*n;
}
- if ((xe==M*dof) && (ye==N)) { /* Assume an edge, not corner */
+ if ((xe==M) && (ye==N)) { /* Assume an edge, not corner */
n8 = rank - m*n +1 - m*n;
n17 = rank - m*n +1;
n26 = rank - m*n +1 + m*n;
@@ -649,30 +662,27 @@
if ((xs==0) && (ys==0) && (ze==P)) { n18 = m*n-1;}
if ((xs==0) && (ye==N) && (zs==0)) { n6 = (size-1)-m*(n-1);}
if ((xs==0) && (ye==N) && (ze==P)) { n24 = m-1;}
- if ((xe==M*dof) && (ys==0) && (zs==0)) { n2 = size-m;}
- if ((xe==M*dof) && (ys==0) && (ze==P)) { n20 = m*n-m;}
- if ((xe==M*dof) && (ye==N) && (zs==0)) { n8 = size-m*n;}
- if ((xe==M*dof) && (ye==N) && (ze==P)) { n26 = 0;}
+ if ((xe==M) && (ys==0) && (zs==0)) { n2 = size-m;}
+ if ((xe==M) && (ys==0) && (ze==P)) { n20 = m*n-m;}
+ if ((xe==M) && (ye==N) && (zs==0)) { n8 = size-m*n;}
+ if ((xe==M) && (ye==N) && (ze==P)) { n26 = 0;}
/* Check for when not X,Y, and Z Periodic */
/* If not X periodic */
- if ((wrap != DMDA_XPERIODIC) && (wrap != DMDA_XYPERIODIC) &&
- (wrap != DMDA_XZPERIODIC) && (wrap != DMDA_XYZPERIODIC)) {
+ if (!DMDAXPeriodic(wrap)) {
if (xs==0) {n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2;}
- if (xe==M*dof) {n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2;}
+ if (xe==M) {n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2;}
}
/* If not Y periodic */
- if ((wrap != DMDA_YPERIODIC) && (wrap != DMDA_XYPERIODIC) &&
- (wrap != DMDA_YZPERIODIC) && (wrap != DMDA_XYZPERIODIC)) {
+ if (!DMDAYPeriodic(wrap)) {
if (ys==0) {n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2;}
if (ye==N) {n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2;}
}
/* If not Z periodic */
- if ((wrap != DMDA_ZPERIODIC) && (wrap != DMDA_XZPERIODIC) &&
- (wrap != DMDA_YZPERIODIC) && (wrap != DMDA_XYZPERIODIC)) {
+ if (!DMDAZPeriodic(wrap)) {
if (zs==0) {n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2;}
if (ze==P) {n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2;}
}
@@ -722,14 +732,13 @@
ierr = PetscLogObjectMemory(da,(Xe-Xs)*(Ye-Ys)*(Ze-Zs)*sizeof(PetscInt));CHKERRQ(ierr);
nn = 0;
-
/* Bottom Level */
for (k=0; k<s_z; k++) {
for (i=1; i<=s_y; i++) {
if (n0 >= 0) { /* left below */
- x_t = lx[n0 % m]*dof;
- y_t = ly[(n0 % (m*n))/m];
- z_t = lz[n0 / (m*n)];
+ x_t = lx[n0 % m];
+ y_t = ly[(n0 % (m*n))/m];
+ z_t = lz[n0 / (m*n)];
s_t = bases[n0] + x_t*y_t*z_t - (s_y-i)*x_t - s_x - (s_z-k-1)*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
}
@@ -741,7 +750,7 @@
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
}
if (n2 >= 0) { /* right below */
- x_t = lx[n2 % m]*dof;
+ x_t = lx[n2 % m];
y_t = ly[(n2 % (m*n))/m];
z_t = lz[n2 / (m*n)];
s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
@@ -751,7 +760,7 @@
for (i=0; i<y; i++) {
if (n3 >= 0) { /* directly left */
- x_t = lx[n3 % m]*dof;
+ x_t = lx[n3 % m];
y_t = y;
z_t = lz[n3 / (m*n)];
s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
@@ -767,7 +776,7 @@
}
if (n5 >= 0) { /* directly right */
- x_t = lx[n5 % m]*dof;
+ x_t = lx[n5 % m];
y_t = y;
z_t = lz[n5 / (m*n)];
s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
@@ -777,7 +786,7 @@
for (i=1; i<=s_y; i++) {
if (n6 >= 0) { /* left above */
- x_t = lx[n6 % m]*dof;
+ x_t = lx[n6 % m];
y_t = ly[(n6 % (m*n))/m];
z_t = lz[n6 / (m*n)];
s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
@@ -791,7 +800,7 @@
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
}
if (n8 >= 0) { /* right above */
- x_t = lx[n8 % m]*dof;
+ x_t = lx[n8 % m];
y_t = ly[(n8 % (m*n))/m];
z_t = lz[n8 / (m*n)];
s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
@@ -804,7 +813,7 @@
for (k=0; k<z; k++) {
for (i=1; i<=s_y; i++) {
if (n9 >= 0) { /* left below */
- x_t = lx[n9 % m]*dof;
+ x_t = lx[n9 % m];
y_t = ly[(n9 % (m*n))/m];
/* z_t = z; */
s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
@@ -818,7 +827,7 @@
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
}
if (n11 >= 0) { /* right below */
- x_t = lx[n11 % m]*dof;
+ x_t = lx[n11 % m];
y_t = ly[(n11 % (m*n))/m];
/* z_t = z; */
s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
@@ -828,7 +837,7 @@
for (i=0; i<y; i++) {
if (n12 >= 0) { /* directly left */
- x_t = lx[n12 % m]*dof;
+ x_t = lx[n12 % m];
y_t = y;
/* z_t = z; */
s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
@@ -840,7 +849,7 @@
for (j=0; j<x; j++) { idx[nn++] = s_t++;}
if (n14 >= 0) { /* directly right */
- x_t = lx[n14 % m]*dof;
+ x_t = lx[n14 % m];
y_t = y;
/* z_t = z; */
s_t = bases[n14] + i*x_t + k*x_t*y_t;
@@ -850,7 +859,7 @@
for (i=1; i<=s_y; i++) {
if (n15 >= 0) { /* left above */
- x_t = lx[n15 % m]*dof;
+ x_t = lx[n15 % m];
y_t = ly[(n15 % (m*n))/m];
/* z_t = z; */
s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t;
@@ -864,7 +873,7 @@
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
}
if (n17 >= 0) { /* right above */
- x_t = lx[n17 % m]*dof;
+ x_t = lx[n17 % m];
y_t = ly[(n17 % (m*n))/m];
/* z_t = z; */
s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t;
@@ -877,7 +886,7 @@
for (k=0; k<s_z; k++) {
for (i=1; i<=s_y; i++) {
if (n18 >= 0) { /* left below */
- x_t = lx[n18 % m]*dof;
+ x_t = lx[n18 % m];
y_t = ly[(n18 % (m*n))/m];
/* z_t = lz[n18 / (m*n)]; */
s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
@@ -891,7 +900,7 @@
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
}
if (n20 >= 0) { /* right below */
- x_t = lx[n20 % m]*dof;
+ x_t = lx[n20 % m];
y_t = ly[(n20 % (m*n))/m];
/* z_t = lz[n20 / (m*n)]; */
s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
@@ -901,7 +910,7 @@
for (i=0; i<y; i++) {
if (n21 >= 0) { /* directly left */
- x_t = lx[n21 % m]*dof;
+ x_t = lx[n21 % m];
y_t = y;
/* z_t = lz[n21 / (m*n)]; */
s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t;
@@ -917,7 +926,7 @@
}
if (n23 >= 0) { /* directly right */
- x_t = lx[n23 % m]*dof;
+ x_t = lx[n23 % m];
y_t = y;
/* z_t = lz[n23 / (m*n)]; */
s_t = bases[n23] + i*x_t + k*x_t*y_t;
@@ -927,7 +936,7 @@
for (i=1; i<=s_y; i++) {
if (n24 >= 0) { /* left above */
- x_t = lx[n24 % m]*dof;
+ x_t = lx[n24 % m];
y_t = ly[(n24 % (m*n))/m];
/* z_t = lz[n24 / (m*n)]; */
s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
@@ -941,56 +950,48 @@
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
}
if (n26 >= 0) { /* right above */
- x_t = lx[n26 % m]*dof;
+ x_t = lx[n26 % m];
y_t = ly[(n26 % (m*n))/m];
/* z_t = lz[n26 / (m*n)]; */
s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
}
}
- }
- base = bases[rank];
- {
- PetscInt nnn = nn/dof,*iidx;
- ierr = PetscMalloc(nnn*sizeof(PetscInt),&iidx);CHKERRQ(ierr);
- for (i=0; i<nnn; i++) {
- iidx[i] = idx[dof*i]/dof;
- }
- ierr = ISCreateBlock(comm,dof,nnn,iidx,PETSC_OWN_POINTER,&from);CHKERRQ(ierr);
}
+
+ ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr);
ierr = PetscLogObjectParent(da,gtol);CHKERRQ(ierr);
+ ierr = ISDestroy(to);CHKERRQ(ierr);
ierr = ISDestroy(from);CHKERRQ(ierr);
- ierr = ISDestroy(to);CHKERRQ(ierr);
- dd->m = m; dd->n = n; dd->p = p;
- dd->xs = xs; dd->xe = xe; dd->ys = ys; dd->ye = ye; dd->zs = zs; dd->ze = ze;
- dd->Xs = Xs; dd->Xe = Xe; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = Zs; dd->Ze = Ze;
+ if (stencil_type == DMDA_STENCIL_STAR) {
+ n0 = sn0; n1 = sn1; n2 = sn2; n3 = sn3; n5 = sn5; n6 = sn6; n7 = sn7;
+ n8 = sn8; n9 = sn9; n11 = sn11; n15 = sn15; n17 = sn17; n18 = sn18;
+ n19 = sn19; n20 = sn20; n21 = sn21; n23 = sn23; n24 = sn24; n25 = sn25;
+ n26 = sn26;
+ }
- ierr = VecDestroy(local);CHKERRQ(ierr);
- ierr = VecDestroy(global);CHKERRQ(ierr);
-
- if (stencil_type == DMDA_STENCIL_STAR) {
+ if ((stencil_type == DMDA_STENCIL_STAR) ||
+ (!DMDAXPeriodic(wrap) && DMDAXGhosted(wrap)) ||
+ (!DMDAYPeriodic(wrap) && DMDAYGhosted(wrap)) ||
+ (!DMDAZPeriodic(wrap) && DMDAZGhosted(wrap))) {
/*
Recompute the local to global mappings, this time keeping the
information about the cross corner processor numbers.
*/
- n0 = sn0; n1 = sn1; n2 = sn2; n3 = sn3; n5 = sn5; n6 = sn6; n7 = sn7;
- n8 = sn8; n9 = sn9; n11 = sn11; n15 = sn15; n17 = sn17; n18 = sn18;
- n19 = sn19; n20 = sn20; n21 = sn21; n23 = sn23; n24 = sn24; n25 = sn25;
- n26 = sn26;
-
nn = 0;
-
/* Bottom Level */
- for (k=0; k<s_z; k++) {
+ for (k=0; k<s_z; k++) {
for (i=1; i<=s_y; i++) {
if (n0 >= 0) { /* left below */
- x_t = lx[n0 % m]*dof;
- y_t = ly[(n0 % (m*n))/m];
- z_t = lz[n0 / (m*n)];
+ x_t = lx[n0 % m];
+ y_t = ly[(n0 % (m*n))/m];
+ z_t = lz[n0 / (m*n)];
s_t = bases[n0] + x_t*y_t*z_t - (s_y-i)*x_t - s_x - (s_z-k-1)*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xs-xs < 0 && Ys-ys < 0 && Zs-zs < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
if (n1 >= 0) { /* directly below */
x_t = x;
@@ -998,23 +999,29 @@
z_t = lz[n1 / (m*n)];
s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
+ } else if (Ys-ys < 0 && Zs-zs < 0) {
+ for (j=0; j<x; j++) { idx[nn++] = -1;}
}
if (n2 >= 0) { /* right below */
- x_t = lx[n2 % m]*dof;
+ x_t = lx[n2 % m];
y_t = ly[(n2 % (m*n))/m];
z_t = lz[n2 / (m*n)];
s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xe-Xe < 0 && Ys-ys < 0 && Zs-zs < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
for (i=0; i<y; i++) {
if (n3 >= 0) { /* directly left */
- x_t = lx[n3 % m]*dof;
+ x_t = lx[n3 % m];
y_t = y;
z_t = lz[n3 / (m*n)];
s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xs-xs < 0 && Zs-zs < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
if (n4 >= 0) { /* middle */
@@ -1023,24 +1030,30 @@
z_t = lz[n4 / (m*n)];
s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
+ } else if (Zs-zs < 0) {
+ for (j=0; j<x; j++) { idx[nn++] = -1;}
}
if (n5 >= 0) { /* directly right */
- x_t = lx[n5 % m]*dof;
+ x_t = lx[n5 % m];
y_t = y;
z_t = lz[n5 / (m*n)];
s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xe-Xe < 0 && Zs-zs < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
for (i=1; i<=s_y; i++) {
if (n6 >= 0) { /* left above */
- x_t = lx[n6 % m]*dof;
+ x_t = lx[n6 % m];
y_t = ly[(n6 % (m*n))/m];
z_t = lz[n6 / (m*n)];
s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xs-xs < 0 && ye-Ye < 0 && Zs-zs < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
if (n7 >= 0) { /* directly above */
x_t = x;
@@ -1048,13 +1061,17 @@
z_t = lz[n7 / (m*n)];
s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
+ } else if (ye-Ye < 0 && Zs-zs < 0) {
+ for (j=0; j<x; j++) { idx[nn++] = -1;}
}
if (n8 >= 0) { /* right above */
- x_t = lx[n8 % m]*dof;
+ x_t = lx[n8 % m];
y_t = ly[(n8 % (m*n))/m];
z_t = lz[n8 / (m*n)];
s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xe-Xe < 0 && ye-Ye < 0 && Zs-zs < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
}
@@ -1063,11 +1080,13 @@
for (k=0; k<z; k++) {
for (i=1; i<=s_y; i++) {
if (n9 >= 0) { /* left below */
- x_t = lx[n9 % m]*dof;
+ x_t = lx[n9 % m];
y_t = ly[(n9 % (m*n))/m];
/* z_t = z; */
s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xs-xs < 0 && Ys-ys < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
if (n10 >= 0) { /* directly below */
x_t = x;
@@ -1075,23 +1094,29 @@
/* z_t = z; */
s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
+ } else if (Ys-ys < 0) {
+ for (j=0; j<x; j++) { idx[nn++] = -1;}
}
if (n11 >= 0) { /* right below */
- x_t = lx[n11 % m]*dof;
+ x_t = lx[n11 % m];
y_t = ly[(n11 % (m*n))/m];
/* z_t = z; */
s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xe-Xe < 0 && Ys-ys < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
for (i=0; i<y; i++) {
if (n12 >= 0) { /* directly left */
- x_t = lx[n12 % m]*dof;
+ x_t = lx[n12 % m];
y_t = y;
/* z_t = z; */
s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xs-xs < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
/* Interior */
@@ -1099,21 +1124,25 @@
for (j=0; j<x; j++) { idx[nn++] = s_t++;}
if (n14 >= 0) { /* directly right */
- x_t = lx[n14 % m]*dof;
+ x_t = lx[n14 % m];
y_t = y;
/* z_t = z; */
s_t = bases[n14] + i*x_t + k*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xe-Xe < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
for (i=1; i<=s_y; i++) {
if (n15 >= 0) { /* left above */
- x_t = lx[n15 % m]*dof;
+ x_t = lx[n15 % m];
y_t = ly[(n15 % (m*n))/m];
/* z_t = z; */
s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xs-xs < 0 && ye-Ye < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
if (n16 >= 0) { /* directly above */
x_t = x;
@@ -1121,13 +1150,17 @@
/* z_t = z; */
s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t;
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
+ } else if (ye-Ye < 0) {
+ for (j=0; j<x; j++) { idx[nn++] = -1;}
}
if (n17 >= 0) { /* right above */
- x_t = lx[n17 % m]*dof;
+ x_t = lx[n17 % m];
y_t = ly[(n17 % (m*n))/m];
/* z_t = z; */
s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xe-Xe < 0 && ye-Ye < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
}
@@ -1136,11 +1169,13 @@
for (k=0; k<s_z; k++) {
for (i=1; i<=s_y; i++) {
if (n18 >= 0) { /* left below */
- x_t = lx[n18 % m]*dof;
+ x_t = lx[n18 % m];
y_t = ly[(n18 % (m*n))/m];
/* z_t = lz[n18 / (m*n)]; */
s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xs-xs < 0 && Ys-ys < 0 && ze-Ze < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
if (n19 >= 0) { /* directly below */
x_t = x;
@@ -1148,23 +1183,29 @@
/* z_t = lz[n19 / (m*n)]; */
s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
+ } else if (Ys-ys < 0 && ze-Ze < 0) {
+ for (j=0; j<x; j++) { idx[nn++] = -1;}
}
if (n20 >= 0) { /* right below */
- x_t = lx[n20 % m]*dof;
+ x_t = lx[n20 % m];
y_t = ly[(n20 % (m*n))/m];
/* z_t = lz[n20 / (m*n)]; */
s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xe-Xe < 0 && Ys-ys < 0 && ze-Ze < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
for (i=0; i<y; i++) {
if (n21 >= 0) { /* directly left */
- x_t = lx[n21 % m]*dof;
+ x_t = lx[n21 % m];
y_t = y;
/* z_t = lz[n21 / (m*n)]; */
s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xs-xs < 0 && ze-Ze < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
if (n22 >= 0) { /* middle */
@@ -1173,24 +1214,30 @@
/* z_t = lz[n22 / (m*n)]; */
s_t = bases[n22] + i*x_t + k*x_t*y_t;
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
+ } else if (ze-Ze < 0) {
+ for (j=0; j<x; j++) { idx[nn++] = -1;}
}
if (n23 >= 0) { /* directly right */
- x_t = lx[n23 % m]*dof;
+ x_t = lx[n23 % m];
y_t = y;
/* z_t = lz[n23 / (m*n)]; */
s_t = bases[n23] + i*x_t + k*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xe-Xe < 0 && ze-Ze < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
for (i=1; i<=s_y; i++) {
if (n24 >= 0) { /* left above */
- x_t = lx[n24 % m]*dof;
- y_t = ly[(n24 % (m*n))/m];
+ x_t = lx[n24 % m];
+ y_t = ly[(n24 % (m*n))/m];
/* z_t = lz[n24 / (m*n)]; */
s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (Xs-xs < 0 && ye-Ye < 0 && ze-Ze < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
if (n25 >= 0) { /* directly above */
x_t = x;
@@ -1198,473 +1245,60 @@
/* z_t = lz[n25 / (m*n)]; */
s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t;
for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
+ } else if (ye-Ye < 0 && ze-Ze < 0) {
+ for (j=0; j<x; j++) { idx[nn++] = -1;}
}
if (n26 >= 0) { /* right above */
- x_t = lx[n26 % m]*dof;
- y_t = ly[(n26 % (m*n))/m];
+ x_t = lx[n26 % m];
+ y_t = ly[(n26 % (m*n))/m];
/* z_t = lz[n26 / (m*n)]; */
s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
+ } else if (xe-Xe < 0 && ye-Ye < 0 && ze-Ze < 0) {
+ for (j=0; j<s_x; j++) { idx[nn++] = -1;}
}
}
- }
- }
- /* redo idx to include "missing" ghost points */
- /* Solve for X,Y, and Z Periodic Case First, Then Modify Solution */
-
- /* Assume Nodes are Internal to the Cube */
-
- n0 = rank - m*n - m - 1;
- n1 = rank - m*n - m;
- n2 = rank - m*n - m + 1;
- n3 = rank - m*n -1;
- n4 = rank - m*n;
- n5 = rank - m*n + 1;
- n6 = rank - m*n + m - 1;
- n7 = rank - m*n + m;
- n8 = rank - m*n + m + 1;
-
- n9 = rank - m - 1;
- n10 = rank - m;
- n11 = rank - m + 1;
- n12 = rank - 1;
- n14 = rank + 1;
- n15 = rank + m - 1;
- n16 = rank + m;
- n17 = rank + m + 1;
-
- n18 = rank + m*n - m - 1;
- n19 = rank + m*n - m;
- n20 = rank + m*n - m + 1;
- n21 = rank + m*n - 1;
- n22 = rank + m*n;
- n23 = rank + m*n + 1;
- n24 = rank + m*n + m - 1;
- n25 = rank + m*n + m;
- n26 = rank + m*n + m + 1;
-
- /* Assume Pieces are on Faces of Cube */
-
- if (xs == 0) { /* First assume not corner or edge */
- n0 = rank -1 - (m*n);
- n3 = rank + m -1 - (m*n);
- n6 = rank + 2*m -1 - (m*n);
- n9 = rank -1;
- n12 = rank + m -1;
- n15 = rank + 2*m -1;
- n18 = rank -1 + (m*n);
- n21 = rank + m -1 + (m*n);
- n24 = rank + 2*m -1 + (m*n);
- }
-
- if (xe == M*dof) { /* First assume not corner or edge */
- n2 = rank -2*m +1 - (m*n);
- n5 = rank - m +1 - (m*n);
- n8 = rank +1 - (m*n);
- n11 = rank -2*m +1;
- n14 = rank - m +1;
- n17 = rank +1;
- n20 = rank -2*m +1 + (m*n);
- n23 = rank - m +1 + (m*n);
- n26 = rank +1 + (m*n);
- }
-
- if (ys==0) { /* First assume not corner or edge */
- n0 = rank + m * (n-1) -1 - (m*n);
- n1 = rank + m * (n-1) - (m*n);
- n2 = rank + m * (n-1) +1 - (m*n);
- n9 = rank + m * (n-1) -1;
- n10 = rank + m * (n-1);
- n11 = rank + m * (n-1) +1;
- n18 = rank + m * (n-1) -1 + (m*n);
- n19 = rank + m * (n-1) + (m*n);
- n20 = rank + m * (n-1) +1 + (m*n);
- }
-
- if (ye == N) { /* First assume not corner or edge */
- n6 = rank - m * (n-1) -1 - (m*n);
- n7 = rank - m * (n-1) - (m*n);
- n8 = rank - m * (n-1) +1 - (m*n);
- n15 = rank - m * (n-1) -1;
- n16 = rank - m * (n-1);
- n17 = rank - m * (n-1) +1;
- n24 = rank - m * (n-1) -1 + (m*n);
- n25 = rank - m * (n-1) + (m*n);
- n26 = rank - m * (n-1) +1 + (m*n);
- }
-
- if (zs == 0) { /* First assume not corner or edge */
- n0 = size - (m*n) + rank - m - 1;
- n1 = size - (m*n) + rank - m;
- n2 = size - (m*n) + rank - m + 1;
- n3 = size - (m*n) + rank - 1;
- n4 = size - (m*n) + rank;
- n5 = size - (m*n) + rank + 1;
- n6 = size - (m*n) + rank + m - 1;
- n7 = size - (m*n) + rank + m ;
- n8 = size - (m*n) + rank + m + 1;
- }
-
- if (ze == P) { /* First assume not corner or edge */
- n18 = (m*n) - (size-rank) - m - 1;
- n19 = (m*n) - (size-rank) - m;
- n20 = (m*n) - (size-rank) - m + 1;
- n21 = (m*n) - (size-rank) - 1;
- n22 = (m*n) - (size-rank);
- n23 = (m*n) - (size-rank) + 1;
- n24 = (m*n) - (size-rank) + m - 1;
- n25 = (m*n) - (size-rank) + m;
- n26 = (m*n) - (size-rank) + m + 1;
- }
-
- if ((xs==0) && (zs==0)) { /* Assume an edge, not corner */
- n0 = size - m*n + rank + m-1 - m;
- n3 = size - m*n + rank + m-1;
- n6 = size - m*n + rank + m-1 + m;
- }
-
- if ((xs==0) && (ze==P)) { /* Assume an edge, not corner */
- n18 = m*n - (size - rank) + m-1 - m;
- n21 = m*n - (size - rank) + m-1;
- n24 = m*n - (size - rank) + m-1 + m;
- }
-
- if ((xs==0) && (ys==0)) { /* Assume an edge, not corner */
- n0 = rank + m*n -1 - m*n;
- n9 = rank + m*n -1;
- n18 = rank + m*n -1 + m*n;
- }
-
- if ((xs==0) && (ye==N)) { /* Assume an edge, not corner */
- n6 = rank - m*(n-1) + m-1 - m*n;
- n15 = rank - m*(n-1) + m-1;
- n24 = rank - m*(n-1) + m-1 + m*n;
- }
-
- if ((xe==M*dof) && (zs==0)) { /* Assume an edge, not corner */
- n2 = size - (m*n-rank) - (m-1) - m;
- n5 = size - (m*n-rank) - (m-1);
- n8 = size - (m*n-rank) - (m-1) + m;
- }
-
- if ((xe==M*dof) && (ze==P)) { /* Assume an edge, not corner */
- n20 = m*n - (size - rank) - (m-1) - m;
- n23 = m*n - (size - rank) - (m-1);
- n26 = m*n - (size - rank) - (m-1) + m;
- }
-
- if ((xe==M*dof) && (ys==0)) { /* Assume an edge, not corner */
- n2 = rank + m*(n-1) - (m-1) - m*n;
- n11 = rank + m*(n-1) - (m-1);
- n20 = rank + m*(n-1) - (m-1) + m*n;
- }
-
- if ((xe==M*dof) && (ye==N)) { /* Assume an edge, not corner */
- n8 = rank - m*n +1 - m*n;
- n17 = rank - m*n +1;
- n26 = rank - m*n +1 + m*n;
- }
-
- if ((ys==0) && (zs==0)) { /* Assume an edge, not corner */
- n0 = size - m + rank -1;
- n1 = size - m + rank;
- n2 = size - m + rank +1;
- }
-
- if ((ys==0) && (ze==P)) { /* Assume an edge, not corner */
- n18 = m*n - (size - rank) + m*(n-1) -1;
- n19 = m*n - (size - rank) + m*(n-1);
- n20 = m*n - (size - rank) + m*(n-1) +1;
- }
-
- if ((ye==N) && (zs==0)) { /* Assume an edge, not corner */
- n6 = size - (m*n-rank) - m * (n-1) -1;
- n7 = size - (m*n-rank) - m * (n-1);
- n8 = size - (m*n-rank) - m * (n-1) +1;
- }
-
- if ((ye==N) && (ze==P)) { /* Assume an edge, not corner */
- n24 = rank - (size-m) -1;
- n25 = rank - (size-m);
- n26 = rank - (size-m) +1;
- }
-
- /* Check for Corners */
- if ((xs==0) && (ys==0) && (zs==0)) { n0 = size -1;}
- if ((xs==0) && (ys==0) && (ze==P)) { n18 = m*n-1;}
- if ((xs==0) && (ye==N) && (zs==0)) { n6 = (size-1)-m*(n-1);}
- if ((xs==0) && (ye==N) && (ze==P)) { n24 = m-1;}
- if ((xe==M*dof) && (ys==0) && (zs==0)) { n2 = size-m;}
- if ((xe==M*dof) && (ys==0) && (ze==P)) { n20 = m*n-m;}
- if ((xe==M*dof) && (ye==N) && (zs==0)) { n8 = size-m*n;}
- if ((xe==M*dof) && (ye==N) && (ze==P)) { n26 = 0;}
-
- /* Check for when not X,Y, and Z Periodic */
-
- /* If not X periodic */
- if (!DMDAXPeriodic(wrap)){
- if (xs==0) {n0 = n3 = n6 = n9 = n12 = n15 = n18 = n21 = n24 = -2;}
- if (xe==M*dof) {n2 = n5 = n8 = n11 = n14 = n17 = n20 = n23 = n26 = -2;}
- }
-
- /* If not Y periodic */
- if (!DMDAYPeriodic(wrap)){
- if (ys==0) {n0 = n1 = n2 = n9 = n10 = n11 = n18 = n19 = n20 = -2;}
- if (ye==N) {n6 = n7 = n8 = n15 = n16 = n17 = n24 = n25 = n26 = -2;}
- }
-
- /* If not Z periodic */
- if (!DMDAZPeriodic(wrap)){
- if (zs==0) {n0 = n1 = n2 = n3 = n4 = n5 = n6 = n7 = n8 = -2;}
- if (ze==P) {n18 = n19 = n20 = n21 = n22 = n23 = n24 = n25 = n26 = -2;}
- }
-
- nn = 0;
-
- /* Bottom Level */
- for (k=0; k<s_z; k++) {
- for (i=1; i<=s_y; i++) {
- if (n0 >= 0) { /* left below */
- x_t = lx[n0 % m]*dof;
- y_t = ly[(n0 % (m*n))/m];
- z_t = lz[n0 / (m*n)];
- s_t = bases[n0] + x_t*y_t*z_t - (s_y-i)*x_t -s_x - (s_z-k-1)*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- if (n1 >= 0) { /* directly below */
- x_t = x;
- y_t = ly[(n1 % (m*n))/m];
- z_t = lz[n1 / (m*n)];
- s_t = bases[n1] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
- for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
- }
- if (n2 >= 0) { /* right below */
- x_t = lx[n2 % m]*dof;
- y_t = ly[(n2 % (m*n))/m];
- z_t = lz[n2 / (m*n)];
- s_t = bases[n2] + x_t*y_t*z_t - (s_y+1-i)*x_t - (s_z-k-1)*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- }
-
- for (i=0; i<y; i++) {
- if (n3 >= 0) { /* directly left */
- x_t = lx[n3 % m]*dof;
- y_t = y;
- z_t = lz[n3 / (m*n)];
- s_t = bases[n3] + (i+1)*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
-
- if (n4 >= 0) { /* middle */
- x_t = x;
- y_t = y;
- z_t = lz[n4 / (m*n)];
- s_t = bases[n4] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
- for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
- }
-
- if (n5 >= 0) { /* directly right */
- x_t = lx[n5 % m]*dof;
- y_t = y;
- z_t = lz[n5 / (m*n)];
- s_t = bases[n5] + i*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- }
-
- for (i=1; i<=s_y; i++) {
- if (n6 >= 0) { /* left above */
- x_t = lx[n6 % m]*dof;
- y_t = ly[(n6 % (m*n))/m];
- z_t = lz[n6 / (m*n)];
- s_t = bases[n6] + i*x_t - s_x + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- if (n7 >= 0) { /* directly above */
- x_t = x;
- y_t = ly[(n7 % (m*n))/m];
- z_t = lz[n7 / (m*n)];
- s_t = bases[n7] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
- for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
- }
- if (n8 >= 0) { /* right above */
- x_t = lx[n8 % m]*dof;
- y_t = ly[(n8 % (m*n))/m];
- z_t = lz[n8 / (m*n)];
- s_t = bases[n8] + (i-1)*x_t + x_t*y_t*z_t - (s_z-k)*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
}
}
-
- /* Middle Level */
- for (k=0; k<z; k++) {
- for (i=1; i<=s_y; i++) {
- if (n9 >= 0) { /* left below */
- x_t = lx[n9 % m]*dof;
- y_t = ly[(n9 % (m*n))/m];
- /* z_t = z; */
- s_t = bases[n9] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- if (n10 >= 0) { /* directly below */
- x_t = x;
- y_t = ly[(n10 % (m*n))/m];
- /* z_t = z; */
- s_t = bases[n10] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
- for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
- }
- if (n11 >= 0) { /* right below */
- x_t = lx[n11 % m]*dof;
- y_t = ly[(n11 % (m*n))/m];
- /* z_t = z; */
- s_t = bases[n11] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- }
-
- for (i=0; i<y; i++) {
- if (n12 >= 0) { /* directly left */
- x_t = lx[n12 % m]*dof;
- y_t = y;
- /* z_t = z; */
- s_t = bases[n12] + (i+1)*x_t - s_x + k*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
-
- /* Interior */
- s_t = bases[rank] + i*x + k*x*y;
- for (j=0; j<x; j++) { idx[nn++] = s_t++;}
-
- if (n14 >= 0) { /* directly right */
- x_t = lx[n14 % m]*dof;
- y_t = y;
- /* z_t = z; */
- s_t = bases[n14] + i*x_t + k*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- }
-
- for (i=1; i<=s_y; i++) {
- if (n15 >= 0) { /* left above */
- x_t = lx[n15 % m]*dof;
- y_t = ly[(n15 % (m*n))/m];
- /* z_t = z; */
- s_t = bases[n15] + i*x_t - s_x + k*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- if (n16 >= 0) { /* directly above */
- x_t = x;
- y_t = ly[(n16 % (m*n))/m];
- /* z_t = z; */
- s_t = bases[n16] + (i-1)*x_t + k*x_t*y_t;
- for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
- }
- if (n17 >= 0) { /* right above */
- x_t = lx[n17 % m]*dof;
- y_t = ly[(n17 % (m*n))/m];
- /* z_t = z; */
- s_t = bases[n17] + (i-1)*x_t + k*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- }
- }
-
- /* Upper Level */
- for (k=0; k<s_z; k++) {
- for (i=1; i<=s_y; i++) {
- if (n18 >= 0) { /* left below */
- x_t = lx[n18 % m]*dof;
- y_t = ly[(n18 % (m*n))/m];
- /* z_t = lz[n18 / (m*n)]; */
- s_t = bases[n18] - (s_y-i)*x_t -s_x + (k+1)*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- if (n19 >= 0) { /* directly below */
- x_t = x;
- y_t = ly[(n19 % (m*n))/m];
- /* z_t = lz[n19 / (m*n)]; */
- s_t = bases[n19] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
- for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
- }
- if (n20 >= 0) { /* right belodof */
- x_t = lx[n20 % m]*dof;
- y_t = ly[(n20 % (m*n))/m];
- /* z_t = lz[n20 / (m*n)]; */
- s_t = bases[n20] - (s_y+1-i)*x_t + (k+1)*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- }
-
- for (i=0; i<y; i++) {
- if (n21 >= 0) { /* directly left */
- x_t = lx[n21 % m]*dof;
- y_t = y;
- /* z_t = lz[n21 / (m*n)]; */
- s_t = bases[n21] + (i+1)*x_t - s_x + k*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
-
- if (n22 >= 0) { /* middle */
- x_t = x;
- y_t = y;
- /* z_t = lz[n22 / (m*n)]; */
- s_t = bases[n22] + i*x_t + k*x_t*y_t;
- for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
- }
-
- if (n23 >= 0) { /* directly right */
- x_t = lx[n23 % m]*dof;
- y_t = y;
- /* z_t = lz[n23 / (m*n)]; */
- s_t = bases[n23] + i*x_t + k*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- }
-
- for (i=1; i<=s_y; i++) {
- if (n24 >= 0) { /* left above */
- x_t = lx[n24 % m]*dof;
- y_t = ly[(n24 % (m*n))/m];
- /* z_t = lz[n24 / (m*n)]; */
- s_t = bases[n24] + i*x_t - s_x + k*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- if (n25 >= 0) { /* directly above */
- x_t = x;
- y_t = ly[(n25 % (m*n))/m];
- /* z_t = lz[n25 / (m*n)]; */
- s_t = bases[n25] + (i-1)*x_t + k*x_t*y_t;
- for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
- }
- if (n26 >= 0) { /* right above */
- x_t = lx[n26 % m]*dof;
- y_t = ly[(n26 % (m*n))/m];
- /* z_t = lz[n26 / (m*n)]; */
- s_t = bases[n26] + (i-1)*x_t + k*x_t*y_t;
- for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
- }
- }
- }
- ierr = PetscFree2(bases,ldims);CHKERRQ(ierr);
- dd->gtol = gtol;
- dd->ltog = ltog;
- dd->idx = idx;
- dd->Nl = nn;
- dd->base = base;
- da->ops->view = DMView_DA_3d;
-
- /*
+ /*
Set the local to global ordering in the global vector, this allows use
of VecSetValuesLocal().
*/
- ierr = ISLocalToGlobalMappingCreate(comm,nn,idx,PETSC_OWN_POINTER,&da->ltogmap);CHKERRQ(ierr);
+ if (nn != (Xe-Xs)*(Ye-Ys)*(Ze-Zs)) {
+ SETERRQ(((PetscObject)da)->comm, PETSC_ERR_SUP,"nn != count_dbg");
+ PetscFunctionReturn(1);
+ }
+ ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_OWN_POINTER,<ogis);CHKERRQ(ierr);
+ ierr = PetscMalloc(nn*dof*sizeof(PetscInt),&idx_cpy);CHKERRQ(ierr);
+ /* ierr = PetscLogObjectMemory(da,nn*dof*sizeof(PetscInt));CHKERRQ(ierr); */
+ ierr = ISGetIndices(ltogis, &idx_full);
+ ierr = PetscMemcpy(idx_cpy,idx_full,nn*dof*sizeof(PetscInt));CHKERRQ(ierr);
+ CHKMEMQ;
+ ierr = ISRestoreIndices(ltogis, &idx_full);
+ ierr = ISLocalToGlobalMappingCreateIS(ltogis,&da->ltogmap);CHKERRQ(ierr);
+ ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr);
+ ierr = ISDestroy(ltogis);CHKERRQ(ierr);
ierr = ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);CHKERRQ(ierr);
ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr);
+ ierr = PetscFree2(bases,ldims);CHKERRQ(ierr);
+ dd->m = m; dd->n = n; dd->p = p;
+ /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
+ dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = zs; dd->ze = ze;
+ dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = Zs; dd->Ze = Ze;
+
+ ierr = VecDestroy(local);CHKERRQ(ierr);
+ ierr = VecDestroy(global);CHKERRQ(ierr);
+
+ dd->gtol = gtol;
+ dd->ltog = ltog;
+ dd->idx = idx_cpy;
+ dd->Nl = nn*dof;
+ dd->base = base;
+ da->ops->view = DMView_DA_3d;
dd->ltol = PETSC_NULL;
dd->ao = PETSC_NULL;
+
PetscFunctionReturn(0);
}
diff -r 1bde78661028 src/dm/impls/da/dacorn.c
--- a/src/dm/impls/da/dacorn.c Thu Mar 10 12:10:19 2011 -0600
+++ b/src/dm/impls/da/dacorn.c Thu Mar 10 13:13:17 2011 -0700
@@ -25,7 +25,7 @@
.keywords: distributed array, get, corners, nodes, local indices, coordinates
-.seealso: DMDAGetGhostCorners(), DMDAGetCoordinates(), DMDASetUniformCoordinates(). DMDAGetGhostedCoordinates(), DMDAGetCoordinateDA()
+.seealso: DMDASetGhostCoordinates(), DMDAGetGhostCorners(), DMDAGetCoordinates(), DMDASetUniformCoordinates(). DMDAGetGhostedCoordinates(), DMDAGetCoordinateDA()
@*/
PetscErrorCode DMDASetCoordinates(DM da,Vec c)
{
@@ -47,6 +47,46 @@
}
#undef __FUNCT__
+#define __FUNCT__ "DMDASetGhostedCoordinates"
+/*@
+ DMDASetGhostedCoordinates - Sets into the DMDA a vector that indicates the
+ coordinates of the local nodes, including ghost nodes.
+
+ Collective on DMDA
+
+ Input Parameter:
++ da - the distributed array
+- c - coordinate vector
+
+ Note:
+ The coordinates of interior ghost points can be set using DMDASetCoordinates()
+ followed by DMDAGetGhostedCoordinates(). This is intended to enable the setting
+ of ghost coordinates outside of the domain.
+
+ Non-ghosted coordinates, if set, are assumed still valid.
+
+ Level: intermediate
+
+.keywords: distributed array, get, corners, nodes, local indices, coordinates
+
+.seealso: DMDASetCoordinates(), DMDAGetGhostCorners(), DMDAGetCoordinates(), DMDASetUniformCoordinates(). DMDAGetGhostedCoordinates(), DMDAGetCoordinateDA()
+@*/
+PetscErrorCode DMDASetGhostedCoordinates(DM da,Vec c)
+{
+ PetscErrorCode ierr;
+ DM_DA *dd = (DM_DA*)da->data;
+
+ PetscFunctionBegin;
+ PetscValidHeaderSpecific(da,DM_CLASSID,1);
+ PetscValidHeaderSpecific(c,VEC_CLASSID,2);
+ ierr = PetscObjectReference((PetscObject)c);CHKERRQ(ierr);
+ if (dd->ghosted_coordinates) {ierr = VecDestroy(dd->ghosted_coordinates);CHKERRQ(ierr);}
+ dd->ghosted_coordinates = c;
+ ierr = VecSetBlockSize(c,dd->dim);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
#define __FUNCT__ "DMDAGetCoordinates"
/*@
DMDAGetCoordinates - Gets the node coordinates associated with a DMDA.
diff -r 1bde78661028 src/dm/impls/da/dainterp.c
--- a/src/dm/impls/da/dainterp.c Thu Mar 10 12:10:19 2011 -0600
+++ b/src/dm/impls/da/dainterp.c Thu Mar 10 13:13:17 2011 -0700
@@ -60,7 +60,7 @@
PetscFunctionBegin;
ierr = DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr);
ierr = DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr);
- if (pt == DMDA_XPERIODIC) {
+ if DMDAXPeriodic(pt) {
ratio = mx/Mx;
if (ratio*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
} else {
@@ -196,7 +196,7 @@
PetscFunctionBegin;
ierr = DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr);
ierr = DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr);
- if (pt == DMDA_XPERIODIC) {
+ if DMDAXPeriodic(pt) {
ratio = mx/Mx;
if (ratio*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
} else {
diff -r 1bde78661028 src/dm/impls/da/gr1.c
--- a/src/dm/impls/da/gr1.c Thu Mar 10 12:10:19 2011 -0600
+++ b/src/dm/impls/da/gr1.c Thu Mar 10 13:13:17 2011 -0700
@@ -44,8 +44,8 @@
ierr = DMDAGetCoordinateDA(da, &cda);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(cda, &xcoor);CHKERRQ(ierr);
if (dim == 1) {
- if (periodic == DMDA_NONPERIODIC) hx = (xmax-xmin)/(M-1);
- else hx = (xmax-xmin)/M;
+ if (DMDAXPeriodic(periodic)) hx = (xmax-xmin)/M;
+ else hx = (xmax-xmin)/(M-1);
ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr);
for (i=0; i<isize; i++) {
coors[i] = xmin + hx*(i+istart);
@@ -201,7 +201,7 @@
ierr = MPI_Send((void*)&array[j+(n-1)*step],1,MPIU_REAL,rank+1,tag1,comm);CHKERRQ(ierr);
ierr = MPI_Send((void*)&xg[n-1],1,MPIU_REAL,rank+1,tag1,comm);CHKERRQ(ierr);
}
- if (!rank && periodic && size > 1) { /* first processor sends first value to last */
+ if (!rank && DMDAXPeriodic(periodic) && size > 1) { /* first processor sends first value to last */
ierr = MPI_Send((void*)&array[j],1,MPIU_REAL,size-1,tag2,comm);CHKERRQ(ierr);
}
@@ -227,7 +227,7 @@
ierr = PetscDrawPoint(draw,xgtmp,tmp,PETSC_DRAW_BLACK);CHKERRQ(ierr);
}
}
- if (rank == size-1 && periodic && size > 1) {
+ if (rank == size-1 && DMDAXPeriodic(periodic) && size > 1) {
ierr = MPI_Recv(&tmp,1,MPIU_REAL,0,tag2,comm,&status);CHKERRQ(ierr);
#if !defined(PETSC_USE_COMPLEX)
ierr = PetscDrawLine(draw,xg[n-2],array[j+step*(n-1)],xg[n-1],tmp,PETSC_DRAW_RED);CHKERRQ(ierr);
diff -r 1bde78661028 src/dm/impls/da/gr2.c
--- a/src/dm/impls/da/gr2.c Thu Mar 10 12:10:19 2011 -0600
+++ b/src/dm/impls/da/gr2.c Thu Mar 10 13:13:17 2011 -0700
@@ -110,7 +110,7 @@
*/
ierr = PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);CHKERRQ(ierr);
if (!xlocal) {
- if (periodic != DMDA_NONPERIODIC || s != 1 || st != DMDA_STENCIL_BOX) {
+ if (!periodic || s != 1 || st != DMDA_STENCIL_BOX) {
/*
if original da is not of stencil width one, or periodic or not a box stencil then
create a special DMDA to handle one level of ghost points for graphics
@@ -135,7 +135,7 @@
ierr = PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);CHKERRQ(ierr);
ierr = PetscObjectDereference((PetscObject)xlocal);CHKERRQ(ierr);
} else {
- if (periodic == DMDA_NONPERIODIC && s == 1 && st == DMDA_STENCIL_BOX) {
+ if (!periodic && s == 1 && st == DMDA_STENCIL_BOX) {
dac = da;
} else {
ierr = PetscObjectQuery((PetscObject)xlocal,"DMDA",(PetscObject*)&dac);CHKERRQ(ierr);
More information about the petsc-dev
mailing list