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,maxnlevels = 100,*refx,*refy,*refz,n,i; 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 ierr = PetscOptionsInt("-da_overlap","Overlap between local grids","DMDASetOverlap",dd->overlap,&dd->overlap,PETSC_NULL);CHKERRQ(ierr); 37 /* Handle DMDA parallel distibution */ 38 ierr = PetscOptionsInt("-da_processors_x","Number of processors in x direction","DMDASetNumProcs",dd->m,&dd->m,PETSC_NULL);CHKERRQ(ierr); 39 if (dd->dim > 1) {ierr = PetscOptionsInt("-da_processors_y","Number of processors in y direction","DMDASetNumProcs",dd->n,&dd->n,PETSC_NULL);CHKERRQ(ierr);} 40 if (dd->dim > 2) {ierr = PetscOptionsInt("-da_processors_z","Number of processors in z direction","DMDASetNumProcs",dd->p,&dd->p,PETSC_NULL);CHKERRQ(ierr);} 41 /* Handle DMDA refinement */ 42 ierr = PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DMDASetRefinementFactor",dd->refine_x,&dd->refine_x,PETSC_NULL);CHKERRQ(ierr); 43 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);} 44 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);} 45 dd->coarsen_x = dd->refine_x; dd->coarsen_y = dd->refine_y; dd->coarsen_z = dd->refine_z; 46 47 /* Get refinement factors, defaults taken from the coarse DMDA */ 48 ierr = PetscMalloc3(maxnlevels,PetscInt,&refx,maxnlevels,PetscInt,&refy,maxnlevels,PetscInt,&refz);CHKERRQ(ierr); 49 ierr = DMDAGetRefinementFactor(da,&refx[0],&refy[0],&refz[0]);CHKERRQ(ierr); 50 for (i=1; i<maxnlevels; i++) { 51 refx[i] = refx[0]; 52 refy[i] = refy[0]; 53 refz[i] = refz[0]; 54 } 55 n = maxnlevels; 56 ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,PETSC_NULL);CHKERRQ(ierr); 57 if (da->levelup - da->leveldown >= 0) dd->refine_x = refx[da->levelup - da->leveldown]; 58 if (da->levelup - da->leveldown >= 1) dd->coarsen_x = refx[da->levelup - da->leveldown - 1]; 59 if (dd->dim > 1) { 60 n = maxnlevels; 61 ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,PETSC_NULL);CHKERRQ(ierr); 62 if (da->levelup - da->leveldown >= 0) dd->refine_y = refy[da->levelup - da->leveldown]; 63 if (da->levelup - da->leveldown >= 1) dd->coarsen_y = refy[da->levelup - da->leveldown - 1]; 64 } 65 if (dd->dim > 2) { 66 n = maxnlevels; 67 ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,PETSC_NULL);CHKERRQ(ierr); 68 if (da->levelup - da->leveldown >= 0) dd->refine_z = refz[da->levelup - da->leveldown]; 69 if (da->levelup - da->leveldown >= 1) dd->coarsen_z = refz[da->levelup - da->leveldown - 1]; 70 } 71 72 if (negativeMNP) {ierr = PetscOptionsInt("-da_refine","Uniformly refine DA one or more times","None",refine,&refine,PETSC_NULL);CHKERRQ(ierr);} 73 ierr = PetscOptionsTail();CHKERRQ(ierr); 74 75 while (refine--) { 76 if (dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){ 77 dd->M = dd->refine_x*dd->M; 78 } else { 79 dd->M = 1 + dd->refine_x*(dd->M - 1); 80 } 81 if (dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){ 82 dd->N = dd->refine_y*dd->N; 83 } else { 84 dd->N = 1 + dd->refine_y*(dd->N - 1); 85 } 86 if (dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){ 87 dd->P = dd->refine_z*dd->P; 88 } else { 89 dd->P = 1 + dd->refine_z*(dd->P - 1); 90 } 91 da->levelup++; 92 if (da->levelup - da->leveldown >= 0) { 93 dd->refine_x = refx[da->levelup - da->leveldown]; 94 dd->refine_y = refy[da->levelup - da->leveldown]; 95 dd->refine_z = refz[da->levelup - da->leveldown]; 96 } 97 if (da->levelup - da->leveldown >= 1) { 98 dd->coarsen_x = refx[da->levelup - da->leveldown - 1]; 99 dd->coarsen_y = refy[da->levelup - da->leveldown - 1]; 100 dd->coarsen_z = refz[da->levelup - da->leveldown - 1]; 101 } 102 } 103 ierr = PetscFree3(refx,refy,refz);CHKERRQ(ierr); 104 PetscFunctionReturn(0); 105 } 106 107 extern PetscErrorCode DMCreateGlobalVector_DA(DM,Vec*); 108 extern PetscErrorCode DMCreateLocalVector_DA(DM,Vec*); 109 extern PetscErrorCode DMGlobalToLocalBegin_DA(DM,Vec,InsertMode,Vec); 110 extern PetscErrorCode DMGlobalToLocalEnd_DA(DM,Vec,InsertMode,Vec); 111 extern PetscErrorCode DMLocalToGlobalBegin_DA(DM,Vec,InsertMode,Vec); 112 extern PetscErrorCode DMLocalToGlobalEnd_DA(DM,Vec,InsertMode,Vec); 113 extern PetscErrorCode DMCreateInterpolation_DA(DM,DM,Mat*,Vec*); 114 extern PetscErrorCode DMCreateColoring_DA(DM,ISColoringType,MatType,ISColoring*); 115 extern PetscErrorCode DMCreateMatrix_DA(DM,MatType,Mat*); 116 extern PetscErrorCode DMCreateCoordinateDM_DA(DM,DM*); 117 extern PetscErrorCode DMRefine_DA(DM,MPI_Comm,DM*); 118 extern PetscErrorCode DMCoarsen_DA(DM,MPI_Comm,DM*); 119 extern PetscErrorCode DMRefineHierarchy_DA(DM,PetscInt,DM[]); 120 extern PetscErrorCode DMCoarsenHierarchy_DA(DM,PetscInt,DM[]); 121 extern PetscErrorCode DMCreateInjection_DA(DM,DM,VecScatter*); 122 extern PetscErrorCode DMCreateAggregates_DA(DM,DM,Mat*); 123 extern PetscErrorCode DMView_DA(DM,PetscViewer); 124 extern PetscErrorCode DMSetUp_DA(DM); 125 extern PetscErrorCode DMDestroy_DA(DM); 126 127 #undef __FUNCT__ 128 #define __FUNCT__ "DMLoad_DA" 129 PetscErrorCode DMLoad_DA(DM da,PetscViewer viewer) 130 { 131 PetscErrorCode ierr; 132 PetscInt dim,m,n,p,dof,swidth; 133 DMDAStencilType stencil; 134 DMDABoundaryType bx,by,bz; 135 PetscInt classid = DM_FILE_CLASSID,subclassid = DMDA_FILE_CLASSID; 136 PetscBool coors; 137 DM dac; 138 Vec c; 139 140 PetscFunctionBegin; 141 ierr = PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr); 142 if (classid != DM_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not DM next in file"); 143 ierr = PetscViewerBinaryRead(viewer,&subclassid,1,PETSC_INT);CHKERRQ(ierr); 144 if (subclassid != DMDA_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not DM DA next in file"); 145 ierr = PetscViewerBinaryRead(viewer,&dim,1,PETSC_INT);CHKERRQ(ierr); 146 ierr = PetscViewerBinaryRead(viewer,&m,1,PETSC_INT);CHKERRQ(ierr); 147 ierr = PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);CHKERRQ(ierr); 148 ierr = PetscViewerBinaryRead(viewer,&p,1,PETSC_INT);CHKERRQ(ierr); 149 ierr = PetscViewerBinaryRead(viewer,&dof,1,PETSC_INT);CHKERRQ(ierr); 150 ierr = PetscViewerBinaryRead(viewer,&swidth,1,PETSC_INT);CHKERRQ(ierr); 151 ierr = PetscViewerBinaryRead(viewer,&bx,1,PETSC_ENUM);CHKERRQ(ierr); 152 ierr = PetscViewerBinaryRead(viewer,&by,1,PETSC_ENUM);CHKERRQ(ierr); 153 ierr = PetscViewerBinaryRead(viewer,&bz,1,PETSC_ENUM);CHKERRQ(ierr); 154 ierr = PetscViewerBinaryRead(viewer,&stencil,1,PETSC_ENUM);CHKERRQ(ierr); 155 156 ierr = DMDASetDim(da, dim);CHKERRQ(ierr); 157 ierr = DMDASetSizes(da, m,n,p);CHKERRQ(ierr); 158 ierr = DMDASetBoundaryType(da, bx, by, bz);CHKERRQ(ierr); 159 ierr = DMDASetDof(da, dof);CHKERRQ(ierr); 160 ierr = DMDASetStencilType(da, stencil);CHKERRQ(ierr); 161 ierr = DMDASetStencilWidth(da, swidth);CHKERRQ(ierr); 162 ierr = DMSetUp(da);CHKERRQ(ierr); 163 ierr = PetscViewerBinaryRead(viewer,&coors,1,PETSC_ENUM);CHKERRQ(ierr); 164 if (coors) { 165 ierr = DMGetCoordinateDM(da,&dac);CHKERRQ(ierr); 166 ierr = DMCreateGlobalVector(dac,&c);CHKERRQ(ierr); 167 ierr = VecLoad(c,viewer);CHKERRQ(ierr); 168 ierr = DMSetCoordinates(da,c);CHKERRQ(ierr); 169 ierr = VecDestroy(&c);CHKERRQ(ierr); 170 } 171 PetscFunctionReturn(0); 172 } 173 174 #undef __FUNCT__ 175 #define __FUNCT__ "DMCreateFieldDecomposition_DA" 176 PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len,char ***namelist, IS **islist, DM** dmlist) 177 { 178 PetscInt i; 179 PetscErrorCode ierr; 180 DM_DA *dd = (DM_DA*)dm->data; 181 PetscInt dof = dd->w; 182 183 PetscFunctionBegin; 184 if (len) *len = dof; 185 if (islist) { 186 Vec v; 187 PetscInt rstart,n; 188 189 ierr = DMGetGlobalVector(dm,&v);CHKERRQ(ierr); 190 ierr = VecGetOwnershipRange(v,&rstart,PETSC_NULL);CHKERRQ(ierr); 191 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 192 ierr = DMRestoreGlobalVector(dm,&v);CHKERRQ(ierr); 193 ierr = PetscMalloc(dof*sizeof(IS),islist);CHKERRQ(ierr); 194 for (i=0; i<dof; i++) { 195 ierr = ISCreateStride(((PetscObject)dm)->comm,n/dof,rstart+i,dof,&(*islist)[i]);CHKERRQ(ierr); 196 } 197 } 198 if (namelist) { 199 ierr = PetscMalloc(dof*sizeof(const char *), namelist);CHKERRQ(ierr); 200 if (dd->fieldname) { 201 for (i=0; i<dof; i++) { 202 ierr = PetscStrallocpy(dd->fieldname[i],&(*namelist)[i]);CHKERRQ(ierr); 203 } 204 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently DMDA must have fieldnames"); 205 } 206 if (dmlist) { 207 DM da; 208 209 ierr = DMDACreate(((PetscObject)dm)->comm, &da);CHKERRQ(ierr); 210 ierr = DMDASetDim(da, dd->dim);CHKERRQ(ierr); 211 ierr = DMDASetSizes(da, dd->M, dd->N, dd->P);CHKERRQ(ierr); 212 ierr = DMDASetNumProcs(da, dd->m, dd->n, dd->p);CHKERRQ(ierr); 213 ierr = DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz);CHKERRQ(ierr); 214 ierr = DMDASetDof(da, 1);CHKERRQ(ierr); 215 ierr = DMDASetStencilType(da, dd->stencil_type);CHKERRQ(ierr); 216 ierr = DMDASetStencilWidth(da, dd->s);CHKERRQ(ierr); 217 ierr = DMSetUp(da);CHKERRQ(ierr); 218 ierr = PetscMalloc(dof*sizeof(DM),dmlist);CHKERRQ(ierr); 219 for (i=0; i<dof-1; i++) {ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr);} 220 for (i=0; i<dof; i++) (*dmlist)[i] = da; 221 } 222 223 PetscFunctionReturn(0); 224 } 225 226 /*MC 227 DMDA = "da" - A DM object that is used to manage data for a structured grid in 1, 2, or 3 dimensions. 228 In the global representation of the vector each process stores a non-overlapping rectangular (or slab in 3d) portion of the grid points. 229 In the local representation these rectangular regions (slabs) are extended in all directions by a stencil width. 230 231 The vectors can be thought of as either cell centered or vertex centered on the mesh. But some variables cannot be cell centered and others 232 vertex centered. 233 234 235 Level: intermediate 236 237 .seealso: DMType, DMCOMPOSITE, DMDACreate(), DMCreate(), DMSetType() 238 M*/ 239 240 241 EXTERN_C_BEGIN 242 #undef __FUNCT__ 243 #define __FUNCT__ "DMCreate_DA" 244 PetscErrorCode DMCreate_DA(DM da) 245 { 246 PetscErrorCode ierr; 247 DM_DA *dd; 248 249 PetscFunctionBegin; 250 PetscValidPointer(da,1); 251 ierr = PetscNewLog(da,DM_DA,&dd);CHKERRQ(ierr); 252 da->data = dd; 253 254 dd->dim = -1; 255 dd->interptype = DMDA_Q1; 256 dd->refine_x = 2; 257 dd->refine_y = 2; 258 dd->refine_z = 2; 259 dd->coarsen_x = 2; 260 dd->coarsen_y = 2; 261 dd->coarsen_z = 2; 262 dd->fieldname = PETSC_NULL; 263 dd->nlocal = -1; 264 dd->Nlocal = -1; 265 dd->M = -1; 266 dd->N = -1; 267 dd->P = -1; 268 dd->m = -1; 269 dd->n = -1; 270 dd->p = -1; 271 dd->w = -1; 272 dd->s = -1; 273 dd->xs = -1; dd->xe = -1; dd->ys = -1; dd->ye = -1; dd->zs = -1; dd->ze = -1; 274 dd->Xs = -1; dd->Xe = -1; dd->Ys = -1; dd->Ye = -1; dd->Zs = -1; dd->Ze = -1; 275 276 dd->overlap = 0; 277 278 dd->gtol = PETSC_NULL; 279 dd->ltog = PETSC_NULL; 280 dd->ltol = PETSC_NULL; 281 dd->ao = PETSC_NULL; 282 dd->base = -1; 283 dd->bx = DMDA_BOUNDARY_NONE; 284 dd->by = DMDA_BOUNDARY_NONE; 285 dd->bz = DMDA_BOUNDARY_NONE; 286 dd->stencil_type = DMDA_STENCIL_BOX; 287 dd->interptype = DMDA_Q1; 288 dd->idx = PETSC_NULL; 289 dd->Nl = -1; 290 dd->lx = PETSC_NULL; 291 dd->ly = PETSC_NULL; 292 dd->lz = PETSC_NULL; 293 294 dd->elementtype = DMDA_ELEMENT_Q1; 295 296 ierr = PetscStrallocpy(VECSTANDARD,(char**)&da->vectype);CHKERRQ(ierr); 297 da->ops->globaltolocalbegin = DMGlobalToLocalBegin_DA; 298 da->ops->globaltolocalend = DMGlobalToLocalEnd_DA; 299 da->ops->localtoglobalbegin = DMLocalToGlobalBegin_DA; 300 da->ops->localtoglobalend = DMLocalToGlobalEnd_DA; 301 da->ops->createglobalvector = DMCreateGlobalVector_DA; 302 da->ops->createlocalvector = DMCreateLocalVector_DA; 303 da->ops->createinterpolation = DMCreateInterpolation_DA; 304 da->ops->getcoloring = DMCreateColoring_DA; 305 da->ops->creatematrix = DMCreateMatrix_DA; 306 da->ops->refine = DMRefine_DA; 307 da->ops->coarsen = DMCoarsen_DA; 308 da->ops->refinehierarchy = DMRefineHierarchy_DA; 309 da->ops->coarsenhierarchy = DMCoarsenHierarchy_DA; 310 da->ops->getinjection = DMCreateInjection_DA; 311 da->ops->getaggregates = DMCreateAggregates_DA; 312 da->ops->destroy = DMDestroy_DA; 313 da->ops->view = 0; 314 da->ops->setfromoptions = DMSetFromOptions_DA; 315 da->ops->setup = DMSetUp_DA; 316 da->ops->load = DMLoad_DA; 317 da->ops->createcoordinatedm = DMCreateCoordinateDM_DA; 318 da->ops->createfielddecomposition = DMCreateFieldDecomposition_DA; 319 PetscFunctionReturn(0); 320 } 321 EXTERN_C_END 322 323 #undef __FUNCT__ 324 #define __FUNCT__ "DMDACreate" 325 /*@ 326 DMDACreate - Creates a DMDA object. 327 328 Collective on MPI_Comm 329 330 Input Parameter: 331 . comm - The communicator for the DMDA object 332 333 Output Parameter: 334 . da - The DMDA object 335 336 Level: advanced 337 338 Developers Note: Since there exists DMDACreate1/2/3d() should this routine even exist? 339 340 .keywords: DMDA, create 341 .seealso: DMDASetSizes(), DMDADuplicate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 342 @*/ 343 PetscErrorCode DMDACreate(MPI_Comm comm, DM *da) 344 { 345 PetscErrorCode ierr; 346 347 PetscFunctionBegin; 348 PetscValidPointer(da,2); 349 ierr = DMCreate(comm,da);CHKERRQ(ierr); 350 ierr = DMSetType(*da,DMDA);CHKERRQ(ierr); 351 PetscFunctionReturn(0); 352 } 353 354 355