From f12cc6abea19f80bacdb5f8939fcf2fed0861fb3 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 6 Sep 2024 14:42:36 +0200 Subject: [PATCH] Fix `angleEdge` at culled mesh boundaries Before this fix, `angleEdge` was simply copied from the base mesh to the culled mesh. With this change, `angleEdge` is flipped by pi radians if either the first cell on the given edge was removed (so the second cell is now the first cell after culling) or if only the second cell existed before (and after culling it is now the first cell). This is because `angleEdge` on a boundary edge always corresponds to a normal vector pointing from the first (and only) cell on the edge to the middle of the edge. --- .../mpas_cell_culler.cpp | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp index 3ac3859a1..082f40b5e 100755 --- a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp +++ b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_cell_culler.cpp @@ -1231,6 +1231,7 @@ int mapAndOutputEdgeFields( const string inputFilename, const string outputFilen double *dvEdgeOld, *dvEdgeNew; double *dcEdgeOld, *dcEdgeNew; double *angleEdgeOld, *angleEdgeNew; + bool *flipAngleEdge; bool hasWeightsOnEdge, hasAngleEdge; @@ -1239,6 +1240,7 @@ int mapAndOutputEdgeFields( const string inputFilename, const string outputFilen cellsOnEdgeNew = new int[nEdgesNew*2]; verticesOnEdgeOld = new int[nEdges*2]; verticesOnEdgeNew = new int[nEdgesNew*2]; + flipAngleEdge = new bool[nEdgesNew]; ncutil::get_var(inputFilename, "cellsOnEdge", cellsOnEdgeOld); ncutil::get_var(inputFilename, "verticesOnEdge", verticesOnEdgeOld); @@ -1269,18 +1271,21 @@ int mapAndOutputEdgeFields( const string inputFilename, const string outputFilen verticesOnEdgeNew[edgeMap.at(iEdge)*2] = vertexMap.at(vertex1) + 1; verticesOnEdgeNew[edgeMap.at(iEdge)*2+1] = vertexMap.at(vertex2) + 1; + flipAngleEdge[edgeMap.at(iEdge)] = false; } else if (cellMap.at(cell2) == -1){ cellsOnEdgeNew[edgeMap.at(iEdge)*2] = cellMap.at(cell1) + 1; cellsOnEdgeNew[edgeMap.at(iEdge)*2+1] = 0; verticesOnEdgeNew[edgeMap.at(iEdge)*2] = vertexMap.at(vertex1) + 1; verticesOnEdgeNew[edgeMap.at(iEdge)*2+1] = vertexMap.at(vertex2) + 1; + flipAngleEdge[edgeMap.at(iEdge)] = false; } else if (cellMap.at(cell1) == -1){ cellsOnEdgeNew[edgeMap.at(iEdge)*2] = cellMap.at(cell2) + 1; cellsOnEdgeNew[edgeMap.at(iEdge)*2+1] = 0; verticesOnEdgeNew[edgeMap.at(iEdge)*2] = vertexMap.at(vertex2) + 1; verticesOnEdgeNew[edgeMap.at(iEdge)*2+1] = vertexMap.at(vertex1) + 1; + flipAngleEdge[edgeMap.at(iEdge)] = true; } } else if(cell2 == -1){ if(cellMap.at(cell1) != -1){ @@ -1292,12 +1297,14 @@ int mapAndOutputEdgeFields( const string inputFilename, const string outputFilen } else { cout << "ERROR: Edge mask is 1, but has no cells." << endl; } + flipAngleEdge[edgeMap.at(iEdge)] = false; } else if(cell1 == -1){ cellsOnEdgeNew[edgeMap.at(iEdge)*2] = cellMap.at(cell2) + 1; cellsOnEdgeNew[edgeMap.at(iEdge)*2+1] = 0; verticesOnEdgeNew[edgeMap.at(iEdge)*2] = vertexMap.at(vertex2) + 1; verticesOnEdgeNew[edgeMap.at(iEdge)*2+1] = vertexMap.at(vertex1) + 1; + flipAngleEdge[edgeMap.at(iEdge)] = true; } else { cout << "ERROR: Edge mask is 1, but has no cells." << endl; } @@ -1441,7 +1448,13 @@ int mapAndOutputEdgeFields( const string inputFilename, const string outputFilen dvEdgeNew[edgeMap.at(iEdge)] = dvEdgeOld[iEdge]; dcEdgeNew[edgeMap.at(iEdge)] = dcEdgeOld[iEdge]; if(hasAngleEdge) { - angleEdgeNew[edgeMap.at(iEdge)] = angleEdgeOld[iEdge]; + if(flipAngleEdge[edgeMap.at(iEdge)]) { + angleEdgeNew[edgeMap.at(iEdge)] = + fmod(angleEdgeOld[iEdge], 2.0*M_PI) - M_PI; + } + else { + angleEdgeNew[edgeMap.at(iEdge)] = angleEdgeOld[iEdge]; + } } } } @@ -1465,6 +1478,7 @@ int mapAndOutputEdgeFields( const string inputFilename, const string outputFilen delete[] dcEdgeNew; delete[] angleEdgeOld; delete[] angleEdgeNew; + delete[] flipAngleEdge; return 0; }/*}}}*/