1 #include <petsc-private/dmmbimpl.h> /*I "petscdm.h" I*/ 2 #include <petsc-private/vecimpl.h> /*I "petscdm.h" I*/ 3 4 #include <petscdmmoab.h> 5 6 static PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool); 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "DMCreateMatrix_Moab" 10 PetscErrorCode DMCreateMatrix_Moab(DM dm, MatType mtype,Mat *J) 11 { 12 PetscErrorCode ierr; 13 ISLocalToGlobalMapping ltog; 14 MatType mltype; 15 PetscInt i,nloc,count,dof,innz,ionz,nsize; 16 DM_Moab *dmmoab = (DM_Moab*)dm->data; 17 moab::Range *range=dmmoab->vowned; 18 PetscInt *gindices,*nnz,*onz; 19 char *tmp=0; 20 moab::ErrorCode merr; 21 22 PetscFunctionBegin; 23 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 24 PetscValidPointer(J,3); 25 nloc = dmmoab->vowned->size() * dmmoab->bs; 26 27 /* create the Matrix and set its type as specified by user */ 28 ierr = MatCreate(dmmoab->pcomm->comm(), J);CHKERRQ(ierr); 29 ierr = MatSetSizes(*J, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 30 ierr = MatSetType(*J, mtype);CHKERRQ(ierr); 31 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 32 33 /* next, need to allocate the non-zero arrays to enable pre-allocation */ 34 ierr = MatGetType(*J, &mltype);CHKERRQ(ierr); /* in case user overrode the default type from command-line, re-check the type */ 35 ierr = PetscStrstr(mltype, "baij", &tmp);CHKERRQ(ierr); 36 nsize = (tmp ? dmmoab->nloc:dmmoab->nloc*dmmoab->bs); 37 38 /* allocate the nnz, onz arrays based on block size and local nodes */ 39 ierr = PetscMalloc(nsize*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 40 ierr = PetscMemzero(nnz,sizeof(PetscInt)*nsize);CHKERRQ(ierr); 41 ierr = PetscMalloc(nsize*sizeof(PetscInt),&onz);CHKERRQ(ierr); 42 ierr = PetscMemzero(onz,sizeof(PetscInt)*nsize);CHKERRQ(ierr); 43 44 /* compute the nonzero pattern based on MOAB connectivity data for local elements */ 45 ierr = DMMoab_Compute_NNZ_From_Connectivity(dm,&innz,nnz,&ionz,onz,(tmp?PETSC_TRUE:PETSC_FALSE));CHKERRQ(ierr); 46 47 if (dmmoab->bs > 1 && tmp) { 48 /* Block matrix created, now set local to global mapping */ 49 ierr = PetscMalloc(sizeof(PetscInt)*range->size()*dmmoab->bs, &gindices);CHKERRQ(ierr); 50 moab::Range::iterator iter; 51 for(iter = range->begin(),count=0; iter != range->end(); iter++,count+=dmmoab->bs) { 52 merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,&(*iter),1,&dof);MBERRNM(merr); 53 for(i=0; i<dmmoab->bs; ++i) 54 gindices[count+i] = (dof)*dmmoab->bs+i; 55 } 56 57 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,range->size(),gindices,PETSC_COPY_VALUES,<og);CHKERRQ(ierr); 58 ierr = MatSetLocalToGlobalMappingBlock(*J,ltog,ltog);CHKERRQ(ierr); 59 60 /* Clean up */ 61 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr); 62 ierr = PetscFree(gindices);CHKERRQ(ierr); 63 64 } else { 65 ierr = MatSetBlockSize(*J, dmmoab->bs);CHKERRQ(ierr); 66 } 67 68 /* set preallocation based on different supported Mat types */ 69 // ierr = MatSeqAIJSetPreallocation(*J, innz, nnz);CHKERRQ(ierr); 70 // ierr = MatMPIAIJSetPreallocation(*J, innz, nnz, ionz, onz);CHKERRQ(ierr); 71 ierr = MatSeqAIJSetPreallocation(*J, innz, 0);CHKERRQ(ierr); 72 ierr = MatMPIAIJSetPreallocation(*J, innz, 0, ionz, 0);CHKERRQ(ierr); 73 ierr = MatSeqBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz);CHKERRQ(ierr); 74 ierr = MatMPIBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz, ionz, onz);CHKERRQ(ierr); 75 76 ierr = MatSetDM(*J, dm);CHKERRQ(ierr); 77 ierr = PetscFree(nnz);CHKERRQ(ierr); 78 ierr = PetscFree(onz);CHKERRQ(ierr); 79 PetscFunctionReturn(0); 80 } 81 82 83 #undef __FUNCT__ 84 #define __FUNCT__ "DMMoab_Compute_NNZ_From_Connectivity" 85 PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm,PetscInt* innz,PetscInt* nnz,PetscInt* ionz,PetscInt* onz,PetscBool isbaij) 86 { 87 PetscErrorCode ierr; 88 PetscInt i,j,f,nloc,vpere,bs,nsize,nghost_found; 89 DM_Moab *dmmoab = (DM_Moab*)dm->data; 90 const moab::EntityHandle *connect; 91 moab::Range *vowned=dmmoab->vowned,*elocal=dmmoab->elocal,*eghost=dmmoab->eghost,adjs,visited; 92 moab::Range::iterator iter; 93 PetscInt dof,doff,ndofs; 94 moab::Tag id_tag; 95 moab::ErrorCode merr; 96 PetscSection section; 97 98 PetscFunctionBegin; 99 bs = dmmoab->bs; 100 nloc = dmmoab->nloc; 101 nsize = (isbaij ? nloc:nloc*bs); 102 103 ierr = DMMoabGetLocalToGlobalTag(dm,&id_tag);CHKERRQ(ierr); 104 105 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 106 107 for(iter = elocal->begin(),i=0; iter != elocal->end(); iter++,i++) { 108 /* Get connectivity information in canonical ordering for the local element */ 109 merr = dmmoab->mbiface->get_connectivity((*iter),connect,vpere);MBERRNM(merr); 110 111 nghost_found=0; 112 /* loop over vertices and update the number of connectivity */ 113 for (j=0;j<vpere;j++) { 114 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 115 if (giter != dmmoab->vghost->end()) nghost_found++; 116 } 117 // PetscPrintf(PETSC_COMM_WORLD, "Element %D, VPERE %D, nGhost %D\n", i,vpere,nghost_found); 118 119 /* loop over vertices and update the number of connectivity */ 120 for (j=0;j<vpere;j++) { 121 122 /* 123 moab::Range::const_iterator giter = visited.find(connect[j]); 124 if (giter != visited.end()) continue; 125 126 visited.insert(connect[j]); 127 */ 128 129 merr = dmmoab->mbiface->get_adjacencies(&connect[j],1,1,true,adjs,moab::Interface::INTERSECT); 130 // PetscPrintf(PETSC_COMM_WORLD, "Point %D, nADJS %D\n", connect[j], adjs.size()); 131 132 merr = dmmoab->mbiface->tag_get_data(id_tag,&connect[j],1,&dof);MBERRNM(merr); 133 ierr = PetscSectionGetDof(section, connect[j], &ndofs); 134 135 for (f=0;f<dmmoab->nfields;f++) { 136 ierr = PetscSectionGetFieldOffset(section, connect[j], f, &doff); 137 138 // nnz[doff]+=(vpere-nghost_found-1); /* leave out self to avoid repeats -> node shared by multiple elements */ 139 nnz[doff]+=adjs.size(); /* leave out self to avoid repeats -> node shared by multiple elements */ 140 onz[doff]+=nghost_found; /* add ghost non-owned nodes */ 141 142 // PetscPrintf(PETSC_COMM_WORLD, "Point %D, Field %D, DOFPP %D \t DOF %D, OFFSET %D : NNZ = %D, ONZ = %D\n", connect[j], f, ndofs, dof, doff, nnz[doff], onz[doff]); 143 adjs.clear(); 144 } 145 } 146 } 147 148 visited.clear(); 149 150 if (innz) *innz=0; 151 if (ionz) *ionz=0; 152 for (i=0;i<nsize;i++) { 153 nnz[i]+=1; /* self count the node */ 154 if (!isbaij) { 155 nnz[i]*=bs; 156 onz[i]*=bs; 157 } 158 // PetscPrintf(PETSC_COMM_WORLD, "Index %D: NNZ = %D, ONZ = %D\n", i, nnz[i], onz[i]); 159 160 if (innz && nnz[i]>*innz) *innz=nnz[i]; 161 if (ionz && onz[i]>*ionz) *ionz=onz[i]; 162 } 163 // PetscPrintf(PETSC_COMM_WORLD, "MAX: NNZ = %D, ONZ = %D\n", *innz, *ionz); 164 PetscFunctionReturn(0); 165 } 166 167 168 #if 0 169 170 #undef __FUNCT__ 171 #define __FUNCT__ "DMMoab_Compute_NNZ_From_Connectivity_OLD" 172 PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity_OLD(DM dm,PetscInt* innz,PetscInt* nnz,PetscInt* ionz,PetscInt* onz,PetscBool isbaij) 173 { 174 PetscErrorCode ierr; 175 PetscInt i,j,k,f,nloc,count,vpere,bs,nsize,nghost_found; 176 DM_Moab *dmmoab = (DM_Moab*)dm->data; 177 const moab::EntityHandle *connect; 178 moab::Range *vowned=dmmoab->vowned,*elocal=dmmoab->elocal,*eghost=dmmoab->eghost,adjs; 179 moab::Range::iterator iter; 180 PetscInt *vertex_ids,firstvtx,dof,doff,ndofs,offset; 181 moab::Tag id_tag; 182 std::vector<moab::EntityHandle> storage; 183 moab::ErrorCode merr; 184 PetscSection section; 185 186 PetscFunctionBegin; 187 bs = dmmoab->bs; 188 nloc = dmmoab->nloc; 189 nsize = (isbaij ? nloc:nloc*bs); 190 191 ierr = DMMoabGetLocalToGlobalTag(dm,&id_tag);CHKERRQ(ierr); 192 193 merr = dmmoab->mbiface->tag_get_data(id_tag,&(*vowned->begin()),1,&firstvtx);MBERRNM(merr); 194 //firstvtx = vertex_ids[0]; 195 196 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 197 198 for(iter = elocal->begin(),i=0; iter != elocal->end(); iter++,i++) { 199 /* Get connectivity information in canonical ordering for the local element */ 200 merr = dmmoab->mbiface->get_connectivity((*iter),connect,vpere,false,&storage);MBERRNM(merr); 201 202 nghost_found=0; 203 /* loop over vertices and update the number of connectivity */ 204 for (j=0;j<vpere;j++) { 205 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 206 if (giter != dmmoab->vghost->end()) nghost_found++; 207 } 208 209 /* loop over vertices and update the number of connectivity */ 210 for (j=0;j<vpere;j++) { 211 212 merr = dmmoab->mbiface->get_adjacencies(&connect[j],1,1,true,adjs,moab::Interface::INTERSECT); 213 PetscPrintf(PETSC_COMM_WORLD, "Point %D, nADJS %D\n", connect[j], adjs.size()); 214 215 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 216 if (giter != dmmoab->vghost->end()) continue; 217 218 merr = dmmoab->mbiface->tag_get_data(id_tag,&connect[j],1,&dof);MBERRNM(merr); 219 ierr = PetscSectionGetDof(section, connect[j], &ndofs); 220 221 for (f=0;f<dmmoab->nfields;f++) { 222 ierr = PetscSectionGetFieldOffset(section, connect[j], f, &doff); 223 224 // nnz[doff]+=ndofs*(vpere-nghost_found-1); /* leave out self to avoid repeats -> node shared by multiple elements */ 225 nnz[doff]+=adjs.size()-nghost_found-1; /* leave out self to avoid repeats -> node shared by multiple elements */ 226 onz[doff]+=nghost_found; /* add ghost non-owned nodes */ 227 228 PetscPrintf(PETSC_COMM_WORLD, "Point %D, Field %D, DOFPP %D \t DOF %D, OFFSET %D : NNZ = %D, ONZ = %D\n", connect[j], f, ndofs, dof, doff, nnz[doff], onz[doff]); 229 adjs.clear(); 230 } 231 232 // /* if block format, then all the block data are local only */ 233 // if (!isbaij) { 234 // offset=(dof-firstvtx)*bs; 235 // for (k=0;k<bs;k++) { 236 // nnz[offset+k]+=vpere-nghost_found-1; /* leave out self to avoid repeats -> node shared by multiple elements */ 237 // onz[offset+k]+=nghost_found; /* found a ghost non-owned node */ 238 // } 239 // } 240 // else { 241 // nnz[(dof-firstvtx)]+=vpere-nghost_found-1; /* leave out self to avoid repeats -> node shared by multiple elements */ 242 // onz[(dof-firstvtx)]+=nghost_found; /* found a ghost non-owned node */ 243 // } 244 245 } 246 } 247 248 /* loop through ghosted elements that contain non-ghosted (locally owned) vertices */ 249 for(iter = eghost->begin(),i=0; iter != eghost->end(); iter++,i++) { 250 /* Get connectivity information in canonical ordering for the local element */ 251 merr = dmmoab->mbiface->get_connectivity((*iter),connect,vpere,false,&storage);MBERRNM(merr); 252 253 nghost_found=0; 254 /* loop over vertices and update the number of connectivity */ 255 for (j=0;j<vpere;j++) { 256 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 257 if (giter != dmmoab->vghost->end()) nghost_found++; 258 } 259 260 if (nghost_found == vpere) continue; /* all vertices are ghosted.. */ 261 262 /* loop over vertices and update the number of connectivity */ 263 for (j=0;j<vpere;j++) { 264 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 265 if (giter != dmmoab->vghost->end()) continue; 266 267 merr = dmmoab->mbiface->tag_get_data(id_tag,&connect[j],1,&dof);MBERRNM(merr); 268 ierr = PetscSectionGetDof(section, connect[j], &ndofs); 269 270 for (f=0;f<dmmoab->nfields;f++) { 271 ierr = PetscSectionGetFieldOffset(section, connect[j], f, &doff); 272 273 onz[doff]+=nghost_found; /* add ghost non-owned nodes */ 274 275 PetscPrintf(PETSC_COMM_WORLD, "Point %D, Field %D, DOF %D, OFFSET %D : ONZ = %D\n", connect[j], f, dof, doff, onz[doff]); 276 277 } 278 279 // if (dof-firstvtx < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Invalid local vertex ordering. Found ghost node with ID less than first vertex: [%i, %i]",dof,firstvtx); 280 // 281 // /* if block format, then all the block data are local only */ 282 // if (!isbaij) { 283 // offset=(dof-firstvtx)*bs; 284 // for (k=0;k<bs;k++) 285 // onz[offset+k]+=nghost_found; 286 // } 287 // else { 288 // onz[(dof-firstvtx)]+=nghost_found; 289 // } 290 } 291 } 292 293 if (innz) *innz=0; 294 if (ionz) *ionz=0; 295 for (i=0;i<nsize;i++) { 296 nnz[i]+=1; /* self count the node */ 297 if (!isbaij) { 298 nnz[i]*=bs; 299 onz[i]*=bs; 300 } 301 PetscPrintf(PETSC_COMM_WORLD, "Index %D: NNZ = %D, ONZ = %D\n", i, nnz[i], onz[i]); 302 303 if (innz && nnz[i]>*innz) *innz=nnz[i]; 304 if (ionz && onz[i]>*ionz) *ionz=onz[i]; 305 } 306 PetscPrintf(PETSC_COMM_WORLD, "MAX: NNZ = %D, ONZ = %D\n", *innz, *ionz); 307 PetscFunctionReturn(0); 308 } 309 310 #endif 311 312