Skeleton for getting geometry. I'm very unsure that this is the right…

… way.

include/dohpfs.h

@@ -70,8 +70,6 @@ EXTERN dErr dFSExpandedToGlobalEnd(dFS,Vec,InsertMode,Vec);
 EXTERN dErr dFSExpandedToDirichlet(dFS,Vec,InsertMode,Vec);
 EXTERN dErr dFSClosureToGlobal(dFS,Vec,Vec,InsertMode,dFSHomogeneousMode);
 EXTERN dErr dFSGlobalToClosure(dFS,Vec,Vec,InsertMode,dFSHomogeneousMode);
-EXTERN dErr dFSGetClosureCoordinates(dFS,Vec*);
-EXTERN dErr dFSGetGeometryVector(dFS,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*);
@@ -96,6 +94,9 @@ EXTERN dErr dFSRotationApplyLocal(dFSRotation,Vec,dFSRotateMode,dFSHomogeneousMo
 EXTERN dErr dFSSetRotation(dFS,dFSRotation);
 EXTERN dErr dFSGetRotation(dFS,dFSRotation*);
 
+EXTERN dErr dFSGetClosureCoordinates(dFS,Vec*);
+EXTERN dErr dFSGetGeometryVector(dFS,Vec*);
+
 /** The Q1 stuff doesn't really belong here, but it is used at the same level of abstraction and I'm too lazy to
 * separate it out yet.  This macro is a rather dirty way to avoid a bunch of code duplication without much runtime cost.
 * When there's time, I should make this more elegant.

src/CMakeLists.txt

@@ -18,6 +18,7 @@ list (APPEND Dohp_SRCS
   fs/interface/fs.c
   fs/interface/fsreg.c
   fs/interface/fsrot.c
+  fs/interface/fsgeom.c
   fs/interface/q1.c
   fs/impls/cont/cont.c)
 

src/fs/interface/fsgeom.c

@@ -0,0 +1,51 @@
+#include "private/fsimpl.h"
+
+/** Get coordinates for every node in closure
+*
+* @param fs Function space
+* @param inx the new vector with block size 3 and the same number of blocks as the closure vector
+**/
+dErr dFSGetClosureCoordinates(dFS fs,Vec *inx)
+{
+  dErr err;
+  dInt n,bs;
+  const VecType type;
+  Vec x;
+
+  dFunctionBegin;
+  dValidHeader(fs,DM_COOKIE,1);
+  dValidPointer(inx,2);
+  *inx = 0;
+  /* Poor man's duplicate with different block size */
+  err = VecGetLocalSize(fs->gc,&n);dCHK(err);
+  err = VecGetBlockSize(fs->gc,&bs);dCHK(err);
+  err = VecGetType(fs->gc,&type);dCHK(err);
+  err = VecCreate(((dObject)fs)->comm,&x);dCHK(err);
+  err = VecSetSizes(x,3*n/bs,PETSC_DETERMINE);dCHK(err);
+  err = VecSetBlockSize(x,3);dCHK(err);
+  err = VecSetType(x,type);dCHK(err);
+  dERROR(1,"not implemented");
+  dFunctionReturn(0);
+}
+
+/** Get displacements of every point in geometry vector.
+*
+* @param fs Function space
+* @param ingeom the new geometry vector with block size 3 and same number of blocks as the closure vector
+*
+* @note I don't yet have a clear view of how this should work.  If we have mesh motion with interesting boundary
+* conditions, we probably need a bonified function space for the geometry.
+*
+* @note It is important to get geometry associated with boundary sets because it will frequently be projected against
+* the boundary.
+**/
+dErr dFSGetGeometryVector(dFS fs,Vec *ingeom)
+{
+
+  dFunctionBegin;
+  dValidHeader(fs,DM_COOKIE,1);
+  dValidPointer(ingeom,2);
+  *ingeom = 0;
+  dERROR(1,"not implemented");
+  dFunctionReturn(0);
+}