Skip to content

Commit

Permalink
Further splitting of dFSBuildSpace_Cont
Browse files Browse the repository at this point in the history
  • Loading branch information
@jedbrown
jedbrown committed Apr 13, 2011
1 parent da0d39f commit e691465
Show file tree
Hide file tree
Showing 5 changed files with 139 additions and 125 deletions.
1 change: 1 addition & 0 deletions include/dohpfsimpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,7 @@ struct _p_dFSRotation {
};

extern dErr dFSCreateLocalToGlobal_Private(dFS fs,dInt n,dInt nc,dInt ngh,dInt *ghidx,dInt rstart);
extern dErr dFSPopulatePartitionedSets_Private(dFS FS,dMeshAdjacency meshadj);
extern dErr dFSBuildSpaceOffsets_Private(dFS fs,dMeshTag indexTag,const dInt inodes[],dInt rstart,dInt crstart,dInt nents,dMeshEH ents[],dInt *ghstart);
extern dErr dFSBuildSpaceVectors_Private(dFS fs,dMeshTag indexTag,const dInt inodes[],dInt rstart,dInt ghents_s,const dMeshEH ghents[]);
extern dErr DMDestroy_dFS(DM);
Expand Down
45 changes: 6 additions & 39 deletions src/fs/impls/cont/cont.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,7 +254,7 @@ dErr dFSBuildSpace_Cont_CreateElemAssemblyMats(dFS fs,const dInt idx[],const dMe
*
* @param fs The space to build
*/
static dErr dFSBuildSpace_Cont(dFS fs)
dErr dFSBuildSpace_Cont(dFS fs)
{
MPI_Comm comm = ((dObject)fs)->comm;
/* \bug The fact that we aren't using our context here indicates that much/all of the logic here could move up into dFS */
Expand All @@ -265,10 +265,9 @@ static dErr dFSBuildSpace_Cont(dFS fs)
iMesh_Instance mi;
dEntTopology *regTopo;
dPolynomialOrder *deg,*regBDeg;
dInt *inodes,*xnodes,nregions,*bstat,ents_a,ents_s,ghents_s,*intdata,*idx;
dInt *inodes,*xnodes,nregions,ents_a,ents_s,ghents_s,*idx;
dInt *xstart,xcnt;
dInt rstart,crstart;
dIInt ierr;
dMeshEH *ents,*ghents;
dErr err;

Expand All @@ -282,39 +281,9 @@ static dErr dFSBuildSpace_Cont(dFS fs)

/* Allocate a workspace that's plenty big, so that we don't have to allocate memory constantly */
ents_a = ma.nents;
err = dMallocA3(ents_a,&ents,ents_a,&intdata,ents_a,&idx);dCHK(err);
err = dMallocA4(ents_a,&ents,ents_a,&idx,ents_a,&deg,ents_a,&inodes);dCHK(err);

/* Partition entities in active set into owned explicit, owned Dirichlet, and ghost */
{
dInt nboundaries,ghstart;
dMeshESH *bdysets;
iMesh_addEntArrToSet(mi,ma.ents,ma.nents,fs->set.explicit,&ierr);dICHK(mi,ierr);
/* Move ghost ents from \a explicitSet to \a ghostSet */
iMesh_getEntitiesRec(mi,fs->set.explicit,dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr);
err = dMeshPartitionOnOwnership(mesh,ents,ents_s,&ghstart);dCHK(err);
iMesh_rmvEntArrFromSet(mi,ents+ghstart,ents_s-ghstart,fs->set.explicit,&ierr);dICHK(mi,ierr);
iMesh_addEntArrToSet(mi,ents+ghstart,ents_s-ghstart,fs->set.ghost,&ierr);dICHK(mi,ierr);
/* Move owned Dirichlet ents from \a explicitSet to \a dirichletSet */
err = dMeshGetNumSubsets(mesh,fs->set.boundaries,0,&nboundaries);dCHK(err);
if (!nboundaries) goto after_boundaries;
err = dMallocA2(nboundaries,&bdysets,nboundaries,&bstat);dCHK(err);
err = dMeshGetSubsets(mesh,fs->set.boundaries,0,bdysets,nboundaries,NULL);dCHK(err);
err = dMeshTagSGetData(mesh,fs->tag.bstatus,bdysets,nboundaries,bstat,nboundaries,dDATA_INT);dCHK(err);
for (int i=0; i<nboundaries; i++) {
if (bstat[i] & dFSBSTATUS_DIRICHLET) {
iMesh_getEntitiesRec(mi,bdysets[i],dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr);
err = dMeshPartitionOnOwnership(mesh,ents,ents_s,&ghstart);dCHK(err);
iMesh_rmvEntArrFromSet(mi,ents,ghstart,fs->set.explicit,&ierr);dICHK(mi,ierr);
iMesh_addEntArrToSet(mi,ents,ghstart,fs->set.dirichlet,&ierr);dICHK(mi,ierr);
}
if (bstat[i] & dFSBSTATUS_WEAK) {
iMesh_getEntitiesRec(mi,bdysets[i],dTYPE_FACE,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr);
iMesh_addEntArrToSet(mi,ents,ents_s,fs->set.weakFace,&ierr);dICHK(mi,ierr);
}
}
err = dFree2(bdysets,bstat);dCHK(err);
}
after_boundaries:
err = dFSPopulatePartitionedSets_Private(fs,&ma);dCHK(err);

err = dMeshPopulateOrderedSet_Private(mesh,fs->set.ordered,fs->set.explicit,fs->set.dirichlet,fs->set.ghost,fs->orderingtype);dCHK(err);
{
Expand All @@ -336,7 +305,6 @@ static dErr dFSBuildSpace_Cont(dFS fs)
if (ents_s != ents_a) dERROR(PETSC_COMM_SELF,PETSC_ERR_PLIB,"wrong set size");

/* Get number of nodes for all entities */
err = dMallocA2(ma.nents,&deg,ma.nents,&inodes);dCHK(err);
err = dMeshTagGetData(mesh,fs->tag.degree,ma.ents,ma.nents,deg,ma.nents,dDATA_INT);dCHK(err);
/* Fill the \a inodes array with the number of interior nodes for each (topology,degree) pair */
err = dJacobiGetNodeCount(fs->jacobi,ma.nents,ma.topo,deg,inodes,NULL);dCHK(err);
Expand Down Expand Up @@ -397,9 +365,8 @@ static dErr dFSBuildSpace_Cont(dFS fs)

err = dFSBuildSpace_Cont_CreateElemAssemblyMats(fs,idx,&ma,deg,&fs->E,&fs->Ep);dCHK(err);

err = dMeshRestoreAdjacency(fs->mesh,fs->set.active,&meshAdj);dCHK(err); /* Any reason to leave this around for longer? */
err = dFree2(deg,inodes);dCHK(err);
err = dFree3(ents,intdata,idx);dCHK(err);
err = dMeshRestoreAdjacency(fs->mesh,fs->set.active,&meshAdj);dCHK(err);
err = dFree4(ents,idx,deg,inodes);dCHK(err);
dFunctionReturn(0);
}

Expand Down
2 changes: 2 additions & 0 deletions src/fs/impls/cont/cont.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,8 @@ typedef struct {
Mat B;
} dFS_Cont;

extern dErr dFSBuildSpace_Cont(dFS fs);

dEXTERN_C_END

#endif /* _CONT_H */
86 changes: 0 additions & 86 deletions src/fs/interface/fs.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -369,92 +369,6 @@ dErr dFSBuildSpace(dFS fs)
dFunctionReturn(0);
}

/* Set offsets (global, closure, local) of first node associated with every entity.
*
* @param indexTag index into inodes
* @param inodes number of interior nodes associated with each entity inodes[idx[i]]
* @param rstart relative start for global explicit dofs
* @param crstart relative start for global closure dofs
* @param nents size of array to hold entities
* @param ents returns with entities in the proper ordering of the FS
* @param ghstart returns starting offset of ghost entities in ents
*/
dErr dFSBuildSpaceOffsets_Private(dFS fs,dMeshTag indexTag,const dInt inodes[],dInt rstart,dInt crstart,dInt nents,dMeshEH ents[],dInt *ghstart)
{
dInt i,scan,nentsExplicit,nentsDirichlet,*idx,*offset;
dMesh mesh;
dErr err;

dFunctionBegin;
err = dFSGetMesh(fs,&mesh);dCHK(err);
err = dMallocA2(nents,&offset,nents,&idx);dCHK(err);
/* We assume that orderedSet contains explicitSet+dirichletSet+ghostSet (in that order) */
err = dMeshGetEnts(mesh,fs->set.ordered,dTYPE_ALL,dTOPO_ALL,ents,nents,NULL);dCHK(err);
err = dMeshGetNumEnts(mesh,fs->set.explicit,dTYPE_ALL,dTOPO_ALL,&nentsExplicit);dCHK(err);
err = dMeshGetNumEnts(mesh,fs->set.dirichlet,dTYPE_ALL,dTOPO_ALL,&nentsDirichlet);dCHK(err);
err = dMeshTagGetData(mesh,indexTag,ents,nents,idx,nents,dDATA_INT);dCHK(err);

/* global offset */
for (i=0,scan=rstart; i<nentsExplicit; scan+=inodes[idx[i++]])
offset[i] = scan;
for ( ; i<nents; i++) offset[i] = -1;
err = dMeshTagSetData(mesh,fs->tag.goffset,ents,nents,offset,nents,dDATA_INT);dCHK(err);

/* global closure offset */
for (i=0,scan=crstart; i<nentsExplicit+nentsDirichlet; scan+=inodes[idx[i++]])
offset[i] = scan;
for ( ; i<nents; i++) offset[i] = -1;
err = dMeshTagSetData(mesh,fs->tag.gcoffset,ents,nents,offset,nents,dDATA_INT);dCHK(err);

/* local index */
for (i=0,scan=0; i<nents; scan+=inodes[idx[i++]])
offset[i] = scan;
err = dMeshTagSetData(mesh,fs->tag.loffset,ents,nents,offset,nents,dDATA_INT);dCHK(err);

/* communicate global and closure offset for ghosts */
err = dMeshTagBcast(mesh,fs->tag.goffset);dCHK(err);
err = dMeshTagBcast(mesh,fs->tag.gcoffset);dCHK(err);

err = dFree2(offset,idx);dCHK(err);
*ghstart = nentsExplicit + nentsDirichlet;
dFunctionReturn(0);
}

dErr dFSBuildSpaceVectors_Private(dFS fs,dMeshTag indexTag,const dInt inodes[],dInt rstart,dInt ghents_s,const dMeshEH ghents[])
{
dErr err;
dInt *gcoffset,*ghidx,*idx;
dMesh mesh;

dFunctionBegin;
err = dMallocA3(ghents_s,&gcoffset,ghents_s,&ghidx,ghents_s,&idx);dCHK(err);
err = dFSGetMesh(fs,&mesh);dCHK(err);
err = dMeshTagGetData(mesh,indexTag,ghents,ghents_s,idx,ghents_s,dDATA_INT);dCHK(err);

/* Retrieve ghost offsets, to create localupdate. */
err = dMeshTagGetData(mesh,fs->tag.gcoffset,ghents,ghents_s,gcoffset,ghents_s,dDATA_INT);dCHK(err);
for (dInt i=0; i<ghents_s; i++) { /* Paranoia: confirm that all ghost entities were updated. */
if (gcoffset[i] < 0) dERROR(PETSC_COMM_SELF,1,"Tag exchange did not work");
}

/* Set ghost indices of every node using \a ghidx, create global vector. */
{
dInt gh=0;
for (dInt i=0; i<ghents_s; i++) {
for (dInt j=0; j<inodes[idx[i]]; j++) ghidx[gh++] = gcoffset[i] + j;
}
if (gh != fs->ngh) dERROR(PETSC_COMM_SELF,1,"Ghost count inconsistent");
err = VecCreateDohp(((dObject)fs)->comm,fs->bs,fs->n,fs->nc,fs->ngh,ghidx,&fs->gvec);dCHK(err);
}
err = dFree3(gcoffset,ghidx,idx);dCHK(err);

/* Create fs->bmapping and fs->mapping */
err = dFSCreateLocalToGlobal_Private(fs,fs->n,fs->nc,fs->ngh,ghidx,rstart);dCHK(err);

err = VecDohpCreateDirichletCache(fs->gvec,&fs->dcache,&fs->dscat);dCHK(err);
dFunctionReturn(0);dCHK(err);
}

dErr dFSCreateExpandedVector(dFS fs,Vec *x)
{
dErr err;
Expand Down
130 changes: 130 additions & 0 deletions src/fs/interface/fsbuild.c
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include <dohpfsimpl.h>
#include <dohpvec.h>
#include <iMesh_extensions.h>

/* Create fs->bmapping and fs->mapping from given layout.
*/
Expand Down Expand Up @@ -49,3 +50,132 @@ dErr dFSCreateLocalToGlobal_Private(dFS fs,dInt n,dInt nc,dInt ngh,dInt *ghidx,d
}
dFunctionReturn(0);
}

/* Partition entities in active set into owned explicit, owned Dirichlet, and ghost */
dErr dFSPopulatePartitionedSets_Private(dFS fs,dMeshAdjacency meshadj)
{
dInt nboundaries,ghstart,ents_a=0,ents_s;
dFSBStatus *bstat;
dMeshESH *bdysets;
dMesh mesh;
iMesh_Instance mi;
dErr err;
dIInt ierr;
dMeshEH *ents;

dFunctionBegin;
err = dFSGetMesh(fs,&mesh);dCHK(err);
err = dMeshGetInstance(mesh,&mi);dCHK(err);
iMesh_addEntArrToSet(mi,meshadj->ents,meshadj->nents,fs->set.explicit,&ierr);dICHK(mi,ierr);
/* Move ghost ents from \a explicitSet to \a ghostSet */
iMesh_getEntitiesRec(mi,fs->set.explicit,dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr);
err = dMeshPartitionOnOwnership(mesh,ents,ents_s,&ghstart);dCHK(err);
iMesh_rmvEntArrFromSet(mi,ents+ghstart,ents_s-ghstart,fs->set.explicit,&ierr);dICHK(mi,ierr);
iMesh_addEntArrToSet(mi,ents+ghstart,ents_s-ghstart,fs->set.ghost,&ierr);dICHK(mi,ierr);
/* Move owned Dirichlet ents from \a explicitSet to \a dirichletSet */
err = dMeshGetNumSubsets(mesh,fs->set.boundaries,0,&nboundaries);dCHK(err);
err = dMallocA2(nboundaries,&bdysets,nboundaries,&bstat);dCHK(err);
err = dMeshGetSubsets(mesh,fs->set.boundaries,0,bdysets,nboundaries,NULL);dCHK(err);
err = dMeshTagSGetData(mesh,fs->tag.bstatus,bdysets,nboundaries,bstat,nboundaries,dDATA_INT);dCHK(err);
for (int i=0; i<nboundaries; i++) {
if (bstat[i] & dFSBSTATUS_DIRICHLET) {
iMesh_getEntitiesRec(mi,bdysets[i],dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr);
err = dMeshPartitionOnOwnership(mesh,ents,ents_s,&ghstart);dCHK(err);
iMesh_rmvEntArrFromSet(mi,ents,ghstart,fs->set.explicit,&ierr);dICHK(mi,ierr);
iMesh_addEntArrToSet(mi,ents,ghstart,fs->set.dirichlet,&ierr);dICHK(mi,ierr);
}
if (bstat[i] & dFSBSTATUS_WEAK) {
iMesh_getEntitiesRec(mi,bdysets[i],dTYPE_FACE,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr);
iMesh_addEntArrToSet(mi,ents,ents_s,fs->set.weakFace,&ierr);dICHK(mi,ierr);
}
}
err = dFree2(bdysets,bstat);dCHK(err);
free(ents);
dFunctionReturn(0);
}

/* Set offsets (global, closure, local) of first node associated with every entity.
*
* @param indexTag index into inodes
* @param inodes number of interior nodes associated with each entity inodes[idx[i]]
* @param rstart relative start for global explicit dofs
* @param crstart relative start for global closure dofs
* @param nents size of array to hold entities
* @param ents returns with entities in the proper ordering of the FS
* @param ghstart returns starting offset of ghost entities in ents
*/
dErr dFSBuildSpaceOffsets_Private(dFS fs,dMeshTag indexTag,const dInt inodes[],dInt rstart,dInt crstart,dInt nents,dMeshEH ents[],dInt *ghstart)
{
dInt i,scan,nentsExplicit,nentsDirichlet,*idx,*offset;
dMesh mesh;
dErr err;

dFunctionBegin;
err = dFSGetMesh(fs,&mesh);dCHK(err);
err = dMallocA2(nents,&offset,nents,&idx);dCHK(err);
/* We assume that orderedSet contains explicitSet+dirichletSet+ghostSet (in that order) */
err = dMeshGetEnts(mesh,fs->set.ordered,dTYPE_ALL,dTOPO_ALL,ents,nents,NULL);dCHK(err);
err = dMeshGetNumEnts(mesh,fs->set.explicit,dTYPE_ALL,dTOPO_ALL,&nentsExplicit);dCHK(err);
err = dMeshGetNumEnts(mesh,fs->set.dirichlet,dTYPE_ALL,dTOPO_ALL,&nentsDirichlet);dCHK(err);
err = dMeshTagGetData(mesh,indexTag,ents,nents,idx,nents,dDATA_INT);dCHK(err);

/* global offset */
for (i=0,scan=rstart; i<nentsExplicit; scan+=inodes[idx[i++]])
offset[i] = scan;
for ( ; i<nents; i++) offset[i] = -1;
err = dMeshTagSetData(mesh,fs->tag.goffset,ents,nents,offset,nents,dDATA_INT);dCHK(err);

/* global closure offset */
for (i=0,scan=crstart; i<nentsExplicit+nentsDirichlet; scan+=inodes[idx[i++]])
offset[i] = scan;
for ( ; i<nents; i++) offset[i] = -1;
err = dMeshTagSetData(mesh,fs->tag.gcoffset,ents,nents,offset,nents,dDATA_INT);dCHK(err);

/* local index */
for (i=0,scan=0; i<nents; scan+=inodes[idx[i++]])
offset[i] = scan;
err = dMeshTagSetData(mesh,fs->tag.loffset,ents,nents,offset,nents,dDATA_INT);dCHK(err);

/* communicate global and closure offset for ghosts */
err = dMeshTagBcast(mesh,fs->tag.goffset);dCHK(err);
err = dMeshTagBcast(mesh,fs->tag.gcoffset);dCHK(err);

err = dFree2(offset,idx);dCHK(err);
*ghstart = nentsExplicit + nentsDirichlet;
dFunctionReturn(0);
}

dErr dFSBuildSpaceVectors_Private(dFS fs,dMeshTag indexTag,const dInt inodes[],dInt rstart,dInt ghents_s,const dMeshEH ghents[])
{
dErr err;
dInt *gcoffset,*ghidx,*idx;
dMesh mesh;

dFunctionBegin;
err = dMallocA3(ghents_s,&gcoffset,ghents_s,&ghidx,ghents_s,&idx);dCHK(err);
err = dFSGetMesh(fs,&mesh);dCHK(err);
err = dMeshTagGetData(mesh,indexTag,ghents,ghents_s,idx,ghents_s,dDATA_INT);dCHK(err);

/* Retrieve ghost offsets, to create localupdate. */
err = dMeshTagGetData(mesh,fs->tag.gcoffset,ghents,ghents_s,gcoffset,ghents_s,dDATA_INT);dCHK(err);
for (dInt i=0; i<ghents_s; i++) { /* Paranoia: confirm that all ghost entities were updated. */
if (gcoffset[i] < 0) dERROR(PETSC_COMM_SELF,1,"Tag exchange did not work");
}

/* Set ghost indices of every node using \a ghidx, create global vector. */
{
dInt gh=0;
for (dInt i=0; i<ghents_s; i++) {
for (dInt j=0; j<inodes[idx[i]]; j++) ghidx[gh++] = gcoffset[i] + j;
}
if (gh != fs->ngh) dERROR(PETSC_COMM_SELF,1,"Ghost count inconsistent");
err = VecCreateDohp(((dObject)fs)->comm,fs->bs,fs->n,fs->nc,fs->ngh,ghidx,&fs->gvec);dCHK(err);
}
err = dFree3(gcoffset,ghidx,idx);dCHK(err);

/* Create fs->bmapping and fs->mapping */
err = dFSCreateLocalToGlobal_Private(fs,fs->n,fs->nc,fs->ngh,ghidx,rstart);dCHK(err);

err = VecDohpCreateDirichletCache(fs->gvec,&fs->dcache,&fs->dscat);dCHK(err);
dFunctionReturn(0);dCHK(err);
}

0 comments on commit e691465

Please sign in to comment.