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 #include <MBTagConventions.hpp> 6 #include <sstream> 7 8 static PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool); 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "DMCreateMatrix_Moab" 12 PetscErrorCode DMCreateMatrix_Moab(DM dm, MatType mtype,Mat *J) 13 { 14 PetscErrorCode ierr; 15 ISLocalToGlobalMapping ltog; 16 MatType mltype; 17 PetscInt i,nloc,count,dof,innz,ionz,nsize; 18 DM_Moab *dmmoab = (DM_Moab*)dm->data; 19 moab::Range *range=dmmoab->vowned; 20 PetscInt *gindices,*nnz,*onz; 21 char *tmp=0; 22 moab::ErrorCode merr; 23 24 PetscFunctionBegin; 25 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 26 PetscValidPointer(J,3); 27 nloc = dmmoab->vowned->size() * dmmoab->bs; 28 // ierr = DMMoabCreateMatrix(dm,PETSC_NULL,dmmoab->bs,dmmoab->vowned,PETSC_TRUE,PETSC_TRUE,J);CHKERRQ(ierr); 29 30 ierr = MatCreate(dmmoab->pcomm->comm(), J);CHKERRQ(ierr); 31 ierr = MatSetSizes(*J, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 32 ierr = MatSetType(*J, mtype);CHKERRQ(ierr); 33 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 34 35 /* allocate the non-zeros properly */ 36 ierr = MatGetType(*J, &mltype);CHKERRQ(ierr); /* in case user overrode the default type from command-line, re-check the type */ 37 ierr = PetscStrstr(mltype, "baij", &tmp);CHKERRQ(ierr); 38 nsize = (tmp ? dmmoab->nloc:dmmoab->nloc*dmmoab->bs); 39 40 /* allocate the nnz, onz arrays based on block size and local nodes */ 41 ierr = PetscMalloc(nsize*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 42 ierr = PetscMemzero(nnz,sizeof(PetscInt)*nsize);CHKERRQ(ierr); 43 ierr = PetscMalloc(nsize*sizeof(PetscInt),&onz);CHKERRQ(ierr); 44 ierr = PetscMemzero(onz,sizeof(PetscInt)*nsize);CHKERRQ(ierr); 45 46 /* compute the nonzero pattern based on MOAB connectivity data for local elements */ 47 ierr = DMMoab_Compute_NNZ_From_Connectivity(dm,&innz,nnz,&ionz,onz,(tmp?PETSC_TRUE:PETSC_FALSE));CHKERRQ(ierr); 48 49 if (dmmoab->bs > 1 && tmp) { 50 // Block matrix created, now set local to global mapping: 51 ierr = PetscMalloc(sizeof(PetscInt)*range->size()*dmmoab->bs, &gindices);CHKERRQ(ierr); 52 moab::Range::iterator iter; 53 for(iter = range->begin(),count=0; iter != range->end(); iter++,count+=dmmoab->bs) { 54 merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,&(*iter),1,&dof);MBERRNM(merr); 55 for(i=0; i<dmmoab->bs; ++i) 56 gindices[count+i] = (dof)*dmmoab->bs+i; 57 } 58 59 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,range->size(),gindices,PETSC_COPY_VALUES,<og);CHKERRQ(ierr); 60 ierr = MatSetLocalToGlobalMappingBlock(*J,ltog,ltog);CHKERRQ(ierr); 61 62 // Clean up: 63 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr); 64 ierr = PetscFree(gindices);CHKERRQ(ierr); 65 66 } else { 67 ierr = MatSetBlockSize(*J, dmmoab->bs);CHKERRQ(ierr); 68 } 69 70 /* set preallocation based on different supported Mat types */ 71 ierr = MatSeqAIJSetPreallocation(*J, innz, nnz);CHKERRQ(ierr); 72 ierr = MatMPIAIJSetPreallocation(*J, innz, nnz, ionz, onz);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,k,nloc,count,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; 92 moab::Range::iterator iter; 93 PetscInt *vertex_ids,firstvtx,dof,offset; 94 moab::Tag id_tag; 95 moab::ErrorCode merr; 96 97 PetscFunctionBegin; 98 bs = dmmoab->bs; 99 nloc = dmmoab->nloc; 100 nsize = (isbaij ? nloc:nloc*bs); 101 102 ierr = DMMoabGetLocalToGlobalTag(dm,&id_tag);CHKERRQ(ierr); 103 merr = dmmoab->mbiface->tag_iterate(id_tag,dmmoab->vlocal->begin(),dmmoab->vlocal->end(), 104 count,reinterpret_cast<void*&>(vertex_ids));MBERRNM(merr); 105 106 merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,&(*vowned->begin()),1,&firstvtx);MBERRNM(merr); 107 108 for(iter = elocal->begin(),i=0; iter != elocal->end(); iter++,i++) { 109 /* Get connectivity information in canonical ordering for the local element */ 110 merr = dmmoab->mbiface->get_connectivity((*iter),connect,vpere);MBERRNM(merr); 111 112 nghost_found=0; 113 /* loop over vertices and update the number of connectivity */ 114 for (j=0;j<vpere;j++) { 115 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 116 if (giter != dmmoab->vghost->end()) nghost_found++; 117 } 118 119 /* loop over vertices and update the number of connectivity */ 120 for (j=0;j<vpere;j++) { 121 merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,&connect[j],1,&dof);MBERRNM(merr); 122 123 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 124 if (giter != dmmoab->vghost->end()) continue; 125 126 /* if block format, then all the block data are local only */ 127 if (!isbaij) { 128 offset=(dof-firstvtx)*bs; 129 for (k=0;k<bs;k++) { 130 nnz[offset+k]+=vpere-1; 131 onz[offset+k]+=nghost_found; /* found a ghost non-owned node */ 132 } 133 } 134 else { 135 nnz[(dof-firstvtx)]+=vpere-1; 136 onz[(dof-firstvtx)]+=nghost_found; /* found a ghost non-owned node */ 137 } 138 } 139 } 140 141 for(iter = eghost->begin(),i=0; iter != eghost->end(); iter++,i++) { 142 /* Get connectivity information in canonical ordering for the local element */ 143 merr = dmmoab->mbiface->get_connectivity((*iter),connect,vpere);MBERRNM(merr); 144 145 nghost_found=0; 146 /* loop over vertices and update the number of connectivity */ 147 for (j=0;j<vpere;j++) { 148 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 149 if (giter != dmmoab->vghost->end()) nghost_found++; 150 } 151 152 /* loop over vertices and update the number of connectivity */ 153 for (j=0;j<vpere;j++) { 154 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 155 if (giter != dmmoab->vghost->end()) continue; 156 157 merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,&connect[j],1,&dof);MBERRNM(merr); 158 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); 159 160 /* if block format, then all the block data are local only */ 161 if (!isbaij) { 162 offset=(dof-firstvtx)*bs; 163 for (k=0;k<bs;k++) 164 onz[offset+k]+=nghost_found; 165 } 166 else { 167 onz[(dof-firstvtx)]+=nghost_found; 168 } 169 } 170 } 171 172 if (innz) *innz=0; 173 if (ionz) *ionz=0; 174 for (i=0;i<nsize;i++) { 175 nnz[i]+=1; /* self count the node */ 176 if (!isbaij) { 177 nnz[i]*=bs; 178 onz[i]*=bs; 179 } 180 if (innz && nnz[i]>*innz) *innz=nnz[i]; 181 if (ionz && onz[i]>*ionz) *ionz=onz[i]; 182 } 183 // PetscPrintf (PETSC_COMM_WORLD, "\n Maximum NNZ = %D and ONZ = %D \n", *innz,*ionz); 184 185 /* Print for debug purposes only */ 186 /* 187 IS nzis,onzis; 188 ierr = ISCreateGeneral(((PetscObject)dm)->comm, nsize, nnz, PETSC_COPY_VALUES, &nzis);CHKERRQ(ierr); 189 ierr = ISCreateGeneral(((PetscObject)dm)->comm, nsize, onz, PETSC_COPY_VALUES, &onzis);CHKERRQ(ierr); 190 ierr = ISView(nzis, 0);CHKERRQ(ierr); 191 ierr = ISView(onzis, 0);CHKERRQ(ierr); 192 ierr = ISDestroy(&onzis);CHKERRQ(ierr); 193 ierr = ISDestroy(&onzis);CHKERRQ(ierr); 194 std::cin.get(); 195 */ 196 197 198 PetscFunctionReturn(0); 199 } 200 201 202