From 54b5567173e6413ad7d0d19ccaa83fd06a7b1d4a Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Tue, 16 Mar 2021 20:27:07 +0100
Subject: [PATCH] Simplex: DataOut::attach_triangulation()

---
 .../numerics/data_out_dof_data.templates.h    |  55 +--
 source/numerics/data_out.cc                   |  12 +-
 tests/simplex/data_out_write_vtk_03.cc        |  72 ++++
 ...rite_vtk_03.with_simplex_support=on.output | 353 ++++++++++++++++++
 4 files changed, 467 insertions(+), 25 deletions(-)
 create mode 100644 tests/simplex/data_out_write_vtk_03.cc
 create mode 100644 tests/simplex/data_out_write_vtk_03.with_simplex_support=on.output

diff --git a/include/deal.II/numerics/data_out_dof_data.templates.h b/include/deal.II/numerics/data_out_dof_data.templates.h
index 4eb6728d67..60d7982c1d 100644
--- a/include/deal.II/numerics/data_out_dof_data.templates.h
+++ b/include/deal.II/numerics/data_out_dof_data.templates.h
@@ -30,6 +30,7 @@
 #include <deal.II/dofs/dof_handler.h>
 
 #include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_pyramid_p.h>
 #include <deal.II/fe/fe_simplex_p.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/fe_wedge_p.h>
@@ -232,14 +233,9 @@ namespace internal
               {
                 const auto reference_cell = (*fe)[i].reference_cell();
 
-                if ((reference_cell == dealii::ReferenceCells::Vertex) ||
-                    (reference_cell == dealii::ReferenceCells::Line) ||
-                    (reference_cell == dealii::ReferenceCells::Quadrilateral) ||
-                    (reference_cell == dealii::ReferenceCells::Hexahedron))
+                if (reference_cell.is_hyper_cube())
                   needs_hypercube_setup |= true;
-                else if ((reference_cell == dealii::ReferenceCells::Triangle) ||
-                         (reference_cell ==
-                          dealii::ReferenceCells::Tetrahedron))
+                else if (reference_cell.is_simplex())
                   needs_simplex_setup |= true;
                 else if (reference_cell == dealii::ReferenceCells::Wedge)
                   needs_wedge_setup |= true;
@@ -322,17 +318,9 @@ namespace internal
                       const auto reference_cell =
                         (*finite_elements[i])[j].reference_cell();
 
-                      if ((reference_cell == dealii::ReferenceCells::Vertex) ||
-                          (reference_cell == dealii::ReferenceCells::Line) ||
-                          (reference_cell ==
-                           dealii::ReferenceCells::Quadrilateral) ||
-                          (reference_cell ==
-                           dealii::ReferenceCells::Hexahedron))
+                      if (reference_cell.is_hyper_cube())
                         quadrature.push_back(*quadrature_hypercube);
-                      else if ((reference_cell ==
-                                dealii::ReferenceCells::Triangle) ||
-                               (reference_cell ==
-                                dealii::ReferenceCells::Tetrahedron))
+                      else if (reference_cell.is_simplex())
                         quadrature.push_back(*quadrature_simplex);
                       else if (reference_cell == dealii::ReferenceCells::Wedge)
                         quadrature.push_back(*quadrature_wedge);
@@ -2213,10 +2201,35 @@ DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::get_fes() const
     }
   if (this->dof_data.empty())
     {
-      finite_elements.resize(1);
-      finite_elements[0] =
-        std::make_shared<dealii::hp::FECollection<dhdim, dhspacedim>>(
-          FE_DGQ<dhdim, dhspacedim>(0));
+      Assert(triangulation != nullptr, ExcNotImplemented());
+
+      const auto &reference_cells = triangulation->get_reference_cells();
+      finite_elements.reserve(reference_cells.size());
+
+      // TODO: below we select linear, discontinuous elements for simplex,
+      // wedge, and pyramid; we would like to select constant functions, but
+      // that is not implemented yet
+      for (const auto &reference_cell : reference_cells)
+        {
+          if (reference_cell.is_hyper_cube())
+            finite_elements.emplace_back(
+              std::make_shared<dealii::hp::FECollection<dhdim, dhspacedim>>(
+                FE_DGQ<dhdim, dhspacedim>(0)));
+          else if (reference_cell.is_simplex())
+            finite_elements.emplace_back(
+              std::make_shared<dealii::hp::FECollection<dhdim, dhspacedim>>(
+                FE_SimplexDGP<dhdim, dhspacedim>(1)));
+          else if (reference_cell == dealii::ReferenceCells::Wedge)
+            finite_elements.emplace_back(
+              std::make_shared<dealii::hp::FECollection<dhdim, dhspacedim>>(
+                FE_WedgeDGP<dhdim, dhspacedim>(1)));
+          else if (reference_cell == dealii::ReferenceCells::Pyramid)
+            finite_elements.emplace_back(
+              std::make_shared<dealii::hp::FECollection<dhdim, dhspacedim>>(
+                FE_PyramidDGP<dhdim, dhspacedim>(1)));
+          else
+            Assert(false, ExcNotImplemented());
+        }
     }
   return finite_elements;
 }
diff --git a/source/numerics/data_out.cc b/source/numerics/data_out.cc
index 610e1b5060..b467243d1f 100644
--- a/source/numerics/data_out.cc
+++ b/source/numerics/data_out.cc
@@ -1084,10 +1084,14 @@ template <int dim, typename DoFHandlerType>
 void
 DataOut<dim, DoFHandlerType>::build_patches(const unsigned int n_subdivisions)
 {
-  build_patches(StaticMappingQ1<DoFHandlerType::dimension,
-                                DoFHandlerType::space_dimension>::mapping,
-                n_subdivisions,
-                no_curved_cells);
+  AssertDimension(this->triangulation->get_reference_cells().size(), 1);
+
+  build_patches(
+    this->triangulation->get_reference_cells()[0]
+      .template get_default_linear_mapping<DoFHandlerType::dimension,
+                                           DoFHandlerType::space_dimension>(),
+    n_subdivisions,
+    no_curved_cells);
 }
 
 
diff --git a/tests/simplex/data_out_write_vtk_03.cc b/tests/simplex/data_out_write_vtk_03.cc
new file mode 100644
index 0000000000..c186d8e072
--- /dev/null
+++ b/tests/simplex/data_out_write_vtk_03.cc
@@ -0,0 +1,72 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test DataOut::write_vtk() for simplex meshes with only a Triangulation
+// attached.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_wedge_p.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim, int spacedim = dim>
+void
+test()
+{
+  Triangulation<dim, spacedim> tria;
+  GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 4 : 2);
+
+  DataOut<dim> data_out;
+
+  data_out.attach_triangulation(tria);
+
+  data_out.build_patches();
+
+#if false
+  std::ofstream output("test." + std::to_string(dim) + ".vtk");
+  data_out.write_vtk(output);
+#else
+  data_out.write_vtk(deallog.get_file_stream());
+#endif
+}
+
+int
+main()
+{
+  initlog();
+
+  test<2>();
+  test<3>();
+}
diff --git a/tests/simplex/data_out_write_vtk_03.with_simplex_support=on.output b/tests/simplex/data_out_write_vtk_03.with_simplex_support=on.output
new file mode 100644
index 0000000000..5375b1394b
--- /dev/null
+++ b/tests/simplex/data_out_write_vtk_03.with_simplex_support=on.output
@@ -0,0 +1,353 @@
+
+# vtk DataFile Version 3.0
+#This file was generated 
+ASCII
+DATASET UNSTRUCTURED_GRID
+
+POINTS 96 double
+0.00000 0.00000 0
+0.250000 0.00000 0
+0.00000 0.250000 0
+0.250000 0.250000 0
+0.00000 0.250000 0
+0.250000 0.00000 0
+0.250000 0.00000 0
+0.500000 0.00000 0
+0.250000 0.250000 0
+0.500000 0.250000 0
+0.250000 0.250000 0
+0.500000 0.00000 0
+0.500000 0.00000 0
+0.750000 0.00000 0
+0.500000 0.250000 0
+0.750000 0.250000 0
+0.500000 0.250000 0
+0.750000 0.00000 0
+0.750000 0.00000 0
+1.00000 0.00000 0
+0.750000 0.250000 0
+1.00000 0.250000 0
+0.750000 0.250000 0
+1.00000 0.00000 0
+0.00000 0.250000 0
+0.250000 0.250000 0
+0.00000 0.500000 0
+0.250000 0.500000 0
+0.00000 0.500000 0
+0.250000 0.250000 0
+0.250000 0.250000 0
+0.500000 0.250000 0
+0.250000 0.500000 0
+0.500000 0.500000 0
+0.250000 0.500000 0
+0.500000 0.250000 0
+0.500000 0.250000 0
+0.750000 0.250000 0
+0.500000 0.500000 0
+0.750000 0.500000 0
+0.500000 0.500000 0
+0.750000 0.250000 0
+0.750000 0.250000 0
+1.00000 0.250000 0
+0.750000 0.500000 0
+1.00000 0.500000 0
+0.750000 0.500000 0
+1.00000 0.250000 0
+0.00000 0.500000 0
+0.250000 0.500000 0
+0.00000 0.750000 0
+0.250000 0.750000 0
+0.00000 0.750000 0
+0.250000 0.500000 0
+0.250000 0.500000 0
+0.500000 0.500000 0
+0.250000 0.750000 0
+0.500000 0.750000 0
+0.250000 0.750000 0
+0.500000 0.500000 0
+0.500000 0.500000 0
+0.750000 0.500000 0
+0.500000 0.750000 0
+0.750000 0.750000 0
+0.500000 0.750000 0
+0.750000 0.500000 0
+0.750000 0.500000 0
+1.00000 0.500000 0
+0.750000 0.750000 0
+1.00000 0.750000 0
+0.750000 0.750000 0
+1.00000 0.500000 0
+0.00000 0.750000 0
+0.250000 0.750000 0
+0.00000 1.00000 0
+0.250000 1.00000 0
+0.00000 1.00000 0
+0.250000 0.750000 0
+0.250000 0.750000 0
+0.500000 0.750000 0
+0.250000 1.00000 0
+0.500000 1.00000 0
+0.250000 1.00000 0
+0.500000 0.750000 0
+0.500000 0.750000 0
+0.750000 0.750000 0
+0.500000 1.00000 0
+0.750000 1.00000 0
+0.500000 1.00000 0
+0.750000 0.750000 0
+0.750000 0.750000 0
+1.00000 0.750000 0
+0.750000 1.00000 0
+1.00000 1.00000 0
+0.750000 1.00000 0
+1.00000 0.750000 0
+
+CELLS 32 128
+	3	0	1	2
+	3	3	4	5
+	3	6	7	8
+	3	9	10	11
+	3	12	13	14
+	3	15	16	17
+	3	18	19	20
+	3	21	22	23
+	3	24	25	26
+	3	27	28	29
+	3	30	31	32
+	3	33	34	35
+	3	36	37	38
+	3	39	40	41
+	3	42	43	44
+	3	45	46	47
+	3	48	49	50
+	3	51	52	53
+	3	54	55	56
+	3	57	58	59
+	3	60	61	62
+	3	63	64	65
+	3	66	67	68
+	3	69	70	71
+	3	72	73	74
+	3	75	76	77
+	3	78	79	80
+	3	81	82	83
+	3	84	85	86
+	3	87	88	89
+	3	90	91	92
+	3	93	94	95
+
+CELL_TYPES 32
+ 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
+POINT_DATA 96
+# vtk DataFile Version 3.0
+#This file was generated 
+ASCII
+DATASET UNSTRUCTURED_GRID
+
+POINTS 160 double
+0.00000 0.00000 0.00000
+0.500000 0.00000 0.00000
+0.00000 0.500000 0.00000
+0.00000 0.00000 0.500000
+0.00000 0.500000 0.00000
+0.500000 0.00000 0.00000
+0.500000 0.500000 0.00000
+0.500000 0.500000 0.500000
+0.500000 0.00000 0.00000
+0.00000 0.00000 0.500000
+0.500000 0.00000 0.500000
+0.500000 0.500000 0.500000
+0.00000 0.500000 0.00000
+0.00000 0.00000 0.500000
+0.500000 0.500000 0.500000
+0.00000 0.500000 0.500000
+0.500000 0.00000 0.00000
+0.00000 0.500000 0.00000
+0.00000 0.00000 0.500000
+0.500000 0.500000 0.500000
+0.500000 0.00000 0.00000
+1.00000 0.00000 0.00000
+1.00000 0.500000 0.00000
+1.00000 0.00000 0.500000
+0.500000 0.00000 0.00000
+1.00000 0.500000 0.00000
+0.500000 0.500000 0.00000
+0.500000 0.500000 0.500000
+0.500000 0.00000 0.00000
+0.500000 0.00000 0.500000
+1.00000 0.00000 0.500000
+0.500000 0.500000 0.500000
+1.00000 0.500000 0.00000
+1.00000 0.00000 0.500000
+1.00000 0.500000 0.500000
+0.500000 0.500000 0.500000
+0.500000 0.00000 0.00000
+1.00000 0.500000 0.00000
+0.500000 0.500000 0.500000
+1.00000 0.00000 0.500000
+0.00000 0.500000 0.00000
+0.500000 0.500000 0.00000
+0.500000 1.00000 0.00000
+0.500000 0.500000 0.500000
+0.00000 0.500000 0.00000
+0.500000 1.00000 0.00000
+0.00000 1.00000 0.00000
+0.00000 1.00000 0.500000
+0.00000 0.500000 0.00000
+0.00000 0.500000 0.500000
+0.500000 0.500000 0.500000
+0.00000 1.00000 0.500000
+0.500000 1.00000 0.00000
+0.500000 0.500000 0.500000
+0.500000 1.00000 0.500000
+0.00000 1.00000 0.500000
+0.00000 0.500000 0.00000
+0.500000 1.00000 0.00000
+0.00000 1.00000 0.500000
+0.500000 0.500000 0.500000
+0.500000 0.500000 0.00000
+1.00000 0.500000 0.00000
+0.500000 1.00000 0.00000
+0.500000 0.500000 0.500000
+0.500000 1.00000 0.00000
+1.00000 0.500000 0.00000
+1.00000 1.00000 0.00000
+1.00000 1.00000 0.500000
+1.00000 0.500000 0.00000
+0.500000 0.500000 0.500000
+1.00000 0.500000 0.500000
+1.00000 1.00000 0.500000
+0.500000 1.00000 0.00000
+0.500000 0.500000 0.500000
+1.00000 1.00000 0.500000
+0.500000 1.00000 0.500000
+1.00000 0.500000 0.00000
+0.500000 1.00000 0.00000
+0.500000 0.500000 0.500000
+1.00000 1.00000 0.500000
+0.00000 0.00000 0.500000
+0.500000 0.00000 0.500000
+0.500000 0.500000 0.500000
+0.500000 0.00000 1.00000
+0.00000 0.00000 0.500000
+0.500000 0.500000 0.500000
+0.00000 0.500000 0.500000
+0.00000 0.500000 1.00000
+0.00000 0.00000 0.500000
+0.00000 0.00000 1.00000
+0.500000 0.00000 1.00000
+0.00000 0.500000 1.00000
+0.500000 0.500000 0.500000
+0.500000 0.00000 1.00000
+0.500000 0.500000 1.00000
+0.00000 0.500000 1.00000
+0.00000 0.00000 0.500000
+0.500000 0.500000 0.500000
+0.00000 0.500000 1.00000
+0.500000 0.00000 1.00000
+0.500000 0.00000 0.500000
+1.00000 0.00000 0.500000
+0.500000 0.500000 0.500000
+0.500000 0.00000 1.00000
+0.500000 0.500000 0.500000
+1.00000 0.00000 0.500000
+1.00000 0.500000 0.500000
+1.00000 0.500000 1.00000
+1.00000 0.00000 0.500000
+0.500000 0.00000 1.00000
+1.00000 0.00000 1.00000
+1.00000 0.500000 1.00000
+0.500000 0.500000 0.500000
+0.500000 0.00000 1.00000
+1.00000 0.500000 1.00000
+0.500000 0.500000 1.00000
+1.00000 0.00000 0.500000
+0.500000 0.500000 0.500000
+0.500000 0.00000 1.00000
+1.00000 0.500000 1.00000
+0.00000 0.500000 0.500000
+0.500000 0.500000 0.500000
+0.00000 1.00000 0.500000
+0.00000 0.500000 1.00000
+0.00000 1.00000 0.500000
+0.500000 0.500000 0.500000
+0.500000 1.00000 0.500000
+0.500000 1.00000 1.00000
+0.500000 0.500000 0.500000
+0.00000 0.500000 1.00000
+0.500000 0.500000 1.00000
+0.500000 1.00000 1.00000
+0.00000 1.00000 0.500000
+0.00000 0.500000 1.00000
+0.500000 1.00000 1.00000
+0.00000 1.00000 1.00000
+0.500000 0.500000 0.500000
+0.00000 1.00000 0.500000
+0.00000 0.500000 1.00000
+0.500000 1.00000 1.00000
+0.500000 0.500000 0.500000
+1.00000 0.500000 0.500000
+1.00000 1.00000 0.500000
+1.00000 0.500000 1.00000
+0.500000 0.500000 0.500000
+1.00000 1.00000 0.500000
+0.500000 1.00000 0.500000
+0.500000 1.00000 1.00000
+0.500000 0.500000 0.500000
+0.500000 0.500000 1.00000
+1.00000 0.500000 1.00000
+0.500000 1.00000 1.00000
+1.00000 1.00000 0.500000
+1.00000 0.500000 1.00000
+1.00000 1.00000 1.00000
+0.500000 1.00000 1.00000
+0.500000 0.500000 0.500000
+1.00000 1.00000 0.500000
+0.500000 1.00000 1.00000
+1.00000 0.500000 1.00000
+
+CELLS 40 200
+	4	0	1	2	3
+	4	4	5	6	7
+	4	8	9	10	11
+	4	12	13	14	15
+	4	16	17	18	19
+	4	20	21	22	23
+	4	24	25	26	27
+	4	28	29	30	31
+	4	32	33	34	35
+	4	36	37	38	39
+	4	40	41	42	43
+	4	44	45	46	47
+	4	48	49	50	51
+	4	52	53	54	55
+	4	56	57	58	59
+	4	60	61	62	63
+	4	64	65	66	67
+	4	68	69	70	71
+	4	72	73	74	75
+	4	76	77	78	79
+	4	80	81	82	83
+	4	84	85	86	87
+	4	88	89	90	91
+	4	92	93	94	95
+	4	96	97	98	99
+	4	100	101	102	103
+	4	104	105	106	107
+	4	108	109	110	111
+	4	112	113	114	115
+	4	116	117	118	119
+	4	120	121	122	123
+	4	124	125	126	127
+	4	128	129	130	131
+	4	132	133	134	135
+	4	136	137	138	139
+	4	140	141	142	143
+	4	144	145	146	147
+	4	148	149	150	151
+	4	152	153	154	155
+	4	156	157	158	159
+
+CELL_TYPES 40
+ 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
+POINT_DATA 160
-- 
2.39.5