Merge branch 'master' of https://github.com/ClimateGlobalChange/tempestremap

Author: paullric

src/ApplyOfflineMap.cpp

@@ -19,6 +19,8 @@
 #include "TempestRemapAPI.h"
 #include "OfflineMap.h"
 #include "netcdfcpp.h"
+#include "GridElements.h"
+#include "FiniteElementTools.h"
 
 #include <fstream>
 
@@ -273,14 +275,79 @@ try {
 		AnnounceStartBlock("Applying first offline map to data");
 	}
 */
+
+	DataArray3D<int> dataGLLNodesIn;
+	DataArray3D<double> dataGLLJacobianIn;
+	
+	DataArray3D<int> dataGLLNodesOut;
+	DataArray3D<double> dataGLLJacobianOut;
+	
+	DataArray3D<int> * pdataGLLNodesIn = NULL;
+			
+	DataArray3D<int> * pdataGLLNodesOut = NULL;
+							
+	Mesh meshOverlap;
+	Mesh meshSource;
+	Mesh meshTarget;
+	
+	Mesh * pmeshSource = NULL;
+	Mesh * pmeshOverlap = NULL;
+		
+		
+	if(optsApply.strInputMesh != ""){
+		
+		//If the source mesh is finite volume, we need the source mesh for local p bounds preservation
+		
+		meshSource.Read(optsApply.strInputMesh);
+		meshSource.ConstructReverseNodeArray();
+		meshSource.ConstructEdgeMap();
+		pmeshSource = &meshSource;
+		
+		
+		//If the source mesh is finite element, we need dataGLLNodes for local bounds preservation
+		if(optsApply.fgll){
+			
+			double dTotalAreaInput = meshSource.CalculateFaceAreas(optsApply.fContainsConcaveFaces);
+			double dNumericalAreaIn = GenerateMetaData(meshSource,optsApply.nPin,false,dataGLLNodesIn,dataGLLJacobianIn);
+			
+			pdataGLLNodesIn = &dataGLLNodesIn;
+			
+		}
+		
+	}
+	
+	if(optsApply.strOverlapMesh != ""){
+		
+		meshOverlap.Read(optsApply.strOverlapMesh);
+		meshOverlap.RemoveZeroEdges();
+		pmeshOverlap = &meshOverlap;
+		
+	}
+	
+	//If the target mesh is finite element, get dataGLLNodesOut
+	if(optsApply.strOutputMesh != ""){
+		
+		meshTarget.Read(optsApply.strOutputMesh);
+		meshTarget.ConstructReverseNodeArray();
+		meshTarget.ConstructEdgeMap();
+		
+		double dTotalAreaOutput = meshTarget.CalculateFaceAreas(optsApply.fContainsConcaveFaces);
+		double dNumericalAreaOut = GenerateMetaData(meshTarget,optsApply.nPout,false,dataGLLNodesOut,dataGLLJacobianOut);
+		
+		pdataGLLNodesOut = &dataGLLNodesOut;
+					
+	}
+		
 	// OfflineMap
 	AnnounceStartBlock("Loading offline map");
 
 	OfflineMap mapRemap;
 	mapRemap.Read(strInputMap);
 	mapRemap.SetFillValueOverrideDbl(optsApply.dFillValueOverride);
 	mapRemap.SetFillValueOverride(static_cast<float>(optsApply.dFillValueOverride));
-	mapRemap.SetEnforcementBounds(optsApply.strEnforceBounds);
+	mapRemap.SetEnforcementBounds(optsApply.strEnforceBounds,pmeshSource,pmeshOverlap,pdataGLLNodesIn,
+								  pdataGLLNodesOut,optsApply.nPin);
+	//mapRemap.SetEnforcementBounds(optsApply.strEnforceBounds);
 
 	AnnounceEndBlock("Done");
 

src/ApplyOfflineMapExe.cpp

@@ -56,11 +56,18 @@ int main(int argc, char** argv) {
 		//CommandLineString(strVariables2, "var2", "");
 		CommandLineString(optsApply.strNColName, "ncol_name", "ncol");
 		CommandLineString(optsApply.strEnforceBounds, "bounds", "");
+		CommandLineString(optsApply.strInputMesh,"mesh_in","");
+		CommandLineString(optsApply.strOutputMesh,"mesh_out","");
+		CommandLineString(optsApply.strOverlapMesh,"mesh_ov","");
+		CommandLineInt(optsApply.nPin, "in_np", 0);
+		CommandLineInt(optsApply.nPout, "out_np", 0);
 		CommandLineBool(optsApply.fOutputDouble, "out_double");
 		CommandLineString(optsApply.strPreserveVariables, "preserve", "");
 		CommandLineBool(optsApply.fPreserveAll, "preserveall");
 		CommandLineDouble(optsApply.dFillValueOverride, "fillvalue", 0.0);
 		CommandLineString(optsApply.strLogDir, "logdir", "");
+		CommandLineBool(optsApply.fContainsConcaveFaces,"concave");
+		CommandLineBool(optsApply.fgll,"gll");
 
 		ParseCommandLine(argc, argv);
 	EndCommandLine(argv)

src/CoordTransforms.h

@@ -337,6 +337,27 @@ inline void StereographicProjectionInv(
 
 ///////////////////////////////////////////////////////////////////////////////
 
+inline void GnomonicProjection(
+	double dLonRad0,
+	double dLatRad0,
+	double dLonRad,
+	double dLatRad,
+	double & dXs,
+	double & dYs
+) {
+	// Forward projection using equations (1)-(3)
+	// https://mathworld.wolfram.com/GnomonicProjection.html
+	
+	double dK = sin(dLatRad0) * sin(dLatRad) + cos(dLatRad0) * cos(dLatRad) * cos(dLonRad - dLonRad0);
+	
+	dXs = cos(dLatRad) * sin(dLonRad - dLonRad0)/dK;
+	
+	dYs = (cos(dLatRad0) * sin(dLatRad) - sin(dLatRad0) * cos(dLatRad) * cos(dLonRad - dLonRad0))/dK;
+		
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
 
 #endif // _COORDTRANSFORMS_H_