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 7 static PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscBool); 8 static PetscErrorCode DMMoab_MatFillMatrixEntries_Private(DM,Mat); 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 ltogb; 16 PetscInt innz,ionz,nlsiz; 17 DM_Moab *dmmoab=(DM_Moab*)dm->data; 18 PetscInt *nnz=0,*onz=0; 19 char *tmp=0; 20 moab::ErrorCode merr; 21 22 PetscFunctionBegin; 23 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 24 PetscValidPointer(J,3); 25 26 /* create the Matrix and set its type as specified by user */ 27 ierr = MatCreate(dmmoab->pcomm->comm(), J);CHKERRQ(ierr); 28 ierr = MatSetSizes(*J, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 29 ierr = MatSetType(*J, mtype);CHKERRQ(ierr); 30 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 31 32 /* next, need to allocate the non-zero arrays to enable pre-allocation */ 33 ierr = MatGetType(*J, &mltype);CHKERRQ(ierr); /* in case user overrode the default type from command-line, re-check the type */ 34 ierr = PetscStrstr(mltype, "baij", &tmp);CHKERRQ(ierr); 35 nsize = (tmp ? dmmoab->nloc:dmmoab->nloc*dmmoab->bs); 36 37 /* allocate the nnz, onz arrays based on block size and local nodes */ 38 ierr = PetscMalloc((nlsiz)*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 39 ierr = PetscMemzero(nnz,sizeof(PetscInt)*(nlsiz));CHKERRQ(ierr); 40 ierr = PetscMalloc(nlsiz*sizeof(PetscInt),&onz);CHKERRQ(ierr); 41 ierr = PetscMemzero(onz,sizeof(PetscInt)*nlsiz);CHKERRQ(ierr); 42 43 /* compute the nonzero pattern based on MOAB connectivity data for local elements */ 44 ierr = DMMoab_Compute_NNZ_From_Connectivity(dm,&innz,nnz,&ionz,onz,(tmp?PETSC_TRUE:PETSC_FALSE));CHKERRQ(ierr); 45 46 /* create the Matrix and set its type as specified by user */ 47 ierr = MatCreate(dmmoab->pcomm->comm(), J);CHKERRQ(ierr); 48 ierr = MatSetSizes(*J, dmmoab->nloc*dmmoab->numFields, dmmoab->nloc*dmmoab->numFields, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 49 ierr = MatSetBlockSize(*J, dmmoab->bs);CHKERRQ(ierr); 50 ierr = MatSetType(*J, mtype);CHKERRQ(ierr); 51 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 52 53 if (!dmmoab->ltog_map) SETERRQ(dmmoab->pcomm->comm(), PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first."); 54 ierr = MatSetLocalToGlobalMapping(*J,dmmoab->ltog_map,dmmoab->ltog_map);CHKERRQ(ierr); 55 ierr = ISLocalToGlobalMappingBlock(dmmoab->ltog_map,dmmoab->bs,<ogb); 56 ierr = MatSetLocalToGlobalMappingBlock(*J,ltogb,ltogb);CHKERRQ(ierr); 57 ierr = ISLocalToGlobalMappingDestroy(<ogb);CHKERRQ(ierr); 58 59 /* set preallocation based on different supported Mat types */ 60 ierr = MatSeqAIJSetPreallocation(*J, innz, nnz);CHKERRQ(ierr); 61 ierr = MatMPIAIJSetPreallocation(*J, innz, nnz, ionz, onz);CHKERRQ(ierr); 62 ierr = MatSeqBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz);CHKERRQ(ierr); 63 ierr = MatMPIBAIJSetPreallocation(*J, dmmoab->bs, innz, nnz, ionz, onz);CHKERRQ(ierr); 64 65 /* clean up temporary memory */ 66 ierr = PetscFree(nnz);CHKERRQ(ierr); 67 ierr = PetscFree(onz);CHKERRQ(ierr); 68 69 /* set up internal matrix data-structures */ 70 ierr = MatSetUp(*J);CHKERRQ(ierr); 71 72 /* set DM reference */ 73 ierr = MatSetDM(*J, dm);CHKERRQ(ierr); 74 75 /* set the correct NNZ pattern by setting matrix entries - make the matrix ready to use */ 76 ierr = DMMoab_MatFillMatrixEntries_Private(dm,*J);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 PetscInt i,f,nloc,vpere,bs,nsize,ivtx,n_nnz,n_onz; 86 DM_Moab *dmmoab = (DM_Moab*)dm->data; 87 const moab::EntityHandle *connect; 88 moab::Range adjs,found,allvlocal,allvghost; 89 moab::Range::iterator iter,jter; 90 std::vector<moab::EntityHandle> storage; 91 moab::EntityHandle vtx; 92 moab::ErrorCode merr; 93 94 PetscFunctionBegin; 95 bs = dmmoab->bs; 96 nloc = dmmoab->nloc; 97 nsize = (isbaij ? nloc:nloc*bs); 98 99 /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */ 100 merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,allvlocal,true);MBERRNM(merr); 101 merr = dmmoab->pcomm->filter_pstatus(allvlocal,PSTATUS_NOT_OWNED,PSTATUS_NOT,-1,&adjs);MBERRNM(merr); 102 allvghost = moab::subtract(allvlocal, adjs); 103 104 /* loop over vertices and update the number of connectivity */ 105 for (j=0;j<vpere;j++) { 106 107 vtx = *iter; 108 adjs.clear(); 109 /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects 110 to the current vertex. We can then decipher if a vertex is ghosted or not and compute the 111 non-zero pattern accordingly. */ 112 merr = dmmoab->mbiface->get_adjacencies(&vtx,1,dmmoab->dim,false,adjs,moab::Interface::INTERSECT); 113 114 /* loop over vertices and update the number of connectivity */ 115 for(jter = adjs.begin(); jter != adjs.end(); jter++) { 116 117 /* Get connectivity information in canonical ordering for the local element */ 118 merr = dmmoab->mbiface->get_connectivity(*jter,connect,vpere,false,&storage);MBERRNM(merr); 119 120 for (i=0; i<vpere; ++i) { 121 if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; 122 if (allvghost.find(connect[i]) != allvghost.end()) n_onz++; 123 else n_nnz++; 124 found.insert(connect[i]); 125 } 126 } 127 128 if (isbaij) { 129 nnz[ivtx]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 130 if (onz) onz[ivtx]=n_onz; /* add ghost non-owned nodes */ 131 } 132 else { 133 for (f=0;f<dmmoab->numFields;f++) { 134 nnz[dmmoab->numFields*ivtx+f]=n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 135 if (onz) onz[dmmoab->numFields*ivtx+f]=n_onz; /* add ghost non-owned nodes */ 136 } 137 } 138 } 139 140 if (innz) *innz=0; 141 if (ionz) *ionz=0; 142 for (i=0;i<nsize;i++) { 143 nnz[i]+=1; /* self count the node */ 144 /* check if we got overzealous */ 145 nnz[i]=(nnz[i]>dmmoab->nloc ? dmmoab->nloc:nnz[i]); 146 if (!isbaij) { 147 nnz[i]*=bs; 148 onz[i]*=bs; 149 } 150 PetscInfo3(dm, "Vertex ID: %D \t NNZ = %D \t ONZ = %D.\n",i,nnz[i],onz[i]); 151 152 if (innz && (nnz[i]>*innz)) *innz=nnz[i]; 153 if ((ionz && onz) && (onz[i]>*ionz)) *ionz=onz[i]; 154 } 155 PetscFunctionReturn(0); 156 } 157 158 159 #undef __FUNCT__ 160 #define __FUNCT__ "DMMoab_MatFillMatrixEntries_Private" 161 PetscErrorCode DMMoab_MatFillMatrixEntries_Private(DM dm, Mat A) 162 { 163 DM_Moab *dmmoab = (DM_Moab*)dm->data; 164 PetscInt nconn = 0,prev_nconn = 0; 165 const moab::EntityHandle *connect; 166 PetscScalar *locala=NULL; 167 PetscInt *dof_indices=NULL; 168 PetscErrorCode ierr; 169 170 PetscFunctionBegin; 171 /* loop over local elements */ 172 for(moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) { 173 const moab::EntityHandle ehandle = *iter; 174 175 /* Get connectivity information: */ 176 ierr = DMMoabGetElementConnectivity(dm, ehandle, &nconn, &connect);CHKERRQ(ierr); 177 178 /* if we have mixed elements or arrays have not been initialized - Allocate now */ 179 if (prev_nconn != nconn) { 180 if (locala) { 181 ierr = PetscFree(locala);CHKERRQ(ierr); 182 ierr = PetscFree(dof_indices);CHKERRQ(ierr); 183 } 184 ierr = PetscMalloc(sizeof(PetscScalar)*nconn*nconn*dmmoab->numFields*dmmoab->numFields,&locala);CHKERRQ(ierr); 185 ierr = PetscMemzero(locala,sizeof(PetscScalar)*nconn*nconn*dmmoab->numFields*dmmoab->numFields);CHKERRQ(ierr); 186 ierr = PetscMalloc(sizeof(PetscInt)*nconn,&dof_indices);CHKERRQ(ierr); 187 PetscInfo2(dm, "Allocating new memory and zeroing out for locala. [Prev: %d, Curr-size: %d].\n",prev_nconn*prev_nconn*dmmoab->numFields*dmmoab->numFields,nconn*nconn*dmmoab->numFields*dmmoab->numFields); 188 prev_nconn=nconn; 189 } 190 191 /* get the global DOF number to appropriately set the element contribution in the RHS vector */ 192 ierr = DMMoabGetDofsBlockedLocal(dm, nconn, connect, dof_indices);CHKERRQ(ierr); 193 194 /* set the values directly into appropriate locations. Can alternately use VecSetValues */ 195 ierr = MatSetValuesBlockedLocal(A, nconn, dof_indices, nconn, dof_indices, locala, INSERT_VALUES);CHKERRQ(ierr); 196 } 197 198 /* clean up memory */ 199 ierr = PetscFree(locala);CHKERRQ(ierr); 200 ierr = PetscFree(dof_indices);CHKERRQ(ierr); 201 202 /* finish assembly */ 203 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 204 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 205 PetscFunctionReturn(0); 206 } 207 208