Skip to content

Commit

Permalink
Bring the common parts of dFSBuildSpace_Cont back in one place instea…
Browse files Browse the repository at this point in the history 
…d fo duplicating in dFSLoadIntoFS_Cont_DHM
  • Loading branch information
@jedbrown
jedbrown committed Apr 23, 2011
1 parent 2429f03 commit 731b053
Show file tree
Hide file tree
Showing 3 changed files with 39 additions and 89 deletions.
1 change: 1 addition & 0 deletions include/dohpfsimpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,6 +151,7 @@ extern dErr dFSCreateLocalToGlobal_Private(dFS fs,dInt n,dInt nc,dInt ngh,dInt *
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 dFSBuildSpaceWithOrderedSet_Private(dFS fs,dMeshAdjacency meshAdj);
extern dErr DMDestroy_dFS(DM);

dEXTERN_C_END
Expand Down
59 changes: 36 additions & 23 deletions src/fs/impls/cont/cont.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -257,34 +257,17 @@ dErr dFSBuildSpace_Cont_CreateElemAssemblyMats(dFS fs,const dInt idx[],const dMe
*/
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 */
dUNUSED dFS_Cont *cont = fs->data;
struct _p_dMeshAdjacency ma;
dMeshAdjacency meshAdj;
dMesh mesh;
iMesh_Instance mi;
dEntTopology *regTopo;
dPolynomialOrder *deg,*regBDeg;
dInt *inodes,*xnodes,nregions,ents_a,ents_s,ghents_s,*idx;
dInt *xstart,xcnt;
dInt rstart,crstart;
dMeshEH *ents,*ghents;
dErr err;
dErr err;
dMeshAdjacency meshAdj;
dMesh mesh;

dFunctionBegin;
dValidHeader(fs,DM_CLASSID,1);
mesh = fs->mesh;
err = dMeshGetInstance(mesh,&mi);dCHK(err);
err = dFSGetMesh(fs,&mesh);dCHK(err);
err = dMeshGetAdjacency(mesh,fs->set.active,&meshAdj);dCHK(err);
err = dMemcpy(&ma,meshAdj,sizeof ma);dCHK(err); /* To have object rather than pointer semantics in this function. */
err = dFSContPropagateDegree(fs,meshAdj);dCHK(err);

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

err = dFSPopulatePartitionedSets_Private(fs,&ma);dCHK(err);
err = dFSPopulatePartitionedSets_Private(fs,meshAdj);dCHK(err);

err = dMeshPopulateOrderedSet_Private(mesh,fs->set.ordered,fs->set.explicit,fs->set.dirichlet,fs->set.ghost,fs->orderingtype);dCHK(err);
{
Expand All @@ -302,6 +285,36 @@ dErr dFSBuildSpace_Cont(dFS fs)
err = dMeshTagCreate(mesh,strbuf,1,dDATA_INT,&fs->tag.orderedsub);dCHK(err);
err = dMeshTagSSetData(mesh,fs->tag.orderedsub,&fs->set.ordered,1,&rank,1,dDATA_INT);dCHK(err);
}
err = dFSBuildSpaceWithOrderedSet_Private(fs,meshAdj);dCHK(err);
err = dMeshRestoreAdjacency(fs->mesh,fs->set.active,&meshAdj);dCHK(err);
dFunctionReturn(0);
}

dErr dFSBuildSpaceWithOrderedSet_Private(dFS fs,dMeshAdjacency meshAdj)
{
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 */
dUNUSED dFS_Cont *cont = fs->data;
struct _p_dMeshAdjacency ma;
dMesh mesh;
iMesh_Instance mi;
dEntTopology *regTopo;
dPolynomialOrder *deg,*regBDeg;
dInt *inodes,*xnodes,nregions,ents_a,ents_s,ghents_s,*idx;
dInt *xstart,xcnt;
dInt rstart,crstart;
dMeshEH *ents,*ghents;
dErr err;

dFunctionBegin;
mesh = fs->mesh;
err = dMeshGetInstance(mesh,&mi);dCHK(err);
err = dMemcpy(&ma,meshAdj,sizeof ma);dCHK(err); /* To have object rather than pointer semantics in this function. */

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

err = dMeshGetEnts(fs->mesh,fs->set.ordered,dTYPE_ALL,dTOPO_ALL,ents,ents_a,&ents_s);dCHK(err);
if (ents_s != ents_a) dERROR(PETSC_COMM_SELF,PETSC_ERR_PLIB,"wrong set size");

Expand Down Expand Up @@ -366,8 +379,8 @@ 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);
err = dFree4(ents,idx,deg,inodes);dCHK(err);
fs->spacebuilt = dTRUE;
dFunctionReturn(0);
}

Expand Down
68 changes: 2 additions & 66 deletions src/fs/impls/cont/contviewdhm.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -366,79 +366,15 @@ dErr dFSLoadIntoFS_Cont_DHM(PetscViewer viewer,const char name[],dFS fs)
herr = H5Dclose(meshobj);dH5CHK(herr,H5Aclose);
// \bug herr = H5Dvlen_reclaim(&fs5);
}
/** @note The FS has the layout, ordering, and boundary status tags set so we are ready to build the function space.
**/

/* @todo Most of this block mirrors dFSBuildSpace_Cont, the common part should be better contained. */
{
dInt ents_a,ents_s,nregions,*inodes,*xnodes,*xstart,*idx,xcnt,rstart,crstart,ghents_s;
dPolynomialOrder *bdeg,*regBDeg;
dEntTopology *regTopo;
dMeshEH *ents,*ghents;
{ // The FS has the layout, ordering, and boundary status tags set so we are ready to build the function space.
dMesh mesh;
dMeshAdjacency meshadj;

err = dFSGetMesh(fs,&mesh);dCHK(err);
err = dMeshGetAdjacency(mesh,fs->set.ordered,&meshadj);dCHK(err);

ents_a = meshadj->nents;
err = dMallocA4(ents_a,&ents,ents_a,&idx,ents_a,&bdeg,ents_a,&inodes);dCHK(err);

err = dFSPopulatePartitionedSets_Private(fs,meshadj);dCHK(err);

/* Ordered set is already populated so no need to call dMeshPopulateOrderedSet_Private() and assign tags afterward.
* This is the key difference between the current code block and dFSBuildSpace_Cont(). */

err = dMeshGetEnts(mesh,fs->set.ordered,dTYPE_ALL,dTOPO_ALL,ents,ents_a,&ents_s);dCHK(err);
if (ents_s != ents_a) dERROR(PETSC_COMM_SELF,PETSC_ERR_PLIB,"wrong set size");

err = dMeshTagGetData(mesh,fs->tag.degree,meshadj->ents,meshadj->nents,bdeg,ents_s,dDATA_INT);dCHK(err);
err = dJacobiGetNodeCount(fs->jacobi,ents_s,meshadj->topo,bdeg,inodes,NULL);dCHK(err);

{
dInt counts[3],rstarts[3];
err = dMeshClassifyCountInt(mesh,meshadj->nents,meshadj->ents,inodes,3,(const dMeshESH[]){fs->set.explicit,fs->set.dirichlet,fs->set.ghost},counts);dCHK(err);
err = MPI_Scan(counts,rstarts,3,MPIU_INT,MPI_SUM,((dObject)fs)->comm);dCHK(err);
for (dInt i=0; i<3; i++) rstarts[i] -= counts[i];
fs->n = counts[0];
fs->nc = counts[0] + counts[1];
fs->ngh = counts[2];
rstart = rstarts[0];
crstart = rstarts[0] + rstarts[1];
}

{
dInt ghstart;
err = dFSBuildSpaceOffsets_Private(fs,meshadj->indexTag,inodes,rstart,crstart,ents_s,ents,&ghstart);dCHK(err);
ghents = ents + ghstart;
ghents_s = ents_s - ghstart;
}

err = dFSBuildSpaceVectors_Private(fs,meshadj->indexTag,inodes,rstart,ghents_s,ghents);dCHK(err);

err = dMeshGetEnts(mesh,fs->set.active,dTYPE_REGION,dTOPO_ALL,ents,ents_a,&ents_s);dCHK(err);
err = dMeshTagGetData(mesh,meshadj->indexTag,ents,ents_s,idx,ents_s,dDATA_INT);dCHK(err);
nregions = ents_s;
err = dMallocA4(nregions+1,&xstart,nregions,&regTopo,nregions,&regBDeg,nregions,&xnodes);dCHK(err);
err = dMeshGetTopo(mesh,ents_s,ents,regTopo);dCHK(err);
err = dMeshTagGetData(mesh,fs->tag.degree,ents,ents_s,regBDeg,nregions,dDATA_INT);dCHK(err);
err = dJacobiGetNodeCount(fs->jacobi,ents_s,regTopo,regBDeg,NULL,xnodes);dCHK(err);

xcnt = xstart[0] = 0;
for (dInt i=0; i<nregions; i++) xstart[i+1] = (xcnt += xnodes[i]);

fs->nelem = nregions;
err = dMallocA(nregions+1,&fs->off);dCHK(err); /* Will be freed by FS */
err = dMemcpy(fs->off,xstart,(nregions+1)*sizeof(xstart[0]));dCHK(err);

err = dFree4(xstart,regTopo,regBDeg,xnodes);dCHK(err);

err = dMeshTagGetData(mesh,fs->tag.degree,meshadj->ents,meshadj->nents,bdeg,meshadj->nents,dDATA_INT);dCHK(err);
err = dFSBuildSpace_Cont_CreateElemAssemblyMats(fs,idx,meshadj,bdeg,&fs->E,&fs->Ep);dCHK(err);

err = dFSBuildSpaceWithOrderedSet_Private(fs,meshadj);dCHK(err);
err = dMeshRestoreAdjacency(mesh,fs->set.ordered,&meshadj);dCHK(err);
err = dFree4(ents,idx,bdeg,inodes);dCHK(err);
fs->spacebuilt = dTRUE;
}

herr = H5Sclose(fsspace);dH5CHK(herr,H5Sclose);
Expand Down

0 comments on commit 731b053

Please sign in to comment.