Vec caas (#103)

* Non-conservative monotone integrated and non-integrated remapping schemes.

* Implementation of CAAS with vector of inputs.

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/FiniteVolumeTools.cpp

@@ -291,6 +291,7 @@ void GetTriangleThatContainsPoint(
 
 				// If this is a new Face, check whether it contains the point
 				if (setAllFaces.find(iNewFace) == setAllFaces.end()) {
+
 					if(DoesTriangleContainPoint(mesh,iNewFace,dX,dY)){
 
 						iFaceFinal = iNewFace;
@@ -377,6 +378,7 @@ void GetFaceThatContainsPoint(
 
 				// If this is a new Face, check whether it contains the point
 				if (setAllFaces.find(iNewFace) == setAllFaces.end()) {
+
 					if(DoesFaceContainPoint(mesh,iNewFace,dX,dY,dZ)){
 
 						iFaceFinal = iNewFace;
@@ -768,6 +770,7 @@ void TriangleLineIntersection(
 
 	_ASSERT(dCoeffs.GetRows() == 3);
 
+
 	//Setup up columns of 3x3 matrix
 	DataArray2D<double> dInterpMat(3,3);
 

src/LinearRemapFV.cpp

@@ -1298,12 +1298,11 @@ void LinearRemapIntegratedGeneralizedBarycentric(
 	// Get SparseMatrix representation of the OfflineMap
 	SparseMatrix<double> & smatMap = mapRemap.GetSparseMatrix();
 
-	//Dual mesh
 	Mesh meshInputDual = meshInput;
 
 	//Construct dual mesh
 	Dual(meshInputDual);
-		
+
 	//Reverse node array
 	meshInputDual.ConstructReverseNodeArray();	
 
@@ -1595,7 +1594,7 @@ void LinearRemapGeneralizedBarycentric(
 
 	//Construct dual mesh
 	Dual(meshInputDual);
-		
+
 	//Reverse node array
 	meshInputDual.ConstructReverseNodeArray();	
 

src/LinearRemapFV.h

@@ -84,7 +84,6 @@ void LinearRemapIntegratedTriangulation(
 ///		Generate the OfflineMap for non-integrated remapping from finite volumes to finite
 ///		volumes using generalized barycentric coordinates using the dual mesh
 ///	</summary>
-
 void LinearRemapIntegratedGeneralizedBarycentric(
 	const Mesh & meshInput,
 	const Mesh & meshOutput,
@@ -98,7 +97,6 @@ void LinearRemapIntegratedGeneralizedBarycentric(
 ///		Generate the OfflineMap for integrated remapping from finite volumes to finite
 ///		volumes using generalized barycentric coordinates using the dual mesh
 ///	</summary>
-
 void LinearRemapGeneralizedBarycentric(
 	const Mesh & meshInput,
 	const Mesh & meshOutput,
@@ -112,7 +110,6 @@ void LinearRemapGeneralizedBarycentric(
 ///		Generate the OfflineMap for non-integrated remapping from finite volumes to finite
 ///		volumes using a triangulation of the source mesh.
 ///	</summary>
-
 void LinearRemapTriangulation(
 	const Mesh & meshInput,
 	const Mesh & meshOutput,