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