Skip to content

Commit

Permalink
Mapping expanded to global and back.
Browse files Browse the repository at this point in the history 
Signed-off-by: Jed Brown <jed@59A2.org>
  • Loading branch information
@jedbrown
jedbrown committed Nov 8, 2008
1 parent b5119d8 commit 1f71a0a
Show file tree
Hide file tree
Showing 7 changed files with 159 additions and 25 deletions.
6 changes: 5 additions & 1 deletion include/dohpfs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,12 @@ EXTERN dErr dFSSetDegree(dFS,dJacobi,dMeshTag);
EXTERN dErr dFSAddBdy(dFS,const char*,dMeshESH,dMeshTag,dBool,PF); /* name, facets, orientation tag, flip orientation?, normal -> constraints */
EXTERN dErr dFSSetFromOptions(dFS);
EXTERN dErr dFSSetType(dFS,const dFSType);
EXTERN dErr dFSCreateLocalVector(dFS,Vec*);
EXTERN dErr dFSCreateExpandedVector(dFS,Vec*);
EXTERN dErr dFSCreateGlobalVector(dFS,Vec*);
EXTERN dErr dFSGlobalToExpandedBegin(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSGlobalToExpandedEnd(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSExpandedToGlobal(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSGetMatrix(dFS,const MatType,Mat*);
EXTERN dErr dFSBuildSpace(dFS);
EXTERN dErr dFSGetBoundaryType(dFS,dInt,const dMeshEH[],dBdyType[]);
EXTERN dErr dFSCreateSubspace(dFS,dInt,dFSBoundary,dFS*);
Expand Down
3 changes: 2 additions & 1 deletion include/dohptype.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,8 @@ typedef unsigned char dEntStatus;

#define dMax(a,b) PetscMax(a,b)
#define dMin(a,b) PetscMin(a,b)
#define dSqr(a) PetscSqr(a)
#define dSqr(a) PetscSqr(a)
#define dAbs(a) PetscAbs(a)

#define dGamma(a) tgamma(a) /* This is defined in math.h as of C99. */

Expand Down
1 change: 1 addition & 0 deletions include/private/fsimpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ struct _p_dFS {
dEFS *efs; /**< Element function space, defined for all entities */
Mat C; /**< full-order constraint matrix (element dofs to local numbering) */
Mat Cp; /**< preconditioning constraint matrix (element dofs to local numbering, as sparse as possible) */
Vec weight; /**< Vector in global space, used to compensate for overcounting after local to global */
void *data;
};

Expand Down
10 changes: 7 additions & 3 deletions src/fs/impls/cont/cont.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,11 @@ static dErr dFSBuildSpace_Cont(dFS fs)
err = dFree(fieldSpec);dCHK(err);

err = SlicedCreate(((dObject)fs)->comm,&fs->sliced);dCHK(err);
err = SlicedSetGhosts(fs->sliced,1,fs->n,fs->nlocal-fs->n,ghostind);dCHK(err);
if (fs->nlocal-fs->n == 0) {
err = SlicedSetGhosts(fs->sliced,1,fs->n,0,NULL);dCHK(err);
} else {
err = SlicedSetGhosts(fs->sliced,1,fs->n,fs->nlocal-fs->n,ghostind);dCHK(err);
}
err = dFree(ghostind);dCHK(err);

/**
Expand Down Expand Up @@ -178,8 +182,8 @@ static dErr dFSBuildSpace_Cont(dFS fs)
err = dJacobiGetConstraintCount(fs->jacobi,nregions,xind,xstart,deg,&ma,nnz,pnnz);dCHK(err);

/* We don't solve systems with these so it will never make sense for them to use a different format */
err = MatCreateSeqAIJ(PETSC_COMM_SELF,fs->m,fs->nlocal,-1,nnz,&fs->C);dCHK(err);
err = MatCreateSeqAIJ(PETSC_COMM_SELF,fs->m,fs->nlocal,-1,pnnz,&fs->Cp);dCHK(err);
err = MatCreateSeqAIJ(PETSC_COMM_SELF,fs->m,fs->nlocal,1,nnz,&fs->C);dCHK(err);
err = MatCreateSeqAIJ(PETSC_COMM_SELF,fs->m,fs->nlocal,1,pnnz,&fs->Cp);dCHK(err);
err = dFree2(nnz,pnnz);dCHK(err);

err = dJacobiAddConstraints(fs->jacobi,nregions,xind,xstart,istart,deg,&ma,fs->C,fs->Cp);dCHK(err);
Expand Down
106 changes: 89 additions & 17 deletions src/fs/interface/fs.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -191,55 +191,127 @@ dErr dFSDestroy(dFS fs)
*/
dErr dFSBuildSpace(dFS fs)
{
dMesh mesh = fs->mesh;
iMesh_Instance mi;
Vec x,g;
dErr err;

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
/**
/*
* FIXME:
*
*/
err = dMeshGetInstance(mesh,&mi);dCHK(err);
if (fs->ops->buildspace) {
err = (*fs->ops->buildspace)(fs);dCHK(err);
}

/* Determine the number of elements in which each dof appears */
err = dFSCreateExpandedVector(fs,&x);dCHK(err);
err = dFSCreateGlobalVector(fs,&g);dCHK(err);
err = VecSet(x,1);dCHK(err);
err = VecZeroEntries(g);dCHK(err);
err = dFSExpandedToGlobal(fs,x,ADD_VALUES,g);dCHK(err);
err = VecGhostUpdateBegin(g,ADD_VALUES,SCATTER_FORWARD);dCHK(err);
err = VecGhostUpdateEnd(g,ADD_VALUES,SCATTER_FORWARD);dCHK(err);
err = VecDestroy(x);dCHK(err);

fs->spacebuilt = true;
dFunctionReturn(0);
}

dErr dFSCreateLocalVector(dFS fs,Vec *U)
dErr dFSCreateExpandedVector(dFS fs,Vec *x)
{
dErr err;

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
dValidPointer(x,2);
err = MatGetVecs(fs->C,NULL,x);dCHK(err);
dFunctionReturn(0);
}

dErr dFSCreateGlobalVector(dFS fs,Vec *g)
{
dErr err;

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
dValidPointer(g,2);
err = SlicedCreateGlobalVector(fs->sliced,g);dCHK(err);
dFunctionReturn(0);
}

dErr dFSGlobalToExpandedBegin(dFS fs,Vec g,InsertMode imode,Vec x)
{
dErr err;

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
dValidPointer(U,2);
*U = 0;
if (fs->sliced) {
err = SlicedCreateLocalVector(fs->sliced,U);dCHK(err);
} else {
dERROR(1,"no sliced");
dValidHeader(g,VEC_COOKIE,2);
dValidHeader(x,VEC_COOKIE,4);
err = VecGhostUpdateBegin(g,imode,SCATTER_FORWARD);dCHK(err);
dFunctionReturn(0);
}

dErr dFSGlobalToExpandedEnd(dFS fs,Vec g,InsertMode imode,Vec x)
{
Vec lform;
dErr err;

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
dValidHeader(g,VEC_COOKIE,2);
dValidHeader(x,VEC_COOKIE,4);
err = VecGhostUpdateEnd(g,imode,SCATTER_FORWARD);dCHK(err);
err = VecGhostGetLocalForm(g,&lform);dCHK(err);
switch (imode) {
case INSERT_VALUES:
err = MatMult(fs->C,lform,x);dCHK(err);
break;
case ADD_VALUES:
err = MatMultAdd(fs->C,lform,x,x);dCHK(err);
break;
default:
dERROR(1,"InsertMode %d not supported",imode);
}
err = VecGhostRestoreLocalForm(g,&lform);dCHK(err);
dFunctionReturn(0);
}

dErr dFSCreateGlobalVector(dFS fs,Vec *U)
dErr dFSExpandedToGlobal(dFS fs,Vec x,InsertMode imode,Vec g)
{
Vec lform;
dErr err;

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
dValidHeader(U,VEC_COOKIE,2);
if (fs->sliced) {
err = SlicedCreateGlobalVector(fs->sliced,U);dCHK(err);
} else {
dERROR(1,"no sliced");
dValidHeader(x,VEC_COOKIE,2);
dValidHeader(g,VEC_COOKIE,4);
err = VecGhostGetLocalForm(g,&lform);dCHK(err);
switch (imode) {
case INSERT_VALUES:
err = MatMultTranspose(fs->C,x,lform);dCHK(err);
break;
case ADD_VALUES:
err = MatMultTransposeAdd(fs->C,x,lform,lform);dCHK(err);
break;
default:
dERROR(1,"InsertMode %d not supported",imode);
}
err = VecGhostRestoreLocalForm(g,&lform);dCHK(err);
dFunctionReturn(0);
}

dErr dFSGetMatrix(dFS fs,const MatType mtype,Mat *J)
{
dErr err;

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
dValidPointer(J,3);
err = SlicedGetMatrix(fs->sliced,mtype,J);dCHK(err);
dFunctionReturn(0);dCHK(err);
}

/**
* Create a function space based on a simpler space
*
Expand Down
10 changes: 7 additions & 3 deletions src/fs/interface/fsreg.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,13 @@ static PetscFList FSList = 0;
*/
static const struct _dFSOps defaultFSOps = { .view = dFSView,
.createglobalvector = dFSCreateGlobalVector,
.createlocalvector = dFSCreateLocalVector,
.destroy = dFSDestroy,
.impldestroy = 0};
.createlocalvector = dFSCreateExpandedVector,
.localtoglobal = dFSExpandedToGlobal,
.globaltolocalbegin = dFSGlobalToExpandedBegin,
.globaltolocalend = dFSGlobalToExpandedEnd,
.getmatrix = dFSGetMatrix,
.destroy = dFSDestroy,
.impldestroy = 0};

dErr dFSCreate(MPI_Comm comm,dFS *infs)
{
Expand Down
48 changes: 48 additions & 0 deletions src/fs/tests/ex1.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,6 +139,52 @@ static dErr examine(dMesh mesh,dMeshTag tag)
dFunctionReturn(0);
}

static dErr useFS(dFS fs)
{
Vec x,y,g;
dScalar *xx;
dReal norm;
dInt nx=-1;
dErr err;

dFunctionBegin;
err = dFSCreateGlobalVector(fs,&g);dCHK(err);
err = dFSCreateExpandedVector(fs,&x);dCHK(err);
err = VecGetSize(x,&nx);dCHK(err);
{
dInt n,rstart,rend;
err = VecGetSize(g,&n);dCHK(err);
err = VecGetOwnershipRange(g,&rstart,&rend);dCHK(err);
err = PetscPrintf(PETSC_COMM_WORLD,"expanded size %d, global size %d, range %d .. %d \n",nx,n,rstart,rend);dCHK(err);
}
err = VecDuplicate(x,&y);dCHK(err);
err = VecGetArray(x,&xx);dCHK(err);
for (dInt i=0; i<nx; i++) {
xx[i] = 1000 + 1.0 * i;
}
err = VecRestoreArray(x,&xx);dCHK(err);
err = VecSet(x,1.0);dCHK(err);
err = dFSExpandedToGlobal(fs,x,INSERT_VALUES,g);dCHK(err);
err = dFSGlobalToExpandedBegin(fs,g,INSERT_VALUES,y);dCHK(err);
err = dFSGlobalToExpandedEnd(fs,g,INSERT_VALUES,y);dCHK(err);
{
dScalar xsum,ysum,gsum;
err = VecSum(x,&xsum);dCHK(err);
err = VecSum(g,&gsum);dCHK(err);
err = VecSum(y,&ysum);dCHK(err);
err = PetscPrintf(PETSC_COMM_WORLD,"|x| = %f, |g| = %f, |y| = %f\n",xsum,gsum,ysum);dCHK(err);
if (dAbs(xsum-gsum) > 1e-14) dERROR(1,"Expanded sum does not match global sum");
/* There are 16 points on the interface between the elements, these points get double-counted in both elements, so
* the expanded sum is larger than the global sum by 32. */
if (dAbs(ysum-gsum-32) > 1e-14) dERROR(1,"Unexpected expanded sum %f != %f + 32",ysum,gsum);
}
if (norm != 0) dERROR(1,"Norm does not agree.");
err = VecDestroy(x);dCHK(err);
err = VecDestroy(y);dCHK(err);
err = VecDestroy(g);dCHK(err);
dFunctionReturn(0);
}

int main(int argc,char *argv[])
{
/* const char pTagName[] = "dohp_partition"; */
Expand Down Expand Up @@ -181,6 +227,8 @@ int main(int argc,char *argv[])
err = dFSSetDegree(fs,jac,dtag);dCHK(err);
err = dFSSetFromOptions(fs);dCHK(err);

err = useFS(fs);dCHK(err);

err = examine(mesh,dtag);dCHK(err);

err = dFSDestroy(fs);dCHK(err);
Expand Down

0 comments on commit 1f71a0a

Please sign in to comment.