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 = MatSeqBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz);CHKERRQ(ierr); 72 ierr = MatMPIBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz, ionz, onz);CHKERRQ(ierr); 73 74 ierr = MatSetDM(*J, dm);CHKERRQ(ierr); 75 ierr = PetscFree(nnz);CHKERRQ(ierr); 76 ierr = PetscFree(onz);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 81 #undef __FUNCT__ 82 #define __FUNCT__ "DMMoab_Compute_NNZ_From_Connectivity" 83 PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm,PetscInt* innz,PetscInt* nnz,PetscInt* ionz,PetscInt* onz,PetscBool isbaij) 84 { 85 PetscErrorCode ierr; 86 PetscInt i,j,k,nloc,count,vpere,bs,nsize,nghost_found; 87 DM_Moab *dmmoab = (DM_Moab*)dm->data; 88 const moab::EntityHandle *connect; 89 moab::Range *vowned=dmmoab->vowned,*elocal=dmmoab->elocal,*eghost=dmmoab->eghost; 90 moab::Range::iterator iter; 91 PetscInt *vertex_ids,firstvtx,dof,offset; 92 moab::Tag id_tag; 93 moab::ErrorCode merr; 94 95 PetscFunctionBegin; 96 bs = dmmoab->bs; 97 nloc = dmmoab->nloc; 98 nsize = (isbaij ? nloc:nloc*bs); 99 100 ierr = DMMoabGetLocalToGlobalTag(dm,&id_tag);CHKERRQ(ierr); 101 merr = dmmoab->mbiface->tag_iterate(id_tag,dmmoab->vlocal->begin(),dmmoab->vlocal->end(), 102 count,reinterpret_cast<void*&>(vertex_ids));MBERRNM(merr); 103 104 merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,&(*vowned->begin()),1,&firstvtx);MBERRNM(merr); 105 106 for(iter = elocal->begin(),i=0; iter != elocal->end(); iter++,i++) { 107 /* Get connectivity information in canonical ordering for the local element */ 108 merr = dmmoab->mbiface->get_connectivity((*iter),connect,vpere);MBERRNM(merr); 109 110 nghost_found=0; 111 /* loop over vertices and update the number of connectivity */ 112 for (j=0;j<vpere;j++) { 113 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 114 if (giter != dmmoab->vghost->end()) nghost_found++; 115 } 116 117 /* loop over vertices and update the number of connectivity */ 118 for (j=0;j<vpere;j++) { 119 merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,&connect[j],1,&dof);MBERRNM(merr); 120 121 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 122 if (giter != dmmoab->vghost->end()) continue; 123 124 /* if block format, then all the block data are local only */ 125 if (!isbaij) { 126 offset=(dof-firstvtx)*bs; 127 for (k=0;k<bs;k++) { 128 nnz[offset+k]+=vpere-1; 129 onz[offset+k]+=nghost_found; /* found a ghost non-owned node */ 130 } 131 } 132 else { 133 nnz[(dof-firstvtx)]+=vpere-1; 134 onz[(dof-firstvtx)]+=nghost_found; /* found a ghost non-owned node */ 135 } 136 } 137 } 138 139 for(iter = eghost->begin(),i=0; iter != eghost->end(); iter++,i++) { 140 /* Get connectivity information in canonical ordering for the local element */ 141 merr = dmmoab->mbiface->get_connectivity((*iter),connect,vpere);MBERRNM(merr); 142 143 nghost_found=0; 144 /* loop over vertices and update the number of connectivity */ 145 for (j=0;j<vpere;j++) { 146 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 147 if (giter != dmmoab->vghost->end()) nghost_found++; 148 } 149 150 /* loop over vertices and update the number of connectivity */ 151 for (j=0;j<vpere;j++) { 152 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 153 if (giter != dmmoab->vghost->end()) continue; 154 155 merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,&connect[j],1,&dof);MBERRNM(merr); 156 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); 157 158 /* if block format, then all the block data are local only */ 159 if (!isbaij) { 160 offset=(dof-firstvtx)*bs; 161 for (k=0;k<bs;k++) 162 onz[offset+k]+=nghost_found; 163 } 164 else { 165 onz[(dof-firstvtx)]+=nghost_found; 166 } 167 } 168 } 169 170 if (innz) *innz=0; 171 if (ionz) *ionz=0; 172 for (i=0;i<nsize;i++) { 173 nnz[i]+=1; /* self count the node */ 174 if (!isbaij) { 175 nnz[i]*=bs; 176 onz[i]*=bs; 177 } 178 if (innz && nnz[i]>*innz) *innz=nnz[i]; 179 if (ionz && onz[i]>*ionz) *ionz=onz[i]; 180 } 181 PetscFunctionReturn(0); 182 } 183 184 185