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 ierr = DataBucketGetSizes(swarm->db,&n,NULL,NULL);CHKERRQ(ierr); 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 ierr = DataBucketGetSizes(swarm->db,&n,NULL,NULL);CHKERRQ(ierr); 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 ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr); 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 ierr = DataFieldRestoreAccess(gfield);CHKERRQ(ierr); 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 PetscErrorCode ierr; 180 181 swarm->field_registration_finalized = PETSC_TRUE; 182 ierr = DataBucketFinalize(swarm->db);CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 186 #undef __FUNCT__ 187 #define __FUNCT__ "DMSwarmSetLocalSizes" 188 PETSC_EXTERN PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer) 189 { 190 DM_Swarm *swarm = (DM_Swarm*)dm->data; 191 PetscErrorCode ierr; 192 193 ierr = DataBucketSetSizes(swarm->db,nlocal,buffer);CHKERRQ(ierr); 194 195 PetscFunctionReturn(0); 196 } 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "DMSwarmRegisterPetscDatatypeField" 200 PETSC_EXTERN PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type) 201 { 202 PetscErrorCode ierr; 203 DM_Swarm *swarm = (DM_Swarm*)dm->data; 204 size_t size; 205 206 if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first"); 207 if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first"); 208 209 if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 210 if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 211 if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 212 if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 213 if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 214 215 switch (type) { 216 case PETSC_CHAR: 217 size = sizeof(PetscChar); 218 break; 219 case PETSC_SHORT: 220 size = sizeof(PetscShort); 221 break; 222 case PETSC_INT: 223 size = sizeof(PetscInt); 224 break; 225 case PETSC_LONG: 226 size = sizeof(Petsc64bitInt); 227 break; 228 case PETSC_FLOAT: 229 size = sizeof(PetscFloat); 230 break; 231 case PETSC_DOUBLE: 232 size = sizeof(PetscReal); 233 break; 234 235 default: 236 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}"); 237 break; 238 } 239 240 /* Load a specific data type into data bucket, specifying textual name and its size in bytes */ 241 ierr = DataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,size,NULL);CHKERRQ(ierr); 242 swarm->db->field[swarm->db->nfields-1]->petsc_type = type; 243 244 PetscFunctionReturn(0); 245 } 246 247 #undef __FUNCT__ 248 #define __FUNCT__ "DMSwarmRegisterUserStructField" 249 PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size) 250 { 251 PetscErrorCode ierr; 252 DM_Swarm *swarm = (DM_Swarm*)dm->data; 253 254 ierr = DataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);CHKERRQ(ierr); 255 swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ; 256 257 PetscFunctionReturn(0); 258 } 259 260 #undef __FUNCT__ 261 #define __FUNCT__ "DMSwarmRegisterUserDatatypeField" 262 PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size) 263 { 264 DM_Swarm *swarm = (DM_Swarm*)dm->data; 265 PetscErrorCode ierr; 266 267 ierr = DataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,size,NULL);CHKERRQ(ierr); 268 swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN; 269 270 PetscFunctionReturn(0); 271 } 272 273 #undef __FUNCT__ 274 #define __FUNCT__ "DMSwarmGetField" 275 PETSC_EXTERN PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data) 276 { 277 DM_Swarm *swarm = (DM_Swarm*)dm->data; 278 DataField gfield; 279 PetscErrorCode ierr; 280 281 ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr); 282 ierr = DataFieldGetAccess(gfield);CHKERRQ(ierr); 283 ierr = DataFieldGetEntries(gfield,data);CHKERRQ(ierr); 284 if (blocksize) {} 285 if (type) { *type = gfield->petsc_type; } 286 287 PetscFunctionReturn(0); 288 } 289 290 #undef __FUNCT__ 291 #define __FUNCT__ "DMSwarmRestoreField" 292 PETSC_EXTERN PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data) 293 { 294 DM_Swarm *swarm = (DM_Swarm*)dm->data; 295 DataField gfield; 296 PetscErrorCode ierr; 297 298 ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr); 299 ierr = DataFieldRestoreAccess(gfield);CHKERRQ(ierr); 300 if (data) *data = NULL; 301 302 PetscFunctionReturn(0); 303 } 304 305 #undef __FUNCT__ 306 #define __FUNCT__ "DMSwarmAddPoint" 307 PETSC_EXTERN PetscErrorCode DMSwarmAddPoint(DM dm) 308 { 309 DM_Swarm *swarm = (DM_Swarm*)dm->data; 310 PetscErrorCode ierr; 311 312 ierr = DataBucketAddPoint(swarm->db);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 #undef __FUNCT__ 317 #define __FUNCT__ "DMSwarmAddNPoints" 318 PETSC_EXTERN PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints) 319 { 320 DM_Swarm *swarm = (DM_Swarm*)dm->data; 321 PetscErrorCode ierr; 322 PetscInt nlocal; 323 324 ierr = DataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);CHKERRQ(ierr); 325 nlocal = nlocal + npoints; 326 ierr = DataBucketSetSizes(swarm->db,nlocal,-1);CHKERRQ(ierr); 327 PetscFunctionReturn(0); 328 } 329 330 #undef __FUNCT__ 331 #define __FUNCT__ "DMSwarmRemovePoint" 332 PETSC_EXTERN PetscErrorCode DMSwarmRemovePoint(DM dm) 333 { 334 DM_Swarm *swarm = (DM_Swarm*)dm->data; 335 PetscErrorCode ierr; 336 337 ierr = DataBucketRemovePoint(swarm->db);CHKERRQ(ierr); 338 PetscFunctionReturn(0); 339 } 340 341 #undef __FUNCT__ 342 #define __FUNCT__ "DMSwarmRemovePointAtIndex" 343 PETSC_EXTERN PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx) 344 { 345 DM_Swarm *swarm = (DM_Swarm*)dm->data; 346 PetscErrorCode ierr; 347 348 ierr = DataBucketRemovePointAtIndex(swarm->db,idx);CHKERRQ(ierr); 349 PetscFunctionReturn(0); 350 } 351 352 #undef __FUNCT__ 353 #define __FUNCT__ "DMDestroy_Swarm" 354 PetscErrorCode DMDestroy_Swarm(DM dm) 355 { 356 DM_Swarm *swarm = (DM_Swarm*)dm->data; 357 PetscErrorCode ierr; 358 359 PetscFunctionBegin; 360 ierr = DataBucketDestroy(&swarm->db);CHKERRQ(ierr); 361 ierr = PetscFree(swarm);CHKERRQ(ierr); 362 PetscFunctionReturn(0); 363 } 364 365 #undef __FUNCT__ 366 #define __FUNCT__ "DMView_Swarm" 367 PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer) 368 { 369 DM_Swarm *swarm = (DM_Swarm*)dm->data; 370 PetscBool iascii,ibinary,ishdf5,isvtk; 371 PetscErrorCode ierr; 372 373 PetscFunctionBegin; 374 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 375 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 376 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 377 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);CHKERRQ(ierr); 378 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr); 379 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 380 if (iascii) { 381 ierr = DataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);CHKERRQ(ierr); 382 } else if (ibinary) { 383 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support"); 384 } else if (ishdf5) { 385 #if defined(PETSC_HAVE_HDF5) 386 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO HDF5 support"); 387 #else 388 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5"); 389 #endif 390 } else if (isvtk) { 391 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support"); 392 } 393 PetscFunctionReturn(0); 394 } 395 396 #undef __FUNCT__ 397 #define __FUNCT__ "DMCreate_Swarm" 398 PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm) 399 { 400 DM_Swarm *swarm; 401 PetscErrorCode ierr; 402 403 PetscFunctionBegin; 404 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 405 ierr = PetscNewLog(dm,&swarm);CHKERRQ(ierr); 406 407 ierr = DataBucketCreate(&swarm->db);CHKERRQ(ierr); 408 swarm->vec_field_set = PETSC_FALSE; 409 410 dm->dim = 0; 411 dm->data = swarm; 412 413 dm->ops->view = DMView_Swarm; 414 dm->ops->load = NULL; 415 dm->ops->setfromoptions = NULL; 416 dm->ops->clone = NULL; 417 dm->ops->setup = NULL; 418 dm->ops->createdefaultsection = NULL; 419 dm->ops->createdefaultconstraints = NULL; 420 dm->ops->createglobalvector = DMCreateGlobalVector_Swarm; 421 dm->ops->createlocalvector = DMCreateLocalVector_Swarm; 422 dm->ops->getlocaltoglobalmapping = NULL; 423 dm->ops->createfieldis = NULL; 424 dm->ops->createcoordinatedm = NULL; 425 dm->ops->getcoloring = NULL; 426 dm->ops->creatematrix = NULL; 427 dm->ops->createinterpolation = NULL; 428 dm->ops->getaggregates = NULL; 429 dm->ops->getinjection = NULL; 430 dm->ops->refine = NULL; 431 dm->ops->coarsen = NULL; 432 dm->ops->refinehierarchy = NULL; 433 dm->ops->coarsenhierarchy = NULL; 434 dm->ops->globaltolocalbegin = NULL; 435 dm->ops->globaltolocalend = NULL; 436 dm->ops->localtoglobalbegin = NULL; 437 dm->ops->localtoglobalend = NULL; 438 dm->ops->destroy = DMDestroy_Swarm; 439 dm->ops->createsubdm = NULL; 440 dm->ops->getdimpoints = NULL; 441 dm->ops->locatepoints = NULL; 442 443 PetscFunctionReturn(0); 444 }