diff --git a/src/ApplyOfflineMap.cpp b/src/ApplyOfflineMap.cpp index c6fad2f..904880b 100644 --- a/src/ApplyOfflineMap.cpp +++ b/src/ApplyOfflineMap.cpp @@ -19,6 +19,8 @@ #include "TempestRemapAPI.h" #include "OfflineMap.h" #include "netcdfcpp.h" +#include "GridElements.h" +#include "FiniteElementTools.h" #include @@ -273,6 +275,69 @@ try { AnnounceStartBlock("Applying first offline map to data"); } */ + + DataArray3D dataGLLNodesIn; + DataArray3D dataGLLJacobianIn; + + DataArray3D dataGLLNodesOut; + DataArray3D dataGLLJacobianOut; + + DataArray3D * pdataGLLNodesIn = NULL; + + DataArray3D * 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"); @@ -280,7 +345,9 @@ try { mapRemap.Read(strInputMap); mapRemap.SetFillValueOverrideDbl(optsApply.dFillValueOverride); mapRemap.SetFillValueOverride(static_cast(optsApply.dFillValueOverride)); - mapRemap.SetEnforcementBounds(optsApply.strEnforceBounds); + mapRemap.SetEnforcementBounds(optsApply.strEnforceBounds,pmeshSource,pmeshOverlap,pdataGLLNodesIn, + pdataGLLNodesOut,optsApply.nPin); + //mapRemap.SetEnforcementBounds(optsApply.strEnforceBounds); AnnounceEndBlock("Done"); diff --git a/src/ApplyOfflineMapExe.cpp b/src/ApplyOfflineMapExe.cpp index 58e94ef..80af8ed 100644 --- a/src/ApplyOfflineMapExe.cpp +++ b/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) diff --git a/src/CoordTransforms.h b/src/CoordTransforms.h index f638a8e..3dda545 100644 --- a/src/CoordTransforms.h +++ b/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_ diff --git a/src/FiniteVolumeTools.cpp b/src/FiniteVolumeTools.cpp index 078d67b..34191ed 100644 --- a/src/FiniteVolumeTools.cpp +++ b/src/FiniteVolumeTools.cpp @@ -231,6 +231,444 @@ void GetAdjacentFaceVectorByEdge( /////////////////////////////////////////////////////////////////////////////// +void GetTriangleThatContainsPoint( + const Mesh & mesh, + int iFaceInitial, + int & iFaceFinal, + double dX, + double dY +) { + // Ensure the ReverseNodeArray has been constructed + if (mesh.edgemap.size() == 0) { + _EXCEPTIONT("EdgeMap is required"); + } + + // First check if point is in initial face + if (DoesTriangleContainPoint(mesh,iFaceInitial,dX,dY)){ + + iFaceFinal = iFaceInitial; + return; + + } + + // Set of all Faces + std::set setAllFaces; + + // Set of current Faces + std::set setCurrentFaces; + + // Insert initial Face + setAllFaces.insert(iFaceInitial); + setCurrentFaces.insert(iFaceInitial); + + while (setAllFaces.size() < mesh.faces.size()) { + + // Set of Faces to examine next + std::set setNextFaces; + + // Loop through all Faces adjacent to Faces in setCurrentFaces + std::set::const_iterator iterCurrentFace = setCurrentFaces.begin(); + for (; iterCurrentFace != setCurrentFaces.end(); iterCurrentFace++) { + + const Face & faceCurrent = mesh.faces[*iterCurrentFace]; + for (int i = 0; i < faceCurrent.edges.size(); i++) { + const FacePair & facepair = + mesh.edgemap.find(faceCurrent.edges[i])->second; + + // New face index + int iNewFace; + if (facepair[0] == *iterCurrentFace) { + iNewFace = facepair[1]; + } else if (facepair[1] == *iterCurrentFace) { + iNewFace = facepair[0]; + } else { + _EXCEPTIONT("Logic error"); + } + + if (iNewFace == InvalidFace) { + continue; + } + + // 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; + return; + + } + else{ + + setAllFaces.insert(iNewFace); + setNextFaces.insert(iNewFace); + + } + } + } + } + + setCurrentFaces = setNextFaces; + } + _EXCEPTIONT("Unable to find a triangle that contains the point"); +} + +/////////////////////////////////////////////////////////////////////////////// + +void GetFaceThatContainsPoint( + const Mesh & mesh, + int iFaceInitial, + int & iFaceFinal, + double dX, + double dY, + double dZ + +){ + + // Ensure the ReverseNodeArray has been constructed + if (mesh.edgemap.size() == 0) { + _EXCEPTIONT("EdgeMap is required"); + } + + // First check if point is in initial face + if (DoesFaceContainPoint(mesh,iFaceInitial,dX,dY,dZ)){ + + iFaceFinal = iFaceInitial; + return; + + } + + // Set of all Faces + std::set setAllFaces; + + // Set of current Faces + std::set setCurrentFaces; + + // Insert initial Face + setAllFaces.insert(iFaceInitial); + setCurrentFaces.insert(iFaceInitial); + + while (setAllFaces.size() < mesh.faces.size()) { + + // Set of Faces to examine next + std::set setNextFaces; + + // Loop through all Faces adjacent to Faces in setCurrentFaces + std::set::const_iterator iterCurrentFace = setCurrentFaces.begin(); + for (; iterCurrentFace != setCurrentFaces.end(); iterCurrentFace++) { + + const Face & faceCurrent = mesh.faces[*iterCurrentFace]; + for (int i = 0; i < faceCurrent.edges.size(); i++) { + const FacePair & facepair = + mesh.edgemap.find(faceCurrent.edges[i])->second; + + // New face index + int iNewFace; + if (facepair[0] == *iterCurrentFace) { + iNewFace = facepair[1]; + } else if (facepair[1] == *iterCurrentFace) { + iNewFace = facepair[0]; + } else { + _EXCEPTIONT("Logic error"); + } + + if (iNewFace == InvalidFace) { + continue; + } + + // 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; + return; + + } + else{ + + setAllFaces.insert(iNewFace); + setNextFaces.insert(iNewFace); + + } + } + } + } + + setCurrentFaces = setNextFaces; + } + _EXCEPTIONT("Unable to find a face that contains the point"); +} + +/////////////////////////////////////////////////////////////////////////////// + +bool DoesFaceContainPoint( + const Mesh & mesh, + int iFace, + double dX, + double dY, + double dZ + +){ + + //Convert sample point to lat/lon coordinates + + double dLonRad0 = 0.0; + double dLatRad0 = 0.0; + + XYZtoRLL_Rad(dX,dY,dZ,dLonRad0,dLatRad0); + + int iEdges = mesh.faces[iFace].edges.size(); + + NodeVector nodesPlane; + + double dCenterX = 0; + double dCenterY = 0; + + // Project face nodes onto plane tangent to the sphere at the sample point whose + // coordinates are dX, dY, dZ + for (int i = 0; i < iEdges; i++){ + + Node nodeCurrent = mesh.nodes[mesh.faces[iFace][i]]; + + double dLonRad = 0.0; + double dLatRad = 0.0; + + //Convert to lat/lon coordinates + + XYZtoRLL_Rad(nodeCurrent.x,nodeCurrent.y,nodeCurrent.z,dLonRad,dLatRad); + + //Project on tangent plane at sample point + + double dGX = 0.0; + double dGY = 0.0; + + GnomonicProjection(dLonRad0,dLatRad0,dLonRad,dLatRad,dGX,dGY); + + Node nodeI(dGX,dGY,0.0); + + nodesPlane.push_back(nodeI); + + //add contributions to planar face centroid + dCenterX += dGX; + dCenterY += dGY; + + } + + //Compute coordinates of planar face centroid and turn it into a Node + + dCenterX /= iEdges; + dCenterY /= iEdges; + + Node nodeCenter(dCenterX,dCenterY,0); + + //loop over all edges + for (int i = 0; i < iEdges; i++){ + + Node nodeI = nodesPlane[i]; + + Node nodeIPlusOne = nodesPlane[(i+1)%iEdges]; + + Node nodeEdge(nodeIPlusOne.x - nodeI.x, + nodeIPlusOne.y - nodeI.y, + 0.0); + + Node nodeEdgeNormal(-nodeEdge.y, nodeEdge.x, 0.0); + + + Node nodeCenterMinusI(nodeCenter.x - nodeI.x, + nodeCenter.y - nodeI.y, + 0.0); + + //This is the projected sampled point (whose coordinates are zero in the gnomonic plane) + //minus nodeI + + Node nodeQ(-nodeI.x,-nodeI.y,0.0); + + //If the signs of the dot products are different, than the sample point is outside the polygon + + if((DotProduct(nodeCenterMinusI,nodeEdgeNormal) > 0 && DotProduct(nodeQ,nodeEdgeNormal) < 0) || + (DotProduct(nodeCenterMinusI,nodeEdgeNormal) < 0 && DotProduct(nodeQ,nodeEdgeNormal) > 0)){ + + return false; + + } + + } + + return true; + +} + +/////////////////////////////////////////////////////////////////////////////// + +bool DoesFaceContainPoint( + const NodeVector & nodesP, + double dX, + double dY, + double dZ + +){ + + //Convert sample point to lat/lon coordinates + + double dLonRad0 = 0.0; + double dLatRad0 = 0.0; + + XYZtoRLL_Rad(dX,dY,dZ,dLonRad0,dLatRad0); + + int iEdges = nodesP.size(); + + NodeVector nodesPlane; + + double dCenterX = 0; + double dCenterY = 0; + + // Project face nodes onto plane tangent to the sphere at the sample point whose + // coordinates are dX, dY, dZ + for (int i = 0; i < iEdges; i++){ + + Node nodeCurrent = nodesP[i]; + + double dLonRad = 0.0; + double dLatRad = 0.0; + + //Convert to lat/lon coordinates + + XYZtoRLL_Rad(nodeCurrent.x,nodeCurrent.y,nodeCurrent.z,dLonRad,dLatRad); + + //Project on tangent plane at sample point + + double dGX = 0.0; + double dGY = 0.0; + + GnomonicProjection(dLonRad0,dLatRad0,dLonRad,dLatRad,dGX,dGY); + + Node nodeI(dGX,dGY,0.0); + + nodesPlane.push_back(nodeI); + + //add contributions to planar face centroid + dCenterX += dGX; + dCenterY += dGY; + + } + + //Compute coordinates of planar face centroid and turn it into a Node + + dCenterX /= iEdges; + dCenterY /= iEdges; + + Node nodeCenter(dCenterX,dCenterY,0); + + //loop over all edges + for (int i = 0; i < iEdges; i++){ + + Node nodeI = nodesPlane[i]; + + Node nodeIPlusOne = nodesPlane[(i+1)%iEdges]; + + Node nodeEdge(nodeIPlusOne.x - nodeI.x, + nodeIPlusOne.y - nodeI.y, + 0.0); + + Node nodeEdgeNormal(-nodeEdge.y, nodeEdge.x, 0.0); + + + Node nodeCenterMinusI(nodeCenter.x - nodeI.x, + nodeCenter.y - nodeI.y, + 0.0); + + //This is the projected sampled point (whose coordinates are zero in the gnomonic plane) + //minus nodeI + + Node nodeQ(-nodeI.x,-nodeI.y,0.0); + + //If the signs of the dot products are different, than the sample point is outside the polygon + + if((DotProduct(nodeCenterMinusI,nodeEdgeNormal) > 0.0 && DotProduct(nodeQ,nodeEdgeNormal) < 0.0) || + (DotProduct(nodeCenterMinusI,nodeEdgeNormal) < 0.0 && DotProduct(nodeQ,nodeEdgeNormal) > 0.0)){ + + return false; + + } + + } + + return true; + +} + +/////////////////////////////////////////////////////////////////////////////// + +void BarycentricCoordinates( + const Mesh & mesh, + int iFace, + double dX, + double dY, + double & dA, + double & dB +){ + + Face face = mesh.faces[iFace]; + + // The input face has to be a triangle + if(face.edges.size() != 3){ + _EXCEPTIONT("The input face must be a triangle"); + } + + Node nodeV1 = mesh.nodes[face[0]]; + Node nodeV2 = mesh.nodes[face[1]]; + Node nodeV3 = mesh.nodes[face[2]]; + + double dX1 = nodeV1.x; + double dY1 = nodeV1.y; + + double dX2 = nodeV2.x; + double dY2 = nodeV2.y; + + double dX3 = nodeV3.x; + double dY3 = nodeV3.y; + + double dDenom = (dY2 - dY3)*(dX1 - dX3) + (dX3 - dX2)*(dY1 - dY3); + + if( abs(dDenom) < ReferenceTolerance ){ + + _EXCEPTIONT("Points are close to colinear"); + + } + + dA = ((dY2 - dY3)*(dX - dX3) + (dX3 - dX2)*(dY - dY3))/dDenom; + + dB = ((dY3 - dY1)*(dX - dX3) + (dX1 - dX3)*(dY - dY3))/dDenom; + +} + +/////////////////////////////////////////////////////////////////////////////// + +bool DoesTriangleContainPoint( + const Mesh & mesh, + int iFace, + double dX, + double dY +){ + + double dA; + double dB; + + BarycentricCoordinates(mesh,iFace,dX,dY,dA,dB); + + if ((0.0 <= dA) && (0.0 <= dB) && (dA + dB <= 1+ReferenceTolerance)){ + return true; + } + else{ + return false; + } + +} + +/////////////////////////////////////////////////////////////////////////////// + void GetAdjacentFaceVectorByNode( const Mesh & mesh, int iFaceInitial, @@ -310,6 +748,252 @@ void GetAdjacentFaceVectorByNode( /////////////////////////////////////////////////////////////////////////////// +void MatVectorMult(const DataArray2D & dMat, + DataArray1D & dRHS, + DataArray1D & dOutput) +{ + + dOutput(0) = (dMat(0,0))*dRHS(0) + (dMat(0,1))*dRHS(1) + (dMat(0,2))*dRHS(2); + dOutput(1) = (dMat(1,0))*dRHS(0) + (dMat(1,1))*dRHS(1) + (dMat(1,2))*dRHS(2); + dOutput(2) = (dMat(2,0))*dRHS(0) + (dMat(2,1))*dRHS(1) + (dMat(2,2))*dRHS(2); + +} + +/////////////////////////////////////////////////////////////////////////////// + +void TriangleLineIntersection( + Node & nodeQ, + NodeVector & nodesP, + DataArray1D & dCoeffs, + double & dCond +) { + + _ASSERT(dCoeffs.GetRows() == 3); + + + //Setup up columns of 3x3 matrix + DataArray2D dInterpMat(3,3); + + dInterpMat(0,0) = nodeQ.x; + dInterpMat(1,0) = nodeQ.y; + dInterpMat(2,0) = nodeQ.z; + + dInterpMat(0,1) = nodesP[1].x - nodesP[0].x; + dInterpMat(1,1) = nodesP[1].y - nodesP[0].y; + dInterpMat(2,1) = nodesP[1].z - nodesP[0].z; + dInterpMat(0,2) = nodesP[2].x - nodesP[0].x; + dInterpMat(1,2) = nodesP[2].y - nodesP[0].y; + dInterpMat(2,2) = nodesP[2].z - nodesP[0].z; + + int m = 3; + int n = 3; + int lda = 3; + int info; + + DataArray1D iPIV; + iPIV.Allocate(3); + + DataArray1D dWork; + dWork.Allocate(3); + + int lWork = 3; + + //Column sums for A and A inverse + DataArray1D dColSumA(3); + + for (int j = 0; j < 3; j++) { + + for (int k = 0; k < 3; k++) { + + dColSumA[j] += fabs(dInterpMat(j,k)); + + } + + } + + DataArray1D dColSumAInv(3); + + dgetrf_(&m, &n, &(dInterpMat(0,0)), + &lda, &(iPIV[0]), &info); + + + dgetri_(&n, &(dInterpMat(0,0)), + &lda, &(iPIV[0]), &(dWork[0]), &lWork, &info); + + + //A inverse column sums + for (int j = 0; j < 3; j++) { + + for (int k = 0; k < 3; k++) { + + dColSumAInv[j] += fabs(dInterpMat(j,k)); + + } + + } + + + //max column sums of A and A inverse + double dMaxColSumA = dColSumA[0]; + + double dMaxColSumAInv = dColSumAInv[0]; + + for (int k = 1; k < 3; k++) { + + if (dColSumA[k] > dMaxColSumA) { + dMaxColSumA = dColSumA[k]; + } + if (dColSumAInv[k] > dMaxColSumAInv) { + dMaxColSumAInv = dColSumAInv[k]; + } + } + + dCond = dMaxColSumAInv * dMaxColSumA; + + if (info < 0) { + _EXCEPTION1("dgetrf_ reports matrix had an illegal value (%i)", info); + } + if (info > 0) { + Announce("WARNING: Singular matrix detected in fit (likely colinear elements)"); + } + + //Set right hand side of linear system + + DataArray1D dRHS(3); + + dRHS(0) = nodeQ.x - nodesP[0].x; + dRHS(1) = nodeQ.y - nodesP[0].y; + dRHS(2) = nodeQ.z - nodesP[0].z; + + //Solve linear system with matrix multiply + + MatVectorMult(dInterpMat, dRHS, dCoeffs); + + +} + +/////////////////////////////////////////////////////////////////////////////// + +void NewtonQuadrilateral( + Node & nodeQ, + NodeVector & nodesP, + DataArray1D & dCoeffs, + bool & fConverged +) { + + _ASSERT(dCoeffs.GetRows() == 3); + + int iMaxIterations = 100; + + //This algorithm is essentially the one that is used in ESMF + + // Four corners of quadrilateral + Node nodeQ0 = nodesP[0]; + Node nodeQ1 = nodesP[1]; + Node nodeQ2 = nodesP[2]; + Node nodeQ3 = nodesP[3]; + + // Jacobian + DataArray2D dJacobian(3,3); + + DataArray1D A(3), B(3), C(3), D(3), E(3), F(3); + + A[0] = nodeQ0.x - nodeQ1.x + nodeQ2.x - nodeQ3.x; + A[1] = nodeQ0.y - nodeQ1.y + nodeQ2.y - nodeQ3.y; + A[2] = nodeQ0.z - nodeQ1.z + nodeQ2.z - nodeQ3.z; + + B[0] = nodeQ1.x - nodeQ0.x; + B[1] = nodeQ1.y - nodeQ0.y; + B[2] = nodeQ1.z - nodeQ0.z; + + C[0] = nodeQ3.x-nodeQ0.x; + C[1] = nodeQ3.y-nodeQ0.y; + C[2] = nodeQ3.z-nodeQ0.z; + + D[0] = nodeQ.x; + D[1] = nodeQ.y; + D[2] = nodeQ.z; + + E[0] = nodeQ0.x - nodeQ.x; + E[1] = nodeQ0.y - nodeQ.y; + E[2] = nodeQ0.z - nodeQ.z; + + for (int i = 0; i < iMaxIterations; i++){ + + double dA = dCoeffs[0]; + double dB = dCoeffs[1]; + double dC = dCoeffs[2]; + + // Calculate Value of function at X + F[0] = dA*dB*A[0] + dA*B[0] + dB*C[0] + dC*D[0] + E[0]; + F[1] = dA*dB*A[1] + dA*B[1] + dB*C[1] + dC*D[1] + E[1]; + F[2] = dA*dB*A[2] + dA*B[2] + dB*C[2] + dC*D[2] + E[2]; + + // If we're close enough to 0.0 then exit + if (F[0]*F[0]+F[1]*F[1]+F[2]*F[2] < 1.0E-15) { + + fConverged = true; + break; + + } + + // Construct Jacobian + dJacobian(0,0) = A[0]*dB + B[0]; + dJacobian(1,0) = A[1]*dB + B[1]; + dJacobian(2,0) = A[2]*dB + B[2]; + + dJacobian(0,1) = A[0]*dA + C[0]; + dJacobian(1,1) = A[1]*dA + C[1]; + dJacobian(2,1) = A[2]*dA + C[2]; + + dJacobian(0,2) = D[0]; + dJacobian(1,2) = D[1]; + dJacobian(2,2) = D[2]; + + //Invert jacobian + int m = 3; + int n = 3; + int lda = 3; + int info; + + DataArray1D iPIV; + iPIV.Allocate(3); + + DataArray1D dWork; + dWork.Allocate(3); + + int lWork = 3; + + //LU decomposition + dgetrf_(&m, &n, &(dJacobian(0,0)), + &lda, &(iPIV[0]), &info); + + //Invert LU decomposition + dgetri_(&n, &(dJacobian(0,0)), + &lda, &(iPIV[0]), &(dWork[0]), &lWork, &info); + + if (info != 0) { + _EXCEPTIONT("Mass matrix inversion error"); + } + + DataArray1D dDeltaX(3); + + //Multiply solution vector by inverse jacobian + MatVectorMult(dJacobian, F, dDeltaX); + + //Update solution + dCoeffs[0] = dCoeffs[0] - dDeltaX[0]; + dCoeffs[1] = dCoeffs[1] - dDeltaX[1]; + dCoeffs[2] = dCoeffs[2] - dDeltaX[2]; + + + } + + +} + +/////////////////////////////////////////////////////////////////////////////// + void BuildIntegrationArray( const Mesh & meshInput, const Mesh & meshOverlap, @@ -1012,4 +1696,3 @@ void InvertFitArray_LeastSquares( } /////////////////////////////////////////////////////////////////////////////// - diff --git a/src/FiniteVolumeTools.h b/src/FiniteVolumeTools.h index 2f873a3..f9556b6 100644 --- a/src/FiniteVolumeTools.h +++ b/src/FiniteVolumeTools.h @@ -115,6 +115,86 @@ void GetAdjacentFaceVectorByEdge( /////////////////////////////////////////////////////////////////////////////// +/// +/// Find the triangle that contains a two-dimensional point +/// +void GetTriangleThatContainsPoint( + const Mesh & mesh, + int iFaceInitial, + int & iFaceFinal, + double dX, + double dY +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Find the triangle that contains a two-dimensional point +/// +void GetFaceThatContainsPoint( + const Mesh & mesh, + int iFaceInitial, + int & iFaceFinal, + double dX, + double dY, + double dZ +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Returns true if a face contains a given point +/// +bool DoesFaceContainPoint( + const Mesh & mesh, + int iFace, + double dX, + double dY, + double dZ + +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Returns true if a face contains a given point +/// +bool DoesFaceContainPoint( + const NodeVector & nodesP, + double dX, + double dY, + double dZ + +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Barycentric coordinates of a two-dimensional point +/// +void BarycentricCoordinates( + const Mesh & mesh, + int iFace, + double dX, + double dY, + double & dA, + double & dB +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Determine if a triangle contains a point +/// +bool DoesTriangleContainPoint( + const Mesh & mesh, + int iFace, + double dX, + double dY +); + +/////////////////////////////////////////////////////////////////////////////// + /// /// Get a vector of adjacent faces using vertex information to determine /// adjacencies. @@ -128,6 +208,41 @@ void GetAdjacentFaceVectorByNode( /////////////////////////////////////////////////////////////////////////////// +/// +/// Multiply a 3-dimensional vector by a 3x3 matrix. +/// +void MatVectorMult( + const DataArray2D & dMat, + DataArray1D & dRHS, + DataArray1D & dOutput +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Intersection of a point and triangle. +/// +void TriangleLineIntersection( + Node & nodeQ, + NodeVector & nodesP, + DataArray1D & dCoeffs, + double & dCond +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Find the coordinates of a point in the basis of a quadrilateral. +/// +void NewtonQuadrilateral( + Node & nodeQ, + NodeVector & nodesP, + DataArray1D & dCoeffs, + bool & fConverged +); + +/////////////////////////////////////////////////////////////////////////////// + /// /// Build the integration array, an operator that integrates a polynomial /// reconstruction (specified as a vector of polynomial coefficients) to diff --git a/src/GenerateICOMesh.cpp b/src/GenerateICOMesh.cpp index cb79a47..670d1f4 100644 --- a/src/GenerateICOMesh.cpp +++ b/src/GenerateICOMesh.cpp @@ -346,120 +346,6 @@ void GenerateIcosahedralQuadGrid( /////////////////////////////////////////////////////////////////////////////// -void Dual( - Mesh & mesh -) { - const int EdgeCountHexagon = 6; - - // Generate ReverseNodeArray - mesh.ConstructReverseNodeArray(); - - // Backup Nodes and Faces - NodeVector nodesOld = mesh.nodes; - FaceVector facesOld = mesh.faces; - - mesh.nodes.clear(); - mesh.faces.clear(); - - // Generate new Node array - for (int i = 0; i < facesOld.size(); i++) { - Node node; - for (int j = 0; j < facesOld[i].edges.size(); j++) { - node.x += nodesOld[facesOld[i][j]].x; - node.y += nodesOld[facesOld[i][j]].y; - node.z += nodesOld[facesOld[i][j]].z; - } - - node.x /= static_cast(facesOld[i].edges.size()); - node.y /= static_cast(facesOld[i].edges.size()); - node.z /= static_cast(facesOld[i].edges.size()); - - double dMag = node.Magnitude(); - - node.x /= dMag; - node.y /= dMag; - node.z /= dMag; - - mesh.nodes.push_back(node); - } - - // Generate new Face array - for (int i = 0; i < nodesOld.size(); i++) { - const int nEdges = mesh.revnodearray[i].size(); - - Face face(EdgeCountHexagon); - Face faceTemp(EdgeCountHexagon); - - int ixNode = 0; - std::set::const_iterator iter = mesh.revnodearray[i].begin(); - for (; iter != mesh.revnodearray[i].end(); iter++) { - faceTemp.SetNode(ixNode, *iter); - ixNode++; - } - - // Reorient Faces - Node nodeCentral = nodesOld[i]; - Node node0 = mesh.nodes[faceTemp[0]] - nodeCentral; - - Node nodeCross = CrossProduct(mesh.nodes[faceTemp[0]], nodeCentral); - - double dNode0Mag = node0.Magnitude(); - - // Determine the angles about the central Node of each Face Node - std::vector dAngles; - dAngles.resize(faceTemp.edges.size()); - dAngles[0] = 0.0; - - for (int j = 1; j < nEdges; j++) { - Node nodeDiff = mesh.nodes[faceTemp[j]] - nodeCentral; - double dNodeDiffMag = nodeDiff.Magnitude(); - - double dSide = DotProduct(nodeCross, nodeDiff); - - double dDotNorm = - DotProduct(node0, nodeDiff) / (dNode0Mag * dNodeDiffMag); - - double dAngle; - if (dDotNorm > 1.0) { - dDotNorm = 1.0; - } - - dAngles[j] = acos(dDotNorm); - - if (dSide > 0.0) { - dAngles[j] = - dAngles[j] + 2.0 * M_PI; - } - } - - // Orient each Face by putting Nodes in order of increasing angle - double dCurrentAngle = 0.0; - face.SetNode(0, faceTemp[0]); - for (int j = 1; j < nEdges; j++) { - int ixNextNode = 1; - double dNextAngle = 2.0 * M_PI; - - for (int k = 1; k < nEdges; k++) { - if ((dAngles[k] > dCurrentAngle) && (dAngles[k] < dNextAngle)) { - ixNextNode = k; - dNextAngle = dAngles[k]; - } - } - - face.SetNode(j, faceTemp[ixNextNode]); - dCurrentAngle = dNextAngle; - } - - // Fill in missing edges - for (int j = nEdges; j < EdgeCountHexagon; j++) { - face.SetNode(j, face[nEdges-1]); - } - - mesh.faces.push_back(face); - } -} - -/////////////////////////////////////////////////////////////////////////////// - extern "C" int GenerateICOMesh(Mesh& mesh, int nResolution, bool fDual, std::string strOutputFile, std::string strOutputFormat) { diff --git a/src/GenerateOfflineMap.cpp b/src/GenerateOfflineMap.cpp index 4a7c0a6..89bcd1f 100644 --- a/src/GenerateOfflineMap.cpp +++ b/src/GenerateOfflineMap.cpp @@ -25,6 +25,8 @@ #include "SparseMatrix.h" #include "STLStringHelper.h" #include "NetCDFUtilities.h" +#include "triangle.h" +#include "FiniteVolumeTools.h" #include "OverlapMesh.h" #include "OfflineMap.h" @@ -437,6 +439,18 @@ try { meshOverlap, (optsAlg.fMonotone)?(1):(optsAlg.nPin), mapRemap); + + //To run the non-conservative monotone remapping schemes, uncomment the following and replace with + //LinearRemapGeneralizedBarycentric, LinearRemapTriangulation, LinearRemapBilinear, + //LinearRemapIntegratedGeneralizedBarycentric, LinearRemapIntegratedTriangulation, + //or LinearRemapIntegratedBilinear + + //LinearRemapIntegratedGeneralizedBarycentric( + //meshSource, + //meshTarget, + //meshOverlap, + //mapRemap); + } // Finite volume input / Finite element output diff --git a/src/GridElements.cpp b/src/GridElements.cpp index 9a46f70..323b7f0 100644 --- a/src/GridElements.cpp +++ b/src/GridElements.cpp @@ -2923,3 +2923,110 @@ void ConvexifyMesh( /////////////////////////////////////////////////////////////////////////////// +void Dual( + Mesh & mesh +) { + + // Generate ReverseNodeArray + mesh.ConstructReverseNodeArray(); + + // Backup Nodes and Faces + NodeVector nodesOld = mesh.nodes; + FaceVector facesOld = mesh.faces; + + mesh.nodes.clear(); + mesh.faces.clear(); + + // Generate new Node array + for (int i = 0; i < facesOld.size(); i++) { + Node node; + for (int j = 0; j < facesOld[i].edges.size(); j++) { + node.x += nodesOld[facesOld[i][j]].x; + node.y += nodesOld[facesOld[i][j]].y; + node.z += nodesOld[facesOld[i][j]].z; + } + + node.x /= static_cast(facesOld[i].edges.size()); + node.y /= static_cast(facesOld[i].edges.size()); + node.z /= static_cast(facesOld[i].edges.size()); + + double dMag = node.Magnitude(); + + node.x /= dMag; + node.y /= dMag; + node.z /= dMag; + + mesh.nodes.push_back(node); + } + + // Generate new Face array + for (int i = 0; i < nodesOld.size(); i++) { + const int nEdges = mesh.revnodearray[i].size(); + + Face face(mesh.revnodearray[i].size()); + Face faceTemp(mesh.revnodearray[i].size()); + + int ixNode = 0; + std::set::const_iterator iter = mesh.revnodearray[i].begin(); + for (; iter != mesh.revnodearray[i].end(); iter++) { + faceTemp.SetNode(ixNode, *iter); + ixNode++; + } + + // Reorient Faces + Node nodeCentral = nodesOld[i]; + Node node0 = mesh.nodes[faceTemp[0]] - nodeCentral; + + Node nodeCross = CrossProduct(mesh.nodes[faceTemp[0]], nodeCentral); + + double dNode0Mag = node0.Magnitude(); + + // Determine the angles about the central Node of each Face Node + std::vector dAngles; + dAngles.resize(faceTemp.edges.size()); + dAngles[0] = 0.0; + + for (int j = 1; j < nEdges; j++) { + Node nodeDiff = mesh.nodes[faceTemp[j]] - nodeCentral; + double dNodeDiffMag = nodeDiff.Magnitude(); + + double dSide = DotProduct(nodeCross, nodeDiff); + + double dDotNorm = + DotProduct(node0, nodeDiff) / (dNode0Mag * dNodeDiffMag); + + double dAngle; + if (dDotNorm > 1.0) { + dDotNorm = 1.0; + } + + dAngles[j] = acos(dDotNorm); + + if (dSide > 0.0) { + dAngles[j] = - dAngles[j] + 2.0 * M_PI; + } + } + + // Orient each Face by putting Nodes in order of increasing angle + double dCurrentAngle = 0.0; + face.SetNode(0, faceTemp[0]); + for (int j = 1; j < nEdges; j++) { + int ixNextNode = 1; + double dNextAngle = 2.0 * M_PI; + + for (int k = 1; k < nEdges; k++) { + if ((dAngles[k] > dCurrentAngle) && (dAngles[k] < dNextAngle)) { + ixNextNode = k; + dNextAngle = dAngles[k]; + } + } + + face.SetNode(j, faceTemp[ixNextNode]); + dCurrentAngle = dNextAngle; + } + + mesh.faces.push_back(face); + } +} + +/////////////////////////////////////////////////////////////////////////////// diff --git a/src/GridElements.h b/src/GridElements.h index 6497649..371bea2 100644 --- a/src/GridElements.h +++ b/src/GridElements.h @@ -1136,5 +1136,11 @@ inline Node InterpolateQuadrilateralNode( /////////////////////////////////////////////////////////////////////////////// +void Dual( + Mesh & mesh +); + +/////////////////////////////////////////////////////////////////////////////// + #endif diff --git a/src/LinearRemapFV.cpp b/src/LinearRemapFV.cpp index 95703e8..12e5dc0 100644 --- a/src/LinearRemapFV.cpp +++ b/src/LinearRemapFV.cpp @@ -20,18 +20,20 @@ #include "GridElements.h" #include "OfflineMap.h" #include "FiniteElementTools.h" -#include "GaussQuadrature.h" #include "GaussLobattoQuadrature.h" #include "TriangularQuadrature.h" #include "MeshUtilitiesFuzzy.h" #include "OverlapMesh.h" +#include "triangle.h" +#include "kdtree.h" +#include "DataArray3D.h" #include "Announce.h" #include "MathHelper.h" -#include "kdtree.h" #include #include +#include /////////////////////////////////////////////////////////////////////////////// @@ -729,6 +731,2039 @@ void LinearRemapFVtoFVInvDist( /////////////////////////////////////////////////////////////////////////////// +void LinearRemapIntegratedTriangulation( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +) { + + + // Verify ReverseNodeArray has been calculated + if (meshInput.edgemap.size() == 0) { + _EXCEPTIONT("EdgeMap has not been calculated for meshInput"); + } + + // Order of triangular quadrature rule + const int TriQuadRuleOrder = 4; + + // Triangular quadrature rule + TriangularQuadratureRule triquadrule(TriQuadRuleOrder); + + const DataArray2D & dG = triquadrule.GetG(); + const DataArray1D & dW = triquadrule.GetW(); + + // Get SparseMatrix representation of the OfflineMap + SparseMatrix & smatMap = mapRemap.GetSparseMatrix(); + + // Array of global indices + std::vector > vecGlobalIndexI(6); + + // Vector storing meshes of each panel + std::vector vecMesh(6); + + // kd-tree for each panel + std::vector vecKDTreePanelI(6); + + // Arrays of panel boundaries in lat/lon coordinates (radians) + DataArray2D dPanelLat(6,2); + DataArray2D dPanelLon(6,2); + + // Width of buffer zone for each panel (degrees) + double dBuff = 30; + + // Define equatorial panel boundaries + for (int i = 0; i < 4; i++){ + + dPanelLon(i,0) = 90*M_PI*i/180; //Left boundary + dPanelLon(i,1) = 90*M_PI*(i+1)/180; //Right Boundary + dPanelLat(i,1) = 45*M_PI/180; //Upper Boundary + dPanelLat(i,0) = -45*M_PI/180; //Lower Boundary + + } + + // Define polar panels + + dPanelLon(4,0) = 0; dPanelLon(4,1) = 2*M_PI; + dPanelLat(4,0) = 45*M_PI/180; dPanelLat(4,1) = M_PI/2; + + dPanelLon(5,0) = 0; dPanelLon(5,1) = 2*M_PI; + dPanelLat(5,1) = -45*M_PI/180; dPanelLat(5,0) = -M_PI/2; + + + // Generate the Delaunay triangulation for each of the 6 panels + for (int i = 0; i < 6; i++){ + + // Array storing the coordinates of the face centroids of panel i + std::vector> dPanelCentroid(2); + + // Coordinates of tangent plane point of tangency for panel i + double dLatTan; + double dLonTan; + + if ((0 <= i) && (i <= 3)){ + + dLatTan = 0; + dLonTan = (i+1)*45*M_PI/180 + 45*i*M_PI/180; + } + else if (i == 4){ + + dLatTan = M_PI/2; + dLonTan = 0; + + } + else{ + + dLonTan = 0; + dLatTan = -M_PI/2; + + } + + // Number of source faces centers for panel i; this is the size + // of in.pointlist + int nNodes = 0; + + for (int j = 0; j < meshInput.faces.size(); j++){ + + // Get coordinates of centroid of current face + + double dLonRad0; + double dLatRad0; + + const Face & face = meshInput.faces[j]; + + Node node = GetFaceCentroid(face,meshInput.nodes); + + node = node.Normalized(); + //node.Normalized(); + + XYZtoRLL_Rad(node.x,node.y,node.z,dLonRad0,dLatRad0); + + // Determine if centroid is in panel i plus buffer zone + bool fPanelContainsCentroid = false; + + if ((0 <= i) && (i <= 3)){ + + if ((dPanelLat(i,0) - dBuff*M_PI/180 <= dLatRad0) && (dLatRad0 <= dPanelLat(i,1)+dBuff*M_PI/180)){ + + if ( i == 0 ){ + + if ( ((2*M_PI - dBuff*M_PI/180 <= dLonRad0) && (dLonRad0 <= 2*M_PI)) || + ((dPanelLon(i,0) <= dLonRad0) && (dLonRad0 <= dPanelLon(i,1) + dBuff*M_PI/180)) ) { + + fPanelContainsCentroid = true; + } + } + + else if ( i == 3 ){ + + if ( ((0 <= dLonRad0) && (dLonRad0 <= 0 + dBuff*M_PI/180)) || + ((dPanelLon(i,0) - dBuff*M_PI/180 <= dLonRad0) && (dLonRad0 <= dPanelLon(i,1) + dBuff*M_PI/180)) ){ + + fPanelContainsCentroid = true; + + } + + } + + else { + + if ( (dPanelLon(i,0) - dBuff*M_PI/180 <= dLonRad0) && (dLonRad0 <= dPanelLon(i,1) + dBuff*M_PI/180)){ + + fPanelContainsCentroid = true; + } + + + } + + } + } + else if (i == 4){ + + if ((dPanelLat(i,0) - dBuff*M_PI/180 <= dLatRad0) && (dLatRad0 <= dPanelLat(i,1))){ + + fPanelContainsCentroid = true; + + } + + } + + else { + + if ((dPanelLat(i,0) <= dLatRad0) && (dLatRad0 <= dPanelLat(i,1) + dBuff*M_PI/180)){ + fPanelContainsCentroid = true; + + } + + } + + if(fPanelContainsCentroid){ + + vecGlobalIndexI[i].push_back(j); + nNodes++; + + // Gnomonic projection of centroid coordinates + double dGX; + double dGY; + GnomonicProjection(dLonTan,dLatTan,dLonRad0,dLatRad0,dGX,dGY); + + // Add Gnomonic coordinates to vector + dPanelCentroid[0].push_back(dGX); + dPanelCentroid[1].push_back(dGY); + + } + } + + // Structures for Delaunay triangulation + + struct triangulateio in, out, vorout; + + in.numberofpoints = nNodes; + in.numberofpointattributes = 0; + in.numberofsegments = 0; + in.numberofholes = 0; + in.numberofregions = 0; + in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL)); + in.segmentlist = (int *) malloc(in.numberofsegments * 2 * sizeof(int));; + in.pointattributelist = (REAL *) NULL; + in.pointmarkerlist = (int *) NULL; + in.trianglelist = (int *) NULL; + in.triangleattributelist = (REAL *) NULL; + in.neighborlist = (int *) NULL; + in.segmentmarkerlist = (int *) NULL; + in.edgelist = (int *) NULL; + in.edgemarkerlist = (int *) NULL; + + // initialize data structure for output triangulation + out.pointlist = (REAL *) NULL; + out.pointattributelist = (REAL *) NULL; + out.pointmarkerlist = (int *) NULL; + out.trianglelist = (int *) NULL; + out.triangleattributelist = (REAL *) NULL; + out.neighborlist = (int *) NULL; + out.segmentlist = (int *) NULL; + out.segmentmarkerlist = (int *) NULL; + out.edgelist = (int *) NULL; + out.edgemarkerlist = (int *) NULL; + + // initialize data structure for output Voronoi diagram (unused) + vorout.pointlist = (REAL *) NULL; + vorout.pointattributelist = (REAL *) NULL; + vorout.edgelist = (int *) NULL; + vorout.normlist = (REAL *) NULL; + + // Add points to in.pointlist + + for (int j = 0; j < nNodes; j++){ + + in.pointlist[2*j+0] = dPanelCentroid[0][j]; + in.pointlist[2*j+1] = dPanelCentroid[1][j]; + + } + + // Compute Delaunay triangulation. Use the 'c' option so that + // the convex hull is included + + char options[256] ="cjzenYQ"; + triangulate(options, &in, &out, &vorout); + + // Convert to mesh object by building face vector + for (int j = 0; j < out.numberoftriangles; j++){ + + Face newFace(3); + + newFace.SetNode(0, out.trianglelist[3*j+0]); + newFace.SetNode(1, out.trianglelist[3*j+1]); + newFace.SetNode(2, out.trianglelist[3*j+2]); + vecMesh[i].faces.push_back(newFace); + + } + + // Build node vector. Note that the Delaunay triangulation preserves + // point indexing. + for (int j = 0; j < out.numberofpoints; j++){ + + Node node(out.pointlist[2*j+0],out.pointlist[2*j+1],0); + vecMesh[i].nodes.push_back(node); + + } + + vecMesh[i].ConstructReverseNodeArray(); + vecMesh[i].RemoveCoincidentNodes(); + vecMesh[i].RemoveZeroEdges(); + vecMesh[i].ConstructEdgeMap(); + + free(in.pointlist); + free(out.pointlist); + free(out.pointattributelist); + free(out.pointmarkerlist); + free(out.trianglelist); + free(out.triangleattributelist); + free(out.neighborlist); + free(out.segmentlist); + free(out.segmentmarkerlist); + free(out.edgelist); + free(out.edgemarkerlist); + + } + + // Construct kd-tree for each panel + + for (int i = 0; i < 6; i++){ + + vecKDTreePanelI[i] = kd_create(3); + + for (int j = 0; j < vecMesh[i].faces.size(); j++){ + + Face face = vecMesh[i].faces[j]; + + Node nodeCenter = GetFaceCentroid(face, vecMesh[i].nodes); + + kd_insert3( + vecKDTreePanelI[i], + nodeCenter.x, + nodeCenter.y, + nodeCenter.z, + (void*)(&(vecMesh[i].faces[j]))); + + } + + } + + // Overlap face index + int ixOverlap = 0; + + // Loop through all source faces + for (int ixFirst = 0; ixFirst < meshInput.faces.size(); ixFirst++) { + + // Output every 1000 overlap elements + if (ixFirst % 1000 == 0) { + Announce("Element %i/%i", ixFirst, meshInput.faces.size()); + } + + // This Face + const Face & faceFirst = meshInput.faces[ixFirst]; + + // Find the set of Faces that overlap faceFirst + int ixOverlapBegin = ixOverlap; + int ixOverlapEnd = ixOverlapBegin; + + for (; ixOverlapEnd < meshOverlap.faces.size(); ixOverlapEnd++) { + if (meshOverlap.vecSourceFaceIx[ixOverlapEnd] != ixFirst) { + break; + } + } + + int nOverlapFaces = ixOverlapEnd - ixOverlapBegin; + + // Loop through all overlap faces associated with this source face + for (int j = 0; j < nOverlapFaces; j++) { + int iTargetFace = meshOverlap.vecTargetFaceIx[ixOverlap + j]; + + const Face & faceOverlap = meshOverlap.faces[ixOverlap + j]; + + int nSubTriangles = faceOverlap.edges.size() - 2; + + // Jacobian at each quadrature point + DataArray2D dQuadPtWeight(nSubTriangles, dW.GetRows()); + DataArray2D dQuadPtNodes(nSubTriangles, dW.GetRows()); + + // Compute quadrature area of each overlap face + double dQuadratureArea = 0.0; + + for (int k = 0; k < nSubTriangles; k++) { + for (int p = 0; p < dW.GetRows(); p++) { + + dQuadPtWeight(k,p) = + CalculateSphericalTriangleJacobianBarycentric( + meshOverlap.nodes[faceOverlap[0]], + meshOverlap.nodes[faceOverlap[k+1]], + meshOverlap.nodes[faceOverlap[k+2]], + dG(p,0), dG(p,1), + &(dQuadPtNodes(k,p))); + + dQuadPtWeight(k,p) *= dW[p]; + + dQuadratureArea += dQuadPtWeight(k,p); + } + } + + // Loop through all sub-triangles of this overlap Face + for (int k = 0; k < nSubTriangles; k++) { + + // Loop through all quadrature nodes on this overlap Face + for (int p = 0; p < dW.GetRows(); p++) { + + // Get quadrature node and pointwise Jacobian + const Node & nodeQ = dQuadPtNodes(k,p); + + // Get lat/lon coordinates of nodeQ + double dLatRad; + double dLonRad; + + XYZtoRLL_Rad(nodeQ.x,nodeQ.y,nodeQ.z,dLonRad,dLatRad); + + // Determine the panel where nodeQ is located + int iPanelIndex; + + for (int i = 0; i < 6; i++){ + + if ((dPanelLat(i,0) <= dLatRad) && (dLatRad <= dPanelLat(i,1)) && + (dPanelLon(i,0) <= dLonRad) && (dLonRad <= dPanelLon(i,1))){ + + iPanelIndex = i; + break; + + } + + } + // Compute Gnomonic projection onto the corresponding plane + + double dGX,dGY; + double dLatTan; + double dLonTan; + + if ((0 <= iPanelIndex) && (iPanelIndex <= 3)){ + + dLatTan = 0; + dLonTan = (iPanelIndex+1)*45*M_PI/180 + 45*iPanelIndex*M_PI/180; + + } + else if (iPanelIndex == 4){ + + dLatTan = M_PI/2; + dLonTan = 0; + + } + else{ + + dLonTan = 0; + dLatTan = -M_PI/2; + + } + + GnomonicProjection(dLonTan,dLatTan,dLonRad,dLatRad,dGX,dGY); + + // Get point closest to dGX and dGY + kdres * kdresTarget = + kd_nearest3( + vecKDTreePanelI[iPanelIndex], + dGX, + dGY, + 0.0); + + Face * pFace = (Face *)(kd_res_item_data(kdresTarget)); + + int iNearestFace = pFace - &(vecMesh[iPanelIndex].faces[0]); + + // Find triangle that contains the Gnomonic projection of nodeQ. + // This is the face whose local index is iFaceFinal + + int iFaceFinal; + + double dA,dB; + + BarycentricCoordinates(vecMesh[iPanelIndex],iNearestFace,dGX,dGY,dA,dB); + + GetTriangleThatContainsPoint(vecMesh[iPanelIndex],iNearestFace,iFaceFinal,dGX,dGY); + + // Get global indices of the vertices of this triangle and + // calculate corresponding centroids + + std::vector vecContributingFaceI(3); + DataArray2D dataContributingCentroids(3,3); + + // Indices on the local mesh of the containing triangle + + int iLocalVertex1 = vecMesh[iPanelIndex].faces[iFaceFinal][0]; + int iLocalVertex2 = vecMesh[iPanelIndex].faces[iFaceFinal][1]; + int iLocalVertex3 = vecMesh[iPanelIndex].faces[iFaceFinal][2]; + + // Indices on the source mesh of the containing triangle + + int iGlobalVertex1 = vecGlobalIndexI[iPanelIndex][iLocalVertex1]; + int iGlobalVertex2 = vecGlobalIndexI[iPanelIndex][iLocalVertex2]; + int iGlobalVertex3 = vecGlobalIndexI[iPanelIndex][iLocalVertex3]; + + vecContributingFaceI[0] = iGlobalVertex1; + vecContributingFaceI[1] = iGlobalVertex2; + vecContributingFaceI[2] = iGlobalVertex3; + + // The centroids of the source mesh are the vertices of the containing triangle + + for (int i = 0; i < 3; i++){ + + Face faceCurrent = meshInput.faces[vecContributingFaceI[i]]; + + Node nodeCenter = GetFaceCentroid(faceCurrent,meshInput.nodes); + + nodeCenter = nodeCenter.Normalized(); + + dataContributingCentroids(i,0) = nodeCenter.x; + dataContributingCentroids(i,1) = nodeCenter.y; + dataContributingCentroids(i,2) = nodeCenter.z; + + + } + + // Vector of areas of subtriangles + std::vector vecSubAreas(3); + + // Area of triangle is obtained by adding up areas of the three + // subtriangles that are formed by the sample point + + double dTriangleArea = 0; + + for (int i = 0; i < 3; i++){ + + Face faceCurrent(3); + + faceCurrent.SetNode(0,0); + faceCurrent.SetNode(1,1); + faceCurrent.SetNode(2,2); + + NodeVector nodesCurrent; + + nodesCurrent.push_back(nodeQ); + + Node node1(dataContributingCentroids((i+1)%3,0),dataContributingCentroids((i+1)%3,1), + dataContributingCentroids((i+1)%3,2)); + + Node node2(dataContributingCentroids((i+2)%3,0),dataContributingCentroids((i+2)%3,1), + dataContributingCentroids((i+2)%3,2)); + + nodesCurrent.push_back(node1); + nodesCurrent.push_back(node2); + + vecSubAreas[i] = CalculateFaceArea(faceCurrent,nodesCurrent); + + dTriangleArea += vecSubAreas[i]; + + } + + // Calculate vector of weights + std::vector vecContributingFaceWeights(3); + + for (int i = 0; i < 3; i++){ + + vecContributingFaceWeights[i] = vecSubAreas[i]/dTriangleArea; + + } + + // Contribution of this quadrature point to the map + for (int i = 0; i < vecContributingFaceI.size(); i++){ + + smatMap(iTargetFace, vecContributingFaceI[i]) += + vecContributingFaceWeights[i] + * dQuadPtWeight(k,p) + * meshOverlap.vecFaceArea[ixOverlap + j] + / dQuadratureArea + / meshOutput.vecFaceArea[iTargetFace]; + + } + } + } + } + + // Increment the current overlap index + ixOverlap += nOverlapFaces; + } + +} + +/////////////////////////////////////////////////////////////////////////////// + +void LinearRemapIntegratedGeneralizedBarycentric( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +) { + + + // Verify ReverseNodeArray has been calculated + if (meshInput.edgemap.size() == 0) { + _EXCEPTIONT("EdgeMap has not been calculated for meshInput"); + } + + // Order of triangular quadrature rule + const int TriQuadRuleOrder = 4; + + // Triangular quadrature rule + TriangularQuadratureRule triquadrule(TriQuadRuleOrder); + + const DataArray2D & dG = triquadrule.GetG(); + const DataArray1D & dW = triquadrule.GetW(); + + // Get SparseMatrix representation of the OfflineMap + SparseMatrix & smatMap = mapRemap.GetSparseMatrix(); + + Mesh meshInputDual = meshInput; + + //Construct dual mesh + Dual(meshInputDual); + + //Reverse node array + meshInputDual.ConstructReverseNodeArray(); + + //Construct edge map + meshInputDual.ConstructEdgeMap(); + + //kd-tree of dual mesh centers + + kdtree * kdTarget = kd_create(3); + + // Vector of centers of the source mesh + for (int i = 0; i < meshInputDual.faces.size(); i++){ + + const Face & face = meshInputDual.faces[i]; + + Node nodeCentroid = GetFaceCentroid(face, meshInputDual.nodes); + + nodeCentroid = nodeCentroid.Normalized(); + + kd_insert3( + kdTarget, + nodeCentroid.x, + nodeCentroid.y, + nodeCentroid.z, + (void*)(&(meshInputDual.faces[i]))); + + } + + + // Overlap face index + int ixOverlap = 0; + + // Loop through all source faces + for (int ixFirst = 0; ixFirst < meshInput.faces.size(); ixFirst++) { + + // Output every 1000 overlap elements + if (ixFirst % 1000 == 0) { + Announce("Element %i/%i", ixFirst, meshInput.faces.size()); + } + + // This Face + const Face & faceFirst = meshInput.faces[ixFirst]; + + // Find the set of Faces that overlap faceFirst + int ixOverlapBegin = ixOverlap; + int ixOverlapEnd = ixOverlapBegin; + + for (; ixOverlapEnd < meshOverlap.faces.size(); ixOverlapEnd++) { + if (meshOverlap.vecSourceFaceIx[ixOverlapEnd] != ixFirst) { + break; + } + } + + int nOverlapFaces = ixOverlapEnd - ixOverlapBegin; + + // Loop through all overlap faces associated with this source face + for (int j = 0; j < nOverlapFaces; j++) { + int iTargetFace = meshOverlap.vecTargetFaceIx[ixOverlap + j]; + + const Face & faceOverlap = meshOverlap.faces[ixOverlap + j]; + + int nSubTriangles = faceOverlap.edges.size() - 2; + + // Jacobian at each quadrature point + DataArray2D dQuadPtWeight(nSubTriangles, dW.GetRows()); + DataArray2D dQuadPtNodes(nSubTriangles, dW.GetRows()); + + // Compute quadrature area of each overlap face + double dQuadratureArea = 0.0; + + for (int k = 0; k < nSubTriangles; k++) { + for (int p = 0; p < dW.GetRows(); p++) { + + dQuadPtWeight(k,p) = + CalculateSphericalTriangleJacobianBarycentric( + meshOverlap.nodes[faceOverlap[0]], + meshOverlap.nodes[faceOverlap[k+1]], + meshOverlap.nodes[faceOverlap[k+2]], + dG(p,0), dG(p,1), + &(dQuadPtNodes(k,p))); + + dQuadPtWeight(k,p) *= dW[p]; + + dQuadratureArea += dQuadPtWeight(k,p); + } + } + + // Loop through all sub-triangles of this overlap Face + for (int k = 0; k < nSubTriangles; k++) { + + // Loop through all quadrature nodes on this overlap Face + for (int p = 0; p < dW.GetRows(); p++) { + + // Get quadrature node and pointwise Jacobian + const Node & nodeQ = dQuadPtNodes(k,p); + + //Get the dual mesh face whose center is nearest to the sample point + + kdres * kdresTarget = + kd_nearest3( + kdTarget, + nodeQ.x, + nodeQ.y, + nodeQ.z); + + Face * pFace = (Face *)(kd_res_item_data(kdresTarget)); + + int iNearestFace = pFace - &(meshInputDual.faces[0]); + + // Find triangle that contains nodeQ. + // This is the face whose local index is iFaceFinal + + int iFaceFinal = 0; + + GetFaceThatContainsPoint(meshInputDual,iNearestFace,iFaceFinal,nodeQ.x,nodeQ.y,nodeQ.z); + + //Pre-compute all triangle areas and subareas + + int iEdges = meshInputDual.faces[iFaceFinal].edges.size(); + + //Subtriangles with the sample point as a vertex (q,m,m+1) where q is the sample point + std::vector vecTriangleSubAreas(iEdges); + + //Subtriangles without sample point as a vertex (m-1,m,m+1) + std::vector vecTriangleAreas(iEdges); + + for (int m = 0; m < iEdges; m++){ + + Face faceTriangleM(3); //Triangle with vertices m-1,m,m+1 + + Face faceSubTriangleM(3); //Triangle with vertices q,m,m+1 + + faceTriangleM.SetNode(0,0); + faceTriangleM.SetNode(1,1); + faceTriangleM.SetNode(2,2); + + faceSubTriangleM.SetNode(0,0); + faceSubTriangleM.SetNode(1,1); + faceSubTriangleM.SetNode(2,2); + + NodeVector nodesTriangleM; + + NodeVector nodesSubTriangleM; + + Node nodeM(meshInputDual.nodes[meshInputDual.faces[iFaceFinal][m]].x, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][m]].y, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][m]].z); + + //c++ doesn't like the modulo operator when the argument is negative + + int iMMinusOne = m - 1; + + if (m == 0){ + + iMMinusOne = iEdges - 1; + + } + + + Node nodeMMinusOne(meshInputDual.nodes[meshInputDual.faces[iFaceFinal][iMMinusOne]].x, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][iMMinusOne]].y, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][iMMinusOne]].z); + + Node nodeMPlusOne(meshInputDual.nodes[meshInputDual.faces[iFaceFinal][(m+1)%iEdges]].x, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][(m+1)%iEdges]].y, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][(m+1)%iEdges]].z); + + nodesTriangleM.push_back(nodeMMinusOne); + nodesTriangleM.push_back(nodeM); + nodesTriangleM.push_back(nodeMPlusOne); + + nodesSubTriangleM.push_back(nodeQ); + nodesSubTriangleM.push_back(nodeM); + nodesSubTriangleM.push_back(nodeMPlusOne); + + double dSubAreaM = CalculateFaceArea(faceSubTriangleM, nodesSubTriangleM); + + vecTriangleSubAreas[m] = dSubAreaM; + + double dTriangleAreaM = CalculateFaceArea(faceTriangleM, nodesTriangleM); + + vecTriangleAreas[m] = dTriangleAreaM; + + + } + + std::vector vecWeights(iEdges,1); + + double dWeightTotal = 0; + + //Double loop over all nodes of the face because each weight depends on all + //triangle subareas + + for (int m = 0; m < iEdges; m++){ + + for (int l = 0; l < iEdges; l++){ + + int iLMinusOne = l - 1; + + if (l == 0){ + + iLMinusOne = iEdges - 1; + + } + + if (l != m && l != iLMinusOne){ + + + vecWeights[m] *= vecTriangleSubAreas[iLMinusOne]*vecTriangleSubAreas[l]; + + } + + } + + vecWeights[m] *= vecTriangleAreas[m]; + + dWeightTotal += vecWeights[m]; + + } + + for (int m = 0; m < iEdges; m++){ + + vecWeights[m] /= dWeightTotal; + + } + + double dFaceArea = CalculateFaceArea(meshInputDual.faces[iFaceFinal],meshInputDual.nodes); + + double dsum = 0; + + for (int s = 0; s < iEdges; s++){ + + dsum += vecTriangleSubAreas[s]; + + } + + // Contribution of this quadrature point to the map + for (int i = 0; i < iEdges; i++){ + + int iContributingFaceI = meshInputDual.faces[iFaceFinal][i]; + + smatMap(iTargetFace, iContributingFaceI) += + vecWeights[i] + * dQuadPtWeight(k,p) + * meshOverlap.vecFaceArea[ixOverlap + j] + / dQuadratureArea + / meshOutput.vecFaceArea[iTargetFace]; + + } + } + } + } + + // Increment the current overlap index + ixOverlap += nOverlapFaces; + } + +} + +/////////////////////////////////////////////////////////////////////////////// + +void LinearRemapGeneralizedBarycentric( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +) { + + + // Verify ReverseNodeArray has been calculated + if (meshInput.edgemap.size() == 0) { + _EXCEPTIONT("EdgeMap has not been calculated for meshInput"); + } + + // Order of triangular quadrature rule + const int TriQuadRuleOrder = 4; + + // Triangular quadrature rule + TriangularQuadratureRule triquadrule(TriQuadRuleOrder); + + const DataArray2D & dG = triquadrule.GetG(); + const DataArray1D & dW = triquadrule.GetW(); + + // Get SparseMatrix representation of the OfflineMap + SparseMatrix & smatMap = mapRemap.GetSparseMatrix(); + + //Dual mesh + Mesh meshInputDual = meshInput; + + //Construct dual mesh + Dual(meshInputDual); + + //Reverse node array + meshInputDual.ConstructReverseNodeArray(); + + //Construct edge map + meshInputDual.ConstructEdgeMap(); + + //kd-tree of dual mesh centers + kdtree * kdTarget = kd_create(3); + + // Vector of centers of the source mesh + for (int i = 0; i < meshInputDual.faces.size(); i++){ + + const Face & face = meshInputDual.faces[i]; + + Node nodeCentroid = GetFaceCentroid(face, meshInputDual.nodes); + + nodeCentroid = nodeCentroid.Normalized(); + + kd_insert3( + kdTarget, + nodeCentroid.x, + nodeCentroid.y, + nodeCentroid.z, + (void*)(&(meshInputDual.faces[i]))); + + } + + + // Loop through all target faces + for (int ixFirst = 0; ixFirst < meshOutput.faces.size(); ixFirst++) { + + // Output every 1000 overlap elements + if (ixFirst % 1000 == 0) { + Announce("Element %i/%i", ixFirst, meshOutput.faces.size()); + } + + // This Face + const Face & faceFirst = meshOutput.faces[ixFirst]; + + // Get node coordinates of each target face center + Node nodeQ = GetFaceCentroid(faceFirst,meshOutput.nodes); + + nodeQ = nodeQ.Normalized(); + + //Get the dual mesh face whose center is nearest to the target face + + kdres * kdresTarget = + kd_nearest3( + kdTarget, + nodeQ.x, + nodeQ.y, + nodeQ.z); + + Face * pFace = (Face *)(kd_res_item_data(kdresTarget)); + + int iNearestFace = pFace - &(meshInputDual.faces[0]); + + // Find triangle that contains nodeQ. + // This is the face whose local index is iFaceFinal + + int iFaceFinal = 0; + + GetFaceThatContainsPoint(meshInputDual,iNearestFace,iFaceFinal,nodeQ.x,nodeQ.y,nodeQ.z); + + //Pre-compute all triangle areas and subareas + + int iEdges = meshInputDual.faces[iFaceFinal].edges.size(); + + //Subtriangles with the sample point as a vertex (q,m,m+1) where q is the sample point + std::vector vecTriangleSubAreas(iEdges); + + //Subtriangles without sample point as a vertex (m-1,m,m+1) + std::vector vecTriangleAreas(iEdges); + + for (int m = 0; m < iEdges; m++){ + + Face faceTriangleM(3); //Triangle with vertices m-1,m,m+1 + + Face faceSubTriangleM(3); //Triangle with vertices q,m,m+1 + + faceTriangleM.SetNode(0,0); + faceTriangleM.SetNode(1,1); + faceTriangleM.SetNode(2,2); + + faceSubTriangleM.SetNode(0,0); + faceSubTriangleM.SetNode(1,1); + faceSubTriangleM.SetNode(2,2); + + NodeVector nodesTriangleM; + + NodeVector nodesSubTriangleM; + + Node nodeM(meshInputDual.nodes[meshInputDual.faces[iFaceFinal][m]].x, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][m]].y, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][m]].z); + + //c++ doesn't like the modulo operator when the argument is negative + + int iMMinusOne = m - 1; + + if (m == 0){ + + iMMinusOne = iEdges - 1; + + } + + + Node nodeMMinusOne(meshInputDual.nodes[meshInputDual.faces[iFaceFinal][iMMinusOne]].x, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][iMMinusOne]].y, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][iMMinusOne]].z); + + Node nodeMPlusOne(meshInputDual.nodes[meshInputDual.faces[iFaceFinal][(m+1)%iEdges]].x, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][(m+1)%iEdges]].y, + meshInputDual.nodes[meshInputDual.faces[iFaceFinal][(m+1)%iEdges]].z); + + nodesTriangleM.push_back(nodeMMinusOne); + nodesTriangleM.push_back(nodeM); + nodesTriangleM.push_back(nodeMPlusOne); + + nodesSubTriangleM.push_back(nodeQ); + nodesSubTriangleM.push_back(nodeM); + nodesSubTriangleM.push_back(nodeMPlusOne); + + double dSubAreaM = CalculateFaceArea(faceSubTriangleM, nodesSubTriangleM); + + vecTriangleSubAreas[m] = dSubAreaM; + + double dTriangleAreaM = CalculateFaceArea(faceTriangleM, nodesTriangleM); + + vecTriangleAreas[m] = dTriangleAreaM; + + + } + + double dfacearea = CalculateFaceArea(meshInputDual.faces[iFaceFinal],meshInputDual.nodes); + + std::vector vecWeights(iEdges,1); + + double dWeightTotal = 0; + + //Double loop over all nodes of the face because each weight depends on all + //triangle subareas + + for (int m = 0; m < iEdges; m++){ + + for (int l = 0; l < iEdges; l++){ + + int iLMinusOne = l - 1; + + if (l == 0){ + + iLMinusOne = iEdges - 1; + + } + + if (l != m && l != iLMinusOne){ + + + vecWeights[m] *= vecTriangleSubAreas[iLMinusOne]*vecTriangleSubAreas[l]; + + } + + } + + vecWeights[m] *= vecTriangleAreas[m]; + + dWeightTotal += vecWeights[m]; + + } + + for (int m = 0; m < iEdges; m++){ + + vecWeights[m] /= dWeightTotal; + + } + + // Contribution of this quadrature point to the map + for (int i = 0; i < iEdges; i++){ + + int iContributingFaceI = meshInputDual.faces[iFaceFinal][i]; + + smatMap(ixFirst, iContributingFaceI) = vecWeights[i]; + + } + } + +} + +/////////////////////////////////////////////////////////////////////////////// + +void LinearRemapTriangulation( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +) { + + // Verify ReverseNodeArray has been calculated + if (meshInput.edgemap.size() == 0) { + _EXCEPTIONT("EdgeMap has not been calculated for meshInput"); + } + + // Order of triangular quadrature rule + const int TriQuadRuleOrder = 4; + + // Triangular quadrature rule + TriangularQuadratureRule triquadrule(TriQuadRuleOrder); + + const DataArray2D & dG = triquadrule.GetG(); + const DataArray1D & dW = triquadrule.GetW(); + + // Get SparseMatrix represntation of the OfflineMap + SparseMatrix & smatMap = mapRemap.GetSparseMatrix(); + + // Array of global indices + std::vector > vecGlobalIndexI(6); + + // Vector storing meshes of each panel + std::vector vecMesh(6); + + // kd-tree for each panel + std::vector vecKDTreePanelI(6); + + // Arrays of panel boundaries in lat/lon coordinates (radians) + DataArray2D dPanelLat(6,2); + DataArray2D dPanelLon(6,2); + + // Width of buffer zone for each panel (degrees) + double dBuff = 30; + + // Define equatorial panel boundaries + for (int i = 0; i < 4; i++){ + + dPanelLon(i,0) = 90*M_PI*i/180; //Left boundary + dPanelLon(i,1) = 90*M_PI*(i+1)/180; //Right Boundary + dPanelLat(i,1) = 45*M_PI/180; //Upper Boundary + dPanelLat(i,0) = -45*M_PI/180; //Lower Boundary + + } + + // Define polar panels + + dPanelLon(4,0) = 0; dPanelLon(4,1) = 2*M_PI; + dPanelLat(4,0) = 45*M_PI/180; dPanelLat(4,1) = M_PI/2; + + dPanelLon(5,0) = 0; dPanelLon(5,1) = 2*M_PI; + dPanelLat(5,1) = -45*M_PI/180; dPanelLat(5,0) = -M_PI/2; + + + // Generate the Delaunay triangulation for each of the 6 panels + for (int i = 0; i < 6; i++){ + + // Array storing the coordinates of the face centroids of panel i + std::vector> dPanelCentroid(2); + + // Coordinates of tangent plane point of tangency for panel i + double dLatTan; + double dLonTan; + + if ((0 <= i) && (i <= 3)){ + + dLatTan = 0; + dLonTan = (i+1)*45*M_PI/180 + 45*i*M_PI/180; + + } + else if (i == 4){ + + dLatTan = M_PI/2; + dLonTan = 0; + + } + else{ + + dLonTan = 0; + dLatTan = -M_PI/2; + + } + + // Number of source faces centers for panel i; this is the size + // of in.pointlist + int nNodes = 0; + + for (int j = 0; j < meshInput.faces.size(); j++){ + + // Get coordinates of centroid of current face + + double dLonRad0; + double dLatRad0; + + const Face & face = meshInput.faces[j]; + + //const Node & node = GetFaceCentroid(face,meshInput.nodes); + Node node = GetFaceCentroid(face,meshInput.nodes); + + node = node.Normalized(); + + XYZtoRLL_Rad(node.x,node.y,node.z,dLonRad0,dLatRad0); + + // Determine if centroid is in panel i plus buffer zone + bool fPanelContainsCentroid = false; + + if ((0 <= i) && (i <= 3)){ + + if ((dPanelLat(i,0) - dBuff*M_PI/180 <= dLatRad0) && (dLatRad0 <= dPanelLat(i,1)+dBuff*M_PI/180)){ + + if ( i == 0 ){ + + if ( ((2*M_PI - dBuff*M_PI/180 <= dLonRad0) && (dLonRad0 <= 2*M_PI)) || + ((dPanelLon(i,0) <= dLonRad0) && (dLonRad0 <= dPanelLon(i,1) + dBuff*M_PI/180)) ) { + + fPanelContainsCentroid = true; + } + } + + else if ( i == 3 ){ + + if ( ((0 <= dLonRad0) && (dLonRad0 <= 0 + dBuff*M_PI/180)) || + ((dPanelLon(i,0) - dBuff*M_PI/180 <= dLonRad0) && (dLonRad0 <= dPanelLon(i,1) + dBuff*M_PI/180)) ){ + + fPanelContainsCentroid = true; + + } + + } + + else { + + if ( (dPanelLon(i,0) - dBuff*M_PI/180 <= dLonRad0) && (dLonRad0 <= dPanelLon(i,1) + dBuff*M_PI/180)){ + + fPanelContainsCentroid = true; + } + + + } + + } + } + else if (i == 4){ + + if ((dPanelLat(i,0) - dBuff*M_PI/180 <= dLatRad0) && (dLatRad0 <= dPanelLat(i,1))){ + + fPanelContainsCentroid = true; + + } + + } + + else { + + if ((dPanelLat(i,0) <= dLatRad0) && (dLatRad0 <= dPanelLat(i,1) + dBuff*M_PI/180)){ + fPanelContainsCentroid = true; + + } + + } + + if(fPanelContainsCentroid){ + + vecGlobalIndexI[i].push_back(j); + nNodes++; + + // Gnomonic projection of centroid coordinates + double dGX; + double dGY; + GnomonicProjection(dLonTan,dLatTan,dLonRad0,dLatRad0,dGX,dGY); + + // Add Gnomonic coordinates to vector + dPanelCentroid[0].push_back(dGX); + dPanelCentroid[1].push_back(dGY); + + } + } + + // Structures for Delaunay triangulation + + struct triangulateio in, out, vorout; + + in.numberofpoints = nNodes; + in.numberofpointattributes = 0; + in.numberofsegments = 0; + in.numberofholes = 0; + in.numberofregions = 0; + in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL)); + in.segmentlist = (int *) malloc(in.numberofsegments * 2 * sizeof(int));; + in.pointattributelist = (REAL *) NULL; + in.pointmarkerlist = (int *) NULL; + in.trianglelist = (int *) NULL; + in.triangleattributelist = (REAL *) NULL; + in.neighborlist = (int *) NULL; + in.segmentmarkerlist = (int *) NULL; + in.edgelist = (int *) NULL; + in.edgemarkerlist = (int *) NULL; + + // initialize data structure for output triangulation + out.pointlist = (REAL *) NULL; + out.pointattributelist = (REAL *) NULL; + out.pointmarkerlist = (int *) NULL; + out.trianglelist = (int *) NULL; + out.triangleattributelist = (REAL *) NULL; + out.neighborlist = (int *) NULL; + out.segmentlist = (int *) NULL; + out.segmentmarkerlist = (int *) NULL; + out.edgelist = (int *) NULL; + out.edgemarkerlist = (int *) NULL; + + // initialize data structure for output Voronoi diagram (unused) + vorout.pointlist = (REAL *) NULL; + vorout.pointattributelist = (REAL *) NULL; + vorout.edgelist = (int *) NULL; + vorout.normlist = (REAL *) NULL; + + // Add points to in.pointlist + + for (int j = 0; j < nNodes; j++){ + + in.pointlist[2*j+0] = dPanelCentroid[0][j]; + in.pointlist[2*j+1] = dPanelCentroid[1][j]; + + } + + // Compute Delaunay triangulation. Use the 'c' option so that + // the convex hull is included + + char options[256] ="cjzenYQ"; + triangulate(options, &in, &out, &vorout); + + // Convert to mesh object by building face vector + for (int j = 0; j < out.numberoftriangles; j++){ + + Face newFace(3); + + newFace.SetNode(0, out.trianglelist[3*j+0]); + newFace.SetNode(1, out.trianglelist[3*j+1]); + newFace.SetNode(2, out.trianglelist[3*j+2]); + vecMesh[i].faces.push_back(newFace); + + } + + // Build node vector. Note that the Delaunay triangulation preserves + // point indexing. + for (int j = 0; j < out.numberofpoints; j++){ + + Node node(out.pointlist[2*j+0],out.pointlist[2*j+1],0); + vecMesh[i].nodes.push_back(node); + + } + + vecMesh[i].ConstructReverseNodeArray(); + vecMesh[i].RemoveCoincidentNodes(); + vecMesh[i].RemoveZeroEdges(); + vecMesh[i].ConstructEdgeMap(); + + free(in.pointlist); + free(out.pointlist); + free(out.pointattributelist); + free(out.pointmarkerlist); + free(out.trianglelist); + free(out.triangleattributelist); + free(out.neighborlist); + free(out.segmentlist); + free(out.segmentmarkerlist); + free(out.edgelist); + free(out.edgemarkerlist); + + } + + // Construct kd-tree for each panel + + for (int i = 0; i < 6; i++){ + + vecKDTreePanelI[i] = kd_create(3); + + for (int j = 0; j < vecMesh[i].faces.size(); j++){ + + Face face = vecMesh[i].faces[j]; + + Node nodeCenter = GetFaceCentroid(face, vecMesh[i].nodes); + + kd_insert3( + vecKDTreePanelI[i], + nodeCenter.x, + nodeCenter.y, + nodeCenter.z, + (void*)(&(vecMesh[i].faces[j]))); + + } + + } + + // Loop through all target faces + for (int ixFirst = 0; ixFirst < meshOutput.faces.size(); ixFirst++) { + + // Output every 1000 overlap elements + if (ixFirst % 1000 == 0) { + Announce("Element %i/%i", ixFirst, meshOutput.faces.size()); + } + + // This Face + const Face & faceFirst = meshOutput.faces[ixFirst]; + + // Get node coordinates of each target face center + Node nodeQ = GetFaceCentroid(faceFirst,meshOutput.nodes); + + nodeQ = nodeQ.Normalized(); + + // Get lat/lon coordinates of nodeQ + double dLatRad; + double dLonRad; + + XYZtoRLL_Rad(nodeQ.x,nodeQ.y,nodeQ.z,dLonRad,dLatRad); + + // Determine the panel where nodeQ is located + int iPanelIndex; + + for (int i = 0; i < 6; i++){ + + if ((dPanelLat(i,0) <= dLatRad) && (dLatRad <= dPanelLat(i,1)) && + (dPanelLon(i,0) <= dLonRad) && (dLonRad <= dPanelLon(i,1))){ + + iPanelIndex = i; + break; + + } + + } + // Compute Gnomonic projection onto the corresponding plane + + double dGX,dGY; + double dLatTan; + double dLonTan; + + if ((0 <= iPanelIndex) && (iPanelIndex <= 3)){ + + dLatTan = 0; + dLonTan = (iPanelIndex+1)*45*M_PI/180 + 45*iPanelIndex*M_PI/180; + + } + else if (iPanelIndex == 4){ + + dLatTan = M_PI/2; + dLonTan = 0; + + } + else{ + + dLonTan = 0; + dLatTan = -M_PI/2; + + } + + GnomonicProjection(dLonTan,dLatTan,dLonRad,dLatRad,dGX,dGY); + + // Get point closest to dGX and dGY + kdres * kdresTarget = + kd_nearest3( + vecKDTreePanelI[iPanelIndex], + dGX, + dGY, + 0.0); + + Face * pFace = (Face *)(kd_res_item_data(kdresTarget)); + + int iNearestFace = pFace - &(vecMesh[iPanelIndex].faces[0]); + + // Find triangle that contains the Gnomonic projection nodeQ. + // This is the face whose local index is iFaceFinal + + int iFaceFinal; + + double dA,dB; + + BarycentricCoordinates(vecMesh[iPanelIndex],iNearestFace,dGX,dGY,dA,dB); + + GetTriangleThatContainsPoint(vecMesh[iPanelIndex],iNearestFace,iFaceFinal,dGX,dGY); + + // Get global indices of the vertices of this triangle and + // calculate corresponding centroids + + std::vector vecContributingFaceI(3); + DataArray2D dataContributingCentroids(3,3); + + // Indices on the local mesh of the containing triangle + + int iLocalVertex1 = vecMesh[iPanelIndex].faces[iFaceFinal][0]; + int iLocalVertex2 = vecMesh[iPanelIndex].faces[iFaceFinal][1]; + int iLocalVertex3 = vecMesh[iPanelIndex].faces[iFaceFinal][2]; + + // Indices on the source mesh of the containing triangle + + int iGlobalVertex1 = vecGlobalIndexI[iPanelIndex][iLocalVertex1]; + int iGlobalVertex2 = vecGlobalIndexI[iPanelIndex][iLocalVertex2]; + int iGlobalVertex3 = vecGlobalIndexI[iPanelIndex][iLocalVertex3]; + + vecContributingFaceI[0] = iGlobalVertex1; + vecContributingFaceI[1] = iGlobalVertex2; + vecContributingFaceI[2] = iGlobalVertex3; + + // The centroids of the source mesh are the vertices of the containing triangle + + for (int i = 0; i < 3; i++){ + + Face faceCurrent = meshInput.faces[vecContributingFaceI[i]]; + + Node nodeCenter = GetFaceCentroid(faceCurrent,meshInput.nodes); + + nodeCenter = nodeCenter.Normalized(); + + dataContributingCentroids(i,0) = nodeCenter.x; + dataContributingCentroids(i,1) = nodeCenter.y; + dataContributingCentroids(i,2) = nodeCenter.z; + + + } + + // Vector of areas of subtriangles + std::vector vecSubAreas(3); + + // Area of triangle is obtained by adding up areas of the three + // subtriangles that are formed by the sample point + + double dTriangleArea = 0; + + for (int i = 0; i < 3; i++){ + + Face faceCurrent(3); + + faceCurrent.SetNode(0,0); + faceCurrent.SetNode(1,1); + faceCurrent.SetNode(2,2); + + NodeVector nodesCurrent; + + nodesCurrent.push_back(nodeQ); + + Node node1(dataContributingCentroids((i+1)%3,0),dataContributingCentroids((i+1)%3,1), + dataContributingCentroids((i+1)%3,2)); + + Node node2(dataContributingCentroids((i+2)%3,0),dataContributingCentroids((i+2)%3,1), + dataContributingCentroids((i+2)%3,2)); + + nodesCurrent.push_back(node1); + nodesCurrent.push_back(node2); + + vecSubAreas[i] = CalculateFaceArea(faceCurrent,nodesCurrent); + + dTriangleArea += vecSubAreas[i]; + + } + + // Calculate vector of weights + std::vector vecContributingFaceWeights(3); + + for (int i = 0; i < 3; i++){ + + vecContributingFaceWeights[i] = vecSubAreas[i]/dTriangleArea; + + } + + // Contribution of this point to the map + for (int i = 0; i < vecContributingFaceI.size(); i++){ + + smatMap(ixFirst, vecContributingFaceI[i]) = vecContributingFaceWeights[i]; + + } + } + + +} + +/////////////////////////////////////////////////////////////////////////////// + +void LinearRemapIntegratedBilinear( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +) { + + // Verify ReverseNodeArray has been calculated + if (meshInput.edgemap.size() == 0) { + _EXCEPTIONT("EdgeMap has not been calculated for meshInput"); + } + + // Order of triangular quadrature rule + const int TriQuadRuleOrder = 4; + + // Triangular quadrature rule + TriangularQuadratureRule triquadrule(TriQuadRuleOrder); + + const DataArray2D & dG = triquadrule.GetG(); + const DataArray1D & dW = triquadrule.GetW(); + + // Get SparseMatrix representation of the OfflineMap + SparseMatrix & smatMap = mapRemap.GetSparseMatrix(); + + Mesh meshInputDual = meshInput; + + //Construct dual mesh + Dual(meshInputDual); + + //Construct edge map + + meshInputDual.ConstructEdgeMap(); + + //kd-tree of dual mesh centers + + kdtree * kdTarget = kd_create(3); + + // Vector of centers of the source mesh + for (int i = 0; i < meshInputDual.faces.size(); i++){ + + const Face & face = meshInputDual.faces[i]; + + Node nodeCentroid = GetFaceCentroid(face, meshInputDual.nodes); + + nodeCentroid = nodeCentroid.Normalized(); + + kd_insert3( + kdTarget, + nodeCentroid.x, + nodeCentroid.y, + nodeCentroid.z, + (void*)(&(meshInputDual.faces[i]))); + + } + + std::vector vecContributingFaceWeights; + + std::vector vecContributingFaceI; + + // Overlap face index + int ixOverlap = 0; + + // Loop through all source faces + for (int ixFirst = 0; ixFirst < meshInput.faces.size(); ixFirst++) { + + // Output every 1000 overlap elements + if (ixFirst % 1000 == 0) { + Announce("Element %i/%i", ixFirst, meshInput.faces.size()); + } + + // This Face + const Face & faceFirst = meshInput.faces[ixFirst]; + + // Find the set of Faces that overlap faceFirst + int ixOverlapBegin = ixOverlap; + int ixOverlapEnd = ixOverlapBegin; + + for (; ixOverlapEnd < meshOverlap.faces.size(); ixOverlapEnd++) { + if (meshOverlap.vecSourceFaceIx[ixOverlapEnd] != ixFirst) { + break; + } + } + + int nOverlapFaces = ixOverlapEnd - ixOverlapBegin; + + // Loop through all overlap faces associated with this source face + for (int j = 0; j < nOverlapFaces; j++) { + int iTargetFace = meshOverlap.vecTargetFaceIx[ixOverlap + j]; + + const Face & faceOverlap = meshOverlap.faces[ixOverlap + j]; + + int nSubTriangles = faceOverlap.edges.size() - 2; + + // Jacobian at each quadrature point + DataArray2D dQuadPtWeight(nSubTriangles, dW.GetRows()); + DataArray2D dQuadPtNodes(nSubTriangles, dW.GetRows()); + + // Compute quadrature area of each overlap face + double dQuadratureArea = 0.0; + + for (int k = 0; k < nSubTriangles; k++) { + for (int p = 0; p < dW.GetRows(); p++) { + + dQuadPtWeight(k,p) = + CalculateSphericalTriangleJacobianBarycentric( + meshOverlap.nodes[faceOverlap[0]], + meshOverlap.nodes[faceOverlap[k+1]], + meshOverlap.nodes[faceOverlap[k+2]], + dG(p,0), dG(p,1), + &(dQuadPtNodes(k,p))); + + dQuadPtWeight(k,p) *= dW[p]; + + dQuadratureArea += dQuadPtWeight(k,p); + } + } + + // Loop through all sub-triangles of this overlap Face + for (int k = 0; k < nSubTriangles; k++) { + + // Loop through all quadrature nodes on this overlap Face + for (int p = 0; p < dW.GetRows(); p++) { + + // Get quadrature node and pointwise Jacobian + Node & nodeQ = dQuadPtNodes(k,p); + + //Get the dual mesh face whose center is nearest to the sample point + + kdres * kdresTarget = + kd_nearest3( + kdTarget, + nodeQ.x, + nodeQ.y, + nodeQ.z); + + Face * pFace = (Face *)(kd_res_item_data(kdresTarget)); + + int iNearestFace = pFace - &(meshInputDual.faces[0]); + + int iFaceFinal = 0; + + //Get the face that contains nodeQ + + GetFaceThatContainsPoint(meshInputDual,iNearestFace,iFaceFinal,nodeQ.x,nodeQ.y,nodeQ.z); + + int iEdges = meshInputDual.faces[iFaceFinal].edges.size(); + + vecContributingFaceWeights.clear(); + + vecContributingFaceI.clear(); + + DataArray1D dCoeffs(3); + + double dCond = 0; + + bool fConverged = false; + + if( iEdges == 3 ){ + + NodeVector nodesP; + + for (int i = 0; i < 3; i++){ + + Node nodeI = meshInputDual.nodes[meshInputDual.faces[iFaceFinal][i]]; + + nodesP.push_back(nodeI); + + } + + TriangleLineIntersection(nodeQ, nodesP, dCoeffs, dCond); + + vecContributingFaceWeights.push_back(1 - dCoeffs[1] - dCoeffs[2]); + vecContributingFaceWeights.push_back(dCoeffs[1]); + vecContributingFaceWeights.push_back(dCoeffs[2]); + + for (int i = 0; i < iEdges; i++){ + + vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][i]); + + } + + } + + else if ( iEdges == 4 ){ + + NodeVector nodesP; + + for (int i = 0; i < 4; i++){ + + Node nodeI = meshInputDual.nodes[meshInputDual.faces[iFaceFinal][i]]; + + nodesP.push_back(nodeI); + + } + + NewtonQuadrilateral(nodeQ, nodesP, dCoeffs, fConverged); + + for (int i = 0; i < iEdges; i++){ + + vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][i]); + + } + + vecContributingFaceWeights.push_back(1.0 - dCoeffs[0] - dCoeffs[1] + dCoeffs[0]*dCoeffs[1]); + vecContributingFaceWeights.push_back(dCoeffs[0]*(1.0 - dCoeffs[1])); + vecContributingFaceWeights.push_back(dCoeffs[0]*dCoeffs[1]); + vecContributingFaceWeights.push_back(dCoeffs[1]*(1.0 - dCoeffs[0])); + + } + + else { + //Loop over the subtrianges until we find one that contains the sample point + + int nSubTriangles = iEdges - 2; + + double dSum; + + for (int k = 0; k < nSubTriangles; k++) { + + Node & node0 = meshInputDual.nodes[meshInputDual.faces[iFaceFinal][0]]; + Node & nodeKPlusOne = meshInputDual.nodes[meshInputDual.faces[iFaceFinal][k+1]]; + Node & nodeKPlusTwo = meshInputDual.nodes[meshInputDual.faces[iFaceFinal][k+2]]; + + NodeVector nodesP; + + nodesP.push_back(node0); + nodesP.push_back(nodeKPlusOne); + nodesP.push_back(nodeKPlusTwo); + + if( DoesFaceContainPoint(nodesP, nodeQ.x, nodeQ.y, nodeQ.z) ){ + + TriangleLineIntersection(nodeQ, nodesP, dCoeffs, dCond); + + vecContributingFaceWeights.push_back(1.0 - dCoeffs[1] - dCoeffs[2]); + vecContributingFaceWeights.push_back(dCoeffs[1]); + vecContributingFaceWeights.push_back(dCoeffs[2]); + + vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][0]); + vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][k+1]); + vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][k+2]); + + break; + + } + + } + + } + + // Contribution of this quadrature point to the map + for (int i = 0; i < vecContributingFaceI.size(); i++){ + + if( vecContributingFaceWeights[i] < 0 || vecContributingFaceWeights[i] > 1 ){ + + _EXCEPTIONT("Non-monotone weight"); + + } + + smatMap(iTargetFace, vecContributingFaceI[i]) += + vecContributingFaceWeights[i] + * dQuadPtWeight(k,p) + * meshOverlap.vecFaceArea[ixOverlap + j] + / dQuadratureArea + / meshOutput.vecFaceArea[iTargetFace]; + } + + } + + } + + } + + // Increment the current overlap index + ixOverlap += nOverlapFaces; + } + +} + +/////////////////////////////////////////////////////////////////////////////// + +void LinearRemapBilinear( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +) { + + // Verify ReverseNodeArray has been calculated + if (meshInput.edgemap.size() == 0) { + _EXCEPTIONT("EdgeMap has not been calculated for meshInput"); + } + + // Order of triangular quadrature rule + const int TriQuadRuleOrder = 4; + + // Triangular quadrature rule + TriangularQuadratureRule triquadrule(TriQuadRuleOrder); + + const DataArray2D & dG = triquadrule.GetG(); + const DataArray1D & dW = triquadrule.GetW(); + + // Get SparseMatrix representation of the OfflineMap + SparseMatrix & smatMap = mapRemap.GetSparseMatrix(); + + Mesh meshInputDual = meshInput; + + //Construct dual mesh + Dual(meshInputDual); + + //Construct edge map + + meshInputDual.ConstructEdgeMap(); + + //kd-tree of dual mesh centers + + kdtree * kdTarget = kd_create(3); + + // Vector of centers of the source mesh + for (int i = 0; i < meshInputDual.faces.size(); i++){ + + const Face & face = meshInputDual.faces[i]; + + Node nodeCentroid = GetFaceCentroid(face, meshInputDual.nodes); + + nodeCentroid = nodeCentroid.Normalized(); + + kd_insert3( + kdTarget, + nodeCentroid.x, + nodeCentroid.y, + nodeCentroid.z, + (void*)(&(meshInputDual.faces[i]))); + + } + + std::vector vecContributingFaceWeights; + + std::vector vecContributingFaceI; + + // Loop through all target faces + for (int ixFirst = 0; ixFirst < meshOutput.faces.size(); ixFirst++) { + + // Output every 1000 overlap elements + if (ixFirst % 1000 == 0) { + Announce("Element %i/%i", ixFirst, meshOutput.faces.size()); + } + + // This Face + const Face & faceFirst = meshOutput.faces[ixFirst]; + + // Get node coordinates of each target face center + Node nodeQ = GetFaceCentroid(faceFirst,meshOutput.nodes); + + nodeQ = nodeQ.Normalized(); + + //Get the dual mesh face whose center is nearest to the target face + + kdres * kdresTarget = + kd_nearest3( + kdTarget, + nodeQ.x, + nodeQ.y, + nodeQ.z); + + Face * pFace = (Face *)(kd_res_item_data(kdresTarget)); + + int iNearestFace = pFace - &(meshInputDual.faces[0]); + + // Find face that contains nodeQ. + // This is the face whose local index is iFaceFinal + + int iFaceFinal = 0; + + GetFaceThatContainsPoint(meshInputDual,iNearestFace,iFaceFinal,nodeQ.x,nodeQ.y,nodeQ.z); + + int iEdges = meshInputDual.faces[iFaceFinal].edges.size(); + + vecContributingFaceWeights.clear(); + + vecContributingFaceI.clear(); + + DataArray1D dCoeffs(3); + + double dCond = 0; + + bool fConverged = false; + + if( iEdges == 3 ){ + + NodeVector nodesP; + + for (int i = 0; i < 3; i++){ + + Node nodeI = meshInputDual.nodes[meshInputDual.faces[iFaceFinal][i]]; + + nodesP.push_back(nodeI); + + } + + TriangleLineIntersection(nodeQ, nodesP, dCoeffs, dCond); + + vecContributingFaceWeights.push_back(1 - dCoeffs[1] - dCoeffs[2]); + vecContributingFaceWeights.push_back(dCoeffs[1]); + vecContributingFaceWeights.push_back(dCoeffs[2]); + + for (int i = 0; i < iEdges; i++){ + + vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][i]); + + } + + } + + else if ( iEdges == 4 ){ + + NodeVector nodesP; + + for (int i = 0; i < 4; i++){ + + Node nodeI = meshInputDual.nodes[meshInputDual.faces[iFaceFinal][i]]; + + nodesP.push_back(nodeI); + + } + + NewtonQuadrilateral(nodeQ, nodesP, dCoeffs, fConverged); + + for (int i = 0; i < iEdges; i++){ + + vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][i]); + + } + + vecContributingFaceWeights.push_back(1.0 - dCoeffs[0] - dCoeffs[1] + dCoeffs[0]*dCoeffs[1]); + vecContributingFaceWeights.push_back(dCoeffs[0]*(1.0 - dCoeffs[1])); + vecContributingFaceWeights.push_back(dCoeffs[0]*dCoeffs[1]); + vecContributingFaceWeights.push_back(dCoeffs[1]*(1.0 - dCoeffs[0])); + + } + + else { + //Loop over the subtrianges until we find one that contains the sample point + + int nSubTriangles = iEdges - 2; + + for (int k = 0; k < nSubTriangles; k++) { + + Node & node0 = meshInputDual.nodes[meshInputDual.faces[iFaceFinal][0]]; + Node & nodeKPlusOne = meshInputDual.nodes[meshInputDual.faces[iFaceFinal][k+1]]; + Node & nodeKPlusTwo = meshInputDual.nodes[meshInputDual.faces[iFaceFinal][k+2]]; + + NodeVector nodesP; + + nodesP.push_back(node0); + nodesP.push_back(nodeKPlusOne); + nodesP.push_back(nodeKPlusTwo); + + if( DoesFaceContainPoint(nodesP, nodeQ.x, nodeQ.y, nodeQ.z) ){ + + TriangleLineIntersection(nodeQ, nodesP, dCoeffs, dCond); + + vecContributingFaceWeights.push_back(1.0 - dCoeffs[1] - dCoeffs[2]); + vecContributingFaceWeights.push_back(dCoeffs[1]); + vecContributingFaceWeights.push_back(dCoeffs[2]); + + vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][0]); + vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][k+1]); + vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][k+2]); + + break; + + } + + } + + } + + // Contribution of each point to the map + for (int i = 0; i < vecContributingFaceI.size(); i++){ + + if( vecContributingFaceWeights[i] < 0 || vecContributingFaceWeights[i] > 1 ){ + + _EXCEPTIONT("Non-monotone weight"); + + } + + int iContributingFaceI = vecContributingFaceI[i]; + + smatMap(ixFirst, iContributingFaceI) = vecContributingFaceWeights[i]; + + } + + } + +} + +/////////////////////////////////////////////////////////////////////////////// + void ForceIntArrayConsistencyConservation( const DataArray1D & vecSourceArea, const DataArray1D & vecTargetArea, diff --git a/src/LinearRemapFV.h b/src/LinearRemapFV.h index 161e627..196a99b 100644 --- a/src/LinearRemapFV.h +++ b/src/LinearRemapFV.h @@ -67,6 +67,84 @@ void LinearRemapFVtoFV( /////////////////////////////////////////////////////////////////////////////// +/// +/// Generate the OfflineMap for integrated remapping from finite volumes to finite +/// volumes using a triangulation of the source mesh. +/// +void LinearRemapIntegratedTriangulation( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Generate the OfflineMap for non-integrated remapping from finite volumes to finite +/// volumes using generalized barycentric coordinates using the dual mesh +/// +void LinearRemapIntegratedGeneralizedBarycentric( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Generate the OfflineMap for integrated remapping from finite volumes to finite +/// volumes using generalized barycentric coordinates using the dual mesh +/// +void LinearRemapGeneralizedBarycentric( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Generate the OfflineMap for non-integrated remapping from finite volumes to finite +/// volumes using a triangulation of the source mesh. +/// +void LinearRemapTriangulation( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Generate the OfflineMap for integrated remapping from finite volumes to finite +/// volumes using the ESMF based bilinear interpolation +/// +void LinearRemapIntegratedBilinear( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Generate the OfflineMap for non-integrated remapping from finite volumes to finite +/// volumes using the ESMF based bilinear interpolation +/// +void LinearRemapBilinear( + const Mesh & meshInput, + const Mesh & meshOutput, + const Mesh & meshOverlap, + OfflineMap & mapRemap +); + +/////////////////////////////////////////////////////////////////////////////// + /// /// Generate the OfflineMap for remapping from finite volumes to finite /// elements using simple sampling of the FV reconstruction. diff --git a/src/OfflineMap.cpp b/src/OfflineMap.cpp index 938f320..d09e379 100644 --- a/src/OfflineMap.cpp +++ b/src/OfflineMap.cpp @@ -21,6 +21,7 @@ #include "GridElements.h" #include "FiniteElementTools.h" #include "STLStringHelper.h" +#include "FiniteVolumeTools.h" #include "Announce.h" #include "Exception.h" @@ -114,12 +115,14 @@ void ParseEnforceBounds( bool fLowerBoundIsFloat = STLStringHelper::IsFloat(it.strLowerBound); bool fUpperBoundIsFloat = STLStringHelper::IsFloat(it.strUpperBound); - if ((it.strLowerBound != "n") && (it.strLowerBound != "l") && (it.strLowerBound != "g") && (!fLowerBoundIsFloat)) { - _EXCEPTION1("Invalid lower bound in bounds string \"%s\", expected \"n\", \"l\", \"g\", or floating point value", + if ((it.strLowerBound != "n") && (it.strLowerBound != "l") && (it.strLowerBound != "g") && (!fLowerBoundIsFloat) && + (it.strLowerBound != "lp")) { + _EXCEPTION1("Invalid lower bound in bounds string \"%s\", expected \"n\", \"l\", \"g\", \"lp\", or floating point value", it.strLowerBound.c_str()); } - if ((it.strUpperBound != "n") && (it.strUpperBound != "l") && (it.strUpperBound != "g") && (!fUpperBoundIsFloat)) { - _EXCEPTION1("Invalid upper bound in bounds string \"%s\", expected \"n\", \"l\", \"g\", or floating point value", + if ((it.strUpperBound != "n") && (it.strUpperBound != "l") && (it.strUpperBound != "g") && (!fUpperBoundIsFloat) && + (it.strUpperBound != "lp")) { + _EXCEPTION1("Invalid upper bound in bounds string \"%s\", expected \"n\", \"l\", \"g\", \"lp\", or floating point value", it.strUpperBound.c_str()); } if (fLowerBoundIsFloat && fUpperBoundIsFloat) { @@ -140,6 +143,65 @@ void ParseEnforceBounds( /////////////////////////////////////////////////////////////////////////////// +void OfflineMap::SetEnforcementBounds( + const std::string & strEnforcementBounds, + Mesh * pmeshSource, + Mesh * pmeshOverlap, + DataArray3D * pdataGLLNodesIn, + DataArray3D * pdataGLLNodesOut, + int iPin +){ + + ParseEnforceBounds(strEnforcementBounds,m_vecEnforcementBounds); + + bool fHasLocalBounds = false; + bool fHasLocalPBounds = false; + + for (int i = 0; i < m_vecEnforcementBounds.size(); i++) { + if ((m_vecEnforcementBounds[i].strLowerBound == "l") || + (m_vecEnforcementBounds[i].strUpperBound == "l") + ) { + fHasLocalBounds = true; + } + } + + for (int i = 0; i < m_vecEnforcementBounds.size(); i++) { + if ((m_vecEnforcementBounds[i].strLowerBound == "lp") || + (m_vecEnforcementBounds[i].strUpperBound == "lp") + ) { + fHasLocalPBounds = true; + } + } + + if (fHasLocalBounds || fHasLocalPBounds) { + if (pmeshOverlap == NULL) { + _EXCEPTIONT("Local bounds enforcement requires the overlap mesh"); + } + if(fHasLocalPBounds){ + if(pmeshSource == NULL){ + _EXCEPTIONT("Local p bounds enforcement requires the source mesh"); + } + } + //if((!fgll) & ((pmeshSource == NULL))){ + // _EXCEPTIONT("Local bounds enforcement with a finite volume source mesh requires source and overlap meshes"); + //} + } + + m_pmeshSource = pmeshSource; + + m_pmeshOverlap = pmeshOverlap; + + m_pdataGLLNodesIn = pdataGLLNodesIn; + + m_pdataGLLNodesOut = pdataGLLNodesOut; + + m_iPin = iPin; + + +} + +/////////////////////////////////////////////////////////////////////////////// + void OfflineMap::InitializeDimensionsFromMeshFile( const std::string & strMeshFile, std::vector & vecDimNames, @@ -441,6 +503,115 @@ void OfflineMap::InitializeTargetDimensions( /////////////////////////////////////////////////////////////////////////////// +void OfflineMap::CAAS( + DataArray1D & dataCorrectedField, + DataArray1D & dataLowerBound, + DataArray1D & dataUpperBound, + double & dMass){ + + DataArray1D dataCorrection(dataLowerBound.GetRows()); + + + for (int i = 0; i < dataLowerBound.GetRows(); i++){ + dataCorrection[i] = fmax(dataLowerBound[i],fmin(dataUpperBound[i],0.0)); + } + + + double dMassL=0.0; + double dMassU=0.0; + + for (int i = 0; i < dataLowerBound.GetRows(); i++){ + dMassL += m_dTargetAreas[i]*dataLowerBound[i]; + dMassU += m_dTargetAreas[i]*dataUpperBound[i]; + } + + double dMassDiff = dMass; + + for (int i = 0; i < dataCorrectedField.GetRows(); i++){ + dMassDiff -= m_dTargetAreas[i]*dataCorrection[i]; + } + + double dLMinusU = abs(dataLowerBound[0]-dataUpperBound[0]); + + for (int i = 0; i < dataLowerBound.GetRows(); i++){ + if(abs(dataLowerBound[i] - dataUpperBound[i]) > dLMinusU){ + + dLMinusU = abs(dataLowerBound[i] - dataUpperBound[i]); + + } + } + + //If the upper and lower bounds are too close together, just clip + if(dMassDiff == 0 || (dLMinusU < 1e-13)){ + + for(int i = 0; i < dataUpperBound.GetRows(); i++){ + + dataCorrectedField[i] += dataCorrection[i]; + + } + + return; + } + else if(dMassL > dMassDiff){ + + Announce("Lower bound mass exceeds target mass by %1.15e: CAAS not applied",dMassL-dMassDiff); + return; + + } + else if(dMassU < dMassDiff){ + + Announce("Target mass exceeds upper bound mass by %1.15e: CAAS not applied",dMassDiff-dMassU); + return; + + } + else{ + + DataArray1D dataMassVec(dataUpperBound.GetRows()); //vector of mass redistribution + + if(dMassDiff > 0.0){ + double dMassCorrectU = 0.0; + for (int i = 0; i < dataUpperBound.GetRows(); i++){ + + dMassCorrectU += m_dTargetAreas[i]*(dataUpperBound[i]-dataCorrection[i]); + + } + + for (int i = 0; i < dataUpperBound.GetRows(); i++){ + + dataMassVec[i] = (dataUpperBound[i]-dataCorrection[i])/dMassCorrectU; + dataCorrection[i]=dataCorrection[i]+dMassDiff*dataMassVec[i]; + + } + } + else{ + double dMassCorrectL = 0.0; + for (int i = 0; i < dataUpperBound.GetRows(); i++){ + + dMassCorrectL += m_dTargetAreas[i]*(dataCorrection[i]-dataLowerBound[i]); + + } + + for (int i = 0; i < dataLowerBound.GetRows(); i++){ + + dataMassVec[i] = (dataCorrection[i]-dataLowerBound[i])/dMassCorrectL; + dataCorrection[i] = dataCorrection[i]+dMassDiff*dataMassVec[i]; + + } + } + + for (int i = 0; i < dataUpperBound.GetRows(); i++){ + + dataCorrectedField[i] += dataCorrection[i]; + + } + + } + + return; +} + +/////////////////////////////////////////////////////////////////////////////// + void OfflineMap::InitializeCoordinatesFromMeshFV( const Mesh & mesh, DataArray1D & dCenterLon, @@ -1864,6 +2035,371 @@ void OfflineMap::Apply( // Apply the offline map to the data m_mapRemap.Apply(dataInDouble, dataOutDouble); + + int iVarIndex = 0; + + bool fCAAS = false; + + for (; iVarIndex < m_vecEnforcementBounds.size(); iVarIndex++){ + if((m_vecEnforcementBounds[iVarIndex].strVariable == vecVariableList[v].c_str()) || + ((m_vecEnforcementBounds.size() == 1) & (m_vecEnforcementBounds[0].strVariable==""))){ + fCAAS = true; + break; + } + } + + //Upper and lower bounds for CAAS + + if(fCAAS){ + + DataArray1D dataLowerBound(nTargetCount); + DataArray1D dataUpperBound(nTargetCount); + DataArray1D x(nTargetCount); + double dMassDiff = dSourceMass; + + for (int i = 0; i < dataLowerBound.GetRows(); i++){ + dMassDiff -= dataOutDouble[i]*m_dTargetAreas[i]; + } + + + if(m_vecEnforcementBounds[iVarIndex].strLowerBound == "g"){ + for(int i = 0; i < nTargetCount; i++){ + dataLowerBound[i] = dSourceMin - dataOutDouble[i]; + } + } + + if(STLStringHelper::IsFloat(m_vecEnforcementBounds[iVarIndex].strLowerBound)){ + for (int i = 0; i < nTargetCount; i++){ + dataLowerBound[i] = std::stod(m_vecEnforcementBounds[iVarIndex].strLowerBound) - dataOutDouble[i]; + } + } + + if(m_vecEnforcementBounds[iVarIndex].strLowerBound == "n"){ + for(int i = 0; i < nTargetCount; i++){ + dataLowerBound[i] = -DBL_MAX; + } + } + + if(m_vecEnforcementBounds[iVarIndex].strUpperBound == "g"){ + for(int i = 0; i < nTargetCount; i++){ + dataUpperBound[i] = dSourceMax - dataOutDouble[i]; + } + } + + if(STLStringHelper::IsFloat(m_vecEnforcementBounds[iVarIndex].strUpperBound)){ + for (int i = 0; i < nTargetCount; i++){ + dataUpperBound[i] = std::stod(m_vecEnforcementBounds[iVarIndex].strUpperBound) - dataOutDouble[i]; + } + } + + if(m_vecEnforcementBounds[iVarIndex].strUpperBound == "n"){ + for (int i = 0; i < nTargetCount; i++){ + dataUpperBound[i] = DBL_MAX; + } + } + + if((m_vecEnforcementBounds[iVarIndex].strLowerBound == "l") || + (m_vecEnforcementBounds[iVarIndex].strLowerBound == "lp") || + (m_vecEnforcementBounds[iVarIndex].strUpperBound == "l") || + (m_vecEnforcementBounds[iVarIndex].strUpperBound == "lp")){ + + double dMinI; + double dMaxI; + + int nTargetFaces; + + if(m_pdataGLLNodesOut != NULL){ + nTargetFaces = m_pdataGLLNodesOut->GetSize(2); + } + else{ + nTargetFaces = nTargetCount; + } + + std::vector vecLocalUpperBound(nTargetCount); + std::vector vecLocalLowerBound(nTargetCount); + + std::vector< std::vector > vecSourceOvTarget(nTargetFaces); + + for (int i=0; ifaces.size(); i++){ + + int ixT = m_pmeshOverlap->vecTargetFaceIx[i]; + int ixS = m_pmeshOverlap->vecSourceFaceIx[i]; + vecSourceOvTarget[ixT].push_back(ixS); + + } + + //FV to FV + if(m_pdataGLLNodesIn == NULL && m_pdataGLLNodesOut == NULL){ + + for (int i = 0; i < nTargetCount; i++){ + + dMaxI=dataInDouble[vecSourceOvTarget[i][0]]; + dMinI=dataInDouble[vecSourceOvTarget[i][0]]; + + //Compute max over interstecting source faces + + for (int j = 0; j < vecSourceOvTarget[i].size(); j++){ + + int k = vecSourceOvTarget[i][j]; + dMaxI = fmax(dMaxI,dataInDouble[k]); + dMinI = fmin(dMinI,dataInDouble[k]); + + } + + if(m_vecEnforcementBounds[iVarIndex].strLowerBound == "lp" || + m_vecEnforcementBounds[iVarIndex].strUpperBound == "lp"){ + + double dMaxIAdj = dMaxI; + double dMinIAdj = dMinI; + + AdjacentFaceVector vecAdjFaces; + + GetAdjacentFaceVectorByEdge(*m_pmeshSource,vecSourceOvTarget[i][0],(m_iPin+1)*(m_iPin+1),vecAdjFaces); + + //Compute max over neighboring faces + + for (int j = 0; j < vecAdjFaces.size(); j++) { + + int k = vecAdjFaces[j].first; + + dMaxIAdj=fmax(dMaxIAdj,dataInDouble[k]); + dMinIAdj=fmin(dMinIAdj,dataInDouble[k]); + + } + + if(m_vecEnforcementBounds[iVarIndex].strLowerBound == "lp"){ + vecLocalLowerBound[i] = dMinIAdj; + } + if(m_vecEnforcementBounds[iVarIndex].strUpperBound == "lp"){ + vecLocalUpperBound[i] = dMaxIAdj; + } + + } + + if(m_vecEnforcementBounds[iVarIndex].strLowerBound == "l"){ + vecLocalLowerBound[i] = dMinI; + } + if(m_vecEnforcementBounds[iVarIndex].strUpperBound == "l"){ + vecLocalUpperBound[i] = dMaxI; + } + + } + + } + + //FE to FE + + else if(m_pdataGLLNodesIn != NULL && m_pdataGLLNodesOut != NULL){ + + std::vector< std::vector > vecMaxI(nTargetCount); + std::vector< std::vector > vecMinI(nTargetCount); + + + for (int i = 0; i < m_pdataGLLNodesOut->GetSize(2); i++){ + + int q = vecSourceOvTarget[i][0]; + + dMaxI = dataInDouble[(*m_pdataGLLNodesIn)[0][0][q]-1]; + dMinI = dataInDouble[(*m_pdataGLLNodesIn)[0][0][q]-1]; + + //Compute maximum over all source elements that intersect target element i + + for (int j = 0; j < vecSourceOvTarget[i].size(); j++){ + + for (int k = 0; k < m_pdataGLLNodesIn->GetSize(0); k++){ + for (int r = 0; r < m_pdataGLLNodesIn->GetSize(1); r++){ + + int m = vecSourceOvTarget[i][j]; + int n = (*m_pdataGLLNodesIn)[k][r][m]-1; + + dMaxI = fmax(dMaxI,dataInDouble[n]); + dMinI = fmin(dMinI,dataInDouble[n]); + + } + } + } + + + for (int l = 0; l < m_pdataGLLNodesOut->GetSize(0); l++){ + for (int s = 0; s < m_pdataGLLNodesOut->GetSize(1); s++){ + + int t = (*m_pdataGLLNodesOut)[l][s][i]-1; + + vecMaxI[t].push_back(dMaxI); + vecMinI[t].push_back(dMinI); + + } + } + + } + + for (int i = 0; i < nTargetCount; i++){ + + double dMaxNodeI = vecMaxI[i][0]; + double dMinNodeI = vecMinI[i][0]; + + + //Determine maximum and minimum for corner and edge nodes + for (int j = 0; j < vecMaxI[i].size(); j++){ + dMaxNodeI = fmax(dMaxNodeI,vecMaxI[i][j]); + dMinNodeI = fmin(dMinNodeI,vecMinI[i][j]); + } + + + if(m_vecEnforcementBounds[iVarIndex].strUpperBound == "l"){ + vecLocalUpperBound[i] = dMaxNodeI; + } + if(m_vecEnforcementBounds[iVarIndex].strLowerBound == "l"){ + vecLocalLowerBound[i] = dMinNodeI; + } + + } + + } + + + //FV to GLL + else if(m_pdataGLLNodesOut != NULL){ + + std::vector< std::vector > vecMaxI(nTargetCount); + std::vector< std::vector > vecMinI(nTargetCount); + + for (int i = 0; i < m_pdataGLLNodesOut->GetSize(2); i++){ + + dMaxI=dataInDouble[vecSourceOvTarget[i][0]]; + dMinI=dataInDouble[vecSourceOvTarget[i][0]]; + + for (int j = 0; j < vecSourceOvTarget[i].size(); j++) { + + int k = vecSourceOvTarget[i][j]; + dMaxI=fmax(dMaxI,dataInDouble[k]); + dMinI=fmin(dMinI,dataInDouble[k]); + + } + + double dMaxIAdj; + double dMinIAdj; + + if(m_vecEnforcementBounds[iVarIndex].strLowerBound == "lp" || + m_vecEnforcementBounds[iVarIndex].strUpperBound == "lp"){ + + dMaxIAdj = dMaxI; + dMinIAdj = dMinI; + + AdjacentFaceVector vecAdjFaces; + + GetAdjacentFaceVectorByEdge(*m_pmeshSource,vecSourceOvTarget[i][0],(m_iPin+1)*(m_iPin+1),vecAdjFaces); + + for (int j = 0; j < vecAdjFaces.size(); j++) { + + int k = vecAdjFaces[j].first; + + dMaxIAdj=fmax(dMaxIAdj,dataInDouble[k]); + dMinIAdj=fmin(dMinIAdj,dataInDouble[k]); + + } + + } + + //For lp bounds, push back the max/min over the intersecting set, and the adjacency set + //For local bounds, push back the max/min over the intersecting set + + for (int l = 0; l < m_pdataGLLNodesOut->GetSize(0); l++){ + for (int s = 0; s < m_pdataGLLNodesOut->GetSize(1); s++){ + + int t = (*m_pdataGLLNodesOut)[l][s][i]-1; + + if(m_vecEnforcementBounds[iVarIndex].strLowerBound == "lp"){ + vecMaxI[t].push_back(dMaxIAdj); + } + else{ + vecMaxI[t].push_back(dMaxI); + } + if(m_vecEnforcementBounds[iVarIndex].strUpperBound == "lp"){ + vecMinI[t].push_back(dMinIAdj); + } + else{ + vecMinI[t].push_back(dMinI); + } + + } + } + } + + for (int i = 0; i < nTargetCount; i++){ + + double dMaxNodeI = vecMaxI[i][0]; + double dMinNodeI = vecMinI[i][0]; + + + //Determine maximum and minimum for corner and edge nodes + for (int j = 0; j < vecMaxI[i].size(); j++){ + dMaxNodeI = fmax(dMaxNodeI,vecMaxI[i][j]); + dMinNodeI = fmin(dMinNodeI,vecMinI[i][j]); + } + + vecLocalUpperBound[i] = dMaxNodeI; + + vecLocalLowerBound[i] = dMinNodeI; + + } + + } + + else{ + + for (int i = 0; i < nTargetCount; i++){ + + int q = vecSourceOvTarget[i][0]; + + dMaxI = dataInDouble[(*m_pdataGLLNodesIn)[0][0][q]-1]; + dMinI = dataInDouble[(*m_pdataGLLNodesIn)[0][0][q]-1]; + + //Compute maximum over all source elements that intersect target element i + + for (int j = 0; j < vecSourceOvTarget[i].size(); j++){ + + for (int k = 0; k < m_pdataGLLNodesIn->GetSize(0); k++){ + for (int r = 0; r < m_pdataGLLNodesIn->GetSize(1); r++){ + + int m = vecSourceOvTarget[i][j]; + int n = (*m_pdataGLLNodesIn)[k][r][m]-1; + + dMaxI = fmax(dMaxI,dataInDouble[n]); + dMinI = fmin(dMinI,dataInDouble[n]); + + } + } + } + + if(m_vecEnforcementBounds[iVarIndex].strLowerBound == "l"){ + vecLocalLowerBound[i] = dMinI; + } + if(m_vecEnforcementBounds[iVarIndex].strUpperBound == "l"){ + vecLocalUpperBound[i] = dMaxI; + } + + } + + } + + for (int i = 0; i < dataLowerBound.GetRows(); i++){ + + if(m_vecEnforcementBounds[iVarIndex].strLowerBound == "l" || + m_vecEnforcementBounds[iVarIndex].strLowerBound == "lp"){ + dataLowerBound[i] = vecLocalLowerBound[i] - dataOutDouble[i]; + } + if(m_vecEnforcementBounds[iVarIndex].strUpperBound == "l" || + m_vecEnforcementBounds[iVarIndex].strUpperBound == "lp"){ + dataUpperBound[i] = vecLocalUpperBound[i] - dataOutDouble[i]; + } + } + + } + + CAAS(dataOutDouble,dataLowerBound,dataUpperBound,dMassDiff); + + } // Announce output mass double dTargetMass = 0.0; diff --git a/src/OfflineMap.h b/src/OfflineMap.h index b0dfb00..8ae74c9 100644 --- a/src/OfflineMap.h +++ b/src/OfflineMap.h @@ -132,6 +132,16 @@ class OfflineMap { ); protected: + + /// + /// Clip and assured sum function. + /// + void CAAS( + DataArray1D & dataCorrectedField, + DataArray1D & dataLowerBound, + DataArray1D & dataUpperBound, + double & dMass); + /// /// Initialize the coordinate arrays for a finite-volume mesh. /// @@ -590,12 +600,20 @@ class OfflineMap { /// Set the bounds enforcement parameters. /// void SetEnforcementBounds( - const std::string & strEnforcementBounds - ) { - ParseEnforceBounds( - strEnforcementBounds, - m_vecEnforcementBounds); - } + const std::string & strEnforcementBounds, + Mesh * pmeshSource, + Mesh * pmeshOverlap, + DataArray3D * pdataGLLNodesIn, + DataArray3D * pdataGLLNodesOut, + int iPin + ); + //void SetEnforcementBounds( + //const std::string & strEnforcementBounds + //) { + //ParseEnforceBounds( + //strEnforcementBounds, + //m_vecEnforcementBounds); + //} /// /// Clear the bounds enforcement parameters. @@ -746,6 +764,31 @@ class OfflineMap { /// Enforcement bounds used in this map. /// EnforceBoundsVector m_vecEnforcementBounds; + + /// + /// Pointer to the source mesh. + /// + Mesh * m_pmeshSource; + + /// + /// Pointer to the overlap mesh. + /// + Mesh * m_pmeshOverlap; + + /// + /// Pointer to the GLL node data for the source mesh. + /// + DataArray3D * m_pdataGLLNodesIn; + + /// + /// Pointer to the GLL node data for the target mesh. + /// + DataArray3D * m_pdataGLLNodesOut; + + /// + /// Order of the reconstruction polynomial for FV mesh. + /// + int m_iPin; }; diff --git a/src/TempestRemapAPI.h b/src/TempestRemapAPI.h index 9933d55..9741102 100644 --- a/src/TempestRemapAPI.h +++ b/src/TempestRemapAPI.h @@ -337,12 +337,20 @@ extern "C" { strOutputDataList(""), strVariables(""), strNColName("ncol"), + strEnforceBounds(""), + strInputMesh(""), + strOutputMesh(""), + strOverlapMesh(""), + nPin(0), + nPout(0), fOutputDouble(false), strOutputFormat("Netcdf4"), strPreserveVariables(""), fPreserveAll(false), dFillValueOverride(0.0), - strLogDir("") + strLogDir(""), + fContainsConcaveFaces(false), + fgll(false) { } public: @@ -379,7 +387,32 @@ extern "C" { /// /// A string describing how bounds enforcement should be performed. /// - std::string strEnforceBounds; + std::string strEnforceBounds; + + /// + /// Input mesh. + /// + std::string strInputMesh; + + /// + /// Output mesh. + /// + std::string strOutputMesh; + + /// + /// Overlap mesh. + /// + std::string strOverlapMesh; + + /// + /// Order of the polynomial on the source mesh. + /// + int nPin; + + /// + /// Order of the polynomial on the target mesh. + /// + int nPout; /// /// Output data using double precision. @@ -410,6 +443,16 @@ extern "C" { /// A directory for writing log files. /// std::string strLogDir; + + /// + /// A mesh contains concave faces. + /// + bool fContainsConcaveFaces; + + /// + /// The source mesh is finite element. + /// + bool fgll; }; ///