1 #ifndef DUNE_MULTIDOMAINGRID_TESTS_OUTPUT_HH
2 #define DUNE_MULTIDOMAINGRID_TESTS_OUTPUT_HH
3
4 #include <iostream>
5 #include <iomanip>
6 #include <sstream>
7 #include <dune/grid/io/file/vtk/vtkwriter.hh>
8 #include <dune/grid/multidomaingrid.hh>
9
10 template<int dim>
11 struct AllLayout
12 {
containsAllLayout13 bool contains(Dune::GeometryType gt) {
14 return true;
15 }
16 };
17
18 template<typename GridView, typename InterfaceIterator>
vtkOut(GridView gv,std::string filename,InterfaceIterator iit,InterfaceIterator iend)19 void vtkOut(GridView gv,std::string filename, InterfaceIterator iit, InterfaceIterator iend) {
20
21 Dune::MultiDomainMCMGMapper<GridView> mapper(gv,[](Dune::GeometryType,int) { return true; });
22
23 std::vector<int> hcid(gv.indexSet().size(0),0);
24 std::vector<int> sdc0(gv.indexSet().size(0),0);
25
26 std::vector<int> sdc1(gv.indexSet().size(0),0);
27 for (const auto & cell : elements(gv)) {
28 typename GridView::IndexSet::IndexType idx = gv.indexSet().index(cell);
29 const typename GridView::Grid::MDGridTraits::template Codim<0>::SubDomainSet& sds = gv.indexSet().subDomains(cell);
30 hcid[idx] = idx;
31 sdc0[idx] = sds.contains(0) ? gv.indexSet().index(0,cell) : -1;
32 sdc1[idx] = sds.contains(1) ? gv.indexSet().index(1,cell) : -1;
33 typename GridView::IndexSet::IndexType tmp;
34 assert((sdc0[idx] > -1) == mapper.contains(0,cell,tmp));
35 assert((sdc1[idx] > -1) == mapper.contains(1,cell,tmp));
36 }
37 std::vector<int> hvid(gv.indexSet().size(2),0);
38 std::vector<int> sdv0(gv.indexSet().size(2),0);
39 std::vector<int> sdv1(gv.indexSet().size(2),0);
40 for (const auto& vertex : vertices(gv)) {
41 typename GridView::IndexSet::IndexType idx = gv.indexSet().index(vertex);
42 const typename GridView::Grid::MDGridTraits::template Codim<2>::SubDomainSet& sds = gv.indexSet().subDomains(vertex);
43 hvid[idx] = idx;
44 sdv0[idx] = sds.contains(0) ? gv.indexSet().index(0,vertex) : -1;
45 sdv1[idx] = sds.contains(1) ? gv.indexSet().index(1,vertex) : -1;
46 typename GridView::IndexSet::IndexType tmp;
47 assert((sdv0[idx] > -1) == mapper.contains(0,vertex,tmp));
48 assert((sdv1[idx] > -1) == mapper.contains(1,vertex,tmp));
49 }
50
51 std::vector<int> borderCells(gv.indexSet().size(0),0);
52 std::vector<int> borderVertices(gv.indexSet().size(2),0);
53
54 for(; iit != iend; ++iit) {
55 gv.grid().subDomain(0).subDomainIntersectionIterator(iit);
56 borderCells[gv.indexSet().index(iit->firstCell())] = 1;
57 borderCells[gv.indexSet().index(iit->secondCell())] = 2;
58 }
59
60 Dune::VTKWriter<GridView> vtkWriter(gv);
61 vtkWriter.addCellData(sdc0,"cell_subdomain0");
62 vtkWriter.addCellData(sdc1,"cell_subdomain1");
63 vtkWriter.addCellData(hcid,"cell_hostIndex");
64 vtkWriter.addCellData(borderCells,"borderCells");
65 vtkWriter.addVertexData(sdv0,"vertex_subdomain0");
66 vtkWriter.addVertexData(sdv1,"vertex_subdomain1");
67 vtkWriter.addVertexData(hvid,"vertex_hostIndex");
68 vtkWriter.write(filename,Dune::VTK::ascii);
69 }
70
71 template<typename GridView>
vtkOut2(GridView gv,std::string filename)72 void vtkOut2(GridView gv,std::string filename) {
73 std::vector<int> cid(gv.indexSet().size(0),0);
74 for (const auto& cell : elements(gv)) {
75 typename GridView::IndexSet::IndexType idx = gv.indexSet().index(cell);
76 cid[idx] = idx;
77 }
78 std::vector<int> vid(gv.indexSet().size(2),0);
79 for (const auto& vertex : vertices(gv)) {
80 typename GridView::IndexSet::IndexType idx = gv.indexSet().index(vertex);
81 vid[idx] = idx;
82 }
83 Dune::VTKWriter<GridView> vtkWriter(gv);
84 vtkWriter.addCellData(cid,"cellIndex");
85 vtkWriter.addVertexData(vid,"vertexIndex");
86 vtkWriter.write(filename,Dune::VTK::ascii);
87 }
88
89
90 template<typename stream>
setup(stream & s,std::string prefix,int counter)91 stream& setup(stream& s, std::string prefix, int counter) {
92 s.str("");
93 s << prefix << "_";
94 if (counter >= 0)
95 s << counter << "_";
96 return s;
97 }
98
99
100 template<typename Grid>
printStatus(Grid & grid,std::string prefix,int counter=-1)101 void printStatus(Grid& grid, std::string prefix, int counter = -1) {
102 std::ostringstream s;
103
104 setup(s,prefix,counter) << "leafGridView";
105 vtkOut(grid.leafGridView(),s.str(),grid.leafSubDomainInterfaceBegin(0,1),grid.leafSubDomainInterfaceEnd(0,1));
106
107 for (int level = 0; level <= grid.maxLevel(); ++level) {
108 setup(s,prefix,counter) << "levelGridView_" << std::setw(2) << std::setfill('0') << level;
109 vtkOut(grid.levelGridView(level),s.str(),grid.levelSubDomainInterfaceBegin(0,1,level),grid.levelSubDomainInterfaceEnd(0,1,level));
110 }
111
112 setup(s,prefix,counter) << "subdomain_0";
113 vtkOut2(grid.subDomain(0).leafGridView(),s.str());
114 setup(s,prefix,counter) << "subdomain_1";
115 vtkOut2(grid.subDomain(1).leafGridView(),s.str());
116 }
117
118
119 #endif // DUNE_MULTIDOMAINGRID_TESTS_OUTPUT_HH
120