Skip to content

Commit

Permalink
Work on handling boundary conditions. [#5] [#10]
Browse files Browse the repository at this point in the history 
* Provide Global, Dirichlet, and Closure vectors.  Global and Dirichlet
  vectors are VecGhost so they have a local form.  Scatters are provided from
  Closure to the others.

* Use MAIJ format for element assembly matrices (allows vector-valued spaces)

* Keep track of boundaries by using ITAPS/MOAB tag conventions instead of
  private data structure (this is useful since the user needs to see the set
  structure if they need to relate to geometric entities).

* Add dFSHomogeneousMode to handle which boundary values should be mapped into
  expanded space.
  • Loading branch information
@jedbrown
jedbrown committed Apr 17, 2009
1 parent 03bbebf commit 69efc49
Show file tree
Hide file tree
Showing 14 changed files with 675 additions and 292 deletions.
24 changes: 17 additions & 7 deletions include/dohpfs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,6 @@ PETSC_EXTERN_CXX_BEGIN

typedef struct _p_dFS *dFS;

typedef struct _p_dFSBoundary *dFSBoundary;

/** User-provided constraint function.
* @param ctx User context
* @param x coordinates of node to generate constraints for (vector of length 3)
Expand All @@ -33,7 +31,16 @@ typedef struct _p_dFSBoundary *dFSBoundary;
* constant for each boundary type). Rationale: setting the number of global dofs separately makes it possible for the
* constraint function to become out of sync with the number of global dofs.
**/
typedef dErr (*dFSBoundaryConstraintFunction)(void *ctx,const dReal x[],const dReal b[],dReal T[],dInt *g);
typedef dErr (*dFSConstraintFunction)(void *ctx,const dReal x[],const dReal b[],dReal T[],dInt *g);

enum { dFSBSTATUS_MASK = 0xffff }; /* Reserved for number of strongly enforced dofs */
typedef enum {
dFSBSTATUS_INTERIOR = 0,
dFSBSTATUS_DIRICHLET = (dFSBSTATUS_MASK+1) << 0, /* Do not represent in global space */
dFSBSTATUS_WEAK = (dFSBSTATUS_MASK+1) << 1, /* A weak form must be evaluated on this element */
} dFSBStatus;

typedef enum {dFS_HOMOGENEOUS, dFS_INHOMOGENEOUS} dFSHomogeneousMode;

#define dFSType char *

Expand All @@ -44,24 +51,27 @@ EXTERN dErr dFSSetMesh(dFS,dMesh,dMeshESH); /* mesh, active set */
EXTERN dErr dFSSetRuleTag(dFS,dJacobi,dMeshTag);
EXTERN dErr dFSSetDegree(dFS,dJacobi,dMeshTag);
EXTERN dErr dFSSetBlockSize(dFS,dInt);
EXTERN dErr dFSRegisterBoundary(dFS,dInt,dTruth,dTruth,dFSBoundaryConstraintFunction,void*);
EXTERN dErr dFSRegisterBoundary(dFS,dInt,dFSBStatus,dFSConstraintFunction,void*);
EXTERN dErr dFSSetFromOptions(dFS);
EXTERN dErr dFSSetType(dFS,const dFSType);
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 dFSCreateDirichletVector(dFS,Vec*);
EXTERN dErr dFSGlobalToExpandedBegin(dFS,Vec,dFSHomogeneousMode,Vec);
EXTERN dErr dFSGlobalToExpandedEnd(dFS,Vec,dFSHomogeneousMode,Vec);
EXTERN dErr dFSExpandedToGlobal(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSExpandedToGlobalBegin(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSExpandedToGlobalEnd(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSExpandedToDirichlet(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSRotateLocalVector(dFS,Vec xloc,Vec save_strong,Vec correct_strong);
EXTERN dErr dFSUnRotateLocalVector(dFS,Vec,Vec,Vec);
EXTERN dErr dFSGetElements(dFS,dInt*,dInt*restrict*,s_dRule*restrict*,s_dEFS*restrict*,dInt*restrict*,dReal(*restrict*)[3]);
EXTERN dErr dFSRestoreElements(dFS,dInt*,dInt*restrict*,s_dRule*restrict*,s_dEFS*restrict*,dInt*restrict*,dReal(*restrict*)[3]);
EXTERN dErr dFSGetWorkspace(dFS,const char[],dReal(*restrict*)[3],dReal(*restrict*)[3][3],dReal*restrict*,dScalar*restrict*,dScalar*restrict*,dScalar*restrict*,dScalar*restrict*);
EXTERN dErr dFSRestoreWorkspace(dFS,const char[],dReal(*restrict*)[3],dReal(*restrict*)[3][3],dReal*restrict*,dScalar*restrict*,dScalar*restrict*,dScalar*restrict*,dScalar*restrict*);
EXTERN dErr dFSMatSetValuesExpanded(dFS,Mat,dInt,const dInt[],dInt,const dInt[],const dScalar[],InsertMode);
EXTERN dErr dFSGetMatrix(dFS,const MatType,Mat*);
EXTERN dErr dFSBuildSpace(dFS);
EXTERN dErr dFSCreateSubspace(dFS,dInt,dFSBoundary,dFS*);

EXTERN dErr dFSDestroy(dFS);
EXTERN dErr dFSView(dFS,PetscViewer);
Expand Down
31 changes: 3 additions & 28 deletions include/dohpjacobi.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,36 +91,11 @@ typedef enum {
*/
typedef struct p_dJacobi *dJacobi;

/**
* This struct is public since it is used to exchange mesh data between the dFS and dJacobi packages.
*
* @example Boundary specs for some boundary conditions:
*
* Dirichlet: dofspernode=0, cperdof does not matter
* Neumann: dofspernode=D, cperdof=1
* Slip: dofspernode=2, cperdof=3
* Normal: dofspernode=1, cperdof=3
* Interior entities: dofspernode=D, cperdof=1
*
*/
struct MeshRegion {
dMeshEH *conn; /**< conn[connoff[i]+j] = j'th vertex of entity [i] */
dInt *connoff;
dInt *adj; /**< adj[adjoff[i]+j] = index of the j'th entity adjacent to entity i */
dInt *adjoff;
dEntTopology *topo;
dEntStatus *status; /**< parallel status */
dInt *deg; /**< deg[3*i..3*i+2] */
dInt *dofspernode; /**< number of dofs per node (0..D), always D for interior nodes */
dInt *cperdof; /**< number of contrstrainst per dof */
dInt nents;
dInt ntype[4],tstart[4];
};

#define dMESHADJACENCY_HAS_CONNECTIVITY 1
typedef struct dMeshAdjacency *dMeshAdjacency;
struct dMeshAdjacency {
dMeshESH set;
dMeshTag indexTag;
dInt nents;
dInt toff[5];
dInt *adjoff,*adjind,*adjperm;
Expand Down Expand Up @@ -165,8 +140,8 @@ EXTERN dErr dEFSApply(dEFS,const dReal[],dInt,const dScalar[],dScalar[restrict],
EXTERN dErr dJacobiPropogateDown(dJacobi,const struct dMeshAdjacency*,dInt[]);
EXTERN dErr dJacobiGetNodeCount(dJacobi,dInt,const dEntTopology[],const dInt[],dInt[],dInt[]);

EXTERN dErr dJacobiGetConstraintCount(dJacobi,dInt,const dInt[],const dInt[],const dInt[],const struct dMeshAdjacency*,dInt[],dInt[]);
EXTERN dErr dJacobiAddConstraints(dJacobi,dInt,const dInt[],const dInt[],const dInt[],const dInt[],const struct dMeshAdjacency*,Mat,Mat);
EXTERN dErr dJacobiGetConstraintCount(dJacobi,dInt,const dInt[],const dInt[],dInt,const dInt[],const dInt[],const struct dMeshAdjacency*,dInt[],dInt[],dInt[]);
EXTERN dErr dJacobiAddConstraints(dJacobi,dInt,const dInt[],const dInt[],dInt,const dInt[],const dInt[],const struct dMeshAdjacency*,Mat,Mat,Mat);

PETSC_EXTERN_CXX_END

Expand Down
12 changes: 6 additions & 6 deletions include/dohpmesh.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,8 @@ EXTERN dErr dMeshDestroyRuleTag(dMesh,dMeshTag);
EXTERN dErr dMeshGetInstance(dMesh,iMesh_Instance*);
EXTERN dErr dMeshGetNumEnts(dMesh,dMeshESH,dEntType,dEntTopology,dInt*);
EXTERN dErr dMeshGetEnts(dMesh,dMeshESH,dEntType,dEntTopology,dMeshEH[],dInt,dInt*);
EXTERN dErr dMeshGetNumSubsets(dMesh,dMeshESH,dInt,dInt*);
EXTERN dErr dMeshGetSubsets(dMesh,dMeshESH,dInt,dMeshESH[],dInt,dInt*);
EXTERN dErr dMeshGetEntsOff(dMesh,dMeshESH,dInt*,dMeshEH**);
EXTERN dErr dMeshGetAdjIndex(dMesh,const dMeshEH[],dInt,const dMeshEH[],dInt,dInt[],dInt*);

Expand All @@ -124,21 +126,19 @@ EXTERN dErr dMeshTagCreate(dMesh mesh,const char[],dInt count,dDataType type,dMe
EXTERN dErr dMeshTagCreateTemp(dMesh mesh,const char[],dInt count,dDataType type,dMeshTag *intag);
EXTERN dErr dMeshTagSetData(dMesh mesh,dMeshTag tag,const dMeshEH ents[],dInt ecount,const void *data,dInt count,dDataType type);
EXTERN dErr dMeshTagGetData(dMesh mesh,dMeshTag tag,const dMeshEH ents[],dInt ecount,void *data,dInt count,dDataType type);
EXTERN dErr dMeshTagSSetData(dMesh mesh,dMeshTag tag,const dMeshESH esets[],dInt ecount,const void *data,dInt count,dDataType type);
EXTERN dErr dMeshTagSGetData(dMesh mesh,dMeshTag tag,const dMeshESH esets[],dInt ecount,void *data,dInt count,dDataType type);
EXTERN dErr dMeshGetTaggedSet(dMesh,dMeshTag,const void*,dMeshESH*);
EXTERN dErr dMeshSetFromOptions(dMesh);
EXTERN dErr dMeshTagBcast(dMesh mesh,dMeshTag tag);
EXTERN dErr dMeshGetStatus(dMesh,dInt,const dMeshEH[],dEntStatus[]);
EXTERN dErr dMeshSetCreate(dMesh,dMeshESH*);
EXTERN dErr dMeshGetStatus(dMesh,const dMeshEH[],dInt,dEntStatus[]);
EXTERN dErr dMeshGetTopo(dMesh,dInt,const dMeshEH[],dEntTopology[]);
EXTERN dErr dMeshGetAdjacency(dMesh,dMeshESH,dMeshAdjacency*);
EXTERN dErr dMeshRestoreAdjacency(dMesh,dMeshESH,dMeshAdjacency*);
EXTERN dErr dMeshGetVertexCoords(dMesh,dInt,const dMeshEH[],dInt**,dReal(**)[3]);
EXTERN dErr dMeshRestoreVertexCoords(dMesh,dInt,const dMeshEH[],dInt**,dReal(**)[3]);

typedef struct _p_dMeshManifold *dMeshManifold;

EXTERN dErr dMeshGetNumSubsets(dMesh,dMeshESH,dInt*);
EXTERN dErr dMeshGetSubsets(dMesh,dMeshESH,dMeshESH[],dInt,dInt*);
EXTERN dErr dMeshPartitionOnOwnership(dMesh,dMeshEH[],dInt,dInt*);

PETSC_EXTERN_CXX_END
#endif
8 changes: 4 additions & 4 deletions include/dohptype.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@
#define _DOHPTYPE_H

#include "dohpconfig.h"
#include "petsc.h"
#include "iMesh.h"
#include <petsc.h>
#include <iMesh_extensions.h>

/**
* These types all have to be exactly the Petsc versions. These typedefs are here just to shorten the names, not to
Expand Down Expand Up @@ -50,8 +50,8 @@ typedef unsigned char dEntStatus;
#define dEntTopoToITAPS(dtopo,itopo) (*(itopo) = (dtopo), 0)
#define dEntTypeToITAPS(dtype,itype) (*(itype) = (dtype), 0)

#define dCHK(err) if (err) {return PetscError(__LINE__,__func__,__FILE__,__SDIR__,(err),0," ");}
#define dERROR(n,...) return PetscError(__LINE__,__func__,__FILE__,__SDIR__,n,1,__VA_ARGS__)
#define dCHK(err) if (err) return PetscError(__LINE__,__func__,__FILE__,__SDIR__,(err),0," ")
#define dERROR(n,...) return PetscError(__LINE__,__func__,__FILE__,__SDIR__,(n),1,__VA_ARGS__)

#define dPrintf PetscPrintf
#define dMemcpy(a,b,c) PetscMemcpy(a,b,c)
Expand Down
4 changes: 1 addition & 3 deletions include/private/dmeshimpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,14 +17,12 @@ struct _p_dMesh {
char *infile,*inoptions;
iMesh_Instance mi;
dMeshESH root,emptyset;
dMeshTag senseTag;
MeshListEH v,e,f,r; /* vertices, edges, faces, vertices */
MeshListEH arf,afe,aev; /* adjacencies region -> face -> edge -> vertex */
MeshListData orf,ofe;
MeshListReal x;
void *data;
dMeshTag manifoldTag,manifoldOrientTag,senseTag;
dMeshManifold *manifoldList;
dInt nManifolds;
};

#endif
52 changes: 34 additions & 18 deletions include/private/fsimpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,12 +8,19 @@ PETSC_EXTERN_CXX_BEGIN

extern PetscLogEvent dLOG_Q1HexComputeQuadrature,dLOG_FSMatSetValuesExpanded;

struct _p_dFSBoundary {
dMeshESH manset;
dTruth strong,expanded;
dFSBoundaryConstraintFunction cfunc;
void *user;
struct _p_dFSBoundary *next;
static inline dInt dFSBStatusStrongCount(dFSBStatus stat) {
return stat & dFSBSTATUS_MASK;
}
static inline dFSBStatus dFSBStatusSetStrongCount(dFSBStatus stat,dInt count) {
return (stat & ~dFSBSTATUS_MASK) /* clear lower bits */ & count;
}
static inline dTruth dFSBStatusValid(dFSBStatus stat) {
return !((stat & dFSBSTATUS_DIRICHLET) && ((stat & dFSBSTATUS_WEAK) || dFSBStatusStrongCount(stat))); /* cannot be Dirichlet and anything else */
}

struct dFSConstraintCtx {
dFSConstraintFunction cfunc;
void *user;
};

typedef struct {
Expand All @@ -40,33 +47,42 @@ struct _p_dFS {
dMesh mesh;
dMeshAdjacency meshAdj;
dMeshTag degreetag,ruletag; /**< tags on regions */
dMeshTag fieldSpecTag; /**< 2-int tag on every entity in active set */
dMeshESH active; /**< regions that will be part of this space */
dMeshESH activeSet; /**< all entities that will be part of this space, weak forms are evaluated on regions in this set */
dQuotient quotient;
dJacobi jacobi;
dFSBoundary bdyStart; /**< Start of linked list of boundary conditions to enforce */
dMeshESH boundaries; /**< Set of boundary conditions to be enforced on this function space */
dMeshTag bdyTag; /**< Tag for each boundary condition */
char bdyTagName[dNAME_LEN]; /**< Usually "NEUMANN_SET" */
dMeshTag bstatusTag; /**< Boundary status tag, every NEUMANN_SET=x set will be tagged */
dMeshTag bdyConstraintTag; /**< User-defined context for enforcing boundary constraints */
dTruth spacebuilt;
dTruth assemblefull; /**< Use full order constraints for assembly */
dInt ruleStrength;
Sliced sliced;
dInt nbdofs;
dInt n,N; /**< length of the owned and global vectors */
dInt nlocal; /**< number of owned+ghost dofs on this process */
dInt rstart; /**< global offset of first owned dof */
Sliced slice; /**< Lower level DM object, manages the global/local aspects with no additional structure */
Sliced dslice; /**< Manages the Dirichlet local/global Dirichlet vectors */
dInt bs; /**< Block size (number of dofs per node) */
dMeshEH *ents; /**< All entities in active set */
dInt m; /**< Number of expanded dofs */
dInt D; /**< Number of dofs per (non-boundary) node */
dInt nelem;
dInt *off; /**< Offset of element dofs in expanded vector */
s_dRule *rule; /**< Integration rule */
s_dEFS *efs; /**< Element function space, defined for all entities */
dInt *vtxoff;
dReal (*vtx)[3];
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) */
dMeshTag goffsetTag; /**< Offset of first node in global vector */
dMeshTag gdoffsetTag; /**< Offset of first node in global Dirichlet vector */
dMeshTag gcoffsetTag; /**< Offset of first node in global closure vector */
dMeshTag loffsetTag; /**< Offset of first node in local vector (split to include both real and Dirichlet vectors, based on dsplit) */
dMeshESH ownedExplicitSet,ghostExplicitSet; /**< Set of all entities that should be represented explicitly */
dMeshESH ownedDirichletSet,ghostDirichletSet; /**< Set of all entities that have full Dirichlet conditions (removed from global system) */
dMeshESH weakFaceSet; /**< Faces on which weak forms need to be evaluated (I'm not sure this is actually needed) */
VecScatter ctod; /**< Scatter from global closure (includes Dirichlet conditions) to local vectors with Dirichlet values */
VecScatter ctog; /**< Scatter from global closure to global vector */
Mat E; /**< full-order element assembly matrix (element nodes to local numbering) */
Mat Ep; /**< preconditioning element assembly matrix (element nodes to local numbering, as sparse as possible) */
Mat Ed; /**< element assembly matrix for Dirichlet nodes */
Vec weight; /**< Vector in global space, used to compensate for overcounting after local to global */
Vec gc; /**< Global closure vector */
Vec d,dl; /**< Global and local Dirichlet vectors */
dInt maxQ;
s_dFSWorkspace workspace[dFS_MAX_WORKSPACES];
void *data;
Expand Down
4 changes: 2 additions & 2 deletions include/private/jacimpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,8 +55,8 @@ struct _dJacobiOps {
dErr (*GetRule)(dJacobi,dInt,const dEntTopology[],const dInt[],dRule);
dErr (*GetEFS)(dJacobi,dInt,const dEntTopology[],const dInt[],dRule,dEFS);
dErr (*GetNodeCount)(dJacobi,dInt,const dEntTopology[],const dInt[],dInt[],dInt[]);
dErr (*GetConstraintCount)(dJacobi,dInt,const dInt[],const dInt[],const dInt[],const struct dMeshAdjacency*,dInt[],dInt[]);
dErr (*AddConstraints)(dJacobi,dInt,const dInt[],const dInt[],const dInt[],const dInt[],const struct dMeshAdjacency*,Mat,Mat);
dErr (*GetConstraintCount)(dJacobi,dInt,const dInt[],const dInt[],dInt,const dInt[],const dInt[],const struct dMeshAdjacency*,dInt[],dInt[],dInt[]);
dErr (*AddConstraints)(dJacobi,dInt,const dInt[],const dInt[],dInt,const dInt[],const dInt[],const struct dMeshAdjacency*,Mat,Mat,Mat);
};

/**
Expand Down

0 comments on commit 69efc49

Please sign in to comment.