Fix second order output in elmer format

doc/element_types.tex

@@ -0,0 +1,282 @@
+\documentclass[convert=pdf2svg]{standalone}
+% \documentclass{article}
+
+\usepackage[T1]{fontenc}
+\usepackage{lmodern}
+\renewcommand{\familydefault}{\sfdefault}
+\usepackage{tikz}
+\usepackage{tikz-3dplot}
+\usetikzlibrary{external}
+\tikzset{external/force remake}
+\tikzset{external/disable dependency files}
+\tikzset{external/aux in dpth={false}}
+% uncomment this to generate a figure for each cell type (and change documentclass to article)
+% \tikzexternalize
+
+\tikzstyle{vertex} = [circle,draw=black,fill=black,scale = 0.5]
+\tdplotsetmaincoords{70}{110}
+
+% cnode(tag,x,y,z,label,label_pos)
+\def\cnode(#1,#2,#3,#4,#5,#6){
+\node (#1) at (#2,#3,#4) [vertex,label=#6:$\mathsf{#5}$] {};
+}
+
+\pagestyle{empty}
+
+\begin{document}
+
+\begin{tabular}{cc}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+SEGMENT &
+SEGMENT3
+\\
+\tikzsetnextfilename{line}
+\begin{tikzpicture}[scale = 2]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,2,below right);
+\draw (n0) -- (n1);
+\end{tikzpicture}
+&
+\tikzsetnextfilename{line3}
+\begin{tikzpicture}[scale = 2]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,2,below right);
+\cnode(n2,1,0,0,3,below right);
+\draw (n0) -- (n2) -- (n1);
+\end{tikzpicture}
+\\[1 em]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+TRIG &
+TRIG6
+\\
+\tikzsetnextfilename{triangle}
+\begin{tikzpicture}[scale = 2]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,2,below right);
+\cnode(n2,0,2,0,3,right);
+\draw (n0) -- (n1) -- (n2) -- (n0);
+\end{tikzpicture}
+&
+\tikzsetnextfilename{triangle6}
+\begin{tikzpicture}[scale = 2]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,2,below right);
+\cnode(n2,0,2,0,3,right);
+\cnode(n3,1,0,0,6,below right);
+\cnode(n4,1,1,0,4,right);
+\cnode(n5,0,1,0,5,below right);
+\draw (n0) -- (n3) -- (n1) -- (n4) -- (n2) -- (n5) -- (n0);
+\end{tikzpicture}
+\\[1 em]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+QUAD &
+QUAD8
+\\
+\tikzsetnextfilename{quad}
+\begin{tikzpicture}[scale = 2]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,2,below right);
+\cnode(n2,2,2,0,3,below right);
+\cnode(n3,0,2,0,4,below right);
+\draw (n0) -- (n1) -- (n2) -- (n3) -- (n0);
+\end{tikzpicture}
+&
+\tikzsetnextfilename{quad8}
+\begin{tikzpicture}[scale = 2]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,2,below right);
+\cnode(n2,2,2,0,3,below right);
+\cnode(n3,0,2,0,4,below right);
+\cnode(n4,1,0,0,5,below right);
+\cnode(n5,2,1,0,8,below right);
+\cnode(n6,1,2,0,6,below right);
+\cnode(n7,0,1,0,7,below right);
+\draw (n0) -- (n4) -- (n1) -- (n5) -- (n2) -- (n6) -- (n3) -- (n7) -- (n0);
+\end{tikzpicture}
+\\[1 em]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+TET &
+TET10
+\\
+\tikzsetnextfilename{tetra}
+\begin{tikzpicture}[scale = 2, tdplot_main_coords]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,3,below right);
+\cnode(n2,0,2,0,2,below right);
+\cnode(n3,0,0,2,4,right);
+\draw (n0) -- (n1) -- (n2) -- (n0);
+\draw (n0) -- (n3);
+\draw (n1) -- (n3);
+\draw (n2) -- (n3);
+\end{tikzpicture}
+&
+\tikzsetnextfilename{tetra10} % VTK
+\begin{tikzpicture}[scale = 2, tdplot_main_coords]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,3,below right);
+\cnode(n2,0,2,0,2,below right);
+\cnode(n3,0,0,2,4,right);
+\cnode(n4,1,0,0,6,below right);
+\cnode(n5,1,1,0,8,below right);
+\cnode(n6,0,1,0,5,below right);
+\cnode(n7,0,0,1,7,below right);
+\cnode(n8,1,0,1,10,below right);
+\cnode(n9,0,1,1,9,right);
+\draw (n0) -- (n4) -- (n1) -- (n5) -- (n2) -- (n6) -- (n0);
+\draw (n0) -- (n7) -- (n3);
+\draw (n1) -- (n8) -- (n3);
+\draw (n2) -- (n9) -- (n3);
+\end{tikzpicture}
+\\[1 em]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+PYRAMID &
+PYRAMID13
+\\
+\tikzsetnextfilename{pyramid}
+\begin{tikzpicture}[scale = 2, tdplot_main_coords]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,4,below right);
+\cnode(n2,2,2,0,3,below right);
+\cnode(n3,0,2,0,2,below right);
+\cnode(n4,1,1,2,5,right);
+\draw (n0) -- (n1) -- (n2) -- (n3) -- (n0);
+\draw (n0) -- (n4);
+\draw (n1) -- (n4);
+\draw (n2) -- (n4);
+\draw (n3) -- (n4);
+\end{tikzpicture}
+&
+\tikzsetnextfilename{pyramid13} % VTK != gmsh
+\begin{tikzpicture}[scale = 2, tdplot_main_coords]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,4,below right);
+\cnode(n2,2,2,0,3,below right);
+\cnode(n3,0,2,0,2,below right);
+\cnode(n4,1,1,2,5,right);
+\cnode(n5,1,0,0,8,below right);
+\cnode(n6,2,1,0,7,below right);
+\cnode(n7,1,2,0,9,below right);
+\cnode(n8,0,1,0,6,below right);
+\cnode(n9,0.5,0.5,1,10,below right);
+\cnode(n10,1.5,0.5,1,13,below right);
+\cnode(n11,1.5,1.5,1,12,below right);
+\cnode(n12,0.5,1.5,1,11,right);
+\draw (n0) -- (n5) -- (n1) -- (n6) -- (n2) -- (n7) -- (n3) -- (n8) -- (n0);
+\draw (n0) -- (n9) -- (n4);
+\draw (n1) -- (n10) -- (n4);
+\draw (n2) -- (n11) -- (n4);
+\draw (n3) -- (n12) -- (n4);
+\end{tikzpicture}
+\\[1 em]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+PRISM &
+PRISM15
+\\
+\tikzsetnextfilename{wedge} % gmsh != VTK
+\begin{tikzpicture}[scale = 2, tdplot_main_coords]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,3,below right);
+\cnode(n2,0,2,0,2,below right);
+\cnode(n3,0,0,2,4,below right);
+\cnode(n4,2,0,2,6,below right);
+\cnode(n5,0,2,2,5,below right);
+\draw (n0) -- (n1) -- (n2) -- (n0);
+\draw (n3) -- (n4) -- (n5) -- (n3);
+\draw (n0) -- (n3);
+\draw (n1) -- (n4);
+\draw (n2) -- (n5);
+\end{tikzpicture}
+&
+\tikzsetnextfilename{wedge15} % VTK != gmsh
+\begin{tikzpicture}[scale = 2, tdplot_main_coords]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,3,below right);
+\cnode(n2,0,2,0,2,below right);
+\cnode(n3,0,0,2,4,below right);
+\cnode(n4,2,0,2,6,below right);
+\cnode(n5,0,2,2,5,below right);
+\cnode(n6,1,0,0,8,below right);
+\cnode(n7,1,1,0,9,below right);
+\cnode(n8,0,1,0,7,below right);
+\cnode(n9,1,0,2,14,below right);
+\cnode(n10,1,1,2,15,below right);
+\cnode(n11,0,1,2,13,below right);
+\cnode(n12,0,0,1,10,below right);
+\cnode(n13,2,0,1,12,below right);
+\cnode(n14,0,2,1,11,below right);
+\draw (n0) -- (n6) -- (n1) -- (n7) -- (n2) -- (n8) -- (n0);
+\draw (n3) -- (n9) -- (n4) -- (n10) -- (n5) -- (n11) -- (n3);
+\draw (n0) -- (n12) -- (n3);
+\draw (n1) -- (n13) -- (n4);
+\draw (n2) -- (n14) -- (n5);
+\end{tikzpicture}
+\\[1 em]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+HEX &
+HEX20
+\\
+\tikzsetnextfilename{hexahedron}
+\begin{tikzpicture}[scale = 2, tdplot_main_coords]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,4,below right);
+\cnode(n2,2,2,0,3,below right);
+\cnode(n3,0,2,0,2,below right);
+\cnode(n4,0,0,2,5,below right);
+\cnode(n5,2,0,2,8,below right);
+\cnode(n6,2,2,2,7,below right);
+\cnode(n7,0,2,2,6,below right);
+\draw (n0) -- (n1) -- (n2) -- (n3) -- (n0);
+\draw (n4) -- (n5) -- (n6) -- (n7) -- (n4);
+\draw (n0) -- (n4);
+\draw (n1) -- (n5);
+\draw (n2) -- (n6);
+\draw (n3) -- (n7);
+\end{tikzpicture}
+&
+\tikzsetnextfilename{hexahedron20} % VTK != gmsh
+\begin{tikzpicture}[scale = 2, tdplot_main_coords]
+\cnode(n0,0,0,0,1,below right);
+\cnode(n1,2,0,0,4,below right);
+\cnode(n2,2,2,0,3,below right);
+\cnode(n3,0,2,0,2,below right);
+\cnode(n4,0,0,2,5,below right);
+\cnode(n5,2,0,2,8,below right);
+\cnode(n6,2,2,2,7,below right);
+\cnode(n7,0,2,2,6,below right);
+\cnode(n8,1,0,0,11,below right);
+\cnode(n9,2,1,0,10,below right);
+\cnode(n10,1,2,0,12,below right);
+\cnode(n11,0,1,0,9,below right);
+\cnode(n12,1,0,2,15,below right);
+\cnode(n13,2,1,2,14,below right);
+\cnode(n14,1,2,2,16,below right);
+\cnode(n15,0,1,2,13,below right);
+\cnode(n16,0,0,1,17,below right);
+\cnode(n17,2,0,1,20,below right);
+\cnode(n18,2,2,1,19,below right);
+\cnode(n19,0,2,1,18,below right);
+\draw (n0) -- (n8) -- (n1) -- (n9) -- (n2) -- (n10) -- (n3) -- (n11) -- (n0);
+\draw (n4) -- (n12) -- (n5) -- (n13) -- (n6) -- (n14) -- (n7) -- (n15) -- (n4);
+\draw (n0) -- (n16) -- (n4);
+\draw (n1) -- (n17) -- (n5);
+\draw (n2) -- (n18) -- (n6);
+\draw (n3) -- (n19) -- (n7);
+\end{tikzpicture}
+
+
+
+\end{tabular}
+
+\end{document}

libsrc/interface/writeelmer.cpp

@@ -24,6 +24,35 @@ void WriteElmerFormat (const Mesh &mesh,
                        const string &filename)
 {
   cout << "write elmer mesh files" << endl;
+
+  std::map<ELEMENT_TYPE, int> tmap;
+  tmap[TRIG] = 303;
+  tmap[TRIG6] = 306;
+  tmap[QUAD] = 404;
+  tmap[QUAD8] = 408;
+  tmap[TET] = 504;
+  tmap[TET10] = 510;
+  tmap[PYRAMID] = 605;
+  tmap[PYRAMID13] = 613;
+  tmap[PRISM] = 706;
+  tmap[PRISM15] = 715;
+  tmap[HEX] = 808;
+  tmap[HEX20] = 820;
+
+  std::map<int, Array<int,int>> pmap;
+  pmap[TRIG]  = {1,2,3};
+  pmap[TRIG6] = {1,2,3, 6,4,5};
+  pmap[QUAD]  = {1,2,3,4};
+  pmap[QUAD8] = {1,2,3,4, 5,8,6,7};
+  pmap[TET]   = {1,2,3,4};
+  pmap[TET10] = {1,2,3,4, 5,8,6,7,9,10};
+  pmap[PYRAMID]={1,2,3,4,5};
+  pmap[PYRAMID13]= {1,2,3,4,5,6,7,8,9,10,11,12,13};
+  pmap[PRISM] = {1,2,3,4,5,6};
+  pmap[PRISM15] = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};
+  pmap[HEX]   = {1,2,3,4,5,6,7,8};
+  pmap[HEX20] = {1,2,3,4,5,8,6,7,8, 9,12,10,11, 17,20,19,18, 13,16,14,15};
+
   int np = mesh.GetNP();
   int ne = mesh.GetNE();
   int nse = mesh.GetNSE();
@@ -50,27 +79,67 @@ void WriteElmerFormat (const Mesh &mesh,
   ofstream outfile_e(str);
   sprintf( str, "%s/mesh.boundary", filename.c_str() );
   ofstream outfile_b(str);
+  sprintf( str, "%s/mesh.names", filename.c_str() );
+  ofstream outfile_names(str);
+
+  for( auto codim : IntRange(0, mesh.GetDimension()-1) )
+  {
+    auto & names = const_cast<Mesh&>(mesh).GetRegionNamesCD(codim);
+
+    for (auto i0 : Range(names) )
+    {
+      if(names[i0] == nullptr)
+        continue;
+      string name = *names[i0];
+      if(name == "" || name == "default")
+        continue;
+      outfile_names << "$" << name << "=" << i0+1 << "\n";
+    }
+  }
+
+  auto get3FacePoints = [](const Element2d & el)
+  {
+      INDEX_3 i3;
+      INDEX_4 i4;
+      auto eltype = el.GetType();
+      switch (eltype)
+      {
+          case TRIG:
+          case TRIG6:
+              i3 = {el[0], el[1], el[2]};
+              i3.Sort();
+              break;
+          case QUAD:
+          case QUAD8:
+              i4 = {el[0], el[1], el[2], el[3]};
+              i4.Sort();
+              i3 = {i4[0], i4[1], i4[2]};
+              break;
+          default:
+             throw Exception("Got invalid type (no face)");
+      }
+      return i3;
+  };
 
   // fill hashtable
 
+  // use lowest three point numbers of lowest-order face to index faces
   INDEX_3_HASHTABLE<int> face2volelement(ne);
 
-  for (i = 1; i <= ne; i++)
+  for (int i = 1; i <= ne; i++)
     {
       const Element & el = mesh.VolumeElement(i);
-      INDEX_3 i3;
+
-      int k, l;
+      // getface not working for second order elements -> reconstruct linear element here
-      for (j = 1; j <= 4; j++)   // loop over faces of tet
+      Element linear_el = el;
+      linear_el.SetNP(el.GetNV()); // GetNV returns 8 for HEX20 for instance
+
+      for (auto j : Range(1,el.GetNFaces()+1))
 	{
-	  l = 0;
+          Element2d face;
-	  for (k = 1; k <= 4; k++)
+          linear_el.GetFace(j, face);
-	    if (k != j)
+	  face2volelement.Set (get3FacePoints(face), i);
-	      {
+          cout << "set " << get3FacePoints(face) << "\tto " << i << endl;
-		l++;
-		i3.I(l) = el.PNum(k);
-	      }
-	  i3.Sort();
-	  face2volelement.Set (i3, i);
 	}
     }
 
@@ -78,10 +147,7 @@ void WriteElmerFormat (const Mesh &mesh,
 //  outfile.setf (ios::fixed, ios::floatfield);
 //  outfile.setf (ios::showpoint);
   
-  outfile_h << np << " " << ne << " " << nse << "\n";
+  std::map<ELEMENT_TYPE, size_t> elcount;
-  outfile_h << "2"     << "\n";
-  outfile_h << "303 " << nse << "\n";
-  outfile_h << "504 " << ne << "\n";
   
   for (i = 1; i <= np; i++)
     {
@@ -97,12 +163,15 @@ void WriteElmerFormat (const Mesh &mesh,
     {
       Element el = mesh.VolumeElement(i);
       if (inverttets) el.Invert();
-      sprintf( str, "5%02d", (int)el.GetNP() );
+      auto eltype = el.GetType();
-      outfile_e << i << " " << el.GetIndex() << " " << str <<  "  ";
+      elcount[eltype]++;
+      outfile_e << i << " " << el.GetIndex() << " " << tmap[eltype] <<  "  ";
+
+      auto & map = pmap[eltype];
       for (j = 1; j <= el.GetNP(); j++)
 	{
 	  outfile_e << " ";
-	  outfile_e << el.PNum(j);
+	  outfile_e << el.PNum(map[j-1]);
 	}
       outfile_e << "\n";
     }
@@ -111,23 +180,29 @@ void WriteElmerFormat (const Mesh &mesh,
     {
       Element2d el = mesh.SurfaceElement(i);
       if (invertsurf) el.Invert();
-      sprintf( str, "3%02d", (int)el.GetNP() );
+      auto eltype = el.GetType();
-      {
+      elcount[eltype]++;
-	  INDEX_3 i3;
+            
-	  for (j = 1; j <= 3; j++) i3.I(j) = el.PNum(j);
+      int elind = face2volelement.Get(get3FacePoints(el));
-	  i3.Sort();
+      cout << "get " << get3FacePoints(el) << "\t " << elind << endl;
 
-	  int elind = face2volelement.Get(i3);
       outfile_b << i << " " << mesh.GetFaceDescriptor(el.GetIndex()).BCProperty() << 
-         " " << elind << " 0 "  << str << "    ";
+         " " << elind << " 0 "  << tmap[eltype] << "    ";
-      }
+
+      auto & map = pmap[el.GetType()];
       for (j = 1; j <= el.GetNP(); j++)
 	{
 	  outfile_b << " ";
-	  outfile_b << el.PNum(j);
+	  outfile_b << el.PNum(map[j-1]);
 	}
       outfile_b << "\n";
     }
+
+  outfile_h << np << " " << ne << " " << nse << "\n";
+  outfile_h << "2"     << "\n";
+
+  for( auto & [eltype,count] : elcount )
+      outfile_h << tmap[eltype] << " " << count << "\n";
 }
 
 }