diff --git a/config/ax_blas.m4 b/config/ax_blas.m4 index 4c84053..97355b5 100644 --- a/config/ax_blas.m4 +++ b/config/ax_blas.m4 @@ -68,11 +68,11 @@ # special exception to the GPL to apply to your modified version as well. AC_DEFUN([AX_BLAS], [ -AC_PREREQ(2.50) +AC_PREREQ([2.69]) ax_blas_ok=no AC_ARG_WITH(blas, - [AC_HELP_STRING([--with-blas=], [use BLAS library ])]) + [AS_HELP_STRING([--with-blas=],[use BLAS library ])]) case $with_blas in yes | "") ;; no) acx_blas_ok=disable ;; @@ -115,7 +115,7 @@ fi # ax_blas_ok to yes, and execute ACTION-IF-FOUND. On failure, set ax_blas_ok # to no and execute ACTION-IF-NOT-FOUND. AC_DEFUN([_AX_BLAS], [ -AC_PREREQ(2.50) +AC_PREREQ([2.69]) ax_blas_ok=no # Get fortran linker names of BLAS functions to check for. diff --git a/configure.ac b/configure.ac index d0f3a1e..4b5c948 100644 --- a/configure.ac +++ b/configure.ac @@ -2,15 +2,15 @@ # Process this file with autoconf to produce a configure script. # Usage: autoreconf -fi # -AC_PREREQ([2.68]) -AC_INIT([TempestRemap], [2.1.2], [paullrich@ucdavis.edu]) +AC_PREREQ([2.69]) +AC_INIT([TempestRemap],[2.1.2],[paullrich@ucdavis.edu]) AC_CONFIG_AUX_DIR([config]) AC_CONFIG_MACRO_DIR([config]) #AC_CONFIG_SRCDIR([test/optimtest.m]) AC_CONFIG_HEADERS([src/TempestConfig.h]) AM_INIT_AUTOMAKE([foreign subdir-objects]) -AM_PROG_LIBTOOL +LT_INIT AM_SILENT_RULES([yes]) LT_INIT([disable-shared]) diff --git a/src/CoordTransforms.h b/src/CoordTransforms.h index 3dda545..286b3b5 100644 --- a/src/CoordTransforms.h +++ b/src/CoordTransforms.h @@ -356,6 +356,26 @@ inline void GnomonicProjection( } +// assume Tx, Tx, Tz is on unit sphere +inline void SimpleGnomonicProjection( + double x, + double y, + double z, + double Tx, + double Ty, + double Tz, + double & Gx, + double & Gy, + double & Gz +) +{ + double dotProduct = x*Tx + y*Ty + z*Tz; + _ASSERT(dotProduct!=0); + Gx = x/dotProduct; + Gy = y/dotProduct; + Gz = z/dotProduct; +} + /////////////////////////////////////////////////////////////////////////////// diff --git a/src/FiniteVolumeTools.cpp b/src/FiniteVolumeTools.cpp index 495118a..e8ccacb 100644 --- a/src/FiniteVolumeTools.cpp +++ b/src/FiniteVolumeTools.cpp @@ -411,19 +411,9 @@ bool DoesFaceContainPoint( ){ - //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; + NodeVector nodes; // Project face nodes onto plane tangent to the sphere at the sample point whose // coordinates are dX, dY, dZ @@ -431,72 +421,10 @@ bool DoesFaceContainPoint( 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; - - } - + nodes.push_back(nodeCurrent); } - return true; + return DoesFaceContainPoint(nodes, dX, dY, dZ); } @@ -510,89 +438,44 @@ bool DoesFaceContainPoint( ){ - //Convert sample point to lat/lon coordinates - - double dLonRad0 = 0.0; - double dLatRad0 = 0.0; - - XYZtoRLL_Rad(dX,dY,dZ,dLonRad0,dLatRad0); + //project nodesP to gnomonic plane tangent at dX, dY, dZ + // dX, dY, dZ is a point on unit sphere int iEdges = nodesP.size(); - NodeVector nodesPlane; + Node T(dX, dY, dZ); - double dCenterX = 0; - double dCenterY = 0; + NodeVector nP; // these nodes will be on plane tangent at dX, dY, dZ (the gnomonic plane) - // Project face nodes onto plane tangent to the sphere at the sample point whose - // coordinates are dX, dY, dZ + // Project face nodes onto plane tangent to the sphere at T + // T is on unit sphere 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 N = nodesP[i]; - Node nodeI(dGX,dGY,0.0); + Node G; - nodesPlane.push_back(nodeI); + SimpleGnomonicProjection(N.x, N.y, N.z, T.x, T.y, T.z, G.x, G.y, G.z); - //add contributions to planar face centroid - dCenterX += dGX; - dCenterY += dGY; + nP.push_back(G); } - //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; - - } - + // first determine orientation of polygon nP; plane is perpendicular to vector T + Node normal1 = CrossProduct(nP[1] - nP[0], nP[2] - nP[1]); + double orientation = DotProduct (T, normal1); + _ASSERT(orientation != 0.); + // assume that the polygon is convex + // also, we already verified that the point T is not close to any of the points nP + for (int i = 0; i< iEdges; i++) + { + Node N = nodesP[i]; + int nextI = (i+1)%iEdges; + Node N1 = nodesP[nextI]; + Node normal = CrossProduct(N - T, N1 - N); + double orient = DotProduct (T, normal); + if (orient * orientation < 0) + return false; } return true; @@ -764,13 +647,11 @@ void MatVectorMult(const DataArray2D & dMat, void TriangleLineIntersection( Node & nodeQ, NodeVector & nodesP, - DataArray1D & dCoeffs, - double & dCond + DataArray1D & dCoeffs ) { - _ASSERT(dCoeffs.GetRows() == 3); + //_ASSERT(dCoeffs.GetRows() == 3); - //Setup up columns of 3x3 matrix DataArray2D dInterpMat(3,3); @@ -813,14 +694,21 @@ void TriangleLineIntersection( DataArray1D dColSumAInv(3); + //LU factorization dgetrf_(&m, &n, &(dInterpMat(0,0)), &lda, &(iPIV[0]), &info); - + //Use LU factorization to find the inverse dgetri_(&n, &(dInterpMat(0,0)), &lda, &(iPIV[0]), &(dWork[0]), &lWork, &info); - + 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)"); + } + //A inverse column sums for (int j = 0; j < 3; j++) { @@ -832,7 +720,6 @@ void TriangleLineIntersection( } - //max column sums of A and A inverse double dMaxColSumA = dColSumA[0]; @@ -848,13 +735,12 @@ void TriangleLineIntersection( } } - 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)"); + //Product of max absolute column sum of A and A inverse for estimate of condition number of A + + double dCond = dMaxColSumAInv * dMaxColSumA; + + if (dCond > FVConditionNumberThreshold) { + Announce("WARNING: Poor conditioning in matrix (%1.15e);", dCond); } //Set right hand side of linear system @@ -869,7 +755,6 @@ void TriangleLineIntersection( MatVectorMult(dInterpMat, dRHS, dCoeffs); - } /////////////////////////////////////////////////////////////////////////////// @@ -877,13 +762,12 @@ void TriangleLineIntersection( void NewtonQuadrilateral( Node & nodeQ, NodeVector & nodesP, - DataArray1D & dCoeffs, - bool & fConverged + DataArray1D & dCoeffs ) { _ASSERT(dCoeffs.GetRows() == 3); - int iMaxIterations = 100; + int iMaxIterations = 150; //This algorithm is essentially the one that is used in ESMF @@ -918,6 +802,12 @@ void NewtonQuadrilateral( E[1] = nodeQ0.y - nodeQ.y; E[2] = nodeQ0.z - nodeQ.z; + DataArray1D iPIV; + iPIV.Allocate(3); + + DataArray1D dWork; + dWork.Allocate(3); + for (int i = 0; i < iMaxIterations; i++){ double dA = dCoeffs[0]; @@ -930,9 +820,9 @@ void NewtonQuadrilateral( 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) { + if (sqrt(F[0]*F[0]+F[1]*F[1]+F[2]*F[2]) < 1e-15 ) { - fConverged = true; + //fConverged = true; break; } @@ -956,13 +846,22 @@ void NewtonQuadrilateral( int lda = 3; int info; - DataArray1D iPIV; - iPIV.Allocate(3); + int lWork = 3; + + //Column sums for A and A inverse + DataArray1D dColSumA(3); + + for (int j = 0; j < 3; j++) { - DataArray1D dWork; - dWork.Allocate(3); + for (int k = 0; k < 3; k++) { + + dColSumA[j] += fabs(dJacobian(j,k)); + + } - int lWork = 3; + } + + DataArray1D dColSumAInv(3); //LU decomposition dgetrf_(&m, &n, &(dJacobian(0,0)), @@ -973,7 +872,41 @@ void NewtonQuadrilateral( &lda, &(iPIV[0]), &(dWork[0]), &lWork, &info); if (info != 0) { - _EXCEPTIONT("Mass matrix inversion error"); + _EXCEPTIONT("Matrix inversion error"); + } + + //A inverse column sums + for (int j = 0; j < 3; j++) { + + for (int k = 0; k < 3; k++) { + + dColSumAInv[j] += fabs(dJacobian(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]; + } + } + + //Product of max absolute column sum of A and A inverse for estimate of condition number of A + + double dCond = dMaxColSumAInv * dMaxColSumA; + + if (dCond > FVConditionNumberThreshold) { + Announce("WARNING: Poor conditioning in matrix (%1.15e);", dCond); } DataArray1D dDeltaX(3); @@ -1698,3 +1631,284 @@ void InvertFitArray_LeastSquares( } /////////////////////////////////////////////////////////////////////////////// + +void GeneralizedBarycentricCoordinates( + Node & nodeQ, + NodeVector & nodesFaceI, + std::vector & vecWeights +) { + + _ASSERT(vecWeights.size() == nodesFaceI.size()); + + int iEdges = nodesFaceI.size(); + + for (int i = 0; i < iEdges; i++){ + + Node nodeI = nodesFaceI[i]; + + Node nodeDiff = nodeI - nodeQ; + + double dMagNodeDiff = nodeDiff.Magnitude(); + + if ( dMagNodeDiff < 1e-8 ){ + + for (int j = 0; j < iEdges; j++){ + + if (j == i){ + + vecWeights[j] = 1; + + } + else { + + vecWeights[j] = 0; + + } + + } + + return; + + } + + } + + + //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); + + //Initialize weights to one + + for (int i = 0; i < iEdges; i++){ + + vecWeights[i] = 1.0; + + } + + 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 = nodesFaceI[m]; + + //c++ doesn't like the modulo operator when the argument is negative + + int iMMinusOne = m - 1; + + if (m == 0){ + + iMMinusOne = iEdges - 1; + + } + + Node nodeMMinusOne = nodesFaceI[iMMinusOne]; + + Node nodeMPlusOne = nodesFaceI[(m+1)%iEdges]; + + 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 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; + + } + + +} + +/////////////////////////////////////////////////////////////////////////////// + +void BilinearWeights( + Node & nodeQ, + NodeVector & nodesFaceI, + Face & faceFaceI, + std::vector & vecWeights, + std::vector & vecContributingFaces +) { + + DataArray1D dCoeffs(3); + + int iEdges = nodesFaceI.size(); + + for (int i = 0; i < iEdges; i++){ + + Node nodeI = nodesFaceI[i]; + + Node nodeDiff = nodeI - nodeQ; + + double dMagNodeDiff = nodeDiff.Magnitude(); + + if ( dMagNodeDiff < 1e-8 ){ + + vecWeights.resize(iEdges,0.0); + + vecContributingFaces.resize(iEdges); + + for (int j = 0; j < iEdges; j++){ + + vecContributingFaces[j] = faceFaceI[j]; + + if (j == i){ + + vecWeights[j] = 1.0; + + } + + } + + return; + + } + + } + + if( iEdges == 3 ){ + + vecWeights.resize(iEdges); + + vecContributingFaces.resize(iEdges); + + TriangleLineIntersection(nodeQ, nodesFaceI, dCoeffs); + + vecWeights[0] = 1 - dCoeffs[1] - dCoeffs[2]; + vecWeights[1] = dCoeffs[1]; + vecWeights[2] = dCoeffs[2]; + + for (int i = 0; i < iEdges; i++){ + + vecContributingFaces[i] = faceFaceI[i]; + + } + + } + + else if ( iEdges == 4 ){ + + vecWeights.resize(iEdges); + + vecContributingFaces.resize(iEdges); + + NewtonQuadrilateral(nodeQ, nodesFaceI, dCoeffs); + + for (int i = 0; i < iEdges; i++){ + + vecContributingFaces[i] = faceFaceI[i]; + + } + + vecWeights[0] = 1.0 - dCoeffs[0] - dCoeffs[1] + dCoeffs[0]*dCoeffs[1]; + vecWeights[1] = dCoeffs[0]*(1.0 - dCoeffs[1]); + vecWeights[2] = dCoeffs[0]*dCoeffs[1]; + vecWeights[3] = dCoeffs[1]*(1.0 - dCoeffs[0]); + + } + + else { + //Loop over the subtriangles until we find one that contains the sample point + + vecWeights.resize(3); + vecContributingFaces.resize(3); + + int nSubTriangles = iEdges - 2; + + for (int k = 0; k < nSubTriangles; k++) { + + Node & node0 = nodesFaceI[0]; + Node & nodeKPlusOne = nodesFaceI[k+1]; + Node & nodeKPlusTwo = nodesFaceI[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); + + vecWeights[0] = 1 - dCoeffs[1] - dCoeffs[2]; + vecWeights[1] = dCoeffs[1]; + vecWeights[2] = dCoeffs[2]; + + vecContributingFaces[0] = faceFaceI[0]; + vecContributingFaces[1] = faceFaceI[k+1]; + vecContributingFaces[2] = faceFaceI[k+2]; + + break; + + } + + } + + } + +} + +/////////////////////////////////////////////////////////////////////////////// diff --git a/src/FiniteVolumeTools.h b/src/FiniteVolumeTools.h index f9556b6..1ce60b0 100644 --- a/src/FiniteVolumeTools.h +++ b/src/FiniteVolumeTools.h @@ -225,8 +225,7 @@ void MatVectorMult( void TriangleLineIntersection( Node & nodeQ, NodeVector & nodesP, - DataArray1D & dCoeffs, - double & dCond + DataArray1D & dCoeffs ); /////////////////////////////////////////////////////////////////////////////// @@ -237,8 +236,7 @@ void TriangleLineIntersection( void NewtonQuadrilateral( Node & nodeQ, NodeVector & nodesP, - DataArray1D & dCoeffs, - bool & fConverged + DataArray1D & dCoeffs ); /////////////////////////////////////////////////////////////////////////////// @@ -303,5 +301,31 @@ void InvertFitArray_LeastSquares( /////////////////////////////////////////////////////////////////////////////// +/// +/// Get the generalized barycentric coordinates for a point. +/// + +void GeneralizedBarycentricCoordinates( + Node & nodeQ, + NodeVector & nodesFaceI, + std::vector & vecWeights +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Get the bilinear weights for a point in a polygon. +/// + +void BilinearWeights( + Node & nodeQ, + NodeVector & nodesFaceI, + Face & faceFaceI, + std::vector & vecWeights, + std::vector & vecContributingFaces +); + +/////////////////////////////////////////////////////////////////////////////// + #endif // _FINITEVOLUMETOOLS_H_ diff --git a/src/GridElements.cpp b/src/GridElements.cpp index 72b185a..dc445ce 100644 --- a/src/GridElements.cpp +++ b/src/GridElements.cpp @@ -3228,3 +3228,143 @@ void Dual( } /////////////////////////////////////////////////////////////////////////////// + +void ConstructLocalDualFace( + const Mesh & mesh, + NodeVector & meshCenters, + int & iNodeX, + Face & faceLocalDual, + NodeVector & nodesFaceLocal +) { + + std::set::const_iterator iterRevNode = mesh.revnodearray[iNodeX].begin(); + + for (; iterRevNode != mesh.revnodearray[iNodeX].end(); iterRevNode++) { + + Node node; + + for (int j = 0; j < mesh.faces[*iterRevNode].edges.size(); j++){ + + node.x += mesh.nodes[mesh.faces[*iterRevNode][j]].x; + node.y += mesh.nodes[mesh.faces[*iterRevNode][j]].y; + node.z += mesh.nodes[mesh.faces[*iterRevNode][j]].z; + + } + + node.x /= static_cast(mesh.faces[*iterRevNode].edges.size()); + node.y /= static_cast(mesh.faces[*iterRevNode].edges.size()); + node.z /= static_cast(mesh.faces[*iterRevNode].edges.size()); + + double dMag = node.Magnitude(); + + node.x /= dMag; + node.y /= dMag; + node.z /= dMag; + + nodesFaceLocal.push_back(node); + + + } + + const int nEdges = mesh.revnodearray[iNodeX].size(); + + + //Face face(mesh.revnodearray[iNodeX].size()); + Face faceTemp(mesh.revnodearray[iNodeX].size()); + + int ixNode = 0; + std::set::const_iterator iter = mesh.revnodearray[iNodeX].begin(); + for (; iter != mesh.revnodearray[iNodeX].end(); iter++) { + faceTemp.SetNode(ixNode, *iter); + ixNode++; + } + + // Reorient Faces + Node nodeCentral = mesh.nodes[iNodeX]; + + Node node0 = meshCenters[faceTemp[0]] - nodeCentral; + + Node nodeCross = CrossProduct(meshCenters[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 = meshCenters[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; + faceLocalDual.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]; + } + } + + faceLocalDual.SetNode(j, faceTemp[ixNextNode]); + dCurrentAngle = dNextAngle; + } + + + +} +/////////////////////////////////////////////////////////////////////////////// + +bool DoesMeshHaveHoles(const Mesh & meshInput) { + + if (meshInput.edgemap.size() == 0) { + _EXCEPTIONT("EdgeMap has not been calculated for meshInput"); + } + + //Loop over source mesh and check + for (int i = 0; i < meshInput.faces.size(); i++){ + + const Face & faceCurrent = meshInput.faces[i]; + + for (int k = 0; k < faceCurrent.edges.size(); k++){ + + const FacePair & facepair = + meshInput.edgemap.find(faceCurrent.edges[k])->second; + + if(facepair[0] == InvalidFace || facepair[1] == InvalidFace){ + + return true; + + } + + } + + } + + return false; + +} + +/////////////////////////////////////////////////////////////////////////////// diff --git a/src/GridElements.h b/src/GridElements.h index a856613..b6b3ca0 100644 --- a/src/GridElements.h +++ b/src/GridElements.h @@ -1142,6 +1142,18 @@ void Dual( /////////////////////////////////////////////////////////////////////////////// +/// +/// Construct the local dual face around a given node +/// +void ConstructLocalDualFace( + const Mesh & mesh, + NodeVector & meshCenters, + int & iNodeX, + Face & faceLocalDual, + NodeVector & nodesFaceLocal +); + +/////////////////////////////////////////////////////////////////////////////// /// /// Find the maximum edge length of an element. @@ -1175,5 +1187,12 @@ Real CalculateFaceAreaTriQuadratureSplit( /////////////////////////////////////////////////////////////////////////////// +/// +/// Determine whether the mesh has holes. +/// +bool DoesMeshHaveHoles(const Mesh & meshInput); + +/////////////////////////////////////////////////////////////////////////////// + #endif diff --git a/src/LinearRemapFV.cpp b/src/LinearRemapFV.cpp index 90ad481..29e3a4c 100644 --- a/src/LinearRemapFV.cpp +++ b/src/LinearRemapFV.cpp @@ -1593,13 +1593,12 @@ void LinearRemapGeneralizedBarycentric( 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; @@ -1611,199 +1610,453 @@ void LinearRemapGeneralizedBarycentric( // Get SparseMatrix representation of the OfflineMap SparseMatrix & smatMap = mapRemap.GetSparseMatrix(); + + //If there are no holes on the source mesh, use the (global) dual of the source. + + if(!DoesMeshHaveHoles(meshInput)){ + + //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 * kdDual = kd_create(3); + + // Vector of centers of the source mesh + for (int i = 0; i < meshInputDual.faces.size(); i++){ - //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++) { + const Face & face = meshInputDual.faces[i]; - // Output every 1000 overlap elements - if (ixFirst % 1000 == 0) { - Announce("Element %i/%i", ixFirst, meshOutput.faces.size()); + Node nodeCentroid = GetFaceCentroid(face, meshInputDual.nodes); + + nodeCentroid = nodeCentroid.Normalized(); + + kd_insert3( + kdDual, + nodeCentroid.x, + nodeCentroid.y, + nodeCentroid.z, + (void*)(&(meshInputDual.faces[i]))); + + } + + // Loop through all target faces + for (int ixFirst = 0; ixFirst < meshOutput.faces.size(); ixFirst++) { - // This Face - const Face & faceFirst = meshOutput.faces[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); + //Get the dual mesh face whose center is nearest to the target face - 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); + kdres * kdresTarget = + kd_nearest3( + kdDual, + nodeQ.x, + nodeQ.y, + nodeQ.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); + Face * pFace = (Face *)(kd_res_item_data(kdresTarget)); - nodesTriangleM.push_back(nodeMMinusOne); - nodesTriangleM.push_back(nodeM); - nodesTriangleM.push_back(nodeMPlusOne); + int iNearestFace = pFace - &(meshInputDual.faces[0]); - nodesSubTriangleM.push_back(nodeQ); - nodesSubTriangleM.push_back(nodeM); - nodesSubTriangleM.push_back(nodeMPlusOne); + // Find triangle that contains nodeQ. + // This is the face whose local index is iFaceFinal - double dSubAreaM = CalculateFaceArea(faceSubTriangleM, nodesSubTriangleM); + int iFaceFinal = 0; + + GetFaceThatContainsPoint(meshInputDual,iNearestFace,iFaceFinal,nodeQ.x,nodeQ.y,nodeQ.z); + + int iEdges = meshInputDual.faces[iFaceFinal].edges.size(); + + //Current dual mesh face + + Face faceDualFace = meshInputDual.faces[iFaceFinal]; + + //Vector of nodes on current dual face + + NodeVector nodesDualFace; + + for (int i = 0; i < iEdges; i++){ + + nodesDualFace.push_back(meshInputDual.nodes[faceDualFace[i]]); + + } + + //Vector of contributing weights + + std::vector vecWeights(iEdges); + + //Compute generalized barycentric coordinates + + GeneralizedBarycentricCoordinates(nodeQ, nodesDualFace, vecWeights); + + //Add weights to map + + for (int i = 0; i < iEdges; i++){ - vecTriangleSubAreas[m] = dSubAreaM; + int iContributingFaceI = meshInputDual.faces[iFaceFinal][i]; - double dTriangleAreaM = CalculateFaceArea(faceTriangleM, nodesTriangleM); + smatMap(ixFirst, iContributingFaceI) = vecWeights[i]; + + } + + } + + } + + //If the source mesh has holes, use local dual faces + + else{ + + //kd-tree of dual mesh centers + kdtree * kdTarget = kd_create(3); + + //Vector of source face centers + NodeVector vecSourceCenterArray; + + //Vector of target face centers + NodeVector vecTargetCenterArray; + + // Vector of centers of the source mesh + for (int i = 0; i < meshInput.faces.size(); i++){ + + const Face & face = meshInput.faces[i]; + + Node nodeCentroid = GetFaceCentroid(face, meshInput.nodes); + + nodeCentroid = nodeCentroid.Normalized(); + + vecSourceCenterArray.push_back(nodeCentroid); + + + } + + for (int i = 0; i < meshOutput.faces.size(); i++){ - vecTriangleAreas[m] = dTriangleAreaM; + const Face & face = meshOutput.faces[i]; + Node nodeCentroid = GetFaceCentroid(face, meshOutput.nodes); + + nodeCentroid = nodeCentroid.Normalized(); + kd_insert3( + kdTarget, + nodeCentroid.x, + nodeCentroid.y, + nodeCentroid.z, + (void*)(&(meshOutput.faces[i]))); + + vecTargetCenterArray.push_back(nodeCentroid); + + } + + std::set setAllTargets; + + // Loop through all source faces + + for (int ixFirst = 0; ixFirst < meshInput.faces.size(); ixFirst++) { + + // Output every 1000 source elements + if (ixFirst % 1000 == 0) { + + Announce("Element %i/%i", ixFirst, meshInput.faces.size()); + + } + + // This Face + const Face & faceFirst = meshInput.faces[ixFirst]; + + //Find the maximum distance from the current source face center to its nodes + + double dMaxNodeDist = 0.0; + + Node faceFirstCenter = vecSourceCenterArray[ixFirst]; + + int iEdges = faceFirst.edges.size(); + + for (int i = 0; i < iEdges; i++){ + + //Great circle distance is just the arccos of dot product + + Node nodeI = meshInput.nodes[faceFirst[i]]; + + double dDistToNodeI = acos(DotProduct(nodeI, faceFirstCenter)); + + if ( dDistToNodeI > dMaxNodeDist ){ + + dMaxNodeDist = dDistToNodeI; + + } + + } + + //Find all the target face centers that are within a given distance of the given source face center + + /* find points closest to the origin and within distance radius */ + //struct *presults; + + struct kdres *presults = kd_nearest_range3(kdTarget,faceFirstCenter.x, faceFirstCenter.y, + faceFirstCenter.z, dMaxNodeDist+.0001); + + + if( kd_res_size(presults) == 0 ){ + + continue; + + } + + //Nodes of the current face + + NodeVector nodesFaceFirst; + + for (int i = 0; i < faceFirst.edges.size(); i++){ + + nodesFaceFirst.push_back(meshInput.nodes[faceFirst[i]]); + + } + + std::vector vecTargetsInFace; + + while( !kd_res_end( presults ) ) { + + + /* get the data and position of the current result item */ + Face * pFace = (Face *)kd_res_item_data( presults ); + + int iNearestTarget = pFace - &(meshOutput.faces[0]); + + if(DoesFaceContainPoint(nodesFaceFirst, vecTargetCenterArray[iNearestTarget].x, + vecTargetCenterArray[iNearestTarget].y,vecTargetCenterArray[iNearestTarget].z)) { + + + if( setAllTargets.find(iNearestTarget) == setAllTargets.end() ){ + + setAllTargets.insert(iNearestTarget); + vecTargetsInFace.push_back(iNearestTarget); + } - - 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; - + + } + + /* go to the next entry */ + kd_res_next( presults ); + + } + + + kd_res_free( presults ); + + std::vector> vecContributingFaceI; + + std::vector> vecContributingWeights; + + + //Loop through all centers of the target faces in the given source face + for (int i = 0; i < vecTargetsInFace.size(); i++){ + + vecContributingFaceI.clear(); + vecContributingWeights.clear(); + + //Current target face center + + Node nodeQ = vecTargetCenterArray[vecTargetsInFace[i]]; + + std::vector vecWeightsCurrentFace(iEdges); + + GeneralizedBarycentricCoordinates(nodeQ, nodesFaceFirst, vecWeightsCurrentFace); + + //loop through all nodes on source face + + for (int k = 0; k < iEdges; k++){ + + //Global index of current node + + int iCurrentNode = meshInput.faces[ixFirst][k]; + + //Number of nodes on local dual test face + + int iEdgesTestK = meshInput.revnodearray[iCurrentNode].size(); + + Face faceLocalDual(iEdgesTestK); + + NodeVector nodesLocalDual; + + + //Current source face node shared by one source face + if ( iEdgesTestK == 1){ + + std::vector vecLocalWeights; + std::vector vecLocalContributingFaces; + + vecLocalWeights.push_back(vecWeightsCurrentFace[k]); + vecLocalContributingFaces.push_back(ixFirst); + + vecContributingWeights.push_back(vecLocalWeights); + vecContributingFaceI.push_back(vecLocalContributingFaces); + + + } + + //Current source node shared by two source faces + + else if (iEdgesTestK == 2){ + + std::set::const_iterator iterRevNodeI = meshInput.revnodearray[iCurrentNode].begin(); + std::set::const_iterator iterRevNodeII = meshInput.revnodearray[iCurrentNode].end(); + + std::vector vecLocalWeights; + std::vector vecLocalContributingFaces; + + vecLocalContributingFaces.push_back(*iterRevNodeI); + vecLocalContributingFaces.push_back(*iterRevNodeII); + + vecLocalWeights.push_back(vecWeightsCurrentFace[k]*0.5); + vecLocalWeights.push_back(vecWeightsCurrentFace[k]*0.5); + + vecContributingFaceI.push_back(vecLocalContributingFaces); + vecContributingWeights.push_back(vecLocalWeights); + + + } + + //Local dual face has at least three edges + + else{ + + //Construct local dual face + + Face faceLocalDual(iEdgesTestK); + + NodeVector nodesLocalDual; + + ConstructLocalDualFace(meshInput,vecSourceCenterArray,iCurrentNode,faceLocalDual,nodesLocalDual); + + //Put the dual face nodes in an ccw oriented vector + + NodeVector nodesLocalDualReordered; + + for (int j = 0; j < nodesLocalDual.size(); j++){ + + nodesLocalDualReordered.push_back(vecSourceCenterArray[faceLocalDual[j]]); + + } + + std::vector vecLocalWeights(iEdgesTestK); + std::vector vecLocalContributingFaces(iEdgesTestK); + + + //If the local dual face around the current node contains the target, then quit + + if (DoesFaceContainPoint(nodesLocalDualReordered, nodeQ.x, nodeQ.y, nodeQ.z)) { + + vecContributingWeights.clear(); + vecContributingFaceI.clear(); + + GeneralizedBarycentricCoordinates(nodeQ, nodesLocalDualReordered, vecLocalWeights); + + for ( int m = 0; m < iEdgesTestK; m++){ + + vecLocalContributingFaces[m] = faceLocalDual[m]; + + } + + vecContributingFaceI.push_back(vecLocalContributingFaces); + vecContributingWeights.push_back(vecLocalWeights); + + break; + + } + + else { + + //If the local dual face around the current node doesn't contain the target, + //use simple average or barycentric coorindates depending on whether the dual face + //contains the current node or not + + Node nodeCurrentNode = meshInput.nodes[iCurrentNode]; + + if (DoesFaceContainPoint(nodesLocalDualReordered, nodeCurrentNode.x, nodeCurrentNode.y, + nodeCurrentNode.z)){ + + //Compute generalized barycentric coordinate of current source face node wrt the local dual face + + GeneralizedBarycentricCoordinates(nodeCurrentNode, nodesLocalDualReordered, vecLocalWeights); + + + for ( int m = 0; m < iEdgesTestK; m++ ){ + + vecLocalContributingFaces[m] = faceLocalDual[m]; + + vecLocalWeights[m] *= vecWeightsCurrentFace[k]; + } - - if (l != m && l != iLMinusOne){ - - - vecWeights[m] *= vecTriangleSubAreas[iLMinusOne]*vecTriangleSubAreas[l]; - + + vecContributingFaceI.push_back(vecLocalContributingFaces); + vecContributingWeights.push_back(vecLocalWeights); + } - + + + else { + + for ( int m = 0; m < iEdgesTestK; m++ ){ + + vecLocalContributingFaces[m] = (faceLocalDual[m]); + + vecLocalWeights[m] = vecWeightsCurrentFace[k]*(1.0/iEdgesTestK); + + } + + vecContributingFaceI.push_back(vecLocalContributingFaces); + vecContributingWeights.push_back(vecLocalWeights); + + } + } - - 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]; - + + } + + //Add contributions to map + + for (int m = 0; m < vecContributingWeights.size(); m++){ + + for (int j = 0; j < vecContributingWeights[m].size(); j++){ + + int iContributingFace = vecContributingFaceI[m][j]; + + smatMap(vecTargetsInFace[i],iContributingFace) += vecContributingWeights[m][j]; + } + } - + + } + + } + } + } /////////////////////////////////////////////////////////////////////////////// @@ -2319,12 +2572,11 @@ void LinearRemapIntegratedBilinear( Dual(meshInputDual); //Construct edge map - meshInputDual.ConstructEdgeMap(); //kd-tree of dual mesh centers - kdtree * kdTarget = kd_create(3); + kdtree * kdTarget = kd_create(3); // Vector of centers of the source mesh for (int i = 0; i < meshInputDual.faces.size(); i++){ @@ -2463,7 +2715,7 @@ void LinearRemapIntegratedBilinear( } - TriangleLineIntersection(nodeQ, nodesP, dCoeffs, dCond); + TriangleLineIntersection(nodeQ, nodesP, dCoeffs); vecContributingFaceWeights.push_back(1 - dCoeffs[1] - dCoeffs[2]); vecContributingFaceWeights.push_back(dCoeffs[1]); @@ -2489,7 +2741,7 @@ void LinearRemapIntegratedBilinear( } - NewtonQuadrilateral(nodeQ, nodesP, dCoeffs, fConverged); + NewtonQuadrilateral(nodeQ, nodesP, dCoeffs); for (int i = 0; i < iEdges; i++){ @@ -2525,7 +2777,7 @@ void LinearRemapIntegratedBilinear( if( DoesFaceContainPoint(nodesP, nodeQ.x, nodeQ.y, nodeQ.z) ){ - TriangleLineIntersection(nodeQ, nodesP, dCoeffs, dCond); + TriangleLineIntersection(nodeQ, nodesP, dCoeffs); vecContributingFaceWeights.push_back(1.0 - dCoeffs[1] - dCoeffs[2]); vecContributingFaceWeights.push_back(dCoeffs[1]); @@ -2574,13 +2826,13 @@ void LinearRemapIntegratedBilinear( /////////////////////////////////////////////////////////////////////////////// -void LinearRemapBilinear( +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"); @@ -2597,198 +2849,493 @@ void LinearRemapBilinear( // Get SparseMatrix representation of the OfflineMap SparseMatrix & smatMap = mapRemap.GetSparseMatrix(); + + //If there are no holes on the source mesh, use the (global) dual of the source. + + if (!DoesMeshHaveHoles(meshInput)) { + + Mesh meshInputDual = meshInput; + + //Construct dual mesh + Dual(meshInputDual); + + //Construct edge map + + meshInputDual.ConstructEdgeMap(); + + //kd-tree of dual mesh centers + kdtree * kdTarget = kd_create(3); + for (int i = 0; i < meshInputDual.faces.size(); i++){ - std::cout << "\nCreating the dual mesh\n"; - 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(); + const Face & face = meshInputDual.faces[i]; - //Get the dual mesh face whose center is nearest to the target face + Node nodeCentroid = GetFaceCentroid(face, meshInputDual.nodes); + + nodeCentroid = nodeCentroid.Normalized(); - kdres * kdresTarget = - kd_nearest3( + kd_insert3( 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); - + 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()); } - - TriangleLineIntersection(nodeQ, nodesP, dCoeffs, dCond); - - vecContributingFaceWeights.push_back(1 - dCoeffs[1] - dCoeffs[2]); - vecContributingFaceWeights.push_back(dCoeffs[1]); - vecContributingFaceWeights.push_back(dCoeffs[2]); - + + // 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(); + + //NodeVector nodesCurrentFace(iEdges); + NodeVector nodesCurrentFace; + for (int i = 0; i < iEdges; i++){ - - vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][i]); - + + nodesCurrentFace.push_back(meshInputDual.nodes[meshInputDual.faces[iFaceFinal][i]]); + } - + + std::vector vecContributingFaceWeights; + + std::vector vecContributingFaceI; + + DataArray1D dCoeffs(3); + + BilinearWeights(nodeQ,nodesCurrentFace,meshInputDual.faces[iFaceFinal],vecContributingFaceWeights,vecContributingFaceI); + + // Contribution of each point to the map + for (int i = 0; i < vecContributingFaceI.size(); i++){ + + if( vecContributingFaceWeights[i] < -1e-12 || vecContributingFaceWeights[i] > 1+1e-12 ){ + std::cout << "\nFound weight value = " << vecContributingFaceWeights[i] << std::endl; + _EXCEPTIONT("Non-monotone weight"); + } + + int iContributingFaceI = vecContributingFaceI[i]; + + smatMap(ixFirst, iContributingFaceI) = vecContributingFaceWeights[i]; + + } + } + + } + + //If the source mesh has holes, use local dual faces + + else { + + //kd-tree of dual mesh centers + kdtree * kdTarget = kd_create(3); + + //Vector of source face centers + NodeVector vecSourceCenterArray; + + //Vector of target face centers + NodeVector vecTargetCenterArray; + + // Vector of centers of the source mesh + for (int i = 0; i < meshInput.faces.size(); i++){ + + const Face & face = meshInput.faces[i]; + + Node nodeCentroid = GetFaceCentroid(face, meshInput.nodes); + + nodeCentroid = nodeCentroid.Normalized(); + + vecSourceCenterArray.push_back(nodeCentroid); + + + } + + for (int i = 0; i < meshOutput.faces.size(); i++){ - else if ( iEdges == 4 ){ - - NodeVector nodesP; - - for (int i = 0; i < 4; i++){ - - Node nodeI = meshInputDual.nodes[meshInputDual.faces[iFaceFinal][i]]; + const Face & face = meshOutput.faces[i]; - nodesP.push_back(nodeI); + Node nodeCentroid = GetFaceCentroid(face, meshOutput.nodes); + + nodeCentroid = nodeCentroid.Normalized(); + kd_insert3( + kdTarget, + nodeCentroid.x, + nodeCentroid.y, + nodeCentroid.z, + (void*)(&(meshOutput.faces[i]))); + + vecTargetCenterArray.push_back(nodeCentroid); + + } + + std::set setAllTargets; + + // Loop through all source faces + for (int ixFirst = 0; ixFirst < meshInput.faces.size(); ixFirst++) { + + // Output every 1000 elements + if (ixFirst % 1000 == 0) { + Announce("Element %i/%i", ixFirst, meshInput.faces.size()); } - - NewtonQuadrilateral(nodeQ, nodesP, dCoeffs, fConverged); - + + // This Face + const Face & faceFirst = meshInput.faces[ixFirst]; + + //Find the maximum distance from the current source face center to its nodes + + double dMaxNodeDist = 0.0; + + Node faceFirstCenter = vecSourceCenterArray[ixFirst]; + + int iEdges = faceFirst.edges.size(); + for (int i = 0; i < iEdges; i++){ - - vecContributingFaceI.push_back(meshInputDual.faces[iFaceFinal][i]); - + + //Great circle distance is just the arccos of dot product + + Node nodeI = meshInput.nodes[faceFirst[i]]; + + double dDistToNodeI = acos(DotProduct(nodeI, faceFirstCenter)); + + if ( dDistToNodeI > dMaxNodeDist ){ + + dMaxNodeDist = dDistToNodeI; + + } + } - - 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; - + + //Find all the target face centers that are within a given distance of the given source face center + + /* find points closest to the origin and within distance radius */ + //struct *presults; + + struct kdres *presults = kd_nearest_range3(kdTarget,faceFirstCenter.x, faceFirstCenter.y, + faceFirstCenter.z, dMaxNodeDist+.0001); + + + if( kd_res_size(presults) == 0 ){ + + continue; + + } + + //Nodes of the current face + + NodeVector nodesFaceFirst; + + for (int i = 0; i < faceFirst.edges.size(); i++){ + + nodesFaceFirst.push_back(meshInput.nodes[faceFirst[i]]); + + } + + //set of target faces that have already been used. This is necessry to avoid double + //counting of target face centers that lie on a source face edge/node + + std::vector vecTargetsInFace; + + while( !kd_res_end( presults ) ) { + + + /* get the data and position of the current result item */ + Face * pFace = (Face *)kd_res_item_data( presults ); + + int iNearestTarget = pFace - &(meshOutput.faces[0]); + + if(DoesFaceContainPoint(nodesFaceFirst, vecTargetCenterArray[iNearestTarget].x, + vecTargetCenterArray[iNearestTarget].y,vecTargetCenterArray[iNearestTarget].z)) { + + if( setAllTargets.find(iNearestTarget) == setAllTargets.end() ){ + + setAllTargets.insert(iNearestTarget); + vecTargetsInFace.push_back(iNearestTarget); + + } + } - + + /* go to the next entry */ + kd_res_next( presults ); + } - - } - - // 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"); - + + kd_res_free( presults ); + + //Loop through all centers of the target faces in the given source face + for (int i = 0; i < vecTargetsInFace.size(); i++){ + + //These are the weights that are input into the spares matrix map + + std::vector> vecContributingFaceI; + + std::vector> vecContributingFaceWeights; + + //These are the temporary weights associated with each source face node. + //They are used if the current target face isn't in any of the local dual faces + + std::vector> vecLocalFacesTemp; + + std::vector> vecLocalWeightsTemp; + + //Current target face center + + Node nodeQ = vecTargetCenterArray[vecTargetsInFace[i]]; + + //Bilinear weights of current target center wrt to nodes of current face + + std::vector vecWeightsCurrentFace(iEdges); + + //Number of nodes we've visited + int iCountK = 0; + + //loop through all nodes on source face + + for (int k = 0; k < iEdges; k++){ + + //Global index of current node + + int iCurrentNode = meshInput.faces[ixFirst][k]; + + //Number of nodes on local dual test face + + int iEdgesTestK = meshInput.revnodearray[iCurrentNode].size(); + + Face faceLocalDual(iEdgesTestK); + + NodeVector nodesLocalDual; + + //Current source face node shared by one source face + if ( iEdgesTestK == 1){ + + std::vector vecLocalWeights; + std::vector vecLocalContributingFaces; + + vecLocalWeights.push_back(1); + vecLocalContributingFaces.push_back(ixFirst); + + vecLocalWeightsTemp.push_back(vecLocalWeights); + vecLocalFacesTemp.push_back(vecLocalContributingFaces); + + iCountK += 1; + + + } + + //Current source node shared by two source faces + + else if (iEdgesTestK == 2){ + + std::set::const_iterator iterRevNodeI = meshInput.revnodearray[iCurrentNode].begin(); + std::set::const_iterator iterRevNodeII = meshInput.revnodearray[iCurrentNode].end(); + + std::vector vecLocalWeights; + std::vector vecLocalContributingFaces; + + vecLocalContributingFaces.push_back(*iterRevNodeI); + vecLocalContributingFaces.push_back(*iterRevNodeII); + + vecLocalWeights.push_back(0.5); + vecLocalWeights.push_back(0.5); + + vecLocalFacesTemp.push_back(vecLocalContributingFaces); + vecLocalWeightsTemp.push_back(vecLocalWeights); + + iCountK += 1; + + + } + + //Local dual face has at least three edges, i.e. current source node shared by more than 3 faces + + else{ + + Face faceLocalDual(iEdgesTestK); + + NodeVector nodesLocalDual; + + ConstructLocalDualFace(meshInput,vecSourceCenterArray,iCurrentNode,faceLocalDual,nodesLocalDual); + + //Put the dual face nodes in an ccw oriented vector + + NodeVector nodesLocalDualReordered; + + for (int j = 0; j < nodesLocalDual.size(); j++){ + + nodesLocalDualReordered.push_back(vecSourceCenterArray[faceLocalDual[j]]); + + } + + std::vector vecLocalWeights; + std::vector vecLocalContributingFaces; + + //If the local dual face around the current node contains the target, then quit + + if (DoesFaceContainPoint(nodesLocalDualReordered, nodeQ.x, nodeQ.y, nodeQ.z)) { + + vecContributingFaceWeights.clear(); + vecContributingFaceI.clear(); + + //Bilinear weights + + BilinearWeights(nodeQ,nodesLocalDualReordered,faceLocalDual,vecLocalWeights,vecLocalContributingFaces); + + vecContributingFaceI.push_back(vecLocalContributingFaces); + vecContributingFaceWeights.push_back(vecLocalWeights); + + break; + + } + + else { + + //Add to count + iCountK += 1; + + //If the local dual face around the current node doesn't contain the target, + //use simple average or bilinear depending on whether the dual face + //contains the current node or not + + Node nodeCurrentNode = meshInput.nodes[iCurrentNode]; + + if (DoesFaceContainPoint(nodesLocalDualReordered, nodeCurrentNode.x, nodeCurrentNode.y, + nodeCurrentNode.z)){ + + //Compute bilinear weights of current source face node wrt the local dual face + + BilinearWeights(nodeCurrentNode,nodesLocalDualReordered,faceLocalDual,vecLocalWeights, + vecLocalContributingFaces); + + + vecLocalFacesTemp.push_back(vecLocalContributingFaces); + vecLocalWeightsTemp.push_back(vecLocalWeights); + + } + + else { + + vecLocalWeights.resize(iEdgesTestK); + vecLocalContributingFaces.resize(iEdgesTestK); + + for ( int m = 0; m < iEdgesTestK; m++ ){ + + vecLocalContributingFaces[m] = (faceLocalDual[m]); + + vecLocalWeights[m] = (1.0/iEdgesTestK); + + } + + vecLocalFacesTemp.push_back(vecLocalContributingFaces); + vecLocalWeightsTemp.push_back(vecLocalWeights); + + + } + + } + + } + + } + + //If we've visited very node of the current source face, we need to do another interpolation + //from these nodes to the current target face center + + if( iCountK == iEdges ){ + + //New vectors to store bilinear weights of current source face wrt current target face center + + std::vector vecLocalFaces; + std::vector vecLocalWeights; + + //New face with local ordering + Face faceCurrentLocal(iEdges); + + for ( int r = 0; r < iEdges; r++ ) { + + faceCurrentLocal.SetNode(r,r); + + } + + //compute bilinear weights + + BilinearWeights(nodeQ, nodesFaceFirst, faceCurrentLocal, vecLocalWeights, vecLocalFaces); + + for (int r = 0; r < vecLocalFaces.size(); r++) { + + for (int q = 0; q < vecLocalWeightsTemp[vecLocalFaces[r]].size(); q++ ){ + + vecLocalWeightsTemp[vecLocalFaces[r]][q] *= vecLocalWeights[r]; + + } + + //Push back to weight and face vectors + vecContributingFaceI.push_back(vecLocalFacesTemp[vecLocalFaces[r]]); + vecContributingFaceWeights.push_back(vecLocalWeightsTemp[vecLocalFaces[r]]); + + } + + } + + for (int m = 0; m < vecContributingFaceWeights.size(); m++){ + + for (int j = 0; j < vecContributingFaceWeights[m].size(); j++){ + + int iContributingFace = vecContributingFaceI[m][j]; + + if( vecContributingFaceWeights[m][j] < -1e-12 || vecContributingFaceWeights[m][j] > 1+1e-12 ){ + std::cout << "\nFound weight value = " << vecContributingFaceWeights[m][j] << std::endl; + _EXCEPTIONT("Non-monotone weight"); + + } + + smatMap(vecTargetsInFace[i],iContributingFace) += vecContributingFaceWeights[m][j]; + + } + + } + } - - int iContributingFaceI = vecContributingFaceI[i]; - - smatMap(ixFirst, iContributingFaceI) = vecContributingFaceWeights[i]; - + } } } + /////////////////////////////////////////////////////////////////////////////// void ForceIntArrayConsistencyConservation(