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