Skip to content

Commit

Permalink
More work on dFS. The boilerplate should be basically finished.
Browse files Browse the repository at this point in the history
A test is started and is leak-free.

Signed-off-by: Jed Brown <jed@59A2.org>
  • Loading branch information
jedbrown committed Sep 11, 2008
1 parent 21795e3 commit 6e5fe96
Show file tree
Hide file tree
Showing 16 changed files with 415 additions and 67 deletions.
9 changes: 6 additions & 3 deletions include/dohpfs.h
Original file line number Diff line number Diff line change
Expand Up @@ -27,12 +27,15 @@ extern PetscCookie dFS_COOKIE;

EXTERN dErr dFSCreate(MPI_Comm,dFS*);
EXTERN dErr dFSSetMesh(dFS,dMesh,dMeshESH,dMeshTag); /* mesh, active set, partition tag */
EXTERN dErr dFSSetQuotient(dFS,dQuotient); /* must be defined at least on the active set and boundary facets */
EXTERN dErr dFSSetDegree(dFS,dMeshTag,dJacobi);
EXTERN dErr dFSSetRuleTag(dFS,dJacobi,dMeshTag);
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 dFSSetUp(dFS);
EXTERN dErr dFSSetFromOptions(dFS);
EXTERN dErr dFSSetType(dFS,const dFSType);
EXTERN dErr dFSSetFromOptions(dFS);
EXTERN dErr dFSCreateLocalVector(dFS,Vec*);
EXTERN dErr dFSCreateGlobalVector(dFS,Vec*);
EXTERN dErr dFSBuildSpace(dFS);

EXTERN dErr dFSDestroy(dFS);
EXTERN dErr dFSView(dFS,PetscViewer);
Expand Down
4 changes: 2 additions & 2 deletions include/dohpjacobi.h
Original file line number Diff line number Diff line change
Expand Up @@ -75,8 +75,8 @@ EXTERN dErr dJacobiRegisterAll(const char[]);
EXTERN dErr dJacobiInitializePackage(const char[]);

EXTERN dErr dJacobiSetDegrees(dJacobi,dInt,dInt);
EXTERN dErr dJacobiGetRule(dJacobi jac,dTopology top,const dInt rsize[],dRule *rule,void **base,dInt *index);
EXTERN dErr dJacobiGetEFS(dJacobi jac,dTopology top,const dInt bsize[],dRule rule,dEFS *efs,void **base,dInt *index);
EXTERN dErr dJacobiGetRule(dJacobi jac,dEntTopology top,const dInt rsize[],dRule *rule,void **base,dInt *index);
EXTERN dErr dJacobiGetEFS(dJacobi jac,dEntTopology top,const dInt bsize[],dRule rule,dEFS *efs,void **base,dInt *index);

EXTERN dErr dRuleView(dRule rule,PetscViewer);
EXTERN dErr dRuleGetSize(dRule rule,dInt *dim,dInt *nnodes);
Expand Down
5 changes: 5 additions & 0 deletions include/dohpmesh.h
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,13 @@
#include <string.h>
#include "dohpbase.h"
#include "dohptype.h"
#include "dohpjacobi.h"

PETSC_EXTERN_CXX_BEGIN

typedef int dMeshInt;
typedef double dMeshReal;
typedef enum iMesh_EntityTopology dMeshEntTopology;
typedef iBase_EntityHandle dMeshEH;
typedef iBase_TagHandle dMeshTag;
typedef iBase_EntitySetHandle dMeshESH;
Expand Down Expand Up @@ -132,6 +134,9 @@ EXTERN dErr dMeshDestroy(dMesh);
EXTERN dErr dMeshView(dMesh,PetscViewer);
EXTERN dErr dMeshRegisterAll(const char path[]);
EXTERN dErr dMeshGetEntSetName(dMesh m,dMeshESH set,char **str);
EXTERN dErr dMeshCreateRuleTagIsotropic(dMesh,dMeshESH,dJacobi,const char*,dInt,dMeshTag*);
EXTERN dErr dMeshDestroyRuleTag(dMesh,dMeshTag);
EXTERN dErr dMeshGetInstance(dMesh,iMesh_Instance*);

PETSC_EXTERN_CXX_END
#endif
6 changes: 4 additions & 2 deletions include/dohptype.h
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,9 @@ typedef PetscErrorCode dErr;
typedef PetscObject dObject;
typedef PetscViewer dViewer;

#define dTopology enum iMesh_EntityTopology

/* #define dEntTopology enum iMesh_EntityTopology */
typedef enum iMesh_EntityTopology dEntTopology;
typedef enum iBase_EntityType dEntType;

#define dCHK(err) CHKERRQ(err);
#define dERROR(n,...) {return PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,n,1,__VA_ARGS__);}
Expand All @@ -32,6 +33,7 @@ typedef PetscViewer dViewer;
#define dValidPointer(a,b) PetscValidPointer(a,b)
#define dMalloc(a,b) PetscMalloc(a,b)
#define dNew(a,b) PetscNew(a,b)
#define dNewLog(a,b,c) PetscNewLog(a,b,c)
#define dFree(a) PetscFree(a)
#define dNewM(n,t,p) (dMalloc((n)*sizeof(t),(p)) || dMemzero(*(p),(n)*sizeof(t)))
#define dMallocM(n,t,p) (dMalloc((n)*sizeof(t),(p)))
Expand Down
11 changes: 8 additions & 3 deletions include/private/fsimpl.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,17 +20,22 @@ struct dFSBoundary {
struct _dFSOps {
DMOPS(dFS)
dErr (*setfromoptions)(dFS);
dErr (*impldestroy)(dFS);
dErr (*buildspace)(dFS);
};

struct _p_dFS {
PETSCHEADER(struct _dFSOps);
dMesh mesh;
dMeshTag partition,degree;
dMeshTag partition,degree,ruletag;
dMeshESH active;
struct dFSBoundary *bdy_start;
struct dFSBoundary *bdylist;
dQuotient quotient;
dJacobi jacobi;
dBool setupcalled;
dBool spacebuilt;
Sliced sliced;
MeshListEH r,f,e,v; /**< region, face, edge, vertex */
PetscInt n,N; /**< length of the local and global vectors */
void *data;
};

Expand Down
6 changes: 3 additions & 3 deletions include/private/jacimpl.h
Original file line number Diff line number Diff line change
Expand Up @@ -67,9 +67,9 @@ struct _dJacobiOps {
dErr (*view)(dJacobi,PetscViewer);
//dErr (*getrule)(dJacobi,dInt,const dInt[],dRule*,dInt*); /**< put a dRule into the output buffer */
//dErr (*getefs)(dJacobi,dInt,const dInt[],const dInt[],dEFS*,dInt*); /**< put a dEFS into the output buffer */
dErr (*getrulesize)(dJacobi,dTopology,dInt*);
dErr (*getrule)(dJacobi jac,dTopology top,const dInt rsize[],dRule *rule,void **base,dInt *index);
dErr (*getefs)(dJacobi jac,dTopology top,const dInt bsize[],dRule rule,dEFS *efs,void **base,dInt *index);
dErr (*getrulesize)(dJacobi,dEntTopology,dInt*);
dErr (*getrule)(dJacobi jac,dEntTopology top,const dInt rsize[],dRule *rule,void **base,dInt *index);
dErr (*getefs)(dJacobi jac,dEntTopology top,const dInt bsize[],dRule rule,dEFS *efs,void **base,dInt *index);
};

/**
Expand Down
8 changes: 4 additions & 4 deletions sandbox/itaps/dohpblock.c
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ dErr dMeshSetQuotient(dMesh m,dQuotient q)
#define __FUNCT__ "main"
int main(int argc, char *argv[])
{
const char *dflt_outfile="dblock.h5m", *outopts="", *pTagName="MATERIAL_SET";
const char *dflt_outfile="dblock.h5m", *outopts="", *pTagName="dohp_partition";
dInt verbosity = 1;
iMesh_Instance mesh;
iBase_EntitySetHandle root;
Expand Down Expand Up @@ -331,8 +331,8 @@ int main(int argc, char *argv[])

/* Create the boundary parent set. */
{
int on_wall(double x[]) {return (x[0] == x0 || x[0] == x1 || x[1] == y0 || x[1] == y1 || x[2] == z0);}
int on_lid(double x[]) {return (x[2] == z1); }
int on_wall(double xx[]) {return (xx[0] == x0 || xx[0] == x1 || xx[1] == y0 || xx[1] == y1 || xx[2] == z0);}
int on_lid(double xx[]) {return (xx[2] == z1); }
typedef int OnBdyFunc(double[]);
const int nbdy = 3;
const char bdyName[3][dNAME_LEN] = {dBDY_ROOT,"WALL","LID"};
Expand Down Expand Up @@ -441,7 +441,7 @@ int main(int argc, char *argv[])
iMesh_getEntities(mesh,bdy.v[i],iBase_ALL_TYPES,iMesh_ALL_TOPOLOGIES,&e.v,&e.a,&e.s,&err);ICHKERRQ(mesh,err);
iMesh_getEntArrTopo(mesh,e.v,e.s,&type.v,&type.a,&type.s,&err);ICHKERRQ(mesh,err);
for (j=0; j<e.s; j++) {
if (type.v[j] > iBase_VERTEX || 1) {
if (type.v[j] > iBase_EDGE || 1) { /* should we set the boundary on all entities or just faces? */
iMesh_setIntData(mesh,e.v[j],bdyNum,i,&err);ICHKERRQ(mesh,err);
}
}
Expand Down
88 changes: 84 additions & 4 deletions src/fs/impls/cont/cont.c
Original file line number Diff line number Diff line change
@@ -1,25 +1,105 @@
#include "cont.h"


static dErr dFSView_Cont(dFS fs,dViewer viewer)
{
struct dFS_Cont *fsc = fs->data;
dBool ascii;
dErr err;

dFunctionBegin;
err = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&ascii);dCHK(err);
if (ascii) {
err = PetscViewerASCIIPrintf(viewer,"dFSView_Cont() nothing here yet\n");dCHK(err);
err = PetscViewerASCIIPrintf(viewer,"Continuous Galerkin function space\n");dCHK(err);
{ /* print aggregate sizes */
PetscMPIInt gm[2],lm[2];
lm[0] = fsc->m; lm[1] = fs->n; /* set local `element' size and `local' size */
err = MPI_Allreduce(lm,gm,2,MPI_INT,MPI_SUM,((dObject)fs)->comm);dCHK(err);
err = PetscViewerASCIIPrintf(viewer,"%d element dofs constrained against %d local dofs assembled to %d global dofs\n",gm[0],gm[1],fs->N);dCHK(err);
}
}
dFunctionReturn(0);
}

/**
* Calculate the sizes of the global and local vectors, create scatter contexts. Assemble the constraint matrix for
* element->global maps.
*
* @param fs
*
* @return
*/
static dErr dFSSetFromOptions_Cont(dFS fs)
{
dMesh mesh = fs->mesh;
struct dFS_Cont *fsc = fs->data;
dBool flg;
dErr err;

dFunctionBegin;
err = PetscOptionsHead("Continuous Galerkin options");dCHK(err);
{
err = PetscOptionsName("-dfs_cont_constraint_matrix","use explicit SeqAIJ constraint matrix for constraints","None",&flg);dCHK(err);
if (flg) { fsc->usecmatrix = true; }
}
err = PetscOptionsTail();dCHK(err);
dFunctionReturn(0);
}

static dErr dFSDestroy_Cont(dFS fs)
{
dErr err;

dFunctionBegin;
err = dFree(fs->data);dCHK(err);
dFunctionReturn(0);
}

static dErr dFSContPropogateDegree(dFS fs)
{
iMesh_Instance mi;
MeshListEH rf=MLZ;
MeshListInt rfo=MLZ;
dErr err;

dFunctionBegin;
err = dMeshGetInstance(fs->mesh,&mi);dCHK(err);
iMesh_getEntArrAdj(mi,fs->r.v,fs->r.s,iBase_FACE,&rf.v,&rf.a,&rf.s,&rfo.v,&rfo.a,&rfo.s,&err);dICHK(mi,err);
/* orient faces with respect to regions, propogate anisotropic degree */
#warning much to do here
MeshListFree(rf); MeshListFree(rfo);
dFunctionReturn(0);
}

static dErr dFSBuildSpace_Cont(dFS fs)
{
dErr err;

dFunctionBegin;
err = dFSContPropogateDegree(fs);dCHK(err);
dFunctionReturn(0);
}

/**
* Create the private structure used by a continuous galerkin function space.
*
* This function does not allocate the constraint matrices.
*
* @param fs the function space
*
* @return err
*/
dErr dFSCreate_Cont(dFS fs)
{
struct dFS_Cont *fsc;
dErr err;

dFunctionBegin;
err = dMalloc(sizeof(dFS_Cont),&fs->data);dCHK(err);
err = dMemzero(fs->data,sizeof(dFS_Cont));dCHK(err);
fs->ops->view = &dFSView_Cont;
err = dNewLog(fs,*fsc,&fsc);dCHK(err);
fs->data = (void*)fsc;
fs->ops->view = dFSView_Cont;
fs->ops->impldestroy = dFSDestroy_Cont;
fs->ops->setfromoptions = dFSSetFromOptions_Cont;
fs->ops->buildspace = dFSBuildSpace_Cont;
dFunctionReturn(0);
}
12 changes: 8 additions & 4 deletions src/fs/impls/cont/cont.h
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,15 @@

PETSC_EXTERN_CXX_BEGIN

dErr dFSCreate_Cont(dFS);
extern dErr dFSCreate_Cont(dFS);

typedef struct {
MeshListEH r,f,e,v;
} dFS_Cont;
struct dFS_Cont {
dBool usecmatrix;
dInt m;
Mat F; /**< facet constraint matrix */
Mat Cprimal; /**< primal constraint matrix (element dofs to local numbering, full order) */
Mat Cdual; /**< dual constraint matrix (element dofs to local numbering, as sparse as possible) */
};

PETSC_EXTERN_CXX_END

Expand Down

0 comments on commit 6e5fe96

Please sign in to comment.