Hi, I'm trying to set up and solve the "reduced system" obtained by omitting all rows and columns associated with zero-prescribed unknowns from the global stiffness matrix K. In Matlab-notation:
u = K(freedofs,freedofs)\F(freedofs) where freedofs are the free unknowns (the purpose is to compare this method with using MatZeroRowsCols to solve the "full system" for problems with many prescribed unknowns (does anyone have benchmark data?)). I've tried various ways to accomplish the same thing in PETSc. The main difficulty seems to be the construction of the IS freedofs (see code sketch below) from the vector NN for use in MatCreateSubMatrix and VecGetSubVector. I've managed to do this but not in such a way that it works for arbitrary partitionings of the mesh (DMDA). I think that something along the lines of the sketch below might be the "correct" way of doing it, but I suspect there is some issue with ghost values which causes freedofs to be slightly incorrect (if imported into Matlab I can see that there are a few duplicate values for example). Any suggestions? Kind regards, Carl-Johan Code sketch: IS freedofs; DMCreateGlobalVector(da, &NN); // Some code to populate NN // NN[i] = 0.0 if DOF i is prescribed // NN[i] = 1.0 if DOF i is free // ... // Get local version of NN DMCreateLocalVector(da, &NNloc); DMGlobalToLocalBegin(da, NN, INSERT_VALUES, NNloc); DMGlobalToLocalEnd(da, NN, INSERT_VALUES, NNloc); VecGetArray(NNloc, &nnlocarr); // Get global indices from local ISLocalToGlobalMapping ltogm; const PetscInt *g_idx; PetscInt locsize; DMGetLocalToGlobalMapping(da, <ogm); VecGetLocalToGlobalMapping(NN, <ogm); ISLocalToGlobalMappingGetIndices(ltogm, &g_idx); ISLocalToGlobalMappingGetSize(ltogm, &locsize); PetscScalar nfreedofs; VecSum(NN, &nfreedofs); PetscInt idx[(PetscInt) nfreedofs]; PetscInt n=0; for (PetscInt i=0; i<locsize; i++) { if (nnlocarr[i] == 1.0) { idx[n++] = g_idx[i]; } } ISCreateGeneral(PETSC_COMM_WORLD,n,idx,PETSC_COPY_VALUES,&freedofs); // Assemble full matrices K and F ... // Extract reduced versions MatCreateSubMatrix(K, freedofs, freedofs, MAT_INITIAL_MATRIX, &Kred); CHKERRQ(ierr); VecGetSubVector(F,freedofs,&Fred); CHKERRQ(ierr);