diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 00000000..6bc99cc9 Binary files /dev/null and b/.DS_Store differ diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 6194c91f..24731755 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -255,13 +255,6 @@ void polympo_setMeshNumVtxs(MPMesh_ptr p_mpmesh, int numVtxs){ p_mesh->setMeshVtxBasedFieldSize(); } -int polympo_getMeshNumVtxs(MPMesh_ptr p_mpmesh) { - checkMPMeshValid(p_mpmesh); //chech vailidity - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - int nVtxs = p_mesh->getNumVertices(); - return nVtxs; -} - void polympo_setMeshNumElms(MPMesh_ptr p_mpmesh, int numElms){ checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; @@ -272,14 +265,10 @@ void polympo_setMeshNumElms(MPMesh_ptr p_mpmesh, int numElms){ p_mesh->setNumElms(numElms); p_mesh->setElm2VtxConn(elm2Vtx); p_mesh->setElm2ElmConn(elm2Elm); + p_mesh->setMeshVtxBasedFieldSize(); } -int polympo_getMeshNumElms(MPMesh_ptr p_mpmesh) { - checkMPMeshValid(p_mpmesh); //chech vailidity - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - int nElms = p_mesh->getNumElements(); - return nElms; -} + void polympo_setMeshTypeGeneralPoly(MPMesh_ptr p_mpmesh){ //chech validity @@ -375,6 +364,20 @@ void polympo_setMeshElm2ElmConn(MPMesh_ptr p_mpmesh, int maxEdges, int nCells, i }); } +int polympo_getMeshNumVtxs(MPMesh_ptr p_mpmesh) { + checkMPMeshValid(p_mpmesh); //chech vailidity + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + int nVtxs = p_mesh->getNumVertices(); + return nVtxs; +} + +int polympo_getMeshNumElms(MPMesh_ptr p_mpmesh) { + checkMPMeshValid(p_mpmesh); //chech vailidity + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + int nElms = p_mesh->getNumElements(); + return nElms; +} + void polympo_setMeshVtxCoords(MPMesh_ptr p_mpmesh, int nVertices, double* xArray, double* yArray, double* zArray){ //chech validity checkMPMeshValid(p_mpmesh); @@ -413,6 +416,38 @@ void polympo_getMeshVtxCoords(MPMesh_ptr p_mpmesh, int nVertices, double* xArray } } +/* +void polympo_setMeshCellCenters(MPMesh_ptr p_mpmesh, int nCells, double* xArray, double* yArray, double* zArray){ + + checkMPMeshValid(p_mpmesh); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + + auto centersArray = p_mesh->getMeshField(); + auto h_centersArray = Kokkos::create_mirror_view(centersArray); + for (int i =0; ip_mesh; + + auto centersArray = p_mesh->getMeshField(); + auto h_centersArray = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),centersArray); + + for (int i=0;i @brief get the number of vertices from the mesh holding by polyMPO - !> @param mpMesh(in) mpMesh object - !> @param numVtxs(return)) the number of vertices - !--------------------------------------------------------------------------- - function polympo_getMeshNumVtxs(mpMesh) result(numVtxs) & - bind(C, NAME = 'polympo_getMeshNumVtxs') - use :: iso_c_binding - type(c_ptr), intent(in), value :: mpMesh - integer(c_int) :: numVtxs - end function - !--------------------------------------------------------------------------- !> @brief set the number of elements of the mesh !> modifies mesh topology polympo_startMeshFill required !> @param mpMesh(in/out) mpMesh object @@ -246,17 +235,6 @@ subroutine polympo_setMeshNumElms(mpMesh,numElms) & integer(c_int), value :: numElms end subroutine !--------------------------------------------------------------------------- - !> @brief get the number of elements from the mesh holding by polyMPO - !> @param mpMesh(in) mpMesh object - !> @param numVtxs(return)) the number of elements - !--------------------------------------------------------------------------- - function polympo_getMeshNumElms(mpMesh) result(numElms) & - bind(C, NAME = 'polympo_getMeshNumElms') - use :: iso_c_binding - type(c_ptr), intent(in), value :: mpMesh - integer(c_int) :: numElms - end function - !--------------------------------------------------------------------------- !> @brief set the polympo mesh element to vertices connectivity !> modifies mesh topology polympo_startMeshFill required !> @param mpmesh(in/out) MPMesh object @@ -285,6 +263,28 @@ subroutine polympo_setMeshElm2ElmConn(mpMesh, maxEdges, nCells, cellsOnCell) & type(c_ptr), intent(in), value :: cellsOnCell end subroutine !--------------------------------------------------------------------------- + !> @brief get the number of vertices from the mesh holding by polyMPO + !> @param mpMesh(in) mpMesh object + !> @param numVtxs(return)) the number of vertices + !--------------------------------------------------------------------------- + function polympo_getMeshNumVtxs(mpMesh) result(numVtxs) & + bind(C, NAME = 'polympo_getMeshNumVtxs') + use :: iso_c_binding + type(c_ptr), intent(in), value :: mpMesh + integer(c_int) :: numVtxs + end function + !--------------------------------------------------------------------------- + !> @brief get the number of elements from the mesh holding by polyMPO + !> @param mpMesh(in) mpMesh object + !> @param numVtxs(return)) the number of elements + !--------------------------------------------------------------------------- + function polympo_getMeshNumElms(mpMesh) result(numElms) & + bind(C, NAME = 'polympo_getMeshNumElms') + use :: iso_c_binding + type(c_ptr), intent(in), value :: mpMesh + integer(c_int) :: numElms + end function + !--------------------------------------------------------------------------- !> @brief set the polympo mesh number of edges per element !> modifies mesh topology polympo_startMeshFill required !> @param mpmesh(in/out) MPMesh object @@ -325,6 +325,32 @@ subroutine polympo_getMeshVtxCoords(mpMesh, nVertices, xArray, yArray, zArray) & type(c_ptr), value :: xArray, yArray, zArray end subroutine !--------------------------------------------------------------------------- + !> @brief set the polympo mesh cell centers + !> @param mpmesh(in/out) MPMesh object + !> @param nVertices(in) length of array in + !> @param x/y/zArray(in) the 1D arrays of cell coordinates + !--------------------------------------------------------------------------- + !subroutine polympo_setMeshCellCenters(mpMesh, nCenters, xArray, yArray, zArray) & + ! bind(C, NAME='polympo_setMeshCellCenters') + ! use :: iso_c_binding + ! type(c_ptr), value :: mpMesh + ! integer(c_int), value:: nVertices + ! type(c_ptr), value :: xArray, yArray, zArray + !end subroutine + !--------------------------------------------------------------------------- + !> @brief get the polympo mesh cell centers + !> @param mpmesh(in/out) MPMesh object + !> @param nVertices(in) length of array in, use for assertion + !> @param x/y/zArray(in/out) the 1D arrays of vertices coordinates + !--------------------------------------------------------------------------- + !subroutine polympo_getMeshCellCenters(mpMesh, nCenters, xArray, yArray, zArray) & + ! bind(C, NAME='polympo_getMeshCellCenters') + ! use :: iso_c_binding + ! type(c_ptr), value :: mpMesh + ! integer(c_int), value :: nVertices + ! type(c_ptr), value :: xArray, yArray, zArray + !end subroutine + !--------------------------------------------------------------------------- !> @brief set the spherical velocity increment mesh array !> from a host array !> @param mpmesh(in/out) MPMesh object diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp index fb409e9f..e55398d1 100644 --- a/src/pmpo_mesh.cpp +++ b/src/pmpo_mesh.cpp @@ -31,4 +31,13 @@ namespace polyMPO{ vtxOnSurfDispIncr_ = DoubleVec2dView(vtxOnSurfDispIncrMapEntry.second,numVtxs_); } + void Mesh::setMeshElmBasedFieldSize(){ + PMT_ALWAYS_ASSERT(meshEdit_); + + auto elmCoordsMapEntry = meshFields2TypeAndString.at(MeshF_ElmCenterXYZ); + PMT_ALWAYS_ASSERT(elmCoordsMapEntry.first == MeshFType_ElmBased); + elmCenterXYZ_ = DoubleVec3dView(elmCoordsMapEntry.second, numElms_); + + } + } // namespace polyMPO diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index 57cda9f6..df890863 100644 --- a/src/pmpo_mesh.hpp +++ b/src/pmpo_mesh.hpp @@ -22,7 +22,8 @@ enum MeshFieldIndex{ MeshF_VtxCoords, MeshF_Vel, MeshF_OnSurfVeloIncr, - MeshF_OnSurfDispIncr + MeshF_OnSurfDispIncr, + MeshF_ElmCenterXYZ }; enum MeshFieldType{ MeshFType_Invalid = -2, @@ -37,7 +38,8 @@ const std::map auto getMeshField(); void setMeshVtxBasedFieldSize(); + + void setMeshElmBasedFieldSize(); void setMeshEdit(bool meshEdit) { meshEdit_ = meshEdit; } //onec MeshType/GeomType is set to valid types, we can't change them anymore @@ -146,6 +152,9 @@ auto Mesh::getMeshField(){ else if constexpr (index==MeshF_OnSurfDispIncr){ return vtxOnSurfDispIncr_; } + else if constexpr (index==MeshF_ElmCenterXYZ){ + return elmCenterXYZ_; + } fprintf(stderr,"Mesh Field Index error!\n"); exit(1); } diff --git a/src/polympo.mod b/src/polympo.mod new file mode 100644 index 00000000..1cbbc66a Binary files /dev/null and b/src/polympo.mod differ diff --git a/test/readMPAS.f90 b/test/readMPAS.f90 index 94ee5d21..aff349e5 100644 --- a/test/readMPAS.f90 +++ b/test/readMPAS.f90 @@ -125,7 +125,8 @@ subroutine readMPASMesh(filename, maxEdges, vertexDegree, & nCells, nVertices, nEdgesOnCell, & onSphere, sphereRadius, & xVertex, yVertex, zVertex, & - verticesOnCell, cellsOnCell) + verticesOnCell, cellsOnCell)!,& + !cellCenters) use :: netcdf use :: iso_c_binding implicit none @@ -142,7 +143,7 @@ subroutine readMPASMesh(filename, maxEdges, vertexDegree, & integer :: ncid, status, nCellsID, nVerticesID, maxEdgesID, vertexDegreeID, & nEdgesOnCellID, xVertexID, yVertexID, zVertexID, & - verticesOnCellID, cellsOnCellID + verticesOnCellID, cellsOnCellID, cellCentersID status = nf90_open(path=trim(filename), mode=nf90_nowrite, ncid=ncid) if (status /= nf90_noerr) then @@ -211,9 +212,10 @@ subroutine readMPASMesh(filename, maxEdges, vertexDegree, & allocate(yVertex(nVertices)) allocate(zVertex(nVertices)) allocate(nEdgesOnCell(nCells)) + !allocate(cellCenters(nCells)) allocate(verticesOnCell(maxEdges, nCells)) allocate(cellsOnCell(maxEdges, nCells)) - + status = nf90_get_att(ncid, nf90_global, "on_a_sphere", onSphere) if (status /= nf90_noerr) then write(0, *) "readMPASMesh: Error on get attribute 'on_a_sphere'" @@ -315,6 +317,13 @@ subroutine readMPASMesh(filename, maxEdges, vertexDegree, & stop end if + !status = nf90_get_var(ncid, cellCentersID, cellCenters) + !if (status /= nf90_noerr) then + ! write(0, *) "readMPASMesh: Error on get var of 'cellCenters'" + ! write(0, *) trim(nf90_strerror(status)) + ! stop + !end if + status = nf90_close(ncid) end subroutine readMPASMesh subroutine setWithMPASMeshByFortran(mpMesh, fileName, n) bind(C, name="setWithMPASMeshByFortran")