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,f,nloc,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,adjs,visited; 90 moab::Range::iterator iter; 91 PetscInt dof,doff,ndofs; 92 moab::Tag id_tag; 93 moab::ErrorCode merr; 94 PetscSection section; 95 96 PetscFunctionBegin; 97 bs = dmmoab->bs; 98 nloc = dmmoab->nloc; 99 nsize = (isbaij ? nloc:nloc*bs); 100 101 ierr = DMMoabGetLocalToGlobalTag(dm,&id_tag);CHKERRQ(ierr); 102 103 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 104 105 for(iter = elocal->begin(),i=0; iter != elocal->end(); iter++,i++) { 106 /* Get connectivity information in canonical ordering for the local element */ 107 merr = dmmoab->mbiface->get_connectivity((*iter),connect,vpere);MBERRNM(merr); 108 109 nghost_found=0; 110 /* loop over vertices and update the number of connectivity */ 111 for (j=0;j<vpere;j++) { 112 moab::Range::const_iterator giter = dmmoab->vghost->find(connect[j]); 113 if (giter != dmmoab->vghost->end()) nghost_found++; 114 } 115 116 /* loop over vertices and update the number of connectivity */ 117 for (j=0;j<vpere;j++) { 118 119 /* 120 moab::Range::const_iterator giter = visited.find(connect[j]); 121 if (giter != visited.end()) continue; 122 123 visited.insert(connect[j]); 124 */ 125 126 merr = dmmoab->mbiface->get_adjacencies(&connect[j],1,1,true,adjs,moab::Interface::INTERSECT); 127 128 merr = dmmoab->mbiface->tag_get_data(id_tag,&connect[j],1,&dof);MBERRNM(merr); 129 ierr = PetscSectionGetDof(section, connect[j], &ndofs); 130 131 for (f=0;f<dmmoab->nfields;f++) { 132 ierr = PetscSectionGetFieldOffset(section, connect[j], f, &doff); 133 134 // nnz[doff]+=(vpere-nghost_found-1); /* leave out self to avoid repeats -> node shared by multiple elements */ 135 nnz[doff]+=adjs.size(); /* leave out self to avoid repeats -> node shared by multiple elements */ 136 onz[doff]+=nghost_found; /* add ghost non-owned nodes */ 137 138 adjs.clear(); 139 } 140 } 141 } 142 143 visited.clear(); 144 145 if (innz) *innz=0; 146 if (ionz) *ionz=0; 147 for (i=0;i<nsize;i++) { 148 nnz[i]+=1; /* self count the node */ 149 if (!isbaij) { 150 nnz[i]*=bs; 151 onz[i]*=bs; 152 } 153 154 if (innz && nnz[i]>*innz) *innz=nnz[i]; 155 if (ionz && onz[i]>*ionz) *ionz=onz[i]; 156 } 157 PetscFunctionReturn(0); 158 } 159 160