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