diff --git a/src/compil.date b/src/compil.date
index 9877436..1ca0c76 100644
--- a/src/compil.date
+++ b/src/compil.date
@@ -1 +1 @@
-#define COMPIL  "2017-07-26 10:59:58"
\ No newline at end of file
+#define COMPIL  "2021-05-20 07:58:25"
\ No newline at end of file
diff --git a/src/dlists.c b/src/dlists.c
index 2ba1d76..3ea3437 100644
--- a/src/dlists.c
+++ b/src/dlists.c
@@ -111,9 +111,12 @@ GLuint listTria(pScene sc,pMesh mesh) {
       else {
         is0 = is1 = is2 = 0;
         if ( mesh->extra->iv ) {
-          if ( pt->v[0] <= mesh->nvn )  is0 = mesh->extra->nv[pt->v[0]];
-          if ( pt->v[1] <= mesh->nvn )  is1 = mesh->extra->nv[pt->v[1]];
-          if ( pt->v[2] <= mesh->nvn )  is2 = mesh->extra->nv[pt->v[2]];
+          is0 = mesh->extra->nv[pt->v[0]];
+          assert ( is0 <= mesh->nvn );
+          is1 = mesh->extra->nv[pt->v[1]];
+          assert ( is1 <= mesh->nvn );
+          is2 = mesh->extra->nv[pt->v[2]];
+          assert ( is2 <= mesh->nvn );
         }
         if ( !is0 && pt->v[0] <= mesh->extra->it )  is0 = mesh->extra->nt[3*(k-1)+1];
         if ( !is1 && pt->v[1] <= mesh->extra->it )  is1 = mesh->extra->nt[3*(k-1)+2];
@@ -248,9 +251,12 @@ GLuint listTriaP2(pScene sc,pMesh mesh) {
       else {
         is0 = is1 = is2 = 0;
         if ( mesh->extra->iv ) {
-          if ( pt->v[0] <= mesh->nvn )  is0 = mesh->extra->nv[pt->v[0]];
-          if ( pt->v[1] <= mesh->nvn )  is1 = mesh->extra->nv[pt->v[1]];
-          if ( pt->v[2] <= mesh->nvn )  is2 = mesh->extra->nv[pt->v[2]];
+          is0 = mesh->extra->nv[pt->v[0]];
+          assert ( is0 <= mesh->nvn );
+          is1 = mesh->extra->nv[pt->v[1]];
+          assert ( is1 <= mesh->nvn );
+          is2 = mesh->extra->nv[pt->v[2]];
+          assert ( is2 <= mesh->nvn );
         }
         if ( !is0 && pt->v[0] <= mesh->extra->it )  is0 = mesh->extra->nt[3*(k-1)+1];
         if ( !is1 && pt->v[1] <= mesh->extra->it )  is1 = mesh->extra->nt[3*(k-1)+2];
@@ -418,14 +424,14 @@ GLuint listQuad(pScene sc,pMesh mesh) {
         }
         is0 = is1 = is2 = is3 = 0;
         if ( mesh->extra->iv ) {
-          if ( pq->v[0] <= mesh->nvn )
-            is0 = mesh->extra->nv[pq->v[0]];
-          if ( pq->v[1] <= mesh->nvn )
-            is1 = mesh->extra->nv[pq->v[1]];
-          if ( pq->v[2] <= mesh->nvn )
-            is2 = mesh->extra->nv[pq->v[2]];
-          if ( pq->v[3] <= mesh->nvn )
-            is3 = mesh->extra->nv[pq->v[3]];
+          is0 = mesh->extra->nv[pq->v[0]];
+          assert ( is0 <= mesh->nvn );
+          is1 = mesh->extra->nv[pq->v[1]];
+          assert ( is1 <= mesh->nvn );
+          is2 = mesh->extra->nv[pq->v[2]];
+          assert ( is2 <= mesh->nvn );
+          is3 = mesh->extra->nv[pq->v[3]];
+          assert ( is3 <= mesh->nvn );
         }
         if ( !is0 && pq->v[0] <= mesh->extra->iq )
           is0 = mesh->extra->nq[4*(k-1)+1];
diff --git a/src/extern.h b/src/extern.h
index cc4793b..5ea2c78 100644
--- a/src/extern.h
+++ b/src/extern.h
@@ -1,5 +1,6 @@
 #ifndef __MEDIT
   extern Canvas   cv;
+  extern float    scale;
   extern ubyte    ddebug,quiet,imprim,option,morphing,animate,saveimg,imgtype;
   extern int      animdep,animfin;
   extern char    *imgtyp[];
diff --git a/src/inout.c b/src/inout.c
index e269f7f..a59296f 100755
--- a/src/inout.c
+++ b/src/inout.c
@@ -3,6 +3,7 @@
 #include "extern.h"
 #include "eigenv.h"
 
+
 int loadMesh(pMesh mesh) {
   pPoint      ppt;
   pEdge       pr;
@@ -461,7 +462,7 @@ int saveMesh(pScene sc,pMesh mesh,char *fileout,ubyte clipon) {
   pTetra     ptt;
   pMaterial  pm;
   float      fp1,fp2,fp3;
-  int        outm,i,k,m,ver,np,nt,nq,ref;
+  int        outm,i,k,m,ver,np,nt,nq,ntt,ref;
 
   ver = GmfFloat;
   if ( !(outm = GmfOpenMesh(fileout,GmfWrite,ver,mesh->dim)) ) {
@@ -517,20 +518,22 @@ int saveMesh(pScene sc,pMesh mesh,char *fileout,ubyte clipon) {
     }
     if ( i == 3 )  nt++;
   }
-  GmfSetKwd(outm,GmfTriangles,nt);
-  for (k=1; k<=mesh->nt; k++) {
-    pt = &mesh->tria[k];
-    if ( !pt->v[0] )  continue;
-    m  = matRef(sc,pt->ref);
-    pm = &sc->material[m];
-    if ( pm->flag )  continue;
-    for (i=0; i<3; i++) {
-      if ( !mesh->point[pt->v[i]].tmp )  break;
+  if ( nt > 0 ) {
+    GmfSetKwd(outm,GmfTriangles,nt);
+    for (k=1; k<=mesh->nt; k++) {
+      pt = &mesh->tria[k];
+      if ( !pt->v[0] )  continue;
+      m  = matRef(sc,pt->ref);
+      pm = &sc->material[m];
+      if ( pm->flag )  continue;
+      for (i=0; i<3; i++) {
+        if ( !mesh->point[pt->v[i]].tmp )  break;
+      }
+      if ( i < 3 )  continue;
+      ref = pt->ref;
+      GmfSetLin(outm,GmfTriangles,mesh->point[pt->v[0]].tmp,mesh->point[pt->v[1]].tmp,
+                mesh->point[pt->v[2]].tmp,ref);
     }
-    if ( i < 3 )  continue;
-    ref = pt->ref;
-    GmfSetLin(outm,GmfTriangles,mesh->point[pt->v[0]].tmp,mesh->point[pt->v[1]].tmp,
-              mesh->point[pt->v[2]].tmp,ref);
   }
 
   /* write quads */
@@ -545,19 +548,21 @@ int saveMesh(pScene sc,pMesh mesh,char *fileout,ubyte clipon) {
       if ( !mesh->point[pq->v[i]].tmp )  break;
     if ( i == 4 )  nq++;
   }
-  GmfSetKwd(outm,GmfQuadrilaterals,nq);
-  for (k=1; k<=mesh->nq; k++) {
-    pq = &mesh->quad[k];
-    if ( !pq->v[0] )  continue;
-    m  = matRef(sc,pq->ref);
-    pm = &sc->material[m];
-    if ( pm->flag )  continue;
-    for (i=0; i<4; i++)
-      if ( !mesh->point[pq->v[i]].tmp )  break;
-    if ( i < 4 )  continue;
-    ref = pq->ref;
-    GmfSetLin(outm,GmfQuadrilaterals,mesh->point[pq->v[0]].tmp,mesh->point[pq->v[1]].tmp,
-              mesh->point[pq->v[2]].tmp,mesh->point[pq->v[3]].tmp,ref);
+  if ( nq > 0 ) {
+    GmfSetKwd(outm,GmfQuadrilaterals,nq);
+    for (k=1; k<=mesh->nq; k++) {
+      pq = &mesh->quad[k];
+      if ( !pq->v[0] )  continue;
+      m  = matRef(sc,pq->ref);
+      pm = &sc->material[m];
+      if ( pm->flag )  continue;
+      for (i=0; i<4; i++)
+        if ( !mesh->point[pq->v[i]].tmp )  break;
+      if ( i < 4 )  continue;
+      ref = pq->ref;
+      GmfSetLin(outm,GmfQuadrilaterals,mesh->point[pq->v[0]].tmp,mesh->point[pq->v[1]].tmp,
+                mesh->point[pq->v[2]].tmp,mesh->point[pq->v[3]].tmp,ref);
+    }
   }
 
   /* write edges */
@@ -578,7 +583,7 @@ int saveMesh(pScene sc,pMesh mesh,char *fileout,ubyte clipon) {
   }
 
   /* write tetrahedra */
-  GmfSetKwd(outm,GmfTetrahedra,mesh->ntet);
+  ntt = 0;
   for (k=1; k<=mesh->ntet; k++) {
     ptt = &mesh->tetra[k];
     if ( !ptt->v[0] )  continue;
@@ -587,19 +592,31 @@ int saveMesh(pScene sc,pMesh mesh,char *fileout,ubyte clipon) {
     if ( pm->flag )  continue;
     for (i=0; i<4; i++)
       if ( !mesh->point[ptt->v[i]].tmp )  break;
-    if ( i < 4 )  continue;
-    ref = ptt->ref;
-    GmfSetLin(outm,GmfTetrahedra,mesh->point[ptt->v[0]].tmp,mesh->point[ptt->v[1]].tmp,
-              mesh->point[ptt->v[2]].tmp,mesh->point[ptt->v[3]].tmp,ref);
+    if ( i == 4 )  ntt++;
+  }
+  if ( ntt > 0 ) {
+    GmfSetKwd(outm,GmfTetrahedra,ntt);
+    for (k=1; k<=mesh->ntet; k++) {
+      ptt = &mesh->tetra[k];
+      if ( !ptt->v[0] )  continue;
+      m  = matRef(sc,ptt->ref);
+      pm = &sc->material[m];
+      if ( pm->flag )  continue;
+      for (i=0; i<4; i++)
+        if ( !mesh->point[ptt->v[i]].tmp )  break;
+      if ( i < 4 )  continue;
+      ref = ptt->ref;
+      GmfSetLin(outm,GmfTetrahedra,mesh->point[ptt->v[0]].tmp,mesh->point[ptt->v[1]].tmp,
+                mesh->point[ptt->v[2]].tmp,mesh->point[ptt->v[3]].tmp,ref);
+    }
   }
 
   /* write hexahedra */
-
   if ( !quiet ) {
     fprintf(stdout,"     TOTAL NUMBER OF VERTICES   %8d\n",np);
     fprintf(stdout,"     TOTAL NUMBER OF TRIANGLES  %8d\n",nt);
     fprintf(stdout,"     TOTAL NUMBER OF QUADS      %8d\n",nq);
-    fprintf(stdout,"     TOTAL NUMBER OF TETRA      %8d\n",mesh->ntet);
+    fprintf(stdout,"     TOTAL NUMBER OF TETRA      %8d\n",ntt);
   }
 
   GmfCloseMesh(outm);
diff --git a/src/material.c b/src/material.c
index 1294bc1..bd264d0 100644
--- a/src/material.c
+++ b/src/material.c
@@ -182,7 +182,7 @@ int matRef(pScene sc,int ref) {
 /* reshape material window */
 void matReshape(int width,int height) {
   glutSetWindow(matwin);
-  glViewport(0,0,Width,Height);
+  glViewport(0,0,scale*Width,scale*Height);
   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
   gluOrtho2D(0,Width,Height,0);
@@ -259,7 +259,7 @@ void matsubReshape(int x,int y) {
   GLfloat    sunpos[4] = {5.,2.,10.0,0.};
 
   glutSetWindow(m_subwin);
-  glViewport(0,0,200,200);
+  glViewport(0,0,scale*200,scale*200);
   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
   gluPerspective(50.0,1.0,0.01f,100.0f);
diff --git a/src/medit.c b/src/medit.c
index c70d1f9..481116e 100644
--- a/src/medit.c
+++ b/src/medit.c
@@ -19,6 +19,7 @@ ubyte     quiet,option,morphing,stereoMode;
 int       menu,amenu,fmenu,femenu,vmenu,mmenu,smenu;
 int       clmenu,cmenu,vwmenu,txmenu,trmenu;
 int       animdep,animfin;
+float     scale=1; // for retina h-dpi screen
 
 
 static void excfun(int sigid) {
diff --git a/src/mlists.c b/src/mlists.c
index ab82e4d..71b4337 100644
--- a/src/mlists.c
+++ b/src/mlists.c
@@ -291,12 +291,12 @@ GLuint listTriaMap(pScene sc,pMesh mesh) {
 
         is0 = is1 = is2 = 0;
         if ( mesh->extra->iv ) {
-          if ( pt->v[0] <= mesh->nvn )
-            is0 = mesh->extra->nv[pt->v[0]];
-          if ( pt->v[1] <= mesh->nvn )
-            is1 = mesh->extra->nv[pt->v[1]];
-          if ( pt->v[2] <= mesh->nvn )
-            is2 = mesh->extra->nv[pt->v[2]];
+          is0 = mesh->extra->nv[pt->v[0]];
+          assert ( is0 <= mesh->nvn );
+          is1 = mesh->extra->nv[pt->v[1]];
+          assert ( is1 <= mesh->nvn );
+          is2 = mesh->extra->nv[pt->v[2]];
+          assert ( is2 <= mesh->nvn );
         }
         if ( !is0 && pt->v[0] <= mesh->extra->it )
           is0 = mesh->extra->nt[3*(k-1)+1];
@@ -481,14 +481,14 @@ GLuint listQuadMap(pScene sc,pMesh mesh) {
       else {
         is0 = is1 = is2 = is3 = 0;
         if ( mesh->extra->iv ) {
-          if ( pq->v[0] <= mesh->nvn )
-            is0 = mesh->extra->nv[pq->v[0]];
-          if ( pq->v[1] <= mesh->nvn )
-            is1 = mesh->extra->nv[pq->v[1]];
-          if ( pq->v[2] <= mesh->nvn )
-            is2 = mesh->extra->nv[pq->v[2]];
-          if ( pq->v[3] <= mesh->nvn )
-            is3 = mesh->extra->nv[pq->v[3]];
+          is0 = mesh->extra->nv[pq->v[0]];
+          assert ( is0 <= mesh->nvn );
+          is1 = mesh->extra->nv[pq->v[1]];
+          assert ( is1 <= mesh->nvn );
+          is2 = mesh->extra->nv[pq->v[2]];
+          assert ( is2 <= mesh->nvn );
+          is3 = mesh->extra->nv[pq->v[3]];
+          assert ( is3 <= mesh->nvn );
         }
         if ( !is0 && pq->v[0] <= mesh->extra->iq )
           is0 = mesh->extra->nq[4*(k-1)+1];
diff --git a/src/parsar.c b/src/parsar.c
index 6a7f17c..fd5401d 100644
--- a/src/parsar.c
+++ b/src/parsar.c
@@ -25,6 +25,7 @@ void usage() {
   fprintf(stdout,"\n** Graphic options\n");
   fprintf(stdout,"-a start stop\t Animation sequence\n");
   fprintf(stdout,"-m  Morphing\n");
+  fprintf(stdout,"-r  Retina h-dpi mode\n");
   fprintf(stdout,"-xv width height\t Visual Schnauzer\n");
   fprintf(stdout,"-stereo\t Stereo mode\n");     
   fprintf(stdout,"\n");
@@ -58,6 +59,8 @@ int parsar(int argc,char *argv[]) {
       option = MORPHING;
     else if ( !strcmp(argv[i],"-iso") )
       option = ISOSURF;
+    else if ( !strcmp(argv[i],"-r") )
+      scale = 2;
     else if ( !strcmp(argv[i],"-stereo") )
       stereoMode = LEFT+RIGHT;
     else if ( !strcmp(argv[i],"-s") )
diff --git a/src/picking.c b/src/picking.c
index 6f8778d..5fd3961 100644
--- a/src/picking.c
+++ b/src/picking.c
@@ -314,6 +314,11 @@ static void drawPoint(pScene sc,pMesh mesh,int k) {
 }
 
 
+static double norm(float *m) {
+  return(sqrt(m[0]*m[0] + m[1]*m[1] + m[2]*m[2]));
+}
+
+
 static void infoData(pScene sc,pMesh mesh,int k,int typel) {
   pSolution  ps;
 
@@ -322,9 +327,10 @@ static void infoData(pScene sc,pMesh mesh,int k,int typel) {
   if ( mesh->nfield == 1 )
     fprintf(stdout,"  Data (scalar): %f\n",ps->bb);
   else if ( mesh->nfield == mesh->dim ) {
-    fprintf(stdout,"  Data (vector): %f %f",ps->m[0],ps->m[1]);
-    if (mesh->dim == 3 )  fprintf(stdout," %f",ps->m[2]); 
-    fprintf(stdout,"\n");
+    if ( mesh->dim == 2 )
+      fprintf(stdout,"  Data (vector): %f  %f  |v|= \n",ps->m[0],ps->m[1],norm(ps->m));
+    else 
+      fprintf(stdout,"  Data (vector): %f %f %f   |v|= %f\n",ps->m[0],ps->m[1],ps->m[2],norm(ps->m)); 
   }
   else if ( mesh->dim == 2 && mesh->nfield == 3 ) {
     fprintf(stdout,"  Data (tensor): %f %f %f\n",
@@ -818,7 +824,7 @@ GLuint pickingScene(pScene sc,int x,int y,int ident) {
 
   /* process color */
   glReadBuffer(GL_BACK_LEFT);
-  glReadPixels(x,viewport[3]-y,1,1,
+  glReadPixels(scale*x,viewport[3]-scale*y,1,1,
                GL_RGBA,GL_UNSIGNED_BYTE,(void *)pixel);
   glDrawBuffer(GL_BACK_LEFT);
 
@@ -859,7 +865,7 @@ GLuint pickingScene(pScene sc,int x,int y,int ident) {
   /* peculiar case: vertex */
   if ( refitem > 0 && ident == LPoint ) {
     if ( mesh->ne ) {
-      refitem = closestPoint(sc,mesh,x,y,refitem,reftype);
+      refitem = closestPoint(sc,mesh,scale*x,scale*y,refitem,reftype);
       reftype = LPoint;
     }
   }
diff --git a/src/scene.c b/src/scene.c
index 13978ff..20b55a2 100644
--- a/src/scene.c
+++ b/src/scene.c
@@ -101,7 +101,7 @@ void reshapeScene(int width,int height) {
   sc->par.xs = width;
   sc->par.ys = height;
 
-  glViewport(0,0,width,height);
+  glViewport(0,0,scale*width,scale*height);
   farclip(GL_TRUE);
 }
 
diff --git a/src/scissor.c b/src/scissor.c
index 2b16204..2e0fd15 100644
--- a/src/scissor.c
+++ b/src/scissor.c
@@ -18,7 +18,7 @@ void scissorScene() {
   p    = sc->persp; 
 
   /* subdivide main window */
-  glViewport(0,0,sc->par.xs,sc->par.ys);
+  glViewport(0,0,scale*sc->par.xs,scale*sc->par.ys);
   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
   gluOrtho2D(0,sc->par.xs,0,sc->par.ys);
@@ -50,8 +50,8 @@ void scissorScene() {
   glEnable(GL_SCISSOR_TEST);
 
   /* draw top right : normal */
-  glViewport(width,height,width,height);
-  glScissor(width,height,width,height);
+  glViewport(scale*width,scale*height,scale*width,scale*height);
+  glScissor(scale*width,scale*height,scale*width,scale*height);
   farclip(1);
 
   glMatrixMode(GL_MODELVIEW);
@@ -64,8 +64,8 @@ void scissorScene() {
   if ( sc->type & S_DECO )  redrawStatusBar(sc);
 
   /* draw top left : top view */
-  glViewport(0,height,width,height);
-  glScissor(0,height,width,height);
+  glViewport(0,scale*height,scale*width,scale*height);
+  glScissor(0,scale*height,scale*width,scale*height);
   farclip(1);
 
   glMatrixMode(GL_MODELVIEW);
@@ -78,8 +78,8 @@ void scissorScene() {
   drawModel(sc);
 
   /* draw bottom left : front */
-  glViewport(0,0,width,height);
-  glScissor(0,0,width,height);
+  glViewport(0,0,scale*width,scale*height);
+  glScissor(0,0,scale*width,scale*height);
   farclip(1);
 
   glMatrixMode(GL_MODELVIEW);
@@ -92,8 +92,8 @@ void scissorScene() {
   drawModel(sc);
 
   /* draw bottom right : right */
-  glViewport(width,0,width,height);
-  glScissor(width,0,width,height);
+  glViewport(scale*width,0,scale*width,scale*height);
+  glScissor(scale*width,0,scale*width,scale*height);
   farclip(1);
 
   glMatrixMode(GL_MODELVIEW);
@@ -107,7 +107,7 @@ void scissorScene() {
 
   glutSwapBuffers();
   glDisable(GL_SCISSOR_TEST);
-  glViewport(0,0,sc->par.xs,sc->par.ys);
+  glViewport(0,0,scale*sc->par.xs,scale*sc->par.ys);
  
   if ( saveimg ) keyFile('H',0,0);
 }