(drafted) Add MG layer for refined mesh

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<mesh_t*, mesh_t*> 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

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++) {

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];
 }

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,18 +51,22 @@ 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);
 
   //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);

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) {

src/mesh/meshPhysicalNodesHex3D.cpp

@@ -33,7 +33,7 @@ void meshPhysicalNodesHex3D(mesh_t *mesh)
   std::vector<dfloat> ym(mesh->Nlocal);
   std::vector<dfloat> 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<dfloat>(mesh->Nlocal);

src/mesh/meshSetup.cpp

@@ -333,10 +333,12 @@ std::pair<mesh_t*, mesh_t*> 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);

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
 

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<double> y(Np);
   std::vector<double> 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());

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);

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