From e25498c2090245049bf99b848cc51559c2d1f2c8 Mon Sep 17 00:00:00 2001 From: Tomer Avgil <132493674+tomeravgil@users.noreply.github.com> Date: Fri, 13 Oct 2023 09:54:49 -0400 Subject: [PATCH 1/3] Added MeshCellCenter Implementation --- .DS_Store | Bin 0 -> 6148 bytes src/pmpo_c.cpp | 57 +++++++++++++++++++++++++++-------- src/pmpo_c.h | 7 +++-- src/pmpo_fortran.f90 | 70 +++++++++++++++++++++++++++++-------------- src/polympo.mod | Bin 0 -> 5341 bytes 5 files changed, 98 insertions(+), 36 deletions(-) create mode 100644 .DS_Store create mode 100644 src/polympo.mod diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..6bc99cc9f0843d5cf5cd4ae0aff7341e5ce7f54f GIT binary patch literal 6148 zcmeHK%}T>S5Z<-brW7Fug&r5Y7L2tN#7l_v1&ruHr6#6mFlI}VnnNk%tS{t~_&m<+ zZop#BB6bFLzxmzGevtiPjB#%f9WiDz#w=)v9F+z^cVlSFBqMShBcF$f48i&c=O*^o z0l&S)a+a}&u>Ai0Nt_ps^U1eb?d_e8)vPcG6ZSGB$lSY4|-wfoEEaO@1k@z7Zf#PVba`e-~_ zt-99!!Qtt}^f`G+<(sCH1KUb=4c71uidoI8pQVXRAHh>+*I9(b05L!e5CfadfH@be z-e%K4t0xACfd&R}e-O|RU4xZIwRJ#;*Jq5k5m7+Lw*;av=o+jvLIi~CQb1kG%@c#` za_|e2=Nhau>T<@_%rK6bxqiHGH9Po)N@v{FNIfw?46HNI(x#2)|2h0JjgS2G60(Q^ zV&I=Kz}pjl;=!WK+4^I7c-9JN_s~!7`Tu0RZzzT>X7FetTf^%=vU={ ObP-U5P)7{>0s~)7Xi0kj literal 0 HcmV?d00001 diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 6194c91f..762671ef 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -255,12 +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); @@ -274,12 +268,7 @@ void polympo_setMeshNumElms(MPMesh_ptr p_mpmesh, int numElms){ p_mesh->setElm2ElmConn(elm2Elm); } -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,36 @@ void polympo_getMeshVtxCoords(MPMesh_ptr p_mpmesh, int nVertices, double* xArray } } +void polympo_setMeshCellCenters(MPMesh_ptr p_mpmesh, int nCenters, 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/polympo.mod b/src/polympo.mod new file mode 100644 index 0000000000000000000000000000000000000000..1cbbc66ac426ad7c7f857e4c51637ccec3514fad GIT binary patch literal 5341 zcmYjUXE+-SyGE^sqBPVd_TIZ_6h(~I4%J$*N9B#YUdPKkLx#7E@P2p&#ZPu7)|CQ2bR%&zbu-voW4-C+Z44)oIuh^-#rCkMSIt z*iYv+^cU-Szxls*X2~QDod>qO-yOh@Q3{DE#~ca{MD@3WhY$Saq3J?`f>FQ5yN`8I z%=ZC-XERsKUr!d{fm|F>sXN0gEsZztLh{Ee&-}zuo`hF3on6&huAHa%htnPOvSB%U z5BvGdwnh8*N45%RROJ8Q4Fxf3!$@2T)|}U@SXfoGPX7@y(lr!pQLoa8#`o=;3`tvI z-|s5J!o^yxbyS40wS#Z*Kzp(aniEOgR>WvxUPhgjh~!Nj%1s^pY6C0>w}f{7q#85$ zntMX>l6u6weL+@+FgMN11@E3qAIVrse!GLX|Ks4Bbp3KBwS3&%^D=S%t~=@BWd2e- zs$|zYepk3Ry~2wjfFL34rk4WWh|V+7&dJR8=aB*bZPX$Em#d&dXl~3mwC0oW8x{Mg zPYBZS3*O>sK>x(!R|*3yLQ@Y)Py5kasj~SCDK9v81Kj$3c@^C044O;x!bFi{r#M^{?$9dHXZ#-FE6)rLpV5MB3T|9W@>tHe_U&3$mOo^yDqaHa zZ7(Pr_}aAjKN3p(PA%ezE9M5+K=xmhf8Avk$)oj24OxEvq+BmU`Bvo7cyF#>aBD^-KtRtt#tR}(7!vG~ zqI(h&jxzz!c5w9*PY%f`(HNW$RmEXd2^W(h{dAbR6~1)A@#cPqTBN#$nEpzX_T8w9 z*}EMmKS*Vq$64t?9JS?VBB8{)x!e~_8=k%u<40u;X}D@*e)9nzD8?!!QJ2?kA(k(% zrjC!A%muhpKXmO@6qrcon71tXHyGYbdpj^=tV5UC@scN$_9~Jm3cwTAW{Ac6olyHmdnI=lCECLo z0#XY(_9rdx2SV(5i)3_zWeS$Bu<5NC zZmV%US#tN0w>5PCo168=(U=p9n}1N?UI9JaV!*>_3ohXuxtHp^&(~NeIS5OWeaPy5 zyfyeJf}jg4x|Sn?CcBs7KTmnO-d3Hu@w7m`*o1T(+B{Blp@9FLd`iq%5kUgl4viC^ zNZgGgzGlr=isPK`W?WG;VMn))RiN-2$SWQx?<1LAl!X< zuE9%yaQM4;LCPM!QYYMPt1btsvBIQMU)gVT(E*j@ zWWtYu2HbmLH?9+zvws_Y&-_39MxJPkc`zVj3HNz_VY}U*`isrnr|C1Q)X|l> zJ#KfW3MIAKIOEzZgt20YB_qP$4ze|sM@2ZxtHY@3@7sPGDAgrz68~6#jm z;JsZ5R~@sAH#CsEGd36otV}a~YRY5V?3uWPxOIDlV+F#nNtZNj@t=|g8zCGwy_0^Yjfl1Lp}{kwjVXqhFByG8<*aQaK^E?NJFG-t4SXqD{_?U`;~SkAdX z+V}7jHpA1gP8HcAJpJaDF#12CIlqhvUP&CEWRvw54F9t83mS{z?1KJT{nKSp@$C74 zPZywx;*p6Gqpm2UiBfJFvXG`7e=-G15q^d#!#lU8pwV^nLmoDU%U%LyGb6QEt$y@{ zCYoszkC>`c)n&DZuZVv-(~#s0`$Q4Axm|Y#4(Zy|FC|RK`9No-udEW8{MDRFc}O7q zdWu>kYd_XA7oUli)|`lNHUu|4-0&xf_?DW~3U3*Y@3Aa1^Xo%F9`5|zV-a5XFM07d zP%9`SPqB~>LmQNlZ(LuR!A*$RE3ea*x?-ldlkJeGV2$dEw8pq8H8A~jWi}kbSx?NF zfaKxF>HH&nuB%`?G=?T}q;DOKvWF{tc`e*Uo4G=Yj#2H>1q43x>jKQMW3CP{0*Q>4 z8Yra!Ic}L~h_{?~jW4}L@~I0Iq}&>Wc^1HW;m)Mgdw_yFf2%~yOUpAq3#^F)DA1*m66jvnd@f(>q~5?sr)4StE8;zh)6Vy4+!1$L3dM|t zVt}iu&S{TdWd*2y85Go3vXodfVqp9nta>DT&;C6%i<=j#&R%GD9NYq> zD7K>}TF)PsLn}rr$?$Pv?Qz{rIUQP>SRBSpYrKDwEdDtyZtSh${F!s=zL+3EYzjOI zMVK=KgT|EwFhq@INM_S#rZ4duzlBA$J*3`~C0`KbZ7q&AiR(@d<8pm;v9jgMtxvhW z$|p3~yzOlb_O+Ezmo4{5CSdK!e!VP4)sD=g{Si?8Fhm3FdG++T;J8$nEJZ9~*qyq< zP{{Eks0nPWK42}bMziNjz`DR)5AZ47h%2mn3G!!_rwsR5p?u`@`~m&gU;z8a3^HVs zOLJwARgZ@CGSo||5R^T@gZNl7TX^1#*e>8pp&G8z^RwwxUL6i_=o(I&_SuUagzd#r zPrZ@-sRU{TQMC)dige7J?&R-6q*httLj9@OuxPEVMeX{JE*H6 zq;O-&b0*VCRk((ZL2svhr2X62S$G8BV)QogY7wPZ8DKWGZBE+a8D|1{q@MC4h`+~h zFRtbJjU>-CZ6gXKdoa6cG|=Djps@q=Sxs#uRXgVqYHOj~;3$0Hwcv6|&cvJ_QB>GD z<^*1n6xE#^+8%W2$ydj?wFbw6GVYzhd=>Wk^@Bw{i|xE9Vhg|Zh@>>9t=se;7|^GB*k zRnS$!FGp^8-4S6z#SXIH+dEpR&Iz$bF%i7i-tIqI`y%-t!i(I3nL3>LQn#ewMfL;< z;u8tc^uzQ^FinOgf8*bSIIm3$8RrI9MgWcqu7i7@onAd&6!r}E*2(R`AR3<7(+&r? z+9-E>ZJuRJNBgQ8K?Ex4>ykr@y;O308ea$Nnd-A;3v$y5=xs+A%NZ6>d~COoz|dm& zTP~4b3qwqArnJUMbrD37m=rt9b|2|pxqNV*aXlL!ZWV3*l#MMn-=?}M#4iO*pn?P( zY(=Zc8RmG8qh(p23TEbtArrP!sHyNUhWI?QtdX&_)rp@RxEASwS=a%8A^6;C1*q*u z26X@X7Xh!gPeX5cn+@&mHikq{~TP-{4YJE44OxEi!)Ti#y4Q zXbO||C}5Di^qc99L5u- zzL+_fOReb0|pv8Lf3-{K*)k4M(fEtYv=r5EcKy(y@>8_cSkc z$}qsNBH20q`xi@j{37Yo(0Yl}(3oO2X)bOdI{SeM&00_VZ@lTvu~ZPBqP9+4zPTod zOxMz;5ZM>r-$L-t!BMhktIcK5t)Uwn0T?l#C;1hjAgs$!p@gAZG><8sSB6|_#X_c3 zm&M0VgMP`hEZONFnx$gjA|PQk426sOcrCm4~QN?=reZTT$y1DU6G%de1fq%*?<@&beC*W47!kU=2zt&dMg`yXMH*P&}_(7m!cR=--%kpLUX zdLCr$2-TgZAhIAd>yh)I)9i=Fi|G9~HZKUo=rL(?v$pmP2O}^^ePLw5U!R~1>50St zfJt-5B5x0<=Kf>Eid^l^!f@67jRZtH*}bso?vrWIvECA;=gpiOi$q6$-pY+%{68}5 zmY=Pa&%bZu=l)uOHU72geS&sE|4TLq`%C`snaJai-~Ts8iLuR@dWxX`jX^8WZ?qib z*nhBKS|c#&FnwWY%{%oeE;fr}P}iEI18cRXFs<+-qa1PqqO5xP5%%8YJH~W1>3rAS zn3x8WhDN~|KB5#Bo*qJ@@=1O zH)l4xhQ5xvS{@`aeokBosOC=TwKo3WTDV`{N|_G`gIYW=9clxZ&Eas@}-y33$k&^M=kY^-2W0heJYPE-INL|>tcR4Ov_MXk^@h_DT%9>_j1C;4hc8z8#~Tq*w{`CVv!iK#j>fVR#Q^*isRRNhx~zVNq%<9*NV-OBc0%)G$S zrsEn+i?;=B%el5u<{W*Q0|WDD2`sR8_&v$sUE*tikvZ`52P(&E=eU2eKmt1ndP3N3JS+=9}|=8D`pBv>2eIRc(; zUha46ZphuZH5fX3H0e1GaKbqp8V>Mk3QdN*$%eGObTux=2ltr>no4KQdD~l3!h2$_ o`=7qHuqwnVJkuC3l&>wdzN{WN{2gkv!?hFaw|C^9)Q^MnKWM>=U;qFB literal 0 HcmV?d00001 From d8a7441f3fd7bea9c9be69701acffd545f3a2217 Mon Sep 17 00:00:00 2001 From: Tomer Avgil <132493674+tomeravgil@users.noreply.github.com> Date: Thu, 19 Oct 2023 11:27:27 -0400 Subject: [PATCH 2/3] Added read files in fortran changed function args --- src/pmpo_c.cpp | 4 ++-- src/pmpo_c.h | 4 ++-- test/readMPAS.f90 | 15 ++++++++++++--- 3 files changed, 16 insertions(+), 7 deletions(-) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 762671ef..b50b72d2 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -416,7 +416,7 @@ void polympo_getMeshVtxCoords(MPMesh_ptr p_mpmesh, int nVertices, double* xArray } } -void polympo_setMeshCellCenters(MPMesh_ptr p_mpmesh, int nCenters, double* xArray, double* yArray, double* zArray){ +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; @@ -431,7 +431,7 @@ void polympo_setMeshCellCenters(MPMesh_ptr p_mpmesh, int nCenters, double* xArra Kokkos::deep_copy(centersArray, h_centersArray); } -void polympo_getMeshCellCenters(MPMesh_ptr p_mpmesh, int nCenters, double* xArray, double* yArray, double* zArray){ +void polympo_getMeshCellCenters(MPMesh_ptr p_mpmesh, int nCells, double* xArray, double* yArray, double* zArray){ checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; diff --git a/src/pmpo_c.h b/src/pmpo_c.h index 1acec339..9adbda27 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -46,8 +46,8 @@ int polympo_getMeshNumVtxs(MPMesh_ptr p_mpmesh); int polympo_getMeshNumElms(MPMesh_ptr p_mpmesh); //Mesh fields -void polympo_setMeshVtxCoords(MPMesh_ptr p_mpmesh, int nVertices, double* xArray, double* yArray, double* zArray); -void polympo_getMeshVtxCoords(MPMesh_ptr p_mpmesh, int nVertices, double* xArray, double* yArray, double* zArray); +void polympo_setMeshVtxCoords(MPMesh_ptr p_mpmesh, int nCells, double* xArray, double* yArray, double* zArray); +void polympo_getMeshVtxCoords(MPMesh_ptr p_mpmesh, int nCells, double* xArray, double* yArray, double* zArray); void polympo_setMeshCellCenters(MPMesh_ptr p_mpmesh, int nCenters, double* xArray, double* yArray, double* zArray); void polympo_getMeshCellCenters(MPMesh_ptr p_mpmesh, int nCenters, double* xArray, double* yArray, double* zArray); void polympo_setMeshOnSurfVeloIncr(MPMesh_ptr p_mpmesh, int nComps, int nVertices, double* array);//vec2d diff --git a/test/readMPAS.f90 b/test/readMPAS.f90 index 94ee5d21..85b7a087 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") From ac27d0d83375385f6e689388d3ea1d8fe8294991 Mon Sep 17 00:00:00 2001 From: Avgil Date: Wed, 25 Oct 2023 13:50:11 -0400 Subject: [PATCH 3/3] added new mesh field elmCenterXYZ --- src/pmpo_c.cpp | 4 +++- src/pmpo_c.h | 4 ++-- src/pmpo_fortran.f90 | 28 ++++++++++++++-------------- src/pmpo_mesh.cpp | 9 +++++++++ src/pmpo_mesh.hpp | 13 +++++++++++-- test/readMPAS.f90 | 18 +++++++++--------- 6 files changed, 48 insertions(+), 28 deletions(-) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index b50b72d2..24731755 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -255,7 +255,6 @@ void polympo_setMeshNumVtxs(MPMesh_ptr p_mpmesh, int numVtxs){ p_mesh->setMeshVtxBasedFieldSize(); } - void polympo_setMeshNumElms(MPMesh_ptr p_mpmesh, int numElms){ checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; @@ -266,6 +265,7 @@ void polympo_setMeshNumElms(MPMesh_ptr p_mpmesh, int numElms){ p_mesh->setNumElms(numElms); p_mesh->setElm2VtxConn(elm2Vtx); p_mesh->setElm2ElmConn(elm2Elm); + p_mesh->setMeshVtxBasedFieldSize(); } @@ -416,6 +416,7 @@ 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); @@ -445,6 +446,7 @@ void polympo_getMeshCellCenters(MPMesh_ptr p_mpmesh, int nCells, double* xArray, zArray[i] = h_centersArray(i,2); } } +*/ void polympo_setMeshOnSurfVeloIncr(MPMesh_ptr p_mpmesh, int nComps, int nVertices, double* array) { //check mpMesh is valid diff --git a/src/pmpo_c.h b/src/pmpo_c.h index 9adbda27..0ff12aaa 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -48,8 +48,8 @@ int polympo_getMeshNumElms(MPMesh_ptr p_mpmesh); //Mesh fields void polympo_setMeshVtxCoords(MPMesh_ptr p_mpmesh, int nCells, double* xArray, double* yArray, double* zArray); void polympo_getMeshVtxCoords(MPMesh_ptr p_mpmesh, int nCells, double* xArray, double* yArray, double* zArray); -void polympo_setMeshCellCenters(MPMesh_ptr p_mpmesh, int nCenters, double* xArray, double* yArray, double* zArray); -void polympo_getMeshCellCenters(MPMesh_ptr p_mpmesh, int nCenters, double* xArray, double* yArray, double* zArray); +//void polympo_setMeshCellCenters(MPMesh_ptr p_mpmesh, int nCenters, double* xArray, double* yArray, double* zArray); +//void polympo_getMeshCellCenters(MPMesh_ptr p_mpmesh, int nCenters, double* xArray, double* yArray, double* zArray); void polympo_setMeshOnSurfVeloIncr(MPMesh_ptr p_mpmesh, int nComps, int nVertices, double* array);//vec2d void polympo_getMeshOnSurfVeloIncr(MPMesh_ptr p_mpmesh, int nComps, int nVertices, double* array);//vec2d void polympo_setMeshOnSurfDispIncr(MPMesh_ptr p_mpmesh, int nComps, int nVertices, double* array);//vec2d diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index bb7c31b4..37c54750 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -330,26 +330,26 @@ subroutine polympo_getMeshVtxCoords(mpMesh, nVertices, xArray, yArray, zArray) & !> @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 + !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 + !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 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/test/readMPAS.f90 b/test/readMPAS.f90 index 85b7a087..aff349e5 100644 --- a/test/readMPAS.f90 +++ b/test/readMPAS.f90 @@ -125,8 +125,8 @@ subroutine readMPASMesh(filename, maxEdges, vertexDegree, & nCells, nVertices, nEdgesOnCell, & onSphere, sphereRadius, & xVertex, yVertex, zVertex, & - verticesOnCell, cellsOnCell, & - cellCenters) + verticesOnCell, cellsOnCell)!,& + !cellCenters) use :: netcdf use :: iso_c_binding implicit none @@ -212,7 +212,7 @@ subroutine readMPASMesh(filename, maxEdges, vertexDegree, & allocate(yVertex(nVertices)) allocate(zVertex(nVertices)) allocate(nEdgesOnCell(nCells)) - allocate(cellCenters(nCells)) + !allocate(cellCenters(nCells)) allocate(verticesOnCell(maxEdges, nCells)) allocate(cellsOnCell(maxEdges, nCells)) @@ -317,12 +317,12 @@ 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_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