1 2 #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMSetFromOptions_DA" 6 PetscErrorCode DMSetFromOptions_DA(DM da) 7 { 8 PetscErrorCode ierr; 9 DM_DA *dd = (DM_DA*)da->data; 10 PetscInt refine = 0; 11 PetscBool negativeMNP = PETSC_FALSE,bM = PETSC_FALSE,bN = PETSC_FALSE, bP = PETSC_FALSE; 12 13 PetscFunctionBegin; 14 PetscValidHeaderSpecific(da,DM_CLASSID,1); 15 16 if (dd->M < 0) { 17 dd->M = -dd->M; 18 bM = PETSC_TRUE; 19 negativeMNP = PETSC_TRUE; 20 } 21 if (dd->dim > 1 && dd->N < 0) { 22 dd->N = -dd->N; 23 bN = PETSC_TRUE; 24 negativeMNP = PETSC_TRUE; 25 } 26 if (dd->dim > 2 && dd->P < 0) { 27 dd->P = -dd->P; 28 bP = PETSC_TRUE; 29 negativeMNP = PETSC_TRUE; 30 } 31 32 ierr = PetscOptionsHead("DMDA Options");CHKERRQ(ierr); 33 if (bM) {ierr = PetscOptionsInt("-da_grid_x","Number of grid points in x direction","DMDASetSizes",dd->M,&dd->M,PETSC_NULL);CHKERRQ(ierr);} 34 if (bN) {ierr = PetscOptionsInt("-da_grid_y","Number of grid points in y direction","DMDASetSizes",dd->N,&dd->N,PETSC_NULL);CHKERRQ(ierr);} 35 if (bP) {ierr = PetscOptionsInt("-da_grid_z","Number of grid points in z direction","DMDASetSizes",dd->P,&dd->P,PETSC_NULL);CHKERRQ(ierr);} 36 /* Handle DMDA parallel distibution */ 37 ierr = PetscOptionsInt("-da_processors_x","Number of processors in x direction","DMDASetNumProcs",dd->m,&dd->m,PETSC_NULL);CHKERRQ(ierr); 38 if (dd->dim > 1) {ierr = PetscOptionsInt("-da_processors_y","Number of processors in y direction","DMDASetNumProcs",dd->n,&dd->n,PETSC_NULL);CHKERRQ(ierr);} 39 if (dd->dim > 2) {ierr = PetscOptionsInt("-da_processors_z","Number of processors in z direction","DMDASetNumProcs",dd->p,&dd->p,PETSC_NULL);CHKERRQ(ierr);} 40 /* Handle DMDA refinement */ 41 ierr = PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DMDASetRefinementFactor",dd->refine_x,&dd->refine_x,PETSC_NULL);CHKERRQ(ierr); 42 if (dd->dim > 1) {ierr = PetscOptionsInt("-da_refine_y","Refinement ratio in y direction","DMDASetRefinementFactor",dd->refine_y,&dd->refine_y,PETSC_NULL);CHKERRQ(ierr);} 43 if (dd->dim > 2) {ierr = PetscOptionsInt("-da_refine_z","Refinement ratio in z direction","DMDASetRefinementFactor",dd->refine_z,&dd->refine_z,PETSC_NULL);CHKERRQ(ierr);} 44 45 if (negativeMNP) {ierr = PetscOptionsInt("-da_refine","Uniformly refine DA one or more times","None",refine,&refine,PETSC_NULL);CHKERRQ(ierr);} 46 ierr = PetscOptionsTail();CHKERRQ(ierr); 47 48 while (refine--) { 49 if (dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){ 50 dd->M = dd->refine_x*dd->M; 51 } else { 52 dd->M = 1 + dd->refine_x*(dd->M - 1); 53 } 54 if (dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){ 55 dd->N = dd->refine_y*dd->N; 56 } else { 57 dd->N = 1 + dd->refine_y*(dd->N - 1); 58 } 59 if (dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){ 60 dd->P = dd->refine_z*dd->P; 61 } else { 62 dd->P = 1 + dd->refine_z*(dd->P - 1); 63 } 64 da->levelup++; 65 } 66 PetscFunctionReturn(0); 67 } 68 69 extern PetscErrorCode DMCreateGlobalVector_DA(DM,Vec*); 70 extern PetscErrorCode DMCreateLocalVector_DA(DM,Vec*); 71 extern PetscErrorCode DMGlobalToLocalBegin_DA(DM,Vec,InsertMode,Vec); 72 extern PetscErrorCode DMGlobalToLocalEnd_DA(DM,Vec,InsertMode,Vec); 73 extern PetscErrorCode DMLocalToGlobalBegin_DA(DM,Vec,InsertMode,Vec); 74 extern PetscErrorCode DMLocalToGlobalEnd_DA(DM,Vec,InsertMode,Vec); 75 extern PetscErrorCode DMCreateInterpolation_DA(DM,DM,Mat*,Vec*); 76 extern PetscErrorCode DMCreateColoring_DA(DM,ISColoringType,const MatType,ISColoring*); 77 extern PetscErrorCode DMCreateMatrix_DA(DM,const MatType,Mat*); 78 extern PetscErrorCode DMRefine_DA(DM,MPI_Comm,DM*); 79 extern PetscErrorCode DMCoarsen_DA(DM,MPI_Comm,DM*); 80 extern PetscErrorCode DMRefineHierarchy_DA(DM,PetscInt,DM[]); 81 extern PetscErrorCode DMCoarsenHierarchy_DA(DM,PetscInt,DM[]); 82 extern PetscErrorCode DMCreateInjection_DA(DM,DM,VecScatter*); 83 extern PetscErrorCode DMCreateAggregates_DA(DM,DM,Mat*); 84 extern PetscErrorCode DMView_DA(DM,PetscViewer); 85 extern PetscErrorCode DMSetUp_DA(DM); 86 extern PetscErrorCode DMDestroy_DA(DM); 87 88 #undef __FUNCT__ 89 #define __FUNCT__ "DMLoad_DA" 90 PetscErrorCode DMLoad_DA(DM da,PetscViewer viewer) 91 { 92 PetscErrorCode ierr; 93 PetscInt dim,m,n,p,dof,swidth; 94 DMDAStencilType stencil; 95 DMDABoundaryType bx,by,bz; 96 PetscInt classid = DM_FILE_CLASSID,subclassid = DMDA_FILE_CLASSID; 97 PetscBool coors; 98 DM dac; 99 Vec c; 100 101 PetscFunctionBegin; 102 ierr = PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr); 103 if (classid != DM_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not DM next in file"); 104 ierr = PetscViewerBinaryRead(viewer,&subclassid,1,PETSC_INT);CHKERRQ(ierr); 105 if (subclassid != DMDA_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not DM DA next in file"); 106 ierr = PetscViewerBinaryRead(viewer,&dim,1,PETSC_INT);CHKERRQ(ierr); 107 ierr = PetscViewerBinaryRead(viewer,&m,1,PETSC_INT);CHKERRQ(ierr); 108 ierr = PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);CHKERRQ(ierr); 109 ierr = PetscViewerBinaryRead(viewer,&p,1,PETSC_INT);CHKERRQ(ierr); 110 ierr = PetscViewerBinaryRead(viewer,&dof,1,PETSC_INT);CHKERRQ(ierr); 111 ierr = PetscViewerBinaryRead(viewer,&swidth,1,PETSC_INT);CHKERRQ(ierr); 112 ierr = PetscViewerBinaryRead(viewer,&bx,1,PETSC_ENUM);CHKERRQ(ierr); 113 ierr = PetscViewerBinaryRead(viewer,&by,1,PETSC_ENUM);CHKERRQ(ierr); 114 ierr = PetscViewerBinaryRead(viewer,&bz,1,PETSC_ENUM);CHKERRQ(ierr); 115 ierr = PetscViewerBinaryRead(viewer,&stencil,1,PETSC_ENUM);CHKERRQ(ierr); 116 117 ierr = DMDASetDim(da, dim);CHKERRQ(ierr); 118 ierr = DMDASetSizes(da, m,n,p);CHKERRQ(ierr); 119 ierr = DMDASetBoundaryType(da, bx, by, bz);CHKERRQ(ierr); 120 ierr = DMDASetDof(da, dof);CHKERRQ(ierr); 121 ierr = DMDASetStencilType(da, stencil);CHKERRQ(ierr); 122 ierr = DMDASetStencilWidth(da, swidth);CHKERRQ(ierr); 123 ierr = DMSetUp(da);CHKERRQ(ierr); 124 ierr = PetscViewerBinaryRead(viewer,&coors,1,PETSC_ENUM);CHKERRQ(ierr); 125 if (coors) { 126 ierr = DMDAGetCoordinateDA(da,&dac);CHKERRQ(ierr); 127 ierr = DMCreateGlobalVector(dac,&c);CHKERRQ(ierr); 128 ierr = VecLoad(c,viewer);CHKERRQ(ierr); 129 ierr = DMDASetCoordinates(da,c);CHKERRQ(ierr); 130 ierr = VecDestroy(&c);CHKERRQ(ierr); 131 } 132 PetscFunctionReturn(0); 133 } 134 135 #undef __FUNCT__ 136 #define __FUNCT__ "DMCreateDecomposition_DA" 137 PetscErrorCode DMCreateDecomposition_DA(DM dm, PetscInt *len,char ***namelist, IS **islist, DM** dmlist) 138 { 139 PetscInt i; 140 PetscErrorCode ierr; 141 DM_DA *dd = (DM_DA*)dm->data; 142 PetscInt dof = dd->w; 143 144 PetscFunctionBegin; 145 if (islist) { 146 Vec v; 147 PetscInt rstart,n; 148 149 ierr = DMGetGlobalVector(dm,&v);CHKERRQ(ierr); 150 ierr = VecGetOwnershipRange(v,&rstart,PETSC_NULL);CHKERRQ(ierr); 151 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 152 ierr = DMRestoreGlobalVector(dm,&v);CHKERRQ(ierr); 153 ierr = PetscMalloc(dof*sizeof(IS),islist);CHKERRQ(ierr); 154 for (i=0; i<dof; i++) { 155 ierr = ISCreateStride(((PetscObject)dm)->comm,n/dof,rstart+i,dof,&(*islist)[i]);CHKERRQ(ierr); 156 } 157 } 158 if (namelist) { 159 ierr = PetscMalloc(dof*sizeof(const char *), namelist);CHKERRQ(ierr); 160 if (dd->fieldname) { 161 for (i=0; i<dof; i++) { 162 ierr = PetscStrallocpy(dd->fieldname[i],&(*namelist)[i]);CHKERRQ(ierr); 163 } 164 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently DMDA must have fieldnames"); 165 } 166 if (dmlist) { 167 DM da; 168 169 ierr = DMDACreate(((PetscObject)dm)->comm, &da);CHKERRQ(ierr); 170 ierr = DMDASetDim(da, dd->dim);CHKERRQ(ierr); 171 ierr = DMDASetSizes(da, dd->M, dd->N, dd->P);CHKERRQ(ierr); 172 ierr = DMDASetNumProcs(da, dd->m, dd->n, dd->p);CHKERRQ(ierr); 173 ierr = DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz);CHKERRQ(ierr); 174 ierr = DMDASetDof(da, 1);CHKERRQ(ierr); 175 ierr = DMDASetStencilType(da, dd->stencil_type);CHKERRQ(ierr); 176 ierr = DMDASetStencilWidth(da, dd->s);CHKERRQ(ierr); 177 ierr = DMSetUp(da);CHKERRQ(ierr); 178 ierr = PetscMalloc(dof*sizeof(DM),dmlist);CHKERRQ(ierr); 179 for (i=0; i<dof-1; i++) {ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr);} 180 for (i=0; i<dof; i++) (*dmlist)[i] = da; 181 } 182 *len = dof; 183 184 PetscFunctionReturn(0); 185 } 186 187 188 EXTERN_C_BEGIN 189 #undef __FUNCT__ 190 #define __FUNCT__ "DMCreate_DA" 191 PetscErrorCode DMCreate_DA(DM da) 192 { 193 PetscErrorCode ierr; 194 DM_DA *dd; 195 196 PetscFunctionBegin; 197 PetscValidPointer(da,1); 198 ierr = PetscNewLog(da,DM_DA,&dd);CHKERRQ(ierr); 199 da->data = dd; 200 201 dd->dim = -1; 202 dd->interptype = DMDA_Q1; 203 dd->refine_x = 2; 204 dd->refine_y = 2; 205 dd->refine_z = 2; 206 dd->fieldname = PETSC_NULL; 207 dd->nlocal = -1; 208 dd->Nlocal = -1; 209 dd->M = -1; 210 dd->N = -1; 211 dd->P = -1; 212 dd->m = -1; 213 dd->n = -1; 214 dd->p = -1; 215 dd->w = -1; 216 dd->s = -1; 217 dd->xs = -1; dd->xe = -1; dd->ys = -1; dd->ye = -1; dd->zs = -1; dd->ze = -1; 218 dd->Xs = -1; dd->Xe = -1; dd->Ys = -1; dd->Ye = -1; dd->Zs = -1; dd->Ze = -1; 219 220 dd->gtol = PETSC_NULL; 221 dd->ltog = PETSC_NULL; 222 dd->ltol = PETSC_NULL; 223 dd->ao = PETSC_NULL; 224 dd->base = -1; 225 dd->bx = DMDA_BOUNDARY_NONE; 226 dd->by = DMDA_BOUNDARY_NONE; 227 dd->bz = DMDA_BOUNDARY_NONE; 228 dd->stencil_type = DMDA_STENCIL_BOX; 229 dd->interptype = DMDA_Q1; 230 dd->idx = PETSC_NULL; 231 dd->Nl = -1; 232 dd->lx = PETSC_NULL; 233 dd->ly = PETSC_NULL; 234 dd->lz = PETSC_NULL; 235 236 dd->elementtype = DMDA_ELEMENT_Q1; 237 238 ierr = PetscStrallocpy(VECSTANDARD,&da->vectype);CHKERRQ(ierr); 239 da->ops->globaltolocalbegin = DMGlobalToLocalBegin_DA; 240 da->ops->globaltolocalend = DMGlobalToLocalEnd_DA; 241 da->ops->localtoglobalbegin = DMLocalToGlobalBegin_DA; 242 da->ops->localtoglobalend = DMLocalToGlobalEnd_DA; 243 da->ops->createglobalvector = DMCreateGlobalVector_DA; 244 da->ops->createlocalvector = DMCreateLocalVector_DA; 245 da->ops->createinterpolation = DMCreateInterpolation_DA; 246 da->ops->getcoloring = DMCreateColoring_DA; 247 da->ops->creatematrix = DMCreateMatrix_DA; 248 da->ops->refine = DMRefine_DA; 249 da->ops->coarsen = DMCoarsen_DA; 250 da->ops->refinehierarchy = DMRefineHierarchy_DA; 251 da->ops->coarsenhierarchy = DMCoarsenHierarchy_DA; 252 da->ops->getinjection = DMCreateInjection_DA; 253 da->ops->getaggregates = DMCreateAggregates_DA; 254 da->ops->destroy = DMDestroy_DA; 255 da->ops->view = 0; 256 da->ops->setfromoptions = DMSetFromOptions_DA; 257 da->ops->setup = DMSetUp_DA; 258 da->ops->load = DMLoad_DA; 259 da->ops->createdecomposition = DMCreateDecomposition_DA; 260 PetscFunctionReturn(0); 261 } 262 EXTERN_C_END 263 264 #undef __FUNCT__ 265 #define __FUNCT__ "DMDACreate" 266 /*@ 267 DMDACreate - Creates a DMDA object. 268 269 Collective on MPI_Comm 270 271 Input Parameter: 272 . comm - The communicator for the DMDA object 273 274 Output Parameter: 275 . da - The DMDA object 276 277 Level: advanced 278 279 Developers Note: Since there exists DMDACreate1/2/3d() should this routine even exist? 280 281 .keywords: DMDA, create 282 .seealso: DMDASetSizes(), DMDADuplicate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 283 @*/ 284 PetscErrorCode DMDACreate(MPI_Comm comm, DM *da) 285 { 286 PetscErrorCode ierr; 287 288 PetscFunctionBegin; 289 PetscValidPointer(da,2); 290 ierr = DMCreate(comm,da);CHKERRQ(ierr); 291 ierr = DMSetType(*da,DMDA);CHKERRQ(ierr); 292 PetscFunctionReturn(0); 293 } 294 295 296