diff --git a/src/mesh/mesh.h b/src/mesh/mesh.h index 5738f24f3..9f0091eb8 100644 --- a/src/mesh/mesh.h +++ b/src/mesh/mesh.h @@ -102,6 +102,9 @@ struct mesh_t { int dim; int Nverts, Nfaces, NfaceVertices; + int uniform = 0; + int isRefine = 0; + int Nbid; int cht; @@ -260,7 +263,7 @@ struct mesh_t { }; std::pair createMesh(MPI_Comm comm, int N, int cubN, bool cht, occa::properties &kernelInfo); -mesh_t *createMeshMG(mesh_t *_mesh, int Nc); +mesh_t *createMeshMG(mesh_t *_mesh, int Nc, int isRefine = 0); occa::properties meshKernelProperties(int N); // serial sort diff --git a/src/mesh/meshBasisHex3D.cpp b/src/mesh/meshBasisHex3D.cpp index dda8e84d3..da28145e9 100644 --- a/src/mesh/meshBasisHex3D.cpp +++ b/src/mesh/meshBasisHex3D.cpp @@ -31,12 +31,16 @@ // ------------------------------------------------------------------------ // HEX 3D NODES // ------------------------------------------------------------------------ -void NodesHex3D(int _N, dfloat *_r, dfloat *_s, dfloat *_t) +void NodesHex3D(int _N, dfloat *_r, dfloat *_s, dfloat *_t, int uniform) { int _Nq = _N + 1; dfloat *r1D = (dfloat *)malloc(_Nq * sizeof(dfloat)); - JacobiGLL(_N, r1D); // Gauss-Legendre-Lobatto nodes + if (unidorm) { + EquispacedNodes1D(_N, r1D); + } else { + JacobiGLL(_N, r1D); // Gauss-Legendre-Lobatto nodes + } // Tensor product for (int k = 0; k < _Nq; k++) { diff --git a/src/mesh/meshGlobalIds.cpp b/src/mesh/meshGlobalIds.cpp index 7ad80c269..09ce9904f 100644 --- a/src/mesh/meshGlobalIds.cpp +++ b/src/mesh/meshGlobalIds.cpp @@ -22,7 +22,7 @@ void meshNekParallelConnectNodes(mesh_t* mesh) dlong localNodeCount = mesh->Np * mesh->Nelements; mesh->globalIds = (hlong*) calloc(localNodeCount, sizeof(hlong)); - hlong ngv = nek::set_glo_num(mesh->N + 1, mesh->cht); + hlong ngv = nek::set_glo_num(mesh->N + 1, mesh->Nelements, mesh->isRefine); for(dlong id = 0; id < localNodeCount; ++id) mesh->globalIds[id] = nekData.glo_num[id]; } diff --git a/src/mesh/meshLoadReferenceNodesHex3D.cpp b/src/mesh/meshLoadReferenceNodesHex3D.cpp index e5d5ce2da..e77c7cc42 100644 --- a/src/mesh/meshLoadReferenceNodesHex3D.cpp +++ b/src/mesh/meshLoadReferenceNodesHex3D.cpp @@ -28,7 +28,7 @@ #include "mesh3D.h" #define NODE_GEN -void meshLoadReferenceNodesHex3D(mesh_t *mesh, int N, int cubN) +void meshLoadReferenceNodesHex3D(mesh_t *mesh, int N, int cubN, int uniform) { mesh->N = N; mesh->Nq = N + 1; @@ -51,7 +51,7 @@ void meshLoadReferenceNodesHex3D(mesh_t *mesh, int N, int cubN) mesh->r = (dfloat*) malloc(mesh->Np * sizeof(dfloat)); mesh->s = (dfloat*) malloc(mesh->Np * sizeof(dfloat)); mesh->t = (dfloat*) malloc(mesh->Np * sizeof(dfloat)); - NodesHex3D(mesh->N, mesh->r, mesh->s, mesh->t); + NodesHex3D(mesh->N, mesh->r, mesh->s, mesh->t, uniform); mesh->faceNodes = (int*) malloc(mesh->Nfaces * mesh->Nfp * sizeof(int)); FaceNodesHex3D(mesh->N, mesh->r, mesh->s, mesh->t, mesh->faceNodes); @@ -59,10 +59,14 @@ void meshLoadReferenceNodesHex3D(mesh_t *mesh, int N, int cubN) //GLL quadrature mesh->gllz = (dfloat*) malloc((mesh->N + 1) * sizeof(dfloat)); mesh->gllw = (dfloat*) malloc((mesh->N + 1) * sizeof(dfloat)); - JacobiGLL(mesh->N, mesh->gllz, mesh->gllw); + if (uniform) { + EquispacedNodes1D(mesh->N, mesh->gllz, mesh->gllw); + } else { + JacobiGLL(mesh->N, mesh->gllz, mesh->gllw); + } mesh->D = (dfloat*) malloc(mesh->Nq * mesh->Nq * sizeof(dfloat)); - Dmatrix1D(mesh->N, mesh->Nq, mesh->gllz, mesh->Nq, mesh->gllz, mesh->D); + Dmatrix1D(mesh->N, mesh->Nq, mesh->gllz, mesh->Nq, mesh->gllz, mesh->D); // TODO uniform? this is not used for now multilevel gMG mesh->DW = (dfloat*) malloc(mesh->Nq * mesh->Nq * sizeof(dfloat)); DWmatrix1D(mesh->N, mesh->D, mesh->DW); diff --git a/src/mesh/meshNekReader.cpp b/src/mesh/meshNekReader.cpp index 57a1eecb6..87e0397d6 100644 --- a/src/mesh/meshNekReader.cpp +++ b/src/mesh/meshNekReader.cpp @@ -34,7 +34,7 @@ void meshNekReaderHex3D(int N, mesh_t *mesh) const int faceMap[] = {1, 2, 3, 4, 0, 5}; // generate element vertex numbering - mesh->Nnodes = nek::set_glo_num(2, mesh->cht); + mesh->Nnodes = nek::set_glo_num(2, mesh->Nelements, mesh->isRefine); mesh->EToV = (hlong *)calloc(mesh->Nelements * mesh->Nverts, sizeof(hlong)); for (int e = 0; e < mesh->Nelements; ++e) { diff --git a/src/mesh/meshPhysicalNodesHex3D.cpp b/src/mesh/meshPhysicalNodesHex3D.cpp index 1c6300555..8f9201e6d 100644 --- a/src/mesh/meshPhysicalNodesHex3D.cpp +++ b/src/mesh/meshPhysicalNodesHex3D.cpp @@ -33,7 +33,7 @@ void meshPhysicalNodesHex3D(mesh_t *mesh) std::vector ym(mesh->Nlocal); std::vector zm(mesh->Nlocal); - nek::xm1N(xm.data(), ym.data(), zm.data(), mesh->N, mesh->Nelements); + nek::xm1N(xm.data(), ym.data(), zm.data(), mesh->N, mesh->Nelements, mesh->uniform); mesh->o_x = platform->device.malloc(mesh->Nlocal); diff --git a/src/mesh/meshSetup.cpp b/src/mesh/meshSetup.cpp index b591fa6e7..6b4a65a2b 100644 --- a/src/mesh/meshSetup.cpp +++ b/src/mesh/meshSetup.cpp @@ -333,10 +333,12 @@ std::pair createMesh(MPI_Comm comm, int N, int cubN, bool cht, return {mesh, meshV}; } -mesh_t *createMeshMG(mesh_t *_mesh, int Nc) +mesh_t *createMeshMG(mesh_t *_mesh, int Nc, int isRefine_) { mesh_t *mesh = new mesh_t(); memcpy(mesh, _mesh, sizeof(mesh_t)); + isRefine = isRefine_; + uniform = (isRefine) ? 1 : 0; const int cubN = 0; meshLoadReferenceNodesHex3D(mesh, Nc, cubN); diff --git a/src/nekInterface/nekInterface.f b/src/nekInterface/nekInterface.f index 97da96ce9..7de5e67a6 100644 --- a/src/nekInterface/nekInterface.f +++ b/src/nekInterface/nekInterface.f @@ -808,22 +808,26 @@ integer function nekf_nbid(isTmsh) return end c----------------------------------------------------------------------- - integer*8 function nekf_set_vert(nx, isTmsh) + integer*8 function nekf_set_vert(nx, nel, isRefine) include 'SIZE' include 'TOTAL' include 'NEKINTF' - integer npts, isTmsh + integer nx, nel, isRefine + common /ivrtx0/ vertex_orig ((2**ldim),lelt) + integer*8 vertex_orig common /ivrtx/ vertex ((2**ldim),lelt) integer*8 vertex integer*8 ngv - nel = nelt - if (isTmsh.eq.0) nel = nelv - call set_vert(glo_num,ngv,nx,nel,vertex,.false.) + if (isRefine.eq.0) then + call set_vert(glo_num,ngv,nx,nel,vertex,.false.) + else ! refined meshes use original vertex + call set_vert(glo_num,ngv,nx,nel,vertex_orig,.false.) + endif nekf_set_vert = ngv diff --git a/src/nekInterface/nekInterfaceAdapter.cpp b/src/nekInterface/nekInterfaceAdapter.cpp index 49a156a1b..3a5e8a769 100644 --- a/src/nekInterface/nekInterfaceAdapter.cpp +++ b/src/nekInterface/nekInterfaceAdapter.cpp @@ -63,6 +63,7 @@ static void (*nek_uic_ptr)(int *); static void (*nek_end_ptr)(void); static void (*nek_restart_ptr)(char *, int *); static void (*nek_map_m_to_n_ptr)(double *a, int *na, double *b, int *nb, int *if3d, double *w, int *nw); +static void (*nek_map2reg_3di_e_ptr)(double *a, int *na, double *b, int *nb); static int (*nek_lglel_ptr)(int *); static void (*nek_bootstrap_ptr)(int *, char *, char *, char *, int, int, int); static void (*nek_setup_ptr)(int *, @@ -526,7 +527,7 @@ void getIC(void) } } -void xm1N(dfloat *_x, dfloat *_y, dfloat *_z, int N, dlong Nelements) +void xm1N(dfloat *_x, dfloat *_y, dfloat *_z, int N, dlong Nelements, int uniform) { const int Np = (N + 1) * (N + 1) * (N + 1); const int nxyz = nekData.nx1 * nekData.nx1 * nekData.nx1; @@ -544,14 +545,27 @@ void xm1N(dfloat *_x, dfloat *_y, dfloat *_z, int N, dlong Nelements) std::vector y(Np); std::vector z(Np); - for (dlong e = 0; e < Nelements; ++e) { - map_m_to_n(x.data(), N + 1, &nekData.xm1[e * nxyz], nekData.nx1); - map_m_to_n(y.data(), N + 1, &nekData.ym1[e * nxyz], nekData.nx1); - map_m_to_n(z.data(), N + 1, &nekData.zm1[e * nxyz], nekData.nx1); - for (int i = 0; i < Np; i++) { - _x[i + e * Np] = x[i]; - _y[i + e * Np] = y[i]; - _z[i + e * Np] = z[i]; + if (uniform) { + for (dlong e = 0; e < Nelements; ++e) { + nek_map2reg_3di_e_ptr(x.data(), N + 1, &nekData.xm1[e * nxyz], nekData.nx1); + nek_map2reg_3di_e_ptr(y.data(), N + 1, &nekData.ym1[e * nxyz], nekData.nx1); + nek_map2reg_3di_e_ptr(z.data(), N + 1, &nekData.zm1[e * nxyz], nekData.nx1); + for (int i = 0; i < Np; i++) { + _x[i + e * Np] = x[i]; + _y[i + e * Np] = y[i]; + _z[i + e * Np] = z[i]; + } + } + } else { + for (dlong e = 0; e < Nelements; ++e) { + map_m_to_n(x.data(), N + 1, &nekData.xm1[e * nxyz], nekData.nx1, uniform); + map_m_to_n(y.data(), N + 1, &nekData.ym1[e * nxyz], nekData.nx1, uniform); + map_m_to_n(z.data(), N + 1, &nekData.zm1[e * nxyz], nekData.nx1, uniform); + for (int i = 0; i < Np; i++) { + _x[i + e * Np] = x[i]; + _y[i + e * Np] = y[i]; + _z[i + e * Np] = z[i]; + } } } } @@ -701,6 +715,8 @@ void set_usr_handles(const char *session_in, int verbose) check_error(dlerror()); nek_map_m_to_n_ptr = (void (*)(double *, int *, double *, int *, int *, double *, int *))dlsym(handle, fname("map_m_to_n")); + nek_map2reg_3di_e_ptr = + (void (*)(double *, int *, double *, int *))dlsym(handle, fname("map2reg_3di_e")); check_error(dlerror()); nek_nbid_ptr = (int (*)(int *))dlsym(handle, fname("nekf_nbid")); check_error(dlerror()); diff --git a/src/nekInterface/nekInterfaceAdapter.hpp b/src/nekInterface/nekInterfaceAdapter.hpp index b39bc58f4..b9e972e14 100644 --- a/src/nekInterface/nekInterfaceAdapter.hpp +++ b/src/nekInterface/nekInterfaceAdapter.hpp @@ -122,7 +122,7 @@ void writeFld(const std::string& filename, bool uniform = false); void finalize(void); -void xm1N(dfloat *x, dfloat *y, dfloat *z, int Nq, dlong Nelements); +void xm1N(dfloat *x, dfloat *y, dfloat *z, int Nq, dlong Nelements, int uniform = 0); long long int localElementIdToGlobal(int _id); int lglel(int e); int setup(int numberActiveFields); @@ -138,7 +138,7 @@ int bcmap(int bid, int ifld); int globalElementIdToRank(long long id); int globalElementIdToLocal(long long id); -long long set_glo_num(int npts, int isTMesh); +long long set_glo_num(int npts, int nel, int isRefine); void bdfCoeff(dfloat *g0, dfloat *coeff, dfloat *dt, int order); void extCoeff(dfloat *coeff, dfloat *dt, int nAB, int nBDF); diff --git a/src/nekInterface/nekRefine.f b/src/nekInterface/nekRefine.f index e449d06f3..e79043c4a 100644 --- a/src/nekInterface/nekRefine.f +++ b/src/nekInterface/nekRefine.f @@ -7,6 +7,7 @@ subroutine usrdat2_oct(ncut) ! interface to oct-refine code include 'TOTAL' parameter(lxyz=lx1*ly1*lz1) + common /c_is1/ glo_num(lxyz*lelt) integer*8 glo_num @@ -277,6 +278,9 @@ subroutine h_refine(glo_num,ncut) $ , pc(lx1*lx1,mxnew),pt(lx1*lx1,mxnew) real pc,pt + common /ivrtx0/ vertex_orig ((2**ldim),lelt) + integer*8 vertex_orig + common /ivrtx/ vertex ((2**ldim),lelt) integer*8 vertex integer*8 ngv @@ -288,6 +292,18 @@ subroutine h_refine(glo_num,ncut) save isym2pre data isym2pre / 1 , 2 , 4 , 3 , 5 , 6 , 8 , 7 / + integer icalld + save icalld + data icalld / 0 / + + if (icalld.eq.0) then + do i = (2**ldim)*nelt + vertex_orig(i,1) = vertex(i) + enddo + if (nio.eq.0) write(6,12) nelt + 13 format('h-refine: backup vertex into vertex_orig',I12) + endif + nblk = ncut**ldim nnew = nblk - 1 lxyc = 2**ldim