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 306 #undef __FUNCT__ 307 #define __FUNCT__ "DMDestroy_Swarm" 308 PetscErrorCode DMDestroy_Swarm(DM dm) 309 { 310 DM_Swarm *swarm = (DM_Swarm*)dm->data; 311 PetscErrorCode ierr; 312 313 PetscFunctionBegin; 314 ierr = DataBucketDestroy(&swarm->db);CHKERRQ(ierr); 315 ierr = PetscFree(swarm);CHKERRQ(ierr); 316 PetscFunctionReturn(0); 317 } 318 319 #undef __FUNCT__ 320 #define __FUNCT__ "DMView_Swarm" 321 PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer) 322 { 323 DM_Swarm *swarm = (DM_Swarm*)dm->data; 324 PetscBool iascii,ibinary,ishdf5,isvtk; 325 PetscErrorCode ierr; 326 327 PetscFunctionBegin; 328 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 329 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 330 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 331 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);CHKERRQ(ierr); 332 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr); 333 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 334 if (iascii) { 335 ierr = DataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);CHKERRQ(ierr); 336 } else if (ibinary) { 337 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support"); 338 } else if (ishdf5) { 339 #if defined(PETSC_HAVE_HDF5) 340 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO HDF5 support"); 341 #else 342 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5"); 343 #endif 344 } else if (isvtk) { 345 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support"); 346 } 347 PetscFunctionReturn(0); 348 } 349 350 #undef __FUNCT__ 351 #define __FUNCT__ "DMCreate_Swarm" 352 PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm) 353 { 354 DM_Swarm *swarm; 355 PetscErrorCode ierr; 356 357 PetscFunctionBegin; 358 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 359 ierr = PetscNewLog(dm,&swarm);CHKERRQ(ierr); 360 361 ierr = DataBucketCreate(&swarm->db);CHKERRQ(ierr); 362 swarm->vec_field_set = PETSC_FALSE; 363 364 dm->dim = 0; 365 dm->data = swarm; 366 367 dm->ops->view = DMView_Swarm; 368 dm->ops->load = NULL; 369 dm->ops->setfromoptions = NULL; 370 dm->ops->clone = NULL; 371 dm->ops->setup = NULL; 372 dm->ops->createdefaultsection = NULL; 373 dm->ops->createdefaultconstraints = NULL; 374 dm->ops->createglobalvector = DMCreateGlobalVector_Swarm; 375 dm->ops->createlocalvector = DMCreateLocalVector_Swarm; 376 dm->ops->getlocaltoglobalmapping = NULL; 377 dm->ops->createfieldis = NULL; 378 dm->ops->createcoordinatedm = NULL; 379 dm->ops->getcoloring = NULL; 380 dm->ops->creatematrix = NULL; 381 dm->ops->createinterpolation = NULL; 382 dm->ops->getaggregates = NULL; 383 dm->ops->getinjection = NULL; 384 dm->ops->refine = NULL; 385 dm->ops->coarsen = NULL; 386 dm->ops->refinehierarchy = NULL; 387 dm->ops->coarsenhierarchy = NULL; 388 dm->ops->globaltolocalbegin = NULL; 389 dm->ops->globaltolocalend = NULL; 390 dm->ops->localtoglobalbegin = NULL; 391 dm->ops->localtoglobalend = NULL; 392 dm->ops->destroy = DMDestroy_Swarm; 393 dm->ops->createsubdm = NULL; 394 dm->ops->getdimpoints = NULL; 395 dm->ops->locatepoints = NULL; 396 397 PetscFunctionReturn(0); 398 }