diff --git a/include/dohpfs.h b/include/dohpfs.h index a869505f..df211642 100644 --- a/include/dohpfs.h +++ b/include/dohpfs.h @@ -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) @@ -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 * @@ -44,16 +51,20 @@ 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*); @@ -61,7 +72,6 @@ EXTERN dErr dFSRestoreWorkspace(dFS,const char[],dReal(*restrict*)[3],dReal(*res 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); diff --git a/include/dohpjacobi.h b/include/dohpjacobi.h index 3de2a016..8947369c 100644 --- a/include/dohpjacobi.h +++ b/include/dohpjacobi.h @@ -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; @@ -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 diff --git a/include/dohpmesh.h b/include/dohpmesh.h index b27bc0be..80a0f8a4 100644 --- a/include/dohpmesh.h +++ b/include/dohpmesh.h @@ -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*); @@ -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 diff --git a/include/dohptype.h b/include/dohptype.h index b666d6a9..2a3f2c56 100644 --- a/include/dohptype.h +++ b/include/dohptype.h @@ -2,8 +2,8 @@ #define _DOHPTYPE_H #include "dohpconfig.h" -#include "petsc.h" -#include "iMesh.h" +#include +#include /** * These types all have to be exactly the Petsc versions. These typedefs are here just to shorten the names, not to @@ -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) diff --git a/include/private/dmeshimpl.h b/include/private/dmeshimpl.h index d9a7bda2..1550e1cd 100644 --- a/include/private/dmeshimpl.h +++ b/include/private/dmeshimpl.h @@ -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 diff --git a/include/private/fsimpl.h b/include/private/fsimpl.h index 81f86785..1a388572 100644 --- a/include/private/fsimpl.h +++ b/include/private/fsimpl.h @@ -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 { @@ -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; diff --git a/include/private/jacimpl.h b/include/private/jacimpl.h index 8da14d7c..ece3e42a 100644 --- a/include/private/jacimpl.h +++ b/include/private/jacimpl.h @@ -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); }; /** diff --git a/src/fs/impls/cont/cont.c b/src/fs/impls/cont/cont.c index a9d356cd..704d6068 100644 --- a/src/fs/impls/cont/cont.c +++ b/src/fs/impls/cont/cont.c @@ -70,143 +70,338 @@ static dErr dFSContPropogateDegree(dFS fs) */ static dErr dFSBuildSpace_Cont(dFS fs) { - MPI_Comm comm = ((dObject)fs)->comm; - /* 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 { dInt start,ndofs; } *fieldSpec; - struct dMeshAdjacency ma; - dMesh mesh; - dMeshTag specTag; - /* MeshListEH v=MLZ,e=MLZ,f=MLZ,r=MLZ; */ - /* MeshListInt in=MLZ,rfo=MLZ,feo=MLZ,rdegree=MLZ; */ - /* dIInt nf; */ - dEntTopology *regTopo; - dInt i,j,gcnt,ghcnt,*inodes,*xnodes,*deg,*rdeg,*ghostind,nregions; - dInt *xind,*xstart,*istart,xcnt,*nnz,*pnnz,*regRDeg,*regBDeg; - dMeshEH *regEnts; - dEntStatus *status; - dErr err; + 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 dMeshAdjacency ma; + dMesh mesh; + iMesh_Instance mi; + dEntTopology *regTopo; + dInt *inodes,*xnodes,*deg,*rdeg,nregions,*bstat,ents_a,ents_s,*intdata,*idx,*ghidx; + dInt *xstart,xcnt,*regRDeg,*regBDeg; + dInt bs,n,ngh,ndirichlet,nghdirichlet,rstart,drstart,crstart,dsplit; + dIInt ierr; + dMeshEH *ents; + dEntStatus *status; + dErr err; dFunctionBegin; dValidHeader(fs,DM_COOKIE,1); + bs = fs->bs; mesh = fs->mesh; - err = dMeshGetAdjacency(mesh,fs->active,&fs->meshAdj);dCHK(err); - err = dMemcpy(&ma,fs->meshAdj,sizeof ma);dCHK(err); /* To have an object rather than pointer semantics in this function. */ + err = dMeshGetInstance(mesh,&mi);dCHK(err); + err = dMeshGetAdjacency(mesh,fs->activeSet,&fs->meshAdj);dCHK(err); + err = dMemcpy(&ma,fs->meshAdj,sizeof ma);dCHK(err); /* To have object rather than pointer semantics in this function. */ err = dFSContPropogateDegree(fs);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,&intdata,ents_a,&idx,ents_a,&ghidx);dCHK(err); + + /* Partition entities in active set into explicitly represented and Dirichlet */ + { + dInt nboundaries; + dMeshESH *bdysets; + iMesh_addEntArrToSet(mi,ma.ents,ma.nents,fs->ownedExplicitSet,&ierr);dICHK(mi,ierr); /* Start with all represented explicitly, pretent they are all owned */ + err = dMeshGetNumSubsets(mesh,fs->boundaries,1,&nboundaries);dCHK(err); + if (!nboundaries) goto after_boundaries; + err = dMallocA2(nboundaries,&bdysets,nboundaries,&bstat);dCHK(err); + err = dMeshGetSubsets(mesh,fs->boundaries,1,bdysets,nboundaries,NULL);dCHK(err); + err = dMeshTagSGetData(mesh,fs->bstatusTag,bdysets,nboundaries,bstat,nboundaries,dDATA_INT);dCHK(err); + for (int i=0; iownedExplicitSet,&ierr);dICHK(mi,ierr); + err = dMeshPartitionOnOwnership(mesh,ents,ents_s,&ghstart);dCHK(err); + iMesh_addEntArrToSet(mi,ents,ghstart,fs->ownedDirichletSet,&ierr);dICHK(mi,ierr); + iMesh_addEntArrToSet(mi,ents+ghstart,ents_s-ghstart,fs->ghostDirichletSet,&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->weakFaceSet,&ierr);dICHK(mi,ierr); + } + } + } + after_boundaries: + + /* Partition explicit entities into owned and ghost */ + { + dInt ghstart; + iMesh_getEntitiesRec(mi,fs->ownedExplicitSet,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->ownedExplicitSet,&ierr);dICHK(mi,ierr); + iMesh_addEntArrToSet(mi,ents+ghstart,ents_s-ghstart,fs->ghostExplicitSet,&ierr);dICHK(mi,ierr); + } + + /* Get number of nodes for all entities, and parallel status */ err = dMallocA5(ma.nents*3,°,ma.nents*3,&rdeg,ma.nents,&inodes,ma.nents,&xnodes,ma.nents,&status);dCHK(err); err = dMeshTagGetData(mesh,fs->degreetag,ma.ents,ma.nents,deg,3*ma.nents,dDATA_INT);dCHK(err); err = dMeshTagGetData(mesh,fs->ruletag,ma.ents,ma.nents,rdeg,3*ma.nents,dDATA_INT);dCHK(err); /* Fill the arrays \a inodes and \a xnodes with the number of interior and expanded nodes for each * (topology,degree) pair */ err = dJacobiGetNodeCount(fs->jacobi,ma.nents,ma.topo,deg,inodes,xnodes);dCHK(err); - err = dMeshGetStatus(mesh,ma.nents,ma.ents,status);dCHK(err); + err = dMeshGetStatus(mesh,ma.ents,ma.nents,status);dCHK(err); + + /* Count the number of nodes in each space (owned, local, owned dirichlet, local dirichlet) */ + n = ngh = ndirichlet = nghdirichlet; + for (int i=0; iownedDirichletSet,ma.ents[i],&isdirichlet,&ierr);dICHK(mi,ierr); + iMesh_isEntContained(mi,fs->ghostDirichletSet,ma.ents[i],&isghostdirichlet,&ierr);dICHK(mi,ierr); + iMesh_isEntContained(mi,fs->ownedExplicitSet,ma.ents[i],&isexplicit,&ierr);dICHK(mi,ierr); + iMesh_isEntContained(mi,fs->ghostExplicitSet,ma.ents[i],&isghostexplicit,&ierr);dICHK(mi,ierr); + if (!!isdirichlet + !!isghostdirichlet + !!isexplicit + !!isghostexplicit != 1) dERROR(1,"should not happen"); + if (isdirichlet) ndirichlet += inodes[i]; + else if (isghostdirichlet) nghdirichlet += inodes[i]; + else if (isexplicit) n += inodes[i]; + else if (isghostexplicit) ngh += inodes[i]; + } + err = MPI_Scan(&n,&rstart,1,MPIU_INT,MPI_SUM,comm);dCHK(err); + rstart -= n; + err = MPI_Scan(&ndirichlet,&drstart,1,MPIU_INT,MPI_SUM,comm);dCHK(err); + drstart -= n; + crstart = rstart + drstart; - fs->n = fs->nlocal = fs->nbdofs = 0; - for (i=0; in += inodes[i]; /* It's a global dof on this process */ + { /* Set global index of first node associated with every explicit entity */ + dInt g = rstart,gh; + iMesh_getEntitiesRec(mi,fs->ownedExplicitSet,dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr); + err = dMeshTagGetData(mesh,ma.indexTag,ents,ents_s,idx,ents_s,dDATA_INT);dCHK(err); + for (dInt i=0; igoffsetTag,ents,ents_s,intdata,ents_s,dDATA_INT);dCHK(err); + iMesh_getEntitiesRec(mi,fs->ghostExplicitSet,dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr); + for (dInt i=0; igoffsetTag,ents,ents_s,intdata,ents_s,dDATA_INT);dCHK(err); + err = dMeshTagBcast(mesh,fs->goffsetTag);dCHK(err); + /* Check that all ghost entities were updated */ + err = dMeshTagGetData(mesh,fs->goffsetTag,ents,ents_s,intdata,ents_s,dDATA_INT);dCHK(err); + for (dInt i=0; inlocal += inodes[i]; /* It's always a local dof */ + /* Set ghost indices of every node using \a ghidx */ + gh = 0; + for (dInt i=0; icomm,&fs->slice);dCHK(err); + err = SlicedSetGhosts(fs->slice,bs,n,gh,ghidx);dCHK(err); } - err = MPI_Scan(&fs->n,&fs->rstart,1,MPI_INT,MPI_SUM,comm);dCHK(err); - fs->rstart -= fs->n; - err = MPI_Allreduce(&fs->n,&fs->N,1,MPI_INT,MPI_SUM,comm);dCHK(err); - - /* Set the global index of the first node associated with every entity (and number of nodes, for a consistency check) */ - err = dMallocA(ma.nents,&fieldSpec);dCHK(err); - for (i=0,gcnt=fs->rstart; iownedDirichletSet,dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr); + err = dMeshTagGetData(mesh,ma.indexTag,ents,ents_s,idx,ents_s,dDATA_INT);dCHK(err); + for (dInt i=0; igdoffsetTag,ents,ents_s,intdata,ents_s,dDATA_INT);dCHK(err); + iMesh_getEntitiesRec(mi,fs->ghostDirichletSet,dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr); + for (dInt i=0; igdoffsetTag,ents,ents_s,intdata,ents_s,dDATA_INT);dCHK(err); + err = dMeshTagBcast(mesh,fs->gdoffsetTag);dCHK(err); + /* Check that all ghost entities were updated */ + err = dMeshTagGetData(mesh,fs->gdoffsetTag,ents,ents_s,intdata,ents_s,dDATA_INT);dCHK(err); + for (dInt i=0; icomm,&fs->dslice);dCHK(err); + err = SlicedSetGhosts(fs->dslice,bs,ndirichlet,gh,ghidx);dCHK(err); + /* Create Dirichlet local and global forms (these are homogeneous for now, the user will have to fill them) */ + err = SlicedCreateGlobalVector(fs->dslice,&fs->d);dCHK(err); + err = VecGhostGetLocalForm(fs->d,&fs->dl);dCHK(err); } - /* Sync the field spec (make unowned entities have global offsets for their nodes) */ - err = dMeshTagCreateTemp(mesh,"field_spec",2,dDATA_INT,&specTag);dCHK(err); - err = dMeshTagSetData(mesh,specTag,ma.ents,ma.nents,fieldSpec,ma.nents*(dInt)sizeof(fieldSpec)/(dInt)sizeof(dInt),dDATA_INT);dCHK(err); - err = dMeshTagBcast(mesh,specTag);dCHK(err); - err = dMeshTagGetData(mesh,specTag,ma.ents,ma.nents,fieldSpec,ma.nents*(dInt)sizeof(fieldSpec)/(dInt)sizeof(dInt),dDATA_INT);dCHK(err); - err = dMeshTagDestroy(mesh,specTag);dCHK(err); - for (i=0; iownedExplicitSet,dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr); + et = ents + ents_s; et_a = ents_a - ents_s; + iMesh_getEntitiesRec(mi,fs->ownedDirichletSet,dTYPE_ALL,dTOPO_ALL,1,&et,&et_a,&et_s,&ierr);dICHK(mi,ierr); + ents_s += et_s; + gc = crstart; + for (dInt i=0; igcoffsetTag,ents,ents_s,intdata,ents_s,dDATA_INT);dCHK(err); + /* Create global closure vector */ + err = VecCreateMPI(comm,(n+ndirichlet)*bs,PETSC_DETERMINE,&fs->gc);dCHK(err); + err = VecSetBlockSize(fs->gc,bs);dCHK(err); } - /* Set the global indices for all ghost dofs */ - err = dMallocA(fs->nlocal-fs->n,&ghostind);dCHK(err); /* global index of each ghosted dof */ - err = dMallocA(ma.nents,&istart);dCHK(err); /* First dof of every entity, needed in this form later */ - for (i=0,ghcnt=0; iownedExplicitSet,dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr); + err = dMeshTagGetData(mesh,ma.indexTag,ents,ents_s,idx,ents_s,dDATA_INT);dCHK(err); + for (dInt i=0; iloffsetTag,ents,ents_s,intdata,ents_s,dDATA_INT);dCHK(err); + + /* ghost explicit */ + iMesh_getEntitiesRec(mi,fs->ghostExplicitSet,dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr); + err = dMeshTagGetData(mesh,ma.indexTag,ents,ents_s,idx,ents_s,dDATA_INT);dCHK(err); + for (dInt i=0; islice,&gf);dCHK(err); + err = VecGhostGetLocalForm(gf,&lf);dCHK(err); + err = VecGetSize(lf,&nl);dCHK(err); + err = VecGhostRestoreLocalForm(gf,&lf);dCHK(err); + err = VecDestroy(gf);dCHK(err); + if (nl != (n+ngh)*bs) dERROR(1,"should not happen"); + if (li*bs != nl) dERROR(1,"Inconsistent sizes, should not happen"); + } +#endif + dsplit = li; + err = dMeshTagSetData(mesh,fs->loffsetTag,ents,ents_s,intdata,ents_s,dDATA_INT);dCHK(err); + + /* owned Dirichlet */ + iMesh_getEntitiesRec(mi,fs->ownedDirichletSet,dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr); + err = dMeshTagGetData(mesh,ma.indexTag,ents,ents_s,idx,ents_s,dDATA_INT);dCHK(err); + for (dInt i=0; iloffsetTag,ents,ents_s,intdata,ents_s,dDATA_INT);dCHK(err); + + /* ghost Dirichlet */ + iMesh_getEntitiesRec(mi,fs->ghostDirichletSet,dTYPE_ALL,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr); + err = dMeshTagGetData(mesh,ma.indexTag,ents,ents_s,idx,ents_s,dDATA_INT);dCHK(err); + for (dInt i=0; idl,&nl);dCHK(err); + if (nl != (ndirichlet+nghdirichlet)*bs) dERROR(1,"should not happen"); + if ((li-dsplit)*bs != nl) dERROR(1,"Inconsistent sizes, should not happen"); } - istart[i] = fieldSpec[i].start; + err = dMeshTagSetData(mesh,fs->loffsetTag,ents,ents_s,intdata,ents_s,dDATA_INT);dCHK(err); } - err = dFree(fieldSpec);dCHK(err); - err = SlicedCreate(((dObject)fs)->comm,&fs->sliced);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); + { + IS from,to; + Vec vec; + /* Create global closure to global scatter (ctog) */ + err = ISCreateStride(comm,n*bs,crstart*bs,1,&from);dCHK(err); + err = ISCreateStride(PETSC_COMM_SELF,n*bs,0,1,&to);dCHK(err); + err = SlicedCreateGlobalVector(fs->slice,&vec);dCHK(err); + err = VecScatterCreate(fs->gc,from,vec,to,&fs->ctog);dCHK(err); + err = ISDestroy(from);dCHK(err); + err = ISDestroy(to);dCHK(err); + err = VecDestroy(vec);dCHK(err); + /* Create global closure to Dirichlet scatter (ctod) */ + err = ISCreateStride(comm,ndirichlet*bs,(crstart+n)*bs,1,&from);dCHK(err); + err = ISCreateStride(PETSC_COMM_SELF,ndirichlet*bs,0,1,&to);dCHK(err); + err = VecScatterCreate(fs->gc,from,fs->d,to,&fs->ctod);dCHK(err); + err = ISDestroy(from);dCHK(err); + err = ISDestroy(to);dCHK(err); } - err = dFree(ghostind);dCHK(err); /** * At this point the local to global mapping is complete. Now we need to assemble the constraint matrices which take - * the local vector to an expanded vector. If the mesh is conforming and there are no strange boundaries (i.e. slip or - * normal) the constraint matrix will be boolean (one unit entry per row) in which case an IS would be sufficient. In - * the general case, there will be some non-conforming elements and some strange boundaries. We assemble a full-order - * constraint matrix and a low-order preconditioning constraint matrix. The full-order matrix will be used for - * residual evaluation and matrix-free Jacobian application. The preconditioning constraint matrix will be used to - * assemble the low-order preconditioner for the Jacobian (or blocks there of). + * the local vector to an expanded vector and the local Dirichlet vector to an expanded. If the mesh is conforming and + * there are no strange boundaries (i.e. slip or normal) the constraint matrix will be boolean (one unit entry per row) + * in which case an IS would be sufficient. In the general case, there will be some non-conforming elements and some + * strange boundaries. We assemble a full-order constraint matrix and a low-order preconditioning constraint matrix. + * The full-order matrix will be used for residual evaluation and matrix-free Jacobian application. The + * preconditioning constraint matrix will be used to assemble the low-order preconditioner for the Jacobian (or blocks + * there of). Even the full-order matrix is cheap to apply, but it's use in preconditioner assembly significantly + * impacts sparsity. * * To generate constraint matrices efficiently, we should preallocate them. We will make the (possibly poor) * assumption that every element with different (must be lower!) order approximation on a downward-adjacent entity will * be constrained against all nodes on the adjacent entity. - * - * First, we define the expanded space which is just all regions in the domain (we're implementing a scalar space - * without boundaries). */ - nregions = ma.toff[dTYPE_REGION+1] - ma.toff[dTYPE_REGION]; - err = dMallocA6(nregions,&xind,nregions+1,&xstart,nregions,®Topo,nregions*3,®RDeg,nregions*3,®BDeg,nregions,®Ents);dCHK(err); - for (i=0,xcnt=0; iactiveSet,dTYPE_REGION,dTOPO_ALL,1,&ents,&ents_a,&ents_s,&ierr);dICHK(mi,ierr); + err = dMeshTagGetData(mesh,ma.indexTag,ents,ents_s,idx,ents_s,dDATA_INT);dCHK(err); + nregions = ents_s; + err = dMallocA4(nregions+1,&xstart,nregions,®Topo,nregions*3,®RDeg,nregions*3,®BDeg);dCHK(err); + xcnt = 0; + for (dInt i=0; iruleStrength); + type = iMesh_TypeFromTopology[regTopo[i]]; + for (dInt j=0; jruleStrength); regBDeg[i*3+j] = deg[ii*3+j]; } + for (dInt j=type; j<3; j++) { + regRDeg[i*3+2] = 1; + regBDeg[i*3+2] = 1; + } xcnt += xnodes[ii]; } - fs->m = xstart[nregions] = xcnt; + xstart[nregions] = xcnt; - err = dMallocA2(xcnt,&nnz,xcnt,&pnnz);dCHK(err); - err = dJacobiGetConstraintCount(fs->jacobi,nregions,xind,xstart,deg,&ma,nnz,pnnz);dCHK(err); + { + dInt *nnz,*pnnz,*dnnz; + Mat E,Ep,Ed; + err = dMallocA3(xcnt,&nnz,xcnt,&pnnz,xcnt,&dnnz);dCHK(err); + err = dMeshTagGetData(mesh,fs->loffsetTag,ma.ents,ma.nents,intdata,ma.nents,dDATA_INT);dCHK(err); + /* To generate element assembly matrices, we need + * \a idx the MeshAdjacency index of every region + * \a xstart offset in expanded vector of first node associated with this region + * \a dsplit integer which splits \a istart into local and local Dirichlet pieces + * \a istart offset in local vectors of first dof associated with each entity (not just regions) + * let \c is=istart[e]. If \c is=dsplit + * then \c is-dsplit is the offset in local Dirichlet vector. The array \a istart is the generic buffer \a intdata + * \a deg integer array of length \c 3*ma.nents which holds the degree of every entity in MeshAdjacency + * \a ma MeshAdjacency (array-based connectivity) + * + * We will create matrices + * \a E full order element assembly matrix + * \a Ep preconditioning element assembly matrix (as sparse as possible) + * \a Ed Dirichlet element assembly matrix + * + * These are preallocated using \a nnz, \a pnnz, \a dnnz respectively. + **/ + err = dJacobiGetConstraintCount(fs->jacobi,nregions,idx,xstart,dsplit,intdata,deg,&ma,nnz,pnnz,dnnz);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 = dFree2(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,xcnt,n+ngh,1,nnz,&E);dCHK(err); + err = MatCreateSeqAIJ(PETSC_COMM_SELF,xcnt,n+ngh,1,pnnz,&Ep);dCHK(err); + if (ndirichlet+nghdirichlet > 0) { + err = MatCreateSeqAIJ(PETSC_COMM_SELF,xcnt,ndirichlet+nghdirichlet,1,dnnz,&Ed);dCHK(err); + } else { + err = MatCreateSeqAIJ(PETSC_COMM_SELF,xcnt,0,0,NULL,&Ed);dCHK(err); + } + err = dFree3(nnz,pnnz,dnnz);dCHK(err); - err = dJacobiAddConstraints(fs->jacobi,nregions,xind,xstart,istart,deg,&ma,fs->C,fs->Cp);dCHK(err); - err = dFree(istart);dCHK(err); - err = dFree5(deg,rdeg,inodes,xnodes,status);dCHK(err); - err = dMeshRestoreAdjacency(mesh,fs->active,&fs->meshAdj);dCHK(err); /* Any reason to leave this around for longer? */ + err = dJacobiAddConstraints(fs->jacobi,nregions,idx,xstart,dsplit,intdata,deg,&ma,E,Ep,Ed);dCHK(err); + err = dFree5(deg,rdeg,inodes,xnodes,status);dCHK(err); + err = dMeshRestoreAdjacency(mesh,fs->activeSet,&fs->meshAdj);dCHK(err); /* Any reason to leave this around for longer? */ - err = MatAssemblyBegin(fs->C,MAT_FINAL_ASSEMBLY);dCHK(err); - err = MatAssemblyBegin(fs->Cp,MAT_FINAL_ASSEMBLY);dCHK(err); - err = MatAssemblyEnd(fs->C,MAT_FINAL_ASSEMBLY);dCHK(err); - err = MatAssemblyEnd(fs->Cp,MAT_FINAL_ASSEMBLY);dCHK(err); + err = MatAssemblyBegin(E,MAT_FINAL_ASSEMBLY);dCHK(err); + err = MatAssemblyBegin(Ep,MAT_FINAL_ASSEMBLY);dCHK(err); + err = MatAssemblyBegin(Ed,MAT_FINAL_ASSEMBLY);dCHK(err); + err = MatAssemblyEnd(E,MAT_FINAL_ASSEMBLY);dCHK(err); + err = MatAssemblyEnd(Ep,MAT_FINAL_ASSEMBLY);dCHK(err); + err = MatAssemblyEnd(Ed,MAT_FINAL_ASSEMBLY);dCHK(err); + + err = MatCreateMAIJ(E,bs,&fs->E);dCHK(err); + err = MatCreateMAIJ(Ep,bs,&fs->Ep);dCHK(err); + err = MatCreateMAIJ(Ed,bs,&fs->Ed);dCHK(err); + + err = MatDestroy(E);dCHK(err); + err = MatDestroy(Ep);dCHK(err); + err = MatDestroy(Ed);dCHK(err); + } /* Get Rule and EFS for domain ents. */ fs->nelem = nregions; @@ -214,8 +409,9 @@ static dErr dFSBuildSpace_Cont(dFS fs) err = dMemcpy(fs->off,xstart,(nregions+1)*sizeof(xstart[0]));dCHK(err); err = dJacobiGetRule(fs->jacobi,nregions,regTopo,regRDeg,fs->rule);dCHK(err); err = dJacobiGetEFS(fs->jacobi,nregions,regTopo,regBDeg,fs->rule,fs->efs);dCHK(err); - err = dMeshGetVertexCoords(mesh,nregions,regEnts,&fs->vtxoff,&fs->vtx);dCHK(err); /* Should be restored by FS on destroy */ - err = dFree6(xind,xstart,regTopo,regRDeg,regBDeg,regEnts);dCHK(err); + err = dMeshGetVertexCoords(mesh,nregions,ents,&fs->vtxoff,&fs->vtx);dCHK(err); /* Should be restored by FS on destroy */ + err = dFree4(xstart,regTopo,regRDeg,regBDeg);dCHK(err); + err = dFree4(ents,intdata,idx,ghidx);dCHK(err); dFunctionReturn(0); } @@ -236,7 +432,7 @@ dErr dFSCreate_Cont(dFS fs) dFunctionBegin; err = dNewLog(fs,*fsc,&fsc);dCHK(err); - fs->D = 1; + fs->bs = 1; fs->data = (void*)fsc; fs->ops->view = dFSView_Cont; fs->ops->impldestroy = dFSDestroy_Cont; diff --git a/src/fs/interface/fs.c b/src/fs/interface/fs.c index 6d9f1ef8..00c84381 100644 --- a/src/fs/interface/fs.c +++ b/src/fs/interface/fs.c @@ -13,8 +13,20 @@ dErr dFSSetMesh(dFS fs,dMesh mesh,dMeshESH active) if (mesh != qmesh) fs->quotient = 0; /* The Quotient is stale */ } fs->mesh = mesh; - fs->active = active; + fs->activeSet = active; err = dMeshGetTag(mesh,fs->bdyTagName,&fs->bdyTag);dCHK(err); + err = dMeshTagCreateTemp(mesh,"boundary_status",1,dDATA_INT,&fs->bstatusTag);dCHK(err); + err = dMeshTagCreateTemp(mesh,"boundary_constraint",sizeof(struct dFSConstraintCtx),dDATA_BYTE,&fs->bdyConstraintTag);dCHK(err); + err = dMeshTagCreateTemp(mesh,"global_offset",1,dDATA_INT,&fs->goffsetTag);dCHK(err); + err = dMeshTagCreateTemp(mesh,"global_dirichlet_offset",1,dDATA_INT,&fs->gdoffsetTag);dCHK(err); + err = dMeshTagCreateTemp(mesh,"local_offset",1,dDATA_INT,&fs->loffsetTag);dCHK(err); + err = dMeshTagCreate(mesh,"global_closure_offset",1,dDATA_INT,&fs->gcoffsetTag);dCHK(err); + err = dMeshSetCreate(mesh,&fs->ownedExplicitSet);dCHK(err); + err = dMeshSetCreate(mesh,&fs->ghostExplicitSet);dCHK(err); + err = dMeshSetCreate(mesh,&fs->ownedDirichletSet);dCHK(err); + err = dMeshSetCreate(mesh,&fs->ghostDirichletSet);dCHK(err); + err = dMeshSetCreate(mesh,&fs->weakFaceSet);dCHK(err); + err = dMeshSetCreate(mesh,&fs->boundaries);dCHK(err); dFunctionReturn(0); } @@ -44,14 +56,14 @@ dErr dFSSetDegree(dFS fs,dJacobi jac,dMeshTag deg) /** Set the block size for a function space. * * @param fs function space -* @param D block size (number of dofs per unconstrained node) +* @param bs block size (number of dofs per node) **/ -dErr dFSSetBlockSize(dFS fs,dInt D) +dErr dFSSetBlockSize(dFS fs,dInt bs) { dFunctionBegin; dValidHeader(fs,DM_COOKIE,1); - fs->D = D; + fs->bs = bs; dFunctionReturn(0); } @@ -60,8 +72,7 @@ dErr dFSSetBlockSize(dFS fs,dInt D) * * @param fs function space object * @param mid Boundary ID, usually the value of the NEUMANN_SET tag -* @param strong If true, the condition will be enforced strongly (all components removed from global system) -* @param expanded If true, entries will be present in expanded space (so that weak forms can be evaluated) +* @param bstat Boundary status, determines how boundary should be represented (weak, global, dirichlet) * @param cfunc constraint function, only makes sense if !strong, but some degrees of freedom must be removed anyway * @param user context for constraint function * @@ -74,21 +85,27 @@ dErr dFSSetBlockSize(dFS fs,dInt D) * declared statically (outside of the function definition) is merely to avoid duplicating information that must be kept * consistent. **/ -dErr dFSRegisterBoundary(dFS fs,dInt mid,dTruth strong,dTruth expanded,dFSBoundaryConstraintFunction cfunc,void *user) +dErr dFSRegisterBoundary(dFS fs,dInt mid,dFSBStatus bstat,dFSConstraintFunction cfunc,void *user) { - dFSBoundary bdy; - dErr err; + dMeshESH bset; + iMesh_Instance mi; + dErr err; + dIInt ierr; dFunctionBegin; dValidHeader(fs,DM_COOKIE,1); - err = dNew(struct _p_dFSBoundary,&bdy);dCHK(err); - err = dMeshGetTaggedSet(fs->mesh,fs->bdyTag,&mid,&bdy->manset);dCHK(err); - bdy->strong = strong; - bdy->expanded = expanded; - bdy->cfunc = cfunc; - bdy->user = user; - bdy->next = fs->bdyStart; - fs->bdyStart = bdy; + err = dFSBStatusValid(bstat);dCHK(err); + if (dFSBStatusStrongCount(bstat) > fs->bs) dERROR(1,"Cannot impose strong conditions on more dofs than the block size"); + err = dMeshGetTaggedSet(fs->mesh,fs->bdyTag,&mid,&bset);dCHK(err); + err = dMeshTagSSetData(fs->mesh,fs->bstatusTag,&bset,1,&bstat,1,dDATA_INT);dCHK(err); + if (cfunc) { + struct dFSConstraintCtx ctx; + ctx.cfunc = cfunc; + ctx.user = user; + err = dMeshTagSSetData(fs->mesh,fs->bdyConstraintTag,&bset,1,&ctx,sizeof(ctx),dDATA_BYTE);dCHK(err); + } + err = dMeshGetInstance(fs->mesh,&mi);dCHK(err); + iMesh_addEntSet(mi,bset,fs->boundaries,&ierr);dICHK(mi,ierr); dFunctionReturn(0); } @@ -119,18 +136,16 @@ dErr dFSView(dFS fs,dViewer viewer) dInt nents[4]; err = PetscViewerASCIIPrintf(viewer,"General information about the mesh topology.\n");dCHK(err); for (dEntType type=dTYPE_VERTEX; typemesh,fs->active,type,dTOPO_ALL,&nents[type]);dCHK(err); + err = dMeshGetNumEnts(fs->mesh,fs->activeSet,type,dTOPO_ALL,&nents[type]);dCHK(err); } err = PetscViewerASCIIPrintf(viewer,"number of vertices=%d edges=%d faces=%d regions=%d\n",nents[0],nents[1],nents[2],nents[3]);dCHK(err); } { /* print aggregate sizes */ PetscMPIInt gm[2],lm[2]; - lm[0] = fs->m; lm[1] = fs->nlocal; /* set local `element' size and `local' size */ + err = MatGetSize(fs->E,&lm[0],&lm[1]);dCHK(err); err = MPI_Reduce(lm,gm,2,MPI_INT,MPI_SUM,0,((dObject)fs)->comm);dCHK(err); err = PetscViewerASCIIPrintf(viewer,"on rank 0, %d/%d element dofs constrained against %d/%d local dofs\n", lm[0],gm[0],lm[1],gm[1]);dCHK(err); - err = PetscViewerASCIIPrintf(viewer,"on rank 0, %d local (%d owned,%d ghosted) out of %d global", - fs->nlocal,fs->n,fs->nlocal-fs->n,fs->N);dCHK(err); } if (fs->ops->view) { err = (*fs->ops->view)(fs,viewer);dCHK(err); @@ -165,11 +180,11 @@ dErr dFSDestroy(dFS fs) default: dERROR(1,"Invalid status %d",w->status); } } - if (fs->sliced) { - err = SlicedDestroy(fs->sliced);dCHK(err); - } - err = MatDestroy(fs->C);dCHK(err); - err = MatDestroy(fs->Cp);dCHK(err); + if (fs->slice) {err = SlicedDestroy(fs->slice);dCHK(err);} + if (fs->dslice) {err = SlicedDestroy(fs->dslice);dCHK(err);} + err = MatDestroy(fs->E);dCHK(err); + err = MatDestroy(fs->Ep);dCHK(err); + err = MatDestroy(fs->Ed);dCHK(err); err = dFree3(fs->rule,fs->efs,fs->off);dCHK(err); err = dMeshRestoreVertexCoords(fs->mesh,fs->nelem,NULL,&fs->vtxoff,&fs->vtx);dCHK(err); err = PetscHeaderDestroy(fs);dCHK(err); @@ -219,7 +234,7 @@ dErr dFSCreateExpandedVector(dFS fs,Vec *x) dFunctionBegin; dValidHeader(fs,DM_COOKIE,1); dValidPointer(x,2); - err = MatGetVecs(fs->C,NULL,x);dCHK(err); + err = MatGetVecs(fs->E,NULL,x);dCHK(err); dFunctionReturn(0); } @@ -230,11 +245,22 @@ dErr dFSCreateGlobalVector(dFS fs,Vec *g) dFunctionBegin; dValidHeader(fs,DM_COOKIE,1); dValidPointer(g,2); - err = SlicedCreateGlobalVector(fs->sliced,g);dCHK(err); + err = SlicedCreateGlobalVector(fs->slice,g);dCHK(err); + dFunctionReturn(0); +} + +dErr dFSCreateDirichletVector(dFS fs,Vec *d) +{ + dErr err; + + dFunctionBegin; + dValidHeader(fs,DM_COOKIE,1); + dValidPointer(d,2); + err = SlicedCreateGlobalVector(fs->dslice,d);dCHK(err); dFunctionReturn(0); } -dErr dFSGlobalToExpandedBegin(dFS dUNUSED fs,Vec g,InsertMode imode,Vec dUNUSED x) +dErr dFSGlobalToExpandedBegin(dFS dUNUSED fs,Vec g,dFSHomogeneousMode dUNUSED hmode,Vec dUNUSED x) { dErr err; @@ -242,11 +268,11 @@ dErr dFSGlobalToExpandedBegin(dFS dUNUSED fs,Vec g,InsertMode imode,Vec dUNUSED dValidHeader(fs,DM_COOKIE,1); dValidHeader(g,VEC_COOKIE,2); dValidHeader(x,VEC_COOKIE,4); - err = VecGhostUpdateBegin(g,imode,SCATTER_FORWARD);dCHK(err); + err = VecGhostUpdateBegin(g,INSERT_VALUES,SCATTER_FORWARD);dCHK(err); dFunctionReturn(0); } -dErr dFSGlobalToExpandedEnd(dFS fs,Vec g,InsertMode imode,Vec x) +dErr dFSGlobalToExpandedEnd(dFS fs,Vec g,dFSHomogeneousMode hmode,Vec x) { Vec lform; dErr err; @@ -255,22 +281,39 @@ dErr dFSGlobalToExpandedEnd(dFS fs,Vec g,InsertMode imode,Vec x) dValidHeader(fs,DM_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); + switch (hmode) { + case dFS_HOMOGENEOUS: + err = VecZeroEntries(x);dCHK(err); break; - case ADD_VALUES: - err = MatMultAdd(fs->C,lform,x,x);dCHK(err); + case dFS_INHOMOGENEOUS: + err = MatMult(fs->Ed,fs->dl,x);dCHK(err); break; default: - dERROR(1,"InsertMode %d not supported",imode); + dERROR(1,"dFSHomogeneousMode %d not supported",hmode); } + err = VecGhostUpdateEnd(g,INSERT_VALUES,SCATTER_FORWARD);dCHK(err); + err = VecGhostGetLocalForm(g,&lform);dCHK(err); + err = MatMultAdd(fs->E,lform,x,x);dCHK(err); err = VecGhostRestoreLocalForm(g,&lform);dCHK(err); dFunctionReturn(0); } +/** Assemble expanded vector into global form. +* +* @param fs Function space +* @param x Expanded vector +* @param imode sum into local form or overwrite it, see note below +* @param g Global vector (must be a VecGhost, as obtained with dFSCreateGlobalVector()), see note 2 below +* +* @note \a imode is treated differently here than in most of PETSc, it never runs a scatter with INSERT_VALUES since +* there are \e always multiple inputs mapped to the same output (at least for continuous spaces). Instead \a imode +* determines whether to clear the local form before assembly (INSERT_VALUES) or whether to just sum into the existing +* vector (ADD_VALUES). +* +* @note Not collective! While the name indicates that it only updates the global vector, it also updates the local +* form. For the collective operation which sums contributions from multiple processes, see dFSExpandedToGlobalBegin(). +* +**/ dErr dFSExpandedToGlobal(dFS fs,Vec x,InsertMode imode,Vec g) { Vec lform; @@ -283,10 +326,10 @@ dErr dFSExpandedToGlobal(dFS fs,Vec x,InsertMode imode,Vec g) err = VecGhostGetLocalForm(g,&lform);dCHK(err); switch (imode) { case INSERT_VALUES: - err = MatMultTranspose(fs->C,x,lform);dCHK(err); + err = MatMultTranspose(fs->E,x,lform);dCHK(err); break; case ADD_VALUES: - err = MatMultTransposeAdd(fs->C,x,lform,lform);dCHK(err); + err = MatMultTransposeAdd(fs->E,x,lform,lform);dCHK(err); break; default: dERROR(1,"InsertMode %d not supported",imode); @@ -295,7 +338,10 @@ dErr dFSExpandedToGlobal(dFS fs,Vec x,InsertMode imode,Vec g) dFunctionReturn(0); } -/** Updates the global values, does \b not broadcast the global values back to the ghosts */ +/** Updates the global values, does \b not broadcast the global values back to the ghosts. +* +* Call VecGhostUpdate{Begin,End}(g,INSERT_VALUES,SCATTER_FORWARD) after this to update ghost values. +**/ dErr dFSExpandedToGlobalBegin(dFS fs,Vec x,InsertMode imode,Vec g) { dErr err; @@ -321,6 +367,27 @@ dErr dFSExpandedToGlobalEnd(dFS dUNUSED fs,Vec dUNUSED x,InsertMode dUNUSED imod dFunctionReturn(0); } +dErr dFSExpandedToDirichlet(dFS fs,Vec x,InsertMode imode,Vec d) +{ + dErr err; + + dFunctionBegin; + dValidHeader(fs,DM_COOKIE,1); + dValidHeader(d,VEC_COOKIE,2); + dValidHeader(x,VEC_COOKIE,4); + switch (imode) { + case INSERT_VALUES: + err = MatMultTranspose(fs->Ed,x,d);dCHK(err); + break; + case ADD_VALUES: + err = MatMultTransposeAdd(fs->Ed,x,d,d);dCHK(err); + break; + default: + dERROR(1,"InsertMode %d not supported",imode); + } + dFunctionReturn(0); +} + dErr dFSGetElements(dFS fs,dInt *n,dInt *restrict*off,s_dRule *restrict*rule,s_dEFS *restrict*efs,dInt *restrict*geomoff,dReal (*restrict*geom)[3]) { @@ -355,10 +422,10 @@ dErr dFSRestoreElements(dFS dUNUSED fs,dInt *n,dInt *restrict*off,s_dRule *restr * @param q Pointer which will hold array of quadrature points, in physical space (not just the local tensor product) * @param jinv Inverse of element Jacobian evaluated at quadrature points, normally just passed on to dEFSApply * @param jw Array to store jacobian determinant times quadrature weights at quadrature points -* @param u first array to hold D values per quadrature point -* @param v second array to hold D values per quadrature point -* @param du first array to hold 3*D values per quadrature point -* @param dv second array to hold 3*D values per quadrature point +* @param u first array to hold \c bs values per quadrature point +* @param v second array to hold \c bs values per quadrature point +* @param du first array to hold \c 3*bs values per quadrature point +* @param dv second array to hold \c 3*bs values per quadrature point */ dErr dFSGetWorkspace(dFS fs,const char name[],dReal (*restrict*q)[3],dReal (*restrict*jinv)[3][3],dReal *restrict*jw,dScalar *restrict*u,dScalar *restrict*v,dScalar *restrict*du,dScalar *restrict*dv) { @@ -379,11 +446,11 @@ dErr dFSGetWorkspace(dFS fs,const char name[],dReal (*restrict*q)[3],dReal (*res } Q = fs->maxQ; for (dInt i=0; iD; + const dInt bs = fs->bs; w = &fs->workspace[i]; switch (w->status) { case 0: /* Not allocated */ - err = dMallocA7(Q,&w->q,Q,&w->jinv,Q,&w->jw,Q*D,&w->u,Q*D,&w->v,Q*D*3,&w->du,Q*D*3,&w->dv);dCHK(err); + err = dMallocA7(Q,&w->q,Q,&w->jinv,Q,&w->jw,Q*bs,&w->u,Q*bs,&w->v,Q*bs*3,&w->du,Q*bs*3,&w->dv);dCHK(err); w->status = 1; case 1: /* Available */ goto found; @@ -450,7 +517,7 @@ dErr dFSGetMatrix(dFS fs,const MatType mtype,Mat *J) dFunctionBegin; dValidHeader(fs,DM_COOKIE,1); dValidPointer(J,3); - err = SlicedGetMatrix(fs->sliced,mtype,J);dCHK(err); + err = SlicedGetMatrix(fs->slice,mtype,J);dCHK(err); dFunctionReturn(0);dCHK(err); } @@ -477,7 +544,7 @@ dErr dFSMatSetValuesExpanded(dFS fs,Mat A,dInt m,const dInt idxm[],dInt n,const dValidPointer(idxn,6); dValidPointer(v,7); err = PetscLogEventBegin(dLOG_FSMatSetValuesExpanded,fs,A,0,0);dCHK(err); - C = (fs->assemblefull) ? fs->C : fs->Cp; + C = (fs->assemblefull) ? fs->E : fs->Ep; err = MatGetRowIJ(C,0,dFALSE,dFALSE,&cn,&ci,&cj,&done);dCHK(err); if (!done) dERROR(1,"Could not get indices"); err = MatGetArray_SeqAIJ(C,&ca);dCHK(err); diff --git a/src/fs/interface/fsreg.c b/src/fs/interface/fsreg.c index 60a996b9..4447c88d 100644 --- a/src/fs/interface/fsreg.c +++ b/src/fs/interface/fsreg.c @@ -8,13 +8,22 @@ static PetscFList FSList = 0; * These default operations are shared with the DM. We are making a two-level inheritance since there may be different * dFS implementations. Certainly both continuous and discontinuous Galerkin are desirable. * -*/ +* Since we use VecGhost, the local and global vectors are really the same. Since there is no way to get a global vector +* from a local vector, the user cannot avoid seeing VecGhost. We do have expanded vectors which are genuinely distinct, +* but I'm very skeptical of the usefulness of identifying them with the "local vector" of a DM. +**/ static const struct _dFSOps defaultFSOps = { .view = dFSView, .createglobalvector = dFSCreateGlobalVector, .createlocalvector = dFSCreateExpandedVector, .localtoglobal = dFSExpandedToGlobal, + /* I think that these don't make sense with the current design. In + * particular, there may be points in the expanded space which are not + * represented in the global vector. In this configuration, an INSERT_MODE + * scatter is not sufficient because it is not surjective. + * .globaltolocalbegin = dFSGlobalToExpandedBegin, .globaltolocalend = dFSGlobalToExpandedEnd, + */ .getmatrix = dFSGetMatrix, .destroy = dFSDestroy, .impldestroy = 0}; diff --git a/src/fs/mesh/interface/mesh.c b/src/fs/mesh/interface/mesh.c index c8491a0c..5aa0003b 100644 --- a/src/fs/mesh/interface/mesh.c +++ b/src/fs/mesh/interface/mesh.c @@ -356,28 +356,39 @@ dErr dMeshTagCreateTemp(dMesh mesh,const char template[],dInt count,dDataType ty dFunctionBegin; dValidHeader(mesh,dMESH_COOKIE,1); dValidPointer(intag,4); - err = PetscSNPrintf(name,sizeof(name),"TEMP_%s_%d",template,unique_id++);dCHK(err); + err = PetscSNPrintf(name,sizeof(name),"__DOHP_%s_%d",template,unique_id++);dCHK(err); err = dMeshTagCreate(mesh,name,count,type,intag);dCHK(err); dFunctionReturn(0); } dErr dMeshTagCreate(dMesh mesh,const char name[],dInt count,dDataType type,dMeshTag *intag) { - dMeshTag tag; - dIInt ierr,itype; - size_t namelen; - dErr err; + dMeshTag tag; + iMesh_Instance mi; + dIInt ierr,itype; + size_t namelen; + dErr err; dFunctionBegin; dValidHeader(mesh,dMESH_COOKIE,1); dValidPointer(intag,4); *intag = 0; - if (count > 0) { - err = dDataTypeToITAPS(type,&itype);dCHK(err); - err = dStrlen(name,&namelen);dCHK(err); - iMesh_createTag(mesh->mi,name,count,itype,&tag,&ierr,(dIInt)namelen);dICHK(mesh->mi,ierr); - *intag = tag; + mi = mesh->mi; + if (count <= 0) dERROR(1,"Tags must have positive count"); + err = dDataTypeToITAPS(type,&itype);dCHK(err); + err = dStrlen(name,&namelen);dCHK(err); + iMesh_getTagHandle(mi,name,&tag,&ierr,(dIInt)namelen); + if (ierr == iBase_TAG_NOT_FOUND) { + iMesh_createTag(mi,name,count,itype,&tag,&ierr,(dIInt)namelen);dICHK(mi,ierr); + } else { + dIInt existing_itype,existing_count; + dICHK(mi,ierr); + iMesh_getTagType(mi,tag,&existing_itype,&ierr);dICHK(mi,ierr); + if (itype != existing_itype) dERROR(1,"tag exists with different type"); + iMesh_getTagSizeValues(mi,tag,&existing_count,&ierr);dICHK(mi,ierr); + if (count != existing_count) dERROR(1,"tag exists with same type but different count"); } + *intag = tag; dFunctionReturn(0); } @@ -420,6 +431,7 @@ dErr dMeshTagGetData(dMesh mesh,dMeshTag tag,const dMeshEH ents[],dInt ecount,vo dValidPointer(ents,3); dValidPointer(data,5); alloc = count * iBase_SizeFromType[type]; + size = 0; /* protect against degenerate case */ iMesh_getArrData(mi,ents,ecount,tag,&dptr,&alloc,&size,&ierr);dICHK(mi,ierr); if (dptr != (char*)data || alloc != count * iBase_SizeFromType[type]) dERROR(1,"Looks like an iMesh inconsistency, the library shouldn't be messing with this"); @@ -427,9 +439,16 @@ dErr dMeshTagGetData(dMesh mesh,dMeshTag tag,const dMeshEH ents[],dInt ecount,vo dFunctionReturn(0); } +/* Might only work for MOAB */ +dErr dMeshTagSSetData(dMesh mesh,dMeshTag tag,const dMeshESH esets[],dInt ecount,const void *data,dInt count,dDataType type) +{ + return dMeshTagSetData(mesh,tag,(dMeshEH*)esets,ecount,data,count,type); +} + +/* Might only work for MOAB */ dErr dMeshTagSGetData(dMesh mesh,dMeshTag tag,const dMeshESH esets[],dInt ecount,void *data,dInt count,dDataType type) { - return dMeshTagGetData(mesh,tag,(dMeshEH*)esets,ecount,data,count,type); /* Maybe only for MOAB */ + return dMeshTagGetData(mesh,tag,(dMeshEH*)esets,ecount,data,count,type); } dErr dMeshGetTaggedSet(dMesh mesh,dMeshTag tag,const void *value,dMeshESH *set) @@ -475,6 +494,31 @@ dErr dMeshGetEnts(dMesh mesh,dMeshESH set,dEntType type,dEntTopology topo,dMeshE dFunctionReturn(0); } +dErr dMeshGetNumSubsets(dMesh mesh,dMeshESH set,dInt hops,dInt *nsubsets) +{ + dIInt ierr,n; + + dFunctionBegin; + dValidHeader(mesh,dMESH_COOKIE,1); + dValidPointer(nsubsets,4); + iMesh_getNumEntSets(mesh->mi,set,hops,&n,&ierr);dICHK(mesh->mi,ierr); + *nsubsets = n; + dFunctionReturn(0); +} + +dErr dMeshGetSubsets(dMesh mesh,dMeshESH set,dInt hops,dMeshESH subsets[],dInt alloc,dInt *size) +{ + dIInt ierr,s; + + dFunctionBegin; + dValidHeader(mesh,dMESH_COOKIE,1); + dValidPointer(subsets,4); + s = 0; + iMesh_getEntSets(mesh->mi,set,hops,&subsets,&alloc,&s,&ierr);dICHK(mesh->mi,ierr); + if (size) *size = s; + dFunctionReturn(0); +} + /** Get entities of every type in a set. * @param mesh mesh object * @param set set to get entities from @@ -514,7 +558,7 @@ dErr dMeshGetEntsOff(dMesh mesh,dMeshESH set,dInt toff[],dMeshEH **inents) * * @return err */ -dErr dMeshGetStatus(dMesh mesh,dInt count,const dMeshEH ents[],dEntStatus status[]) +dErr dMeshGetStatus(dMesh mesh,const dMeshEH ents[],dInt count,dEntStatus status[]) { dMPIInt size; dMeshTag tag; @@ -561,6 +605,18 @@ dErr dMeshTagBcast(dMesh m,dMeshTag tag) dFunctionReturn(0); } +dErr dMeshSetCreate(dMesh mesh,dMeshESH *inset) +{ + dIInt ierr; + + dFunctionBegin; + dValidHeader(mesh,dMESH_COOKIE,1); + dValidPointer(inset,2); + *inset = 0; + iMesh_createEntSet(mesh->mi,0,inset,&ierr);dICHK(mesh->mi,ierr); + dFunctionReturn(0); +} + dErr dMeshLoad(dMesh mesh) { iMesh_Instance mi = mesh->mi; @@ -606,13 +662,6 @@ dErr dMeshLoad(dMesh mesh) } else dICHK(mi,ierr); iMesh_getEntSetEHData(mi,root,emptysettag,(iBase_EntityHandle*)&mesh->emptyset,&ierr);dICHK(mi,ierr); - /* Get the manifold tag, create if it doesn't exist */ - iMesh_getTagHandle(mi,dTAG_MANIFOLD_ID,&mesh->manifoldTag,&ierr,sizeof(dTAG_MANIFOLD_ID)); - if (ierr == iBase_TAG_NOT_FOUND) { - err = dInfo(mesh,"Missing tag \"%s\". Creating but this may cause problems",dTAG_MANIFOLD_ID);dCHK(err); - iMesh_createTag(mi,dTAG_MANIFOLD_ID,1,iBase_INTEGER,&mesh->manifoldTag,&ierr,sizeof(dTAG_MANIFOLD_ID));dICHK(mi,ierr); - } else dICHK(mi,ierr); - /* Get all entities of each type. */ iMesh_getEntities(mi,root,iBase_REGION,iMesh_ALL_TOPOLOGIES,&mesh->r.v,&mesh->r.a,&mesh->r.s,&ierr);dICHK(mi,ierr); iMesh_getEntities(mi,root,iBase_FACE,iMesh_ALL_TOPOLOGIES,&mesh->f.v,&mesh->f.a,&mesh->f.s,&ierr);dICHK(mi,ierr); @@ -945,7 +994,6 @@ dErr dMeshGetAdjacency(dMesh mesh,dMeshESH set,dMeshAdjacency *inadj) { iMesh_Instance mi = mesh->mi; struct dMeshAdjacency ma; - dMeshTag indexTag; dMeshEH *adj,*conn; dEntType type; dInt i,rank,cnt,nadj,*connoff,tnents,*eind; @@ -1002,8 +1050,8 @@ dErr dMeshGetAdjacency(dMesh mesh,dMeshESH set,dMeshAdjacency *inadj) } /* Create the tag to hold indices and set it with strictly increasing values */ - err = dMeshTagCreateTemp(mesh,"index",1,dDATA_INT,&indexTag);dCHK(err); - err = dMeshTagSetData(mesh,indexTag,ma.ents,ma.nents,eind,ma.nents,dDATA_INT);dCHK(err); + err = dMeshTagCreateTemp(mesh,"index",1,dDATA_INT,&ma.indexTag);dCHK(err); + err = dMeshTagSetData(mesh,ma.indexTag,ma.ents,ma.nents,eind,ma.nents,dDATA_INT);dCHK(err); err = dFree(eind);dCHK(err); /* @@ -1052,9 +1100,8 @@ dErr dMeshGetAdjacency(dMesh mesh,dMeshESH set,dMeshAdjacency *inadj) } } - err = dMeshTagGetData(mesh,indexTag,adj,nadj,ma.adjind,nadj,dDATA_INT);dCHK(err); /* get indices of adj ents */ + err = dMeshTagGetData(mesh,ma.indexTag,adj,nadj,ma.adjind,nadj,dDATA_INT);dCHK(err); /* get indices of adj ents */ err = dFree(adj);dCHK(err); - err = dMeshTagDestroy(mesh,indexTag);dCHK(err); /* Determine permutation of adjacent entities */ err = dMeshAdjacencyPermutations_Private(&ma,connoff,conn);dCHK(err); @@ -1079,6 +1126,7 @@ dErr dMeshRestoreAdjacency(dMesh dUNUSED mesh,dMeshESH set,dMeshAdjacency *inma) ma = *inma; *inma = 0; if (set != ma->set) dERROR(1,"Adjacency for the wrong set"); + err = dMeshTagDestroy(mesh,ma->indexTag);dCHK(err); err = dFree3(ma->ents,ma->adjoff,ma->topo);dCHK(err); err = dFree2(ma->adjind,ma->adjperm);dCHK(err); #if defined(dMESHADJACENCY_HAS_CONNECTIVITY) @@ -1149,3 +1197,26 @@ dErr dMeshRestoreVertexCoords(dMesh dUNUSED mesh,dUNUSED dInt n,const dUNUSED dM err = dFree2(*inx,*inxoff);dCHK(err); dFunctionReturn(0); } + + +dErr dMeshPartitionOnOwnership(dMesh mesh,dMeshEH ents[],dInt n,dInt *ghstart) +{ + dEntStatus *status; + dMeshEH *pents; + dInt owned; + dErr err; + + dFunctionBegin; + err = dMallocA2(n,&status,n,&pents);dCHK(err); + err = dMeshGetStatus(mesh,ents,n,status);dCHK(err); + owned = 0; + for (dInt i=0; ifs = fs; + err = dFSCreateDirichletVector(fs,&elp->d);dCHK(err); err = dFSCreateExpandedVector(fs,&elp->x);dCHK(err); err = VecDuplicate(elp->x,&elp->y);dCHK(err); @@ -291,8 +292,8 @@ static dErr EllipFunction(SNES dUNUSED snes,Vec gx,Vec gy,void *ctx) dErr err; dFunctionBegin; - err = dFSGlobalToExpandedBegin(fs,gx,INSERT_VALUES,elp->x);dCHK(err); - err = dFSGlobalToExpandedEnd(fs,gx,INSERT_VALUES,elp->x);dCHK(err); + err = dFSGlobalToExpandedBegin(fs,gx,dFS_INHOMOGENEOUS,elp->x);dCHK(err); + err = dFSGlobalToExpandedEnd(fs,gx,dFS_INHOMOGENEOUS,elp->x);dCHK(err); err = VecGetArray(elp->x,&x);dCHK(err); err = VecZeroEntries(elp->y);dCHK(err); err = VecGetArray(elp->y,&y);dCHK(err); @@ -313,6 +314,20 @@ static dErr EllipFunction(SNES dUNUSED snes,Vec gx,Vec gy,void *ctx) } err = dFSRestoreWorkspace(fs,__func__,&q,&jinv,&jw,&u,&v,&du,&dv);dCHK(err); err = dFSRestoreElements(fs,&n,&off,&rule,&efs,&geomoff,&geom);dCHK(err); +#if 0 + if (0) { + dMeshESH *sets; + dInt nsets; + err = dFSGetBoundarySets(fs,100,&nsets,&sets);dCHK(err); + for (dInt s=0; sx,&x);dCHK(err); err = VecRestoreArray(elp->y,&y);dCHK(err); err = VecZeroEntries(gy);dCHK(err); /* Necessary? */ @@ -335,8 +350,8 @@ static dErr EllipShellMatMult(Mat J,Vec gx,Vec gy) dFunctionBegin; err = MatShellGetContext(J,(void**)&elp);dCHK(err); fs = elp->fs; - err = dFSGlobalToExpandedBegin(fs,gx,INSERT_VALUES,elp->x);dCHK(err); - err = dFSGlobalToExpandedEnd(fs,gx,INSERT_VALUES,elp->x);dCHK(err); + err = dFSGlobalToExpandedBegin(fs,gx,dFS_HOMOGENEOUS,elp->x);dCHK(err); + err = dFSGlobalToExpandedEnd(fs,gx,dFS_HOMOGENEOUS,elp->x);dCHK(err); err = VecGetArray(elp->x,&x);dCHK(err); err = VecZeroEntries(elp->y);dCHK(err); err = VecGetArray(elp->y,&y);dCHK(err); @@ -456,8 +471,8 @@ static dErr EllipErrorNorms(Ellip elp,Vec gx,dReal errorNorms[static 3],dReal ge dFunctionBegin; err = dMemzero(errorNorms,3*sizeof(errorNorms));dCHK(err); err = dMemzero(gerrorNorms,3*sizeof(gerrorNorms));dCHK(err); - err = dFSGlobalToExpandedBegin(fs,gx,INSERT_VALUES,elp->x);dCHK(err); - err = dFSGlobalToExpandedEnd(fs,gx,INSERT_VALUES,elp->x);dCHK(err); + err = dFSGlobalToExpandedBegin(fs,gx,dFS_INHOMOGENEOUS,elp->x);dCHK(err); + err = dFSGlobalToExpandedEnd(fs,gx,dFS_INHOMOGENEOUS,elp->x);dCHK(err); err = VecGetArray(elp->x,&x);dCHK(err); err = dFSGetElements(fs,&n,&off,&rule,&efs,&geomoff,&geom);dCHK(err); err = dFSGetWorkspace(fs,__func__,&q,&jinv,&jw,&u,NULL,(dReal**)&du,NULL);dCHK(err); diff --git a/src/jacobi/impls/tensor/tensor.c b/src/jacobi/impls/tensor/tensor.c index 714ffeab..0fe03636 100644 --- a/src/jacobi/impls/tensor/tensor.c +++ b/src/jacobi/impls/tensor/tensor.c @@ -224,30 +224,24 @@ static dErr dJacobiGetNodeCount_Tensor(dUNUSED dJacobi jac,dInt count,const dEnt dFunctionReturn(0); } -static dErr dJacobiGetConstraintCount_Tensor(dUNUSED dJacobi jac,dInt nx,const dInt xi[],const dUNUSED dInt xs[],const dInt deg[],const struct dMeshAdjacency *ma,dInt nnz[],dInt pnnz[]) +static dErr dJacobiGetConstraintCount_Tensor(dUNUSED dJacobi jac,dInt nx,const dInt xi[],const dUNUSED dInt xs[],dInt dUNUSED dsplit, + const dInt dUNUSED is[],const dInt dUNUSED deg[],const struct dMeshAdjacency *ma, + dInt nnz[],dInt pnnz[],dInt dnnz[]) { - dInt i,j; dFunctionBegin; - for (i=0; itopo[i]) { - case dTOPO_HEX: - for (dInt i0=1; i0topo[i]); - } - } else { - for (j=xs[i]; jtopo[ei]) { + case dTOPO_HEX: + /* \todo Check whether adjacent entities are explicit or Dirichlet when determining the number of entries */ + /* \todo Handle nonconforming elements */ + for (dInt j=xs[i]; jtopo[i]); } } dFunctionReturn(0); @@ -296,18 +290,33 @@ static inline dErr dGeomPermQuadIndex(dInt perm,const dInt dim[],const dInt ij[2 dFunctionReturn(0); } -static dErr dJacobiAddConstraints_Tensor(dJacobi dUNUSED jac,dInt nx,const dInt xi[],const dInt xs[],const dInt is[],const dInt deg[],const struct dMeshAdjacency *ma,Mat C,Mat Cp) +/** Chooses whether to insert the value in the element assembly matrices or in the Dirichlet assembly matrix */ +static dErr PrivateMatSetValue(Mat E,Mat Ep,Mat Ed,dInt dsplit,dInt row,dInt col,MatScalar v,InsertMode imode) +{ + dErr err; + + dFunctionBegin; + if (col < dsplit) { + err = MatSetValue(E,row,col,v,imode);dCHK(err); + err = MatSetValue(Ep,row,col,v,imode);dCHK(err); + } else { + err = MatSetValue(Ed,row,col-dsplit,v,imode);dCHK(err); + } + dFunctionReturn(0); +} + +static dErr dJacobiAddConstraints_Tensor(dJacobi dUNUSED jac,dInt nx,const dInt xi[],const dInt xs[],dInt dsplit,const dInt is[],const dInt deg[],const struct dMeshAdjacency *ma,Mat matE,Mat matEp,Mat matEd) { dInt elem,i,j,k,l,e[12],eP[12],v[8]; dInt nrow,ncol,irow[10],icol[30]; - dScalar interp[30]; /* FIXME: check bounds */ + dScalar interp[30]; /* \todo implement nonconforming elements so that this array is actually used (instead of just the first entry) */ const dInt *f,*fP; const dInt *aI=ma->adjind,*aO=ma->adjoff,*aP=ma->adjperm; dErr err; dFunctionBegin; for (elem=0; elemtopo[ei]) { @@ -345,8 +354,7 @@ static dErr dJacobiAddConstraints_Tensor(dJacobi dUNUSED jac,dInt nx,const dInt for (i=0; i<8; i++) { /* Set vertices, always conforming until we have h-nonconforming meshes */ const dInt T[8][3] = {{0,0,0},{d0-1,0,0},{d0-1,d1-1,0},{0,d1-1,0}, {0,0,d2-1},{d0-1,0,d2-1},{d0-1,d1-1,d2-1},{0,d1-1,d2-1}}; - err = MatSetValue(C, xs[elem]+(T[i][0]*d1+T[i][1])*d2+T[i][2],is[v[i]],1,INSERT_VALUES);dCHK(err); - err = MatSetValue(Cp,xs[elem]+(T[i][0]*d1+T[i][1])*d2+T[i][2],is[v[i]],1,INSERT_VALUES);dCHK(err); + err = PrivateMatSetValue(matE,matEp,matEd,dsplit,xs[elem]+(T[i][0]*d1+T[i][1])*d2+T[i][2],is[v[i]],1,INSERT_VALUES);dCHK(err); } for (i=0; i<12; i++) { /* Set edges */ const struct {dInt start[3],incd,inci,end;} E[12] = { /* How to traverse the interior of the edge in forward order */ @@ -363,8 +371,7 @@ static dErr dJacobiAddConstraints_Tensor(dJacobi dUNUSED jac,dInt nx,const dInt case 1: icol[ncol++] = is[e[i]] - (j-(end-inci))/inci; break; /* traverse the edge in reverse */ } interp[0] = 1; - err = MatSetValues(C,nrow,irow,ncol,icol,interp,INSERT_VALUES);dCHK(err); - err = MatSetValues(Cp,nrow,irow,ncol,icol,interp,INSERT_VALUES);dCHK(err); + err = PrivateMatSetValue(matE,matEp,matEd,dsplit,irow[0],icol[0],interp[0],INSERT_VALUES);dCHK(err); } } for (i=0; i<6; i++) { /* Faces */ @@ -385,8 +392,7 @@ static dErr dJacobiAddConstraints_Tensor(dJacobi dUNUSED jac,dInt nx,const dInt err = dGeomPermQuadIndex(fP[i],faceDim,facejk,&faceIndex);dCHK(err); icol[ncol++] = is[f[i]] + faceIndex; interp[0] = 1; - err = MatSetValues(C, nrow,irow,ncol,icol,interp,INSERT_VALUES);dCHK(err); - err = MatSetValues(Cp,nrow,irow,ncol,icol,interp,INSERT_VALUES);dCHK(err); + err = PrivateMatSetValue(matE,matEp,matEd,dsplit,irow[0],icol[0],interp[0],INSERT_VALUES);dCHK(err); } } } @@ -397,8 +403,8 @@ static dErr dJacobiAddConstraints_Tensor(dJacobi dUNUSED jac,dInt nx,const dInt irow[0] = xs[elem] + (j*d1+k)*d2+l; icol[0] = is[ei] + ((j-1)*(d1-2)+(k-1))*(d2-2)+(l-1); interp[0] = 1; - err = MatSetValues(C, 1,irow,1,icol,interp,INSERT_VALUES);dCHK(err); - err = MatSetValues(Cp,1,irow,1,icol,interp,INSERT_VALUES);dCHK(err); + if (icol[0] >= dsplit) dERROR(1,"Should not happen"); + err = PrivateMatSetValue(matE,matEp,matEd,dsplit,irow[0],icol[0],interp[0],INSERT_VALUES);dCHK(err); } } } diff --git a/src/jacobi/interface/jacobi.c b/src/jacobi/interface/jacobi.c index 0d9d2385..bf83aa68 100644 --- a/src/jacobi/interface/jacobi.c +++ b/src/jacobi/interface/jacobi.c @@ -493,7 +493,21 @@ dErr dEFSApply(dEFS efs,const dReal mapdata[],dInt dofs,const dScalar in[],dScal dFunctionReturn(0); } -dErr dJacobiGetConstraintCount(dJacobi jac,dInt nx,const dInt xi[],const dInt xs[],const dInt deg[],const struct dMeshAdjacency *ma,dInt nnz[],dInt pnnz[]) +/** Get constraint counts for element assembly matrices. +* +* @param jac Jacobi object +* @param nx Number of entities in expanded vector +* @param xi Index of expanded entity in full adjacency (length \c nx) +* @param xs Offset in expanded vector of first node associated with each expanded entity (length \c nx+1) +* @param dsplit Local offset beyond which \a is refers to the Dirichlet vector +* @param is Starting index of interior nodes for each entity in \a ma +* @param deg Basis degree for each entity in \a ma +* @param ma MeshAdjacency object +* @param nnz Number of nonzeros per row of element assembly matrix +* @param pnnz Number of nonzeros per row of preconditioning element assembly matrix +* @param dnnz Numbor of nonzeros per row of dirichlet element assembly matrix +*/ +dErr dJacobiGetConstraintCount(dJacobi jac,dInt nx,const dInt xi[],const dInt xs[],dInt dsplit,const dInt is[],const dInt deg[],const struct dMeshAdjacency *ma,dInt nnz[],dInt pnnz[],dInt dnnz[]) { dErr err; @@ -501,15 +515,19 @@ dErr dJacobiGetConstraintCount(dJacobi jac,dInt nx,const dInt xi[],const dInt xs dValidHeader(jac,dJACOBI_COOKIE,1); dValidPointer(xi,3); dValidPointer(xs,4); - dValidPointer(deg,5); - dValidPointer(ma,6); - dValidPointer(nnz,7); - dValidPointer(pnnz,8); - err = (*jac->ops->GetConstraintCount)(jac,nx,xi,xs,deg,ma,nnz,pnnz);dCHK(err); + dValidPointer(is,6); + dValidPointer(deg,7); + dValidPointer(ma,8); + dValidPointer(nnz,9); + dValidPointer(pnnz,10); + dValidPointer(dnnz,11); + err = (*jac->ops->GetConstraintCount)(jac,nx,xi,xs,dsplit,is,deg,ma,nnz,pnnz,dnnz);dCHK(err); dFunctionReturn(0); } -dErr dJacobiAddConstraints(dJacobi jac,dInt nx,const dInt xi[],const dInt xs[],const dInt is[],const dInt deg[],const struct dMeshAdjacency *ma,Mat C,Mat Cp) +/** Actually assemble the element assembly matrices, see dJacobiGetConstraintCount() +*/ +dErr dJacobiAddConstraints(dJacobi jac,dInt nx,const dInt xi[],const dInt xs[],dInt dsplit,const dInt is[],const dInt deg[],const struct dMeshAdjacency *ma,Mat E,Mat Ep,Mat Ed) { dErr err; @@ -517,10 +535,12 @@ dErr dJacobiAddConstraints(dJacobi jac,dInt nx,const dInt xi[],const dInt xs[],c dValidHeader(jac,dJACOBI_COOKIE,1); dValidPointer(xi,3); dValidPointer(xs,4); - dValidPointer(deg,5); - dValidPointer(ma,6); - dValidHeader(C,MAT_COOKIE,7); - dValidHeader(Cp,MAT_COOKIE,8); - err = (*jac->ops->AddConstraints)(jac,nx,xi,xs,is,deg,ma,C,Cp);dCHK(err); + dValidPointer(is,6); + dValidPointer(deg,7); + dValidPointer(ma,8); + dValidHeader(E,MAT_COOKIE,9); + dValidHeader(Ep,MAT_COOKIE,10); + dValidHeader(Ep,MAT_COOKIE,11); + err = (*jac->ops->AddConstraints)(jac,nx,xi,xs,dsplit,is,deg,ma,E,Ep,Ed);dCHK(err); dFunctionReturn(0); }