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 152 #undef __FUNCT__ 153 #define __FUNCT__ "DMMoabGetFieldDof" 154 PetscErrorCode DMMoabGetFieldDof(DM dm,moab::EntityHandle point,PetscInt field,PetscInt* dof) 155 { 156 DM_Moab *dmmoab; 157 158 PetscFunctionBegin; 159 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 160 dmmoab = (DM_Moab*)(dm)->data; 161 162 *dof=dmmoab->gidmap[(PetscInt)point]; 163 PetscFunctionReturn(0); 164 } 165 166 167 #undef __FUNCT__ 168 #define __FUNCT__ "DMMoabGetFieldDofs" 169 PetscErrorCode DMMoabGetFieldDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt field,PetscInt* dof) 170 { 171 PetscInt i; 172 PetscErrorCode ierr; 173 DM_Moab *dmmoab; 174 175 PetscFunctionBegin; 176 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 177 PetscValidPointer(points,2); 178 dmmoab = (DM_Moab*)(dm)->data; 179 180 if (!dof) { 181 ierr = PetscMalloc(sizeof(PetscInt)*npoints, &dof);CHKERRQ(ierr); 182 } 183 184 /* first get the local indices */ 185 if (dmmoab->bs > 1) { 186 for (i=0; i<npoints; ++i) 187 dof[i] = dmmoab->gidmap[(PetscInt)points[i]]*dmmoab->numFields+field; 188 } 189 else { 190 /* assume all fields have equal distribution */ 191 for (i=0; i<npoints; ++i) 192 dof[i] = dmmoab->gidmap[(PetscInt)points[i]]+field*dmmoab->n; 193 } 194 PetscFunctionReturn(0); 195 } 196 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "DMMoabGetFieldDofsLocal" 200 PetscErrorCode DMMoabGetFieldDofsLocal(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt field,PetscInt* dof) 201 { 202 PetscInt i,offset; 203 PetscErrorCode ierr; 204 DM_Moab *dmmoab; 205 206 PetscFunctionBegin; 207 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 208 PetscValidPointer(points,2); 209 dmmoab = (DM_Moab*)(dm)->data; 210 211 if (!dof) { 212 ierr = PetscMalloc(sizeof(PetscInt)*npoints, &dof);CHKERRQ(ierr); 213 } 214 215 if (dmmoab->bs > 1) { 216 for (i=0; i<npoints; ++i) 217 dof[i] = dmmoab->lidmap[(PetscInt)points[i]]*dmmoab->numFields+field; 218 } 219 else { 220 offset = field*dmmoab->n; /* assume all fields have equal distribution */ 221 for (i=0; i<npoints; ++i) 222 dof[i] = dmmoab->lidmap[(PetscInt)points[i]]+offset; 223 } 224 PetscFunctionReturn(0); 225 } 226 227 228 #undef __FUNCT__ 229 #define __FUNCT__ "DMMoabGetDofs" 230 PetscErrorCode DMMoabGetDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof) 231 { 232 PetscInt i,field,offset; 233 PetscErrorCode ierr; 234 DM_Moab *dmmoab; 235 236 PetscFunctionBegin; 237 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 238 PetscValidPointer(points,2); 239 dmmoab = (DM_Moab*)(dm)->data; 240 241 if (!dof) { 242 ierr = PetscMalloc(sizeof(PetscInt)*dmmoab->numFields*npoints, &dof);CHKERRQ(ierr); 243 } 244 245 if (dmmoab->bs > 1) { 246 for (field=0; field<dmmoab->numFields; ++field) { 247 for (i=0; i<npoints; ++i) 248 dof[i*dmmoab->numFields+field] = dmmoab->gidmap[(PetscInt)points[i]]*dmmoab->numFields+field; 249 } 250 } 251 else { 252 for (field=0; field<dmmoab->numFields; ++field) { 253 offset = field*dmmoab->n; /* assume all fields have equal distribution */ 254 for (i=0; i<npoints; ++i) 255 dof[i*dmmoab->numFields+field] = dmmoab->gidmap[(PetscInt)points[i]]+offset; 256 } 257 } 258 PetscFunctionReturn(0); 259 } 260 261 262 #undef __FUNCT__ 263 #define __FUNCT__ "DMMoabGetDofsLocal" 264 PetscErrorCode DMMoabGetDofsLocal(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof) 265 { 266 PetscInt i,field,offset; 267 PetscErrorCode ierr; 268 DM_Moab *dmmoab; 269 270 PetscFunctionBegin; 271 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 272 PetscValidPointer(points,2); 273 dmmoab = (DM_Moab*)(dm)->data; 274 275 if (!dof) { 276 ierr = PetscMalloc(sizeof(PetscInt)*dmmoab->numFields*npoints, &dof);CHKERRQ(ierr); 277 } 278 279 if (dmmoab->bs > 1) { 280 for (field=0; field<dmmoab->numFields; ++field) { 281 for (i=0; i<npoints; ++i) 282 dof[i*dmmoab->numFields+field] = dmmoab->lidmap[(PetscInt)points[i]]*dmmoab->numFields+field; 283 } 284 } 285 else { 286 for (field=0; field<dmmoab->numFields; ++field) { 287 offset = field*dmmoab->n; /* assume all fields have equal distribution */ 288 for (i=0; i<npoints; ++i) 289 dof[i*dmmoab->numFields+field] = dmmoab->lidmap[(PetscInt)points[i]]+offset; 290 } 291 } 292 PetscFunctionReturn(0); 293 } 294 295 296 #undef __FUNCT__ 297 #define __FUNCT__ "DMMoabGetDofsBlocked" 298 PetscErrorCode DMMoabGetDofsBlocked(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof) 299 { 300 PetscInt i; 301 DM_Moab *dmmoab; 302 PetscErrorCode ierr; 303 304 PetscFunctionBegin; 305 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 306 PetscValidPointer(points,2); 307 dmmoab = (DM_Moab*)(dm)->data; 308 309 if (!dof) { 310 ierr = PetscMalloc(sizeof(PetscInt)*npoints, &dof);CHKERRQ(ierr); 311 } 312 313 for (i=0; i<npoints; ++i) { 314 dof[i]=dmmoab->gidmap[(PetscInt)points[i]]; 315 } 316 PetscFunctionReturn(0); 317 } 318 319 320 #undef __FUNCT__ 321 #define __FUNCT__ "DMMoabGetDofsBlockedLocal" 322 PetscErrorCode DMMoabGetDofsBlockedLocal(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof) 323 { 324 PetscInt i; 325 DM_Moab *dmmoab; 326 PetscErrorCode ierr; 327 328 PetscFunctionBegin; 329 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 330 PetscValidPointer(points,2); 331 dmmoab = (DM_Moab*)(dm)->data; 332 333 if (!dof) { 334 ierr = PetscMalloc(sizeof(PetscInt)*npoints, &dof);CHKERRQ(ierr); 335 } 336 337 for (i=0; i<npoints; ++i) 338 dof[i] = dmmoab->lidmap[(PetscInt)points[i]]; 339 PetscFunctionReturn(0); 340 } 341 342 343 #undef __FUNCT__ 344 #define __FUNCT__ "DMMoabGetVertexDofsBlocked" 345 PetscErrorCode DMMoabGetVertexDofsBlocked(DM dm,PetscInt** dof) 346 { 347 DM_Moab *dmmoab; 348 349 PetscFunctionBegin; 350 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 351 dmmoab = (DM_Moab*)(dm)->data; 352 353 *dof = dmmoab->gidmap; 354 PetscFunctionReturn(0); 355 } 356 357 358 #undef __FUNCT__ 359 #define __FUNCT__ "DMMoabGetVertexDofsBlockedLocal" 360 PetscErrorCode DMMoabGetVertexDofsBlockedLocal(DM dm,PetscInt** dof) 361 { 362 DM_Moab *dmmoab; 363 364 PetscFunctionBegin; 365 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 366 PetscValidPointer(dof,2); 367 dmmoab = (DM_Moab*)(dm)->data; 368 369 *dof = dmmoab->lidmap; 370 PetscFunctionReturn(0); 371 } 372 373 374 375