Cleanup Abaqus export, implement 1d and 2d elements

libsrc/interface/writeabaqus.cpp

@@ -15,15 +15,103 @@
 namespace netgen
 {
 
+using std::vector;
 
+struct AbaqusElementType
+{
+  const char * name;
+  const vector<int> permutation;
 
+  AbaqusElementType(const char * name, const vector<int> & permutation)
+    : name(name), permutation(permutation)
+  {}
+};
+
+static inline const AbaqusElementType & GetAbaqusType(int dim, int num_nodes)
+{
+  // maps num_nodes to AbaqusElementType for each dimension
+  typedef std::map<int, AbaqusElementType> AbaqusElementTypes;
+  static const std::map<int, AbaqusElementType> abaqus_eltypes[3] =
+  {
+    // 1D
+    AbaqusElementTypes{
+      {2, AbaqusElementType{"T2D2", vector{0,1}}},
+    },
+    // 2D
+    AbaqusElementTypes{
+      {3, AbaqusElementType{"CPS3", vector{0,1,2}}},
+    },
+    // 3D
+    AbaqusElementTypes{
+      {4, AbaqusElementType{"C3D4", vector{0,1,3,2}}},
+      {10, AbaqusElementType{"C3D10", vector{0,1,3,2,4,8,6,5,7,9}}},
+    }
+  };
+
+  const auto & eltypes = abaqus_eltypes[dim-1];
+  if (eltypes.count(num_nodes) > 0)
+    return eltypes.at(num_nodes);
+  else
+    throw Exception("unsupported " + ToString(dim)+"d Element type with " + ToString(num_nodes) + " nodes");
+}
+
+static void WritePoints ( const Mesh & mesh, ostream & out )
+{
+  out << "*Node" << endl;
+  for(auto pi : mesh.Points().Range() )
+  {
+    out << pi+1-PointIndex::BASE << ", ";
+    auto p = mesh[pi];
+    out << p[0] << ", " << p[1] << ", " << p[2] << '\n';
+  }
+}
+
+template <typename ElIndex>
+static void WriteElement(ostream & out, const Mesh& mesh, ElIndex ei, const vector<int> & permutation, int & el_counter)
+{
+  el_counter++;
+  auto el = mesh[ei];
+  out << el_counter;
+  for(auto i : Range(el.PNums()))
+    out << ", " << el[permutation[i]]+1-PointIndex::BASE;
+  out << '\n';
+}
+
+template <typename ElIndex, typename Elements>
+static void WriteElements ( ostream & out, const Mesh & mesh, int dim, const Elements & el_range, int & el_counter)
+{
+  // map index, num_nodes to elements
+  std::map<std::tuple<int, int>, Array<ElIndex>> elset_map;
+
+  for(auto ei : el_range)
+    {
+      const auto & el = mesh[ei];
+      int index = 0;
+      if constexpr(std::is_same_v<ElIndex,SegmentIndex>)
+        index = el.edgenr;
+      else
+        index = el.GetIndex();
+      elset_map[{index, el.GetNP()}].Append(ei);
+    }
+
+  for(auto & [key, elems] : elset_map)
+    {
+      auto [index, num_nodes] = key;
+      auto name = mesh.GetRegionName(elems[0]);
+      if (name == "") name = "default";
+      PrintMessage (5, index, ": ", name);
+      const auto & eltype = GetAbaqusType(dim, num_nodes) ;
+      out << "*Element, type=" << eltype.name << ", ELSET=" << name << endl;
+      for(auto ei : elems)
+        WriteElement(out, mesh, ei, eltype.permutation, el_counter);
+    }
+}
 
 void WriteAbaqusFormat (const Mesh & mesh,
 			const filesystem::path & filename)
 
 {
-      
+  PrintMessage (1, "Write Abaqus Mesh");
-  cout << "\nWrite Abaqus Volume Mesh" << endl;
 
   ofstream outfile (filename);
 
@@ -32,96 +120,17 @@ void WriteAbaqusFormat (const Mesh & mesh,
 
   outfile.precision(8);
 
-  outfile << "*Node" << endl;
+  int element_counter = 0;
-
+  WritePoints(mesh, outfile);
-  int np = mesh.GetNP();
+  if(mesh.GetDimension() < 3)
-  int ne = mesh.GetNE();
+    WriteElements<SegmentIndex>(outfile, mesh, 1, mesh.LineSegments().Range(), element_counter);
-  int i, j, k;
+  WriteElements<SurfaceElementIndex>(outfile, mesh, 2, mesh.SurfaceElements().Range(), element_counter);
-
+  WriteElements<ElementIndex>(outfile, mesh, 3, mesh.VolumeElements().Range(), element_counter);
-  for (i = 1; i <= np; i++)
-    {
-      outfile << i << ", ";
-      outfile << mesh.Point(i)(0) << ", ";
-      outfile << mesh.Point(i)(1) << ", ";
-      outfile << mesh.Point(i)(2) << "\n";
-    }
-
-  int elemcnt = 0; //element counter
-  int finished = 0;
-  int indcnt = 1; //index counter
-
-  while (!finished)
-    {
-      int actcnt = 0;
-      const Element & el1 = mesh.VolumeElement(1);
-      int non = el1.GetNP();
-      if (non == 4)
-	{
-	  outfile << "*Element, type=C3D4, ELSET=PART" << indcnt << endl;
-	} 
-      else if (non == 10)
-	{
-	  outfile << "*Element, type=C3D10, ELSET=PART" << indcnt << endl;
-	} 
-      else
-	{
-	  cout << "unsupported Element type!!!" << endl;	  
-	}
-
-      for (i = 1; i <= ne; i++)
-	{
-	  const Element & el = mesh.VolumeElement(i);
-	      
-	  if (el.GetIndex() == indcnt)
-	    {
-	      actcnt++;
-	      if (el.GetNP() != non) 
-		{
-		  cout << "different element-types in a subdomain are not possible!!!" << endl;
-		  continue;
-		}
-		  
-	      elemcnt++;
-	      outfile << elemcnt << ", ";
-	      if (non == 4)
-		{
-		  outfile << el.PNum(1) << ", ";
-		  outfile << el.PNum(2) << ", ";
-		  outfile << el.PNum(4) << ", ";
-		  outfile << el.PNum(3) << "\n";
-		}
-	      else if (non == 10)
-		{
-		  outfile << el.PNum(1) << ", ";
-		  outfile << el.PNum(2) << ", ";
-		  outfile << el.PNum(4) << ", ";
-		  outfile << el.PNum(3) << ", ";
-		  outfile << el.PNum(5) << ", ";
-		  outfile << el.PNum(9) << ", ";
-		  outfile << el.PNum(7) << ", " << "\n";
-		  outfile << el.PNum(6) << ", ";
-		  outfile << el.PNum(8) << ", ";
-		  outfile << el.PNum(10) << "\n";
-		}
-	      else
-		{
-		  cout << "unsupported Element type!!!" << endl;
-		  for (j = 1; j <= el.GetNP(); j++)
-		    {
-		      outfile << el.PNum(j);
-		      if (j != el.GetNP()) outfile << ", ";
-		    }
-		  outfile << "\n";
-		}
-	    }
-	}	  
-      indcnt++;
-      if (elemcnt == ne) {finished = 1; cout << "all elements found by Index!" << endl;}
-      if (actcnt == 0) {finished = 1;}
-    }
 
+  // Write identifications (untested!)
   if (mesh.GetIdentifications().GetMaxNr())
     {
+      const auto np = mesh.GetNP();
       // periodic identification, implementation for
       // Helmut J. Boehm, TU Vienna
 	  
@@ -138,27 +147,27 @@ void WriteAbaqusFormat (const Mesh & mesh,
       NgArray<INDEX_2> pairs;
       NgBitArray master(np), help(np);
       master.Set();
-      for (i = 1; i <= 3; i++)
+      for (int i = 1; i <= 3; i++)
 	{
 	  mesh.GetIdentifications().GetPairs (i, pairs);
 	  help.Clear();
-	  for (j = 1; j <= pairs.Size(); j++)
+	  for (int j = 1; j <= pairs.Size(); j++)
 	    {
 	      help.Set (pairs.Get(j).I1());
 	    }
 	  master.And (help);
 	}
-      for (i = 1; i <= np; i++)
+      for (int i = 1; i <= np; i++)
 	if (master.Test(i))
 	  masternode = i;
 
       cout << "masternode = " << masternode << " = "
 	   << mesh.Point(masternode) << endl;
       NgArray<int> minions(3);
-      for (i = 1; i <= 3; i++)
+      for (int i = 1; i <= 3; i++)
 	{
 	  mesh.GetIdentifications().GetPairs (i, pairs);
-	  for (j = 1; j <= pairs.Size(); j++)
+	  for (int j = 1; j <= pairs.Size(); j++)
 	    {
 	      if (pairs.Get(j).I1() == masternode)
 		minions.Elem(i) = pairs.Get(j).I2();
@@ -179,12 +188,12 @@ void WriteAbaqusFormat (const Mesh & mesh,
 	      << "**POINT_fixed\n"
 	      << "**\n"
 	      << "*BOUNDARY, OP=NEW\n";
-      for (j = 1; j <= 3; j++)
+      for (int j = 1; j <= 3; j++)
 	outfile << masternode << ", " << j << ",,    0.\n";
 
       outfile << "**\n"
 	      << "*BOUNDARY, OP=NEW\n";
-      for (j = 1; j <= 3; j++)
+      for (int j = 1; j <= 3; j++)
 	{
 	  Vec3d v(mesh.Point(masternode), mesh.Point(minions.Get(j)));
 	  double vlen = v.Length();
@@ -203,18 +212,18 @@ void WriteAbaqusFormat (const Mesh & mesh,
 
       NgBitArray eliminated(np);
       eliminated.Clear();
-      for (i = 1; i <= mesh.GetIdentifications().GetMaxNr(); i++)
+      for (int i = 1; i <= mesh.GetIdentifications().GetMaxNr(); i++)
 	{
 	  mesh.GetIdentifications().GetPairs (i, pairs);
 	  if (!pairs.Size())
 	    continue;
 	      
-	  for (j = 1; j <= pairs.Size(); j++)
+	  for (int j = 1; j <= pairs.Size(); j++)
 	    if (pairs.Get(j).I1() != masternode && 
 		!eliminated.Test(pairs.Get(j).I2()))
 	      {
 		eliminated.Set (pairs.Get(j).I2());
-		for (k = 1; k <= 3; k++)
+		for (int k = 1; k <= 3; k++)
 		  {
 		    mpc << "4" << "\n";
 		    mpc << pairs.Get(j).I2() << "," << k << ", -1.0, ";
@@ -227,7 +236,7 @@ void WriteAbaqusFormat (const Mesh & mesh,
     }
 
 
-  cout << "done" << endl;
+  PrintMessage(1, "done");
 }
 
 static RegisterUserFormat reg_abaqus ("Abaqus Format", {".mesh"}, nullopt, WriteAbaqusFormat);