[petsc-users] Extracting reduced stiffness matrix

[Carl-Johan Thore]

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, &ltogm); 
VecGetLocalToGlobalMapping(NN, &ltogm); 
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);