1 #include <petsc-private/dmmbimpl.h> /*I "petscdm.h" I*/ 2 3 #include <petscdmmoab.h> 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "DMMoabSetFieldVector" 7 PetscErrorCode DMMoabSetFieldVector(DM dm, PetscInt ifield, Vec fvec) 8 { 9 DM_Moab *dmmoab; 10 moab::Tag vtag,ntag; 11 const PetscScalar *varray; 12 PetscScalar *farray; 13 moab::ErrorCode merr; 14 PetscErrorCode ierr; 15 std::string tag_name; 16 17 PetscFunctionBegin; 18 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 19 dmmoab = (DM_Moab*)(dm)->data; 20 21 if ((ifield < 0) || (ifield >= dmmoab->numFields)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The field %d should be positive and less than %d.", ifield, dmmoab->numFields); 22 23 /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */ 24 merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield],1,moab::MB_TYPE_DOUBLE,ntag, 25 moab::MB_TAG_DENSE|moab::MB_TAG_CREAT);MBERRNM(merr); 26 27 ierr = DMMoabGetVecTag(fvec,&vtag);CHKERRQ(ierr); 28 29 merr = dmmoab->mbiface->tag_get_name(vtag, tag_name); 30 if (!tag_name.length() && merr !=moab::MB_SUCCESS) { 31 ierr = VecGetArrayRead(fvec,&varray);CHKERRQ(ierr); 32 /* use the entity handle and the Dof index to set the right value */ 33 merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)varray);MBERRNM(merr); 34 ierr = VecRestoreArrayRead(fvec,&varray);CHKERRQ(ierr); 35 } 36 else { 37 ierr = PetscMalloc(dmmoab->nloc*sizeof(PetscScalar),&farray);CHKERRQ(ierr); 38 /* we are using a MOAB Vec - directly copy the tag data to new one */ 39 merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)farray);MBERRNM(merr); 40 merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray);MBERRNM(merr); 41 /* make sure the parallel exchange for ghosts are done appropriately */ 42 ierr = PetscFree(farray);CHKERRQ(ierr); 43 } 44 merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vowned);MBERRNM(merr); 45 PetscFunctionReturn(0); 46 } 47 48 49 #undef __FUNCT__ 50 #define __FUNCT__ "DMMoabSetGlobalFieldVector" 51 PetscErrorCode DMMoabSetGlobalFieldVector(DM dm, Vec fvec) 52 { 53 DM_Moab *dmmoab; 54 moab::Tag vtag,ntag; 55 const PetscScalar *varray; 56 PetscScalar *farray; 57 moab::ErrorCode merr; 58 PetscErrorCode ierr; 59 PetscInt i,ifield; 60 std::string tag_name; 61 moab::Range::iterator iter; 62 63 PetscFunctionBegin; 64 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 65 dmmoab = (DM_Moab*)(dm)->data; 66 67 /* get the Tag corresponding to the global vector - possible that there is no tag associated.. */ 68 ierr = DMMoabGetVecTag(fvec,&vtag);CHKERRQ(ierr); 69 merr = dmmoab->mbiface->tag_get_name(vtag, tag_name); 70 ierr = PetscMalloc(dmmoab->nloc*sizeof(PetscScalar),&farray);CHKERRQ(ierr); 71 if (!tag_name.length() && merr !=moab::MB_SUCCESS) { 72 /* not a MOAB vector - use VecGetSubVector to get the parts as needed */ 73 74 ierr = VecGetArrayRead(fvec,&varray);CHKERRQ(ierr); 75 for (ifield=0; ifield<dmmoab->numFields; ++ifield) { 76 77 /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */ 78 merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield],1,moab::MB_TYPE_DOUBLE,ntag, 79 moab::MB_TAG_DENSE|moab::MB_TAG_CREAT);MBERRNM(merr); 80 81 for(i=0;i<dmmoab->nloc;i++) { 82 if (dmmoab->bs == 1) 83 farray[i]=varray[ifield*dmmoab->nloc+i]; 84 else 85 farray[i]=varray[i*dmmoab->numFields+ifield]; 86 } 87 88 /* use the entity handle and the Dof index to set the right value */ 89 merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray);MBERRNM(merr); 90 } 91 ierr = VecRestoreArrayRead(fvec,&varray);CHKERRQ(ierr); 92 } 93 else { 94 ierr = PetscMalloc(dmmoab->nloc*dmmoab->bs*sizeof(PetscScalar),&varray);CHKERRQ(ierr); 95 96 /* we are using a MOAB Vec - directly copy the tag data to new one */ 97 merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)varray);MBERRNM(merr); 98 for (ifield=0; ifield<dmmoab->numFields; ++ifield) { 99 100 /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */ 101 merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield],1,moab::MB_TYPE_DOUBLE,ntag, 102 moab::MB_TAG_DENSE|moab::MB_TAG_CREAT);MBERRNM(merr); 103 104 /* we are using a MOAB Vec - directly copy the tag data to new one */ 105 for(i=0; i < dmmoab->nloc; i++) { 106 farray[i] = varray[i*dmmoab->bs+ifield]; 107 } 108 109 merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray);MBERRNM(merr); 110 /* make sure the parallel exchange for ghosts are done appropriately */ 111 merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal);MBERRNM(merr); 112 } 113 ierr = PetscFree(varray);CHKERRQ(ierr); 114 } 115 ierr = PetscFree(farray);CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "DMMoabSetFields" 122 PetscErrorCode DMMoabSetFields(DM dm,PetscInt numFields,const char** fields) 123 { 124 PetscErrorCode ierr; 125 PetscInt i; 126 DM_Moab *dmmoab; 127 128 PetscFunctionBegin; 129 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 130 dmmoab = (DM_Moab*)(dm)->data; 131 132 /* first deallocate the existing field structure */ 133 if (dmmoab->fieldNames) { 134 for(i=0; i<dmmoab->numFields; i++) { 135 ierr = PetscFree(dmmoab->fieldNames[i]);CHKERRQ(ierr); 136 } 137 ierr = PetscFree(dmmoab->fieldNames);CHKERRQ(ierr); 138 } 139 140 /* now re-allocate and assign field names */ 141 dmmoab->numFields = numFields; 142 ierr = PetscMalloc(sizeof(char*)*numFields,&dmmoab->fieldNames);CHKERRQ(ierr); 143 if (fields) { 144 for(i=0; i<dmmoab->numFields; i++) { 145 ierr = PetscStrallocpy(fields[i], (char**) &dmmoab->fieldNames[i]);CHKERRQ(ierr); 146 } 147 } 148 PetscFunctionReturn(0); 149 } 150 151 #undef __FUNCT__ 152 #define __FUNCT__ "DMMoabSetFieldName" 153 /*@C 154 DMMoabSetFieldName - Sets the name of a field in the DM 155 156 Not Collective 157 158 Input Parameters: 159 + dm - the DM object 160 . field - the field number 161 - fieldName - the field name 162 163 Level: developer 164 Note: Needs to be called after DMMoabSetFields with correct numFields 165 166 .seealso: DMMoabSetFields() 167 @*/ 168 PetscErrorCode DMMoabSetFieldName(DM dm, PetscInt field, const char fieldName[]) 169 { 170 PetscErrorCode ierr; 171 DM_Moab *dmmoab; 172 173 PetscFunctionBegin; 174 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 175 PetscValidCharPointer(fieldName,3); 176 dmmoab = (DM_Moab*)(dm)->data; 177 178 if ((field < 0) || (field >= dmmoab->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "DM field %d should be in [%d, %d)", field, 0, dmmoab->numFields); 179 ierr = PetscFree(dmmoab->fieldNames[field]);CHKERRQ(ierr); 180 ierr = PetscStrallocpy(fieldName, (char**) &dmmoab->fieldNames[field]);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 185 #undef __FUNCT__ 186 #define __FUNCT__ "DMMoabGetFieldDof" 187 PetscErrorCode DMMoabGetFieldDof(DM dm,moab::EntityHandle point,PetscInt field,PetscInt* dof) 188 { 189 DM_Moab *dmmoab; 190 191 PetscFunctionBegin; 192 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 193 dmmoab = (DM_Moab*)(dm)->data; 194 195 *dof=dmmoab->gidmap[(PetscInt)point]; 196 PetscFunctionReturn(0); 197 } 198 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "DMMoabGetFieldDofs" 202 PetscErrorCode DMMoabGetFieldDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt field,PetscInt* dof) 203 { 204 PetscInt i; 205 PetscErrorCode ierr; 206 DM_Moab *dmmoab; 207 208 PetscFunctionBegin; 209 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 210 PetscValidPointer(points,2); 211 dmmoab = (DM_Moab*)(dm)->data; 212 213 if (!dof) { 214 ierr = PetscMalloc(sizeof(PetscInt)*npoints, &dof);CHKERRQ(ierr); 215 } 216 217 /* first get the local indices */ 218 if (dmmoab->bs > 1) { 219 for (i=0; i<npoints; ++i) 220 dof[i] = dmmoab->gidmap[(PetscInt)points[i]]*dmmoab->numFields+field; 221 } 222 else { 223 /* assume all fields have equal distribution */ 224 for (i=0; i<npoints; ++i) 225 dof[i] = dmmoab->gidmap[(PetscInt)points[i]]+field*dmmoab->n; 226 } 227 PetscFunctionReturn(0); 228 } 229 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "DMMoabGetFieldDofsLocal" 233 PetscErrorCode DMMoabGetFieldDofsLocal(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt field,PetscInt* dof) 234 { 235 PetscInt i,offset; 236 PetscErrorCode ierr; 237 DM_Moab *dmmoab; 238 239 PetscFunctionBegin; 240 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 241 PetscValidPointer(points,2); 242 dmmoab = (DM_Moab*)(dm)->data; 243 244 if (!dof) { 245 ierr = PetscMalloc(sizeof(PetscInt)*npoints, &dof);CHKERRQ(ierr); 246 } 247 248 if (dmmoab->bs > 1) { 249 for (i=0; i<npoints; ++i) 250 dof[i] = dmmoab->lidmap[(PetscInt)points[i]]*dmmoab->numFields+field; 251 } 252 else { 253 offset = field*dmmoab->n; /* assume all fields have equal distribution */ 254 for (i=0; i<npoints; ++i) 255 dof[i] = dmmoab->lidmap[(PetscInt)points[i]]+offset; 256 } 257 PetscFunctionReturn(0); 258 } 259 260 261 #undef __FUNCT__ 262 #define __FUNCT__ "DMMoabGetDofs" 263 PetscErrorCode DMMoabGetDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof) 264 { 265 PetscInt i,field,offset; 266 PetscErrorCode ierr; 267 DM_Moab *dmmoab; 268 269 PetscFunctionBegin; 270 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 271 PetscValidPointer(points,2); 272 dmmoab = (DM_Moab*)(dm)->data; 273 274 if (!dof) { 275 ierr = PetscMalloc(sizeof(PetscInt)*dmmoab->numFields*npoints, &dof);CHKERRQ(ierr); 276 } 277 278 if (dmmoab->bs > 1) { 279 for (field=0; field<dmmoab->numFields; ++field) { 280 for (i=0; i<npoints; ++i) 281 dof[i*dmmoab->numFields+field] = dmmoab->gidmap[(PetscInt)points[i]]*dmmoab->numFields+field; 282 } 283 } 284 else { 285 for (field=0; field<dmmoab->numFields; ++field) { 286 offset = field*dmmoab->n; /* assume all fields have equal distribution */ 287 for (i=0; i<npoints; ++i) 288 dof[i*dmmoab->numFields+field] = dmmoab->gidmap[(PetscInt)points[i]]+offset; 289 } 290 } 291 PetscFunctionReturn(0); 292 } 293 294 295 #undef __FUNCT__ 296 #define __FUNCT__ "DMMoabGetDofsLocal" 297 PetscErrorCode DMMoabGetDofsLocal(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof) 298 { 299 PetscInt i,field,offset; 300 PetscErrorCode ierr; 301 DM_Moab *dmmoab; 302 303 PetscFunctionBegin; 304 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 305 PetscValidPointer(points,2); 306 dmmoab = (DM_Moab*)(dm)->data; 307 308 if (!dof) { 309 ierr = PetscMalloc(sizeof(PetscInt)*dmmoab->numFields*npoints, &dof);CHKERRQ(ierr); 310 } 311 312 if (dmmoab->bs > 1) { 313 for (field=0; field<dmmoab->numFields; ++field) { 314 for (i=0; i<npoints; ++i) 315 dof[i*dmmoab->numFields+field] = dmmoab->lidmap[(PetscInt)points[i]]*dmmoab->numFields+field; 316 } 317 } 318 else { 319 for (field=0; field<dmmoab->numFields; ++field) { 320 offset = field*dmmoab->n; /* assume all fields have equal distribution */ 321 for (i=0; i<npoints; ++i) 322 dof[i*dmmoab->numFields+field] = dmmoab->lidmap[(PetscInt)points[i]]+offset; 323 } 324 } 325 PetscFunctionReturn(0); 326 } 327 328 329 #undef __FUNCT__ 330 #define __FUNCT__ "DMMoabGetDofsBlocked" 331 PetscErrorCode DMMoabGetDofsBlocked(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof) 332 { 333 PetscInt i; 334 DM_Moab *dmmoab; 335 PetscErrorCode ierr; 336 337 PetscFunctionBegin; 338 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 339 PetscValidPointer(points,2); 340 dmmoab = (DM_Moab*)(dm)->data; 341 342 if (!dof) { 343 ierr = PetscMalloc(sizeof(PetscInt)*npoints, &dof);CHKERRQ(ierr); 344 } 345 346 for (i=0; i<npoints; ++i) { 347 dof[i]=dmmoab->gidmap[(PetscInt)points[i]]; 348 } 349 PetscFunctionReturn(0); 350 } 351 352 353 #undef __FUNCT__ 354 #define __FUNCT__ "DMMoabGetDofsBlockedLocal" 355 PetscErrorCode DMMoabGetDofsBlockedLocal(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof) 356 { 357 PetscInt i; 358 DM_Moab *dmmoab; 359 PetscErrorCode ierr; 360 361 PetscFunctionBegin; 362 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 363 PetscValidPointer(points,2); 364 dmmoab = (DM_Moab*)(dm)->data; 365 366 if (!dof) { 367 ierr = PetscMalloc(sizeof(PetscInt)*npoints, &dof);CHKERRQ(ierr); 368 } 369 370 for (i=0; i<npoints; ++i) 371 dof[i] = dmmoab->lidmap[(PetscInt)points[i]]; 372 PetscFunctionReturn(0); 373 } 374 375 376 #undef __FUNCT__ 377 #define __FUNCT__ "DMMoabGetVertexDofsBlocked" 378 PetscErrorCode DMMoabGetVertexDofsBlocked(DM dm,PetscInt** dof) 379 { 380 DM_Moab *dmmoab; 381 382 PetscFunctionBegin; 383 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 384 dmmoab = (DM_Moab*)(dm)->data; 385 386 *dof = dmmoab->gidmap; 387 PetscFunctionReturn(0); 388 } 389 390 391 #undef __FUNCT__ 392 #define __FUNCT__ "DMMoabGetVertexDofsBlockedLocal" 393 PetscErrorCode DMMoabGetVertexDofsBlockedLocal(DM dm,PetscInt** dof) 394 { 395 DM_Moab *dmmoab; 396 397 PetscFunctionBegin; 398 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 399 PetscValidPointer(dof,2); 400 dmmoab = (DM_Moab*)(dm)->data; 401 402 *dof = dmmoab->lidmap; 403 PetscFunctionReturn(0); 404 } 405 406 407 408