1 2 #define PETSCDM_DLL 3 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/ 4 #include "data_bucket.h" 5 6 //typedef PetscErrorCode (*swarm_project)(DM,DM,Vec) DMSwarmProjectMethod; /* swarm, geometry, result */ 7 8 //typedef enum { PROJECT_DMDA_AQ1=0, PROJECT_DMDA_P0 } DMSwarmDMDAProjectionType; 9 10 #if 0 11 12 /* Defines what the local space will be */ 13 PetscErrorCode DMSwarmSetOverlap(void) 14 { 15 16 PetscFunctionReturn(0); 17 } 18 19 20 /* coordinates */ 21 /* 22 DMGetCoordinateDM returns self 23 DMGetCoordinates and DMGetCoordinatesLocal return same thing 24 Local view could be used to define overlapping information 25 */ 26 27 #endif 28 29 #undef __FUNCT__ 30 #define __FUNCT__ "DMSwarmVectorDefineField" 31 PETSC_EXTERN PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[]) 32 { 33 DM_Swarm *swarm = (DM_Swarm*)dm->data; 34 PetscErrorCode ierr; 35 PetscInt bs,n; 36 PetscScalar *array; 37 PetscDataType type; 38 39 DataBucketGetSizes(swarm->db,&n,NULL,NULL); 40 ierr = DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);CHKERRQ(ierr); 41 42 /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */ 43 if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL"); 44 45 PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname); 46 swarm->vec_field_set = PETSC_TRUE; 47 swarm->vec_field_bs = 1;//bs; 48 swarm->vec_field_nlocal = n; 49 50 PetscFunctionReturn(0); 51 } 52 53 /* requires DMSwarmDefineFieldVector has been called */ 54 #undef __FUNCT__ 55 #define __FUNCT__ "DMCreateGlobalVector_Swarm" 56 PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec) 57 { 58 DM_Swarm *swarm = (DM_Swarm*)dm->data; 59 PetscErrorCode ierr; 60 Vec x; 61 char name[PETSC_MAX_PATH_LEN]; 62 63 if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first"); 64 PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name); 65 ierr = VecCreate(PetscObjectComm((PetscObject)dm),&x);CHKERRQ(ierr); 66 ierr = PetscObjectSetName((PetscObject)x,name);CHKERRQ(ierr); 67 ierr = VecSetSizes(x,swarm->vec_field_nlocal,PETSC_DETERMINE);CHKERRQ(ierr); 68 ierr = VecSetBlockSize(x,swarm->vec_field_bs);CHKERRQ(ierr); 69 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 70 *vec = x; 71 72 PetscFunctionReturn(0); 73 } 74 75 /* requires DMSwarmDefineFieldVector has been called */ 76 #undef __FUNCT__ 77 #define __FUNCT__ "DMCreateLocalVector_Swarm" 78 PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec) 79 { 80 DM_Swarm *swarm = (DM_Swarm*)dm->data; 81 PetscErrorCode ierr; 82 Vec x; 83 84 char name[PETSC_MAX_PATH_LEN]; 85 86 if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first"); 87 PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name); 88 ierr = VecCreate(PETSC_COMM_SELF,&x);CHKERRQ(ierr); 89 ierr = PetscObjectSetName((PetscObject)x,name);CHKERRQ(ierr); 90 ierr = VecSetSizes(x,swarm->vec_field_nlocal,swarm->vec_field_nlocal);CHKERRQ(ierr); 91 ierr = VecSetBlockSize(x,swarm->vec_field_bs);CHKERRQ(ierr); 92 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 93 *vec = x; 94 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "DMSwarmCreateGlobalVectorFromField" 100 PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec) 101 { 102 DM_Swarm *swarm = (DM_Swarm*)dm->data; 103 PetscErrorCode ierr; 104 PetscInt bs,n; 105 PetscScalar *array; 106 Vec x; 107 PetscDataType type; 108 char name[PETSC_MAX_PATH_LEN]; 109 110 DataBucketGetSizes(swarm->db,&n,NULL,NULL); 111 ierr = DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);CHKERRQ(ierr); 112 113 /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */ 114 if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL"); 115 116 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)dm),1,n,array,&x);CHKERRQ(ierr); 117 118 PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarm_VecFieldInPlace_%s",fieldname); 119 ierr = PetscObjectComposeFunction((PetscObject)x,name,DMSwarmDestroyGlobalVectorFromField);CHKERRQ(ierr); 120 121 /* Set guard */ 122 123 *vec = x; 124 PetscFunctionReturn(0); 125 } 126 127 #undef __FUNCT__ 128 #define __FUNCT__ "DMSwarmDestroyGlobalVectorFromField" 129 PETSC_EXTERN PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec) 130 { 131 DM_Swarm *swarm = (DM_Swarm*)dm->data; 132 PetscErrorCode ierr; 133 DataField gfield; 134 char name[PETSC_MAX_PATH_LEN]; 135 void (*fptr)(void); 136 137 /* get data field */ 138 DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield); 139 140 /* check vector is an inplace array */ 141 PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarm_VecFieldInPlace_%s",fieldname); 142 ierr = PetscObjectQueryFunction((PetscObject)(*vec),name,&fptr);CHKERRQ(ierr); 143 if (!fptr) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Vector being destroyed was not created from DMSwarm field(%s)",fieldname); 144 145 /* restore data field */ 146 DataFieldRestoreAccess(gfield); 147 148 ierr = VecDestroy(vec);CHKERRQ(ierr); 149 150 PetscFunctionReturn(0); 151 } 152 153 /* 154 PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm,const PetscInt nf,const char *fieldnames[],Vec *vec) 155 { 156 PetscFunctionReturn(0); 157 } 158 159 PETSC_EXTERN PetscErrorCode DMSwarmRestoreGlobalVectorFromFields(DM dm,Vec *vec) 160 { 161 PetscFunctionReturn(0); 162 } 163 */ 164 165 #undef __FUNCT__ 166 #define __FUNCT__ "DMSwarmInitializeFieldRegister" 167 PETSC_EXTERN PetscErrorCode DMSwarmInitializeFieldRegister(DM dm) 168 { 169 DM_Swarm *swarm = (DM_Swarm*)dm->data; 170 swarm->field_registration_initialized = PETSC_TRUE; 171 PetscFunctionReturn(0); 172 } 173 174 #undef __FUNCT__ 175 #define __FUNCT__ "DMSwarmFinalizeFieldRegister" 176 PETSC_EXTERN PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm) 177 { 178 DM_Swarm *swarm = (DM_Swarm*)dm->data; 179 swarm->field_registration_finalized = PETSC_TRUE; 180 DataBucketFinalize(swarm->db); 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "DMSwarmSetLocalSizes" 186 PETSC_EXTERN PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer) 187 { 188 DM_Swarm *swarm = (DM_Swarm*)dm->data; 189 190 DataBucketSetSizes(swarm->db,nlocal,buffer); 191 192 PetscFunctionReturn(0); 193 } 194 195 #undef __FUNCT__ 196 #define __FUNCT__ "DMSwarmRegisterPetscDatatypeField" 197 PETSC_EXTERN PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type) 198 { 199 DM_Swarm *swarm = (DM_Swarm*)dm->data; 200 size_t size; 201 202 if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first"); 203 if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first"); 204 205 if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 206 if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 207 if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 208 if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 209 if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 210 211 switch (type) { 212 case PETSC_CHAR: 213 size = sizeof(PetscChar); 214 break; 215 case PETSC_SHORT: 216 size = sizeof(PetscShort); 217 break; 218 case PETSC_INT: 219 size = sizeof(PetscInt); 220 break; 221 case PETSC_LONG: 222 size = sizeof(Petsc64bitInt); 223 break; 224 case PETSC_FLOAT: 225 size = sizeof(PetscFloat); 226 break; 227 case PETSC_DOUBLE: 228 size = sizeof(PetscReal); 229 break; 230 231 default: 232 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 233 break; 234 } 235 236 /* Load a specific data type into data bucket, specifying textual name and its size in bytes */ 237 DataBucketRegisterField(swarm->db,fieldname,size,NULL); 238 swarm->db->field[swarm->db->nfields-1]->petsc_type = type; 239 240 PetscFunctionReturn(0); 241 } 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "DMSwarmRegisterUserStructField" 245 PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size) 246 { 247 DM_Swarm *swarm = (DM_Swarm*)dm->data; 248 249 DataBucketRegisterField(swarm->db,fieldname,size,NULL); 250 swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ; 251 252 PetscFunctionReturn(0); 253 } 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "DMSwarmRegisterUserDatatypeField" 257 PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size) 258 { 259 DM_Swarm *swarm = (DM_Swarm*)dm->data; 260 261 DataBucketRegisterField(swarm->db,fieldname,size,NULL); 262 swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN; 263 264 PetscFunctionReturn(0); 265 } 266 267 #undef __FUNCT__ 268 #define __FUNCT__ "DMSwarmGetField" 269 PETSC_EXTERN PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data) 270 { 271 DM_Swarm *swarm = (DM_Swarm*)dm->data; 272 DataField gfield; 273 274 DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield); 275 DataFieldGetAccess(gfield); 276 DataFieldGetEntries(gfield,data); 277 if (blocksize) {} 278 if (type) { *type = gfield->petsc_type; } 279 280 PetscFunctionReturn(0); 281 } 282 283 #undef __FUNCT__ 284 #define __FUNCT__ "DMSwarmRestoreField" 285 PETSC_EXTERN PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data) 286 { 287 DM_Swarm *swarm = (DM_Swarm*)dm->data; 288 DataField gfield; 289 290 DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield); 291 DataFieldRestoreAccess(gfield); 292 if (data) *data = NULL; 293 294 PetscFunctionReturn(0); 295 } 296 297 298 #undef __FUNCT__ 299 #define __FUNCT__ "DMDestroy_Swarm" 300 PetscErrorCode DMDestroy_Swarm(DM dm) 301 { 302 DM_Swarm *swarm = (DM_Swarm*)dm->data; 303 PetscErrorCode ierr; 304 305 PetscFunctionBegin; 306 DataBucketDestroy(&swarm->db); 307 ierr = PetscFree(swarm);CHKERRQ(ierr); 308 PetscFunctionReturn(0); 309 } 310 311 #undef __FUNCT__ 312 #define __FUNCT__ "DMView_Swarm" 313 PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer) 314 { 315 DM_Swarm *swarm = (DM_Swarm*)dm->data; 316 PetscBool iascii,ibinary,ishdf5,isvtk; 317 PetscErrorCode ierr; 318 319 PetscFunctionBegin; 320 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 321 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 322 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 323 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);CHKERRQ(ierr); 324 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr); 325 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 326 if (iascii) { 327 DataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT); 328 } else if (ibinary) { 329 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support"); 330 } else if (ishdf5) { 331 #if defined(PETSC_HAVE_HDF5) 332 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO HDF5 support"); 333 #else 334 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5"); 335 #endif 336 } else if (isvtk) { 337 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support"); 338 } 339 PetscFunctionReturn(0); 340 } 341 342 #undef __FUNCT__ 343 #define __FUNCT__ "DMCreate_Swarm" 344 PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm) 345 { 346 DM_Swarm *swarm; 347 PetscErrorCode ierr; 348 349 PetscFunctionBegin; 350 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 351 ierr = PetscNewLog(dm,&swarm);CHKERRQ(ierr); 352 353 DataBucketCreate(&swarm->db); 354 swarm->vec_field_set = PETSC_FALSE; 355 356 dm->dim = 0; 357 dm->data = swarm; 358 359 dm->ops->view = DMView_Swarm; 360 dm->ops->load = NULL; 361 dm->ops->setfromoptions = NULL; 362 dm->ops->clone = NULL; 363 dm->ops->setup = NULL; 364 dm->ops->createdefaultsection = NULL; 365 dm->ops->createdefaultconstraints = NULL; 366 dm->ops->createglobalvector = DMCreateGlobalVector_Swarm; 367 dm->ops->createlocalvector = DMCreateLocalVector_Swarm; 368 dm->ops->getlocaltoglobalmapping = NULL; 369 dm->ops->createfieldis = NULL; 370 dm->ops->createcoordinatedm = NULL; 371 dm->ops->getcoloring = NULL; 372 dm->ops->creatematrix = NULL; 373 dm->ops->createinterpolation = NULL; 374 dm->ops->getaggregates = NULL; 375 dm->ops->getinjection = NULL; 376 dm->ops->refine = NULL; 377 dm->ops->coarsen = NULL; 378 dm->ops->refinehierarchy = NULL; 379 dm->ops->coarsenhierarchy = NULL; 380 dm->ops->globaltolocalbegin = NULL; 381 dm->ops->globaltolocalend = NULL; 382 dm->ops->localtoglobalbegin = NULL; 383 dm->ops->localtoglobalend = NULL; 384 dm->ops->destroy = DMDestroy_Swarm; 385 dm->ops->createsubdm = NULL; 386 dm->ops->getdimpoints = NULL; 387 dm->ops->locatepoints = NULL; 388 389 PetscFunctionReturn(0); 390 }