Find GDAL and set up the basic projections in stokesjako.c

CMakeLists.txt

@@ -6,6 +6,7 @@ list (APPEND CMAKE_MODULE_PATH "${Dohp_SOURCE_DIR}/CMake")
 # Normally PETSc is built with MPI, if not, use CC=mpicc, etc
 find_package (PETSc REQUIRED)
 find_package (ITAPS COMPONENTS MESH GEOM REL)
+find_package (GDAL)
 find_path (MEMCHECK_HEADER valgrind/memcheck.h)
 if (MEMCHECK_HEADER)
   set (dUSE_VALGRIND TRUE)

src/fs/tests/CMakeLists.txt

@@ -13,7 +13,10 @@ dohp_link_executable (ellip ellip.c)
 dohp_link_executable (bu bu.c)
 dohp_link_executable (cunit cunit.c)
 dohp_link_executable (elast elast.c elastexact.c)
-dohp_link_executable (stokes stokes.c stokesjako.c stokesexact.c)
+
+# Stokes needs to link GDAL too so dohp_link_executable does not work
+add_executable (stokes stokes.c stokesjako.c stokesexact.c)
+target_link_libraries (stokes dohp ${GDAL_LIBRARY})
 
 dohp_add_test (fs-mf-projection-v1 1 fs-ex1 -snes_mf -ksp_type minres -const_BDeg 5 -snes_monitor_short -ksp_converged_reason -dfs_ordering_type natural -proj_version 1)
 #dohp_add_test (fs-mf-projection-v2 1 fs-ex1 -snes_mf -ksp_type minres -const_BDeg 5 -snes_monitor_short -ksp_converged_reason -dfs_ordering_type natural -proj_version 2)

src/fs/tests/stokes.c

@@ -184,8 +184,6 @@ static dErr StokesSetFromOptions(Stokes stk)
     }
     err = PetscOptionsList("-stokes_case","Which sort of case to run","",StokesCaseList,scasename,scasename,sizeof(scasename),NULL);dCHK(err);
   } err = PetscOptionsEnd();dCHK(err);
-  err = StokesCaseSetType(stk->scase,scasename);dCHK(err);
-  err = StokesCaseSetFromOptions(stk->scase);dCHK(err);
 
   err = dMeshCreate(stk->comm,&mesh);dCHK(err);
   err = dMeshSetInFile(mesh,"dblock.h5m",NULL);dCHK(err);
@@ -264,6 +262,9 @@ static dErr StokesSetFromOptions(Stokes stk)
   }
   err = dJacobiDestroy(&jac);dCHK(err);
   err = dMeshDestroy(&mesh);dCHK(err);
+
+  err = StokesCaseSetType(stk->scase,scasename);dCHK(err);
+  err = StokesCaseSetFromOptions(stk->scase);dCHK(err);
   dFunctionReturn(0);
 }
 

src/fs/tests/stokesimpl.h

@@ -19,6 +19,7 @@ struct _n_StokesCase {
   dErr (*destroy)(StokesCase);
   struct StokesRheology rheo;
   dReal gravity;
+  dReal bbox[3][2];
   dBool reality; // The "solution" is just a guess or boundary conditions
   void *data;
 };

src/fs/tests/stokesjako.c

@@ -1,4 +1,10 @@
 #include "stokesimpl.h"
+#include <ogr_srs_api.h>
+
+#define dOCHK(err) do {                                                 \
+    if (PetscUnlikely((err) != OGRERR_NONE))                             \
+      dERROR(PETSC_COMM_SELF,PETSC_ERR_LIB,"OGR library routine failed with error code %d",err); \
+  } while (0)
 
 // Trivial gravity model
 static void StokesCaseSolution_Gravity(StokesCase dUNUSED scase,const dReal dUNUSED x[3],dScalar u[],dScalar du[],dScalar *p,dScalar dp[])
@@ -26,15 +32,19 @@ static dErr StokesCaseCreate_Gravity(StokesCase scase)
 
 // A real implementation
 typedef struct {
-  int placeholder;
+  OGRSpatialReferenceH utmref; // UTM zone 22N, the coordinates with Roman's surface elevation
+  OGRSpatialReferenceH llref;  // Longitude-Latitude, for the CReSIS bed elevation
+  OGRSpatialReferenceH ianref; // Ian Joughin's surface velocity
+  OGRCoordinateTransformationH fromll;  // Convert from LonLat to current
+  OGRCoordinateTransformationH fromian; // Convert from Ian's projection to current
 } StokesCase_Jako;
 
 static dErr StokesCaseSolution_Jako(StokesCase scase,const dReal x[3],dScalar u[],dScalar du[],dScalar *p,dScalar dp[])
 {                               /* Defines inhomogeneous Dirichlet boundary conditions */
-  StokesCase_Jako *jako = scase->data;
+  dUNUSED StokesCase_Jako *jako = scase->data;
 
   dFunctionBegin;
-  if (x[0] > 580000. || jako->placeholder) {
+  if (x[0] > 580000.) {
     u[0] = -200;
     u[1] = -100;
     u[2] = 0;
@@ -54,6 +64,54 @@ static void StokesCaseSolution_Jako_Void(StokesCase scase,const dReal x[3],dScal
   dErr err;
   err = StokesCaseSolution_Jako(scase,x,u,du,p,dp);CHKERRV(err);
 }
+static dErr JakoViewWKT(OGRSpatialReferenceH ref,const char *name,PetscViewer viewer)
+{
+  dErr err;
+  OGRErr oerr;
+  char *wkt;
+
+  dFunctionBegin;
+  oerr = OSRExportToPrettyWkt(ref,&wkt,0);dOCHK(oerr);
+  err = PetscViewerASCIIPrintf(viewer,"WKT %s: %s\n\n",name,wkt);dCHK(err);
+  OGRFree(wkt);
+  dFunctionReturn(0);
+}
+static dErr StokesCaseView_Jako(StokesCase scase,PetscViewer viewer)
+{
+  StokesCase_Jako *jako = scase->data;
+  dErr err;
+
+  dFunctionBegin;
+  err = PetscViewerASCIIPrintf(viewer,"StokesCase: Jakobshavn\n");dCHK(err);
+  err = JakoViewWKT(jako->utmref,"UTM",viewer);dCHK(err);
+  err = JakoViewWKT(jako->llref,"LonLat",viewer);dCHK(err);
+  err = JakoViewWKT(jako->ianref,"Ian",viewer);dCHK(err);
+  dFunctionReturn(0);
+}
+static dErr StokesCaseSetUp_Jako(StokesCase scase)
+{
+  StokesCase_Jako *jako = scase->data;
+  dErr err;
+  OGRErr oerr;
+
+  dFunctionBegin;
+  jako->utmref = OSRNewSpatialReference(NULL);
+  oerr = OSRSetProjCS(jako->utmref,"UTM 22N (WGS84)");dOCHK(oerr);
+  oerr = OSRSetWellKnownGeogCS(jako->utmref,"WGS84");dOCHK(oerr);
+  oerr = OSRSetUTM(jako->utmref,22,1);dOCHK(oerr);
+
+  jako->llref  = OSRNewSpatialReference(NULL);
+  oerr = OSRSetWellKnownGeogCS(jako->llref,"WGS84");dOCHK(oerr);
+
+  jako->ianref = OSRNewSpatialReference(NULL);
+  oerr = OSRSetProjCS(jako->ianref,"Stereographic Greenland (WGS84)");dOCHK(oerr);
+  oerr = OSRSetWellKnownGeogCS(jako->ianref,"WGS84");dOCHK(oerr);
+  oerr = OSRSetStereographic(jako->ianref,70.0,-45.0,100.0,-217.75e3,-2302.0e3);dOCHK(oerr);
+
+  err = StokesCaseView_Jako(scase,PETSC_VIEWER_STDOUT_WORLD);dCHK(err);
+  dFunctionReturn(0);
+}
+
 static dErr StokesCaseSetFromOptions_Jako(StokesCase scase)
 {
   char fname[256] = "unknown";
@@ -64,6 +122,7 @@ static dErr StokesCaseSetFromOptions_Jako(StokesCase scase)
     err = PetscOptionsString("-jako_surface_velocity","File to read surface velocity from (assume same projection as model, e.g. UTM)","",fname,fname,sizeof(fname),&flg);dCHK(err);
     if (!flg) dERROR(scase->comm,PETSC_ERR_USER,"User must provide surface velocity file with -jako_surface_velocity FILENAME");
   } err = PetscOptionsTail();dCHK(err);
+  err = StokesCaseSetUp_Jako(scase);dCHK(err);
   dFunctionReturn(0);
 }
 static dErr StokesCaseDestroy_Jako(StokesCase scase)
@@ -72,7 +131,8 @@ static dErr StokesCaseDestroy_Jako(StokesCase scase)
   dErr err;
 
   dFunctionBegin;
-  jako->placeholder = 0;
+  OSRDestroySpatialReference(jako->utmref); jako->utmref = NULL;
+  OSRDestroySpatialReference(jako->llref);  jako->llref  = NULL;
   err = dFree(scase->data);dCHK(err);
   dFunctionReturn(0);
 }