Line data Source code
1 : /*
2 : * Copyright (c) 2011-2022: G-CSC, Goethe University Frankfurt
3 : * Author: Lukas Larisch
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 :
34 : #include "lib_disc/domain.h"
35 : #include "lib_disc/function_spaces/grid_function.h"
36 : #include "lib_disc/function_spaces/grid_function_user_data.h"
37 :
38 : #include "lib_disc/io/vtkoutput.h" ////IteratorProvider
39 :
40 : namespace ug{
41 :
42 : template <typename TDomain, typename TAlgebra>
43 0 : class GridPointsOrdering
44 : {
45 : typedef GridFunction<TDomain, TAlgebra> TGridFunction;
46 : typedef GridFunctionNumberData<GridFunction<TDomain, TAlgebra> > TGridFunctionNumberData;
47 :
48 : typedef ug::Attachment<int> AVrtIndex;
49 :
50 : public:
51 0 : GridPointsOrdering(SmartPtr<TGridFunction> spGridFct, const char* name)
52 : {
53 :
54 0 : m_u = spGridFct->clone_without_values();
55 0 : m_name = name;
56 0 : m_n = 0;
57 0 : m_numVert = 0;
58 0 : m_numElem = 0;
59 0 : m_numConn = 0;
60 :
61 0 : create_vtkoutput_ordering();
62 :
63 0 : std::vector<size_t> index(m_numVert);
64 :
65 : size_t k = 0;
66 :
67 : typedef typename IteratorProvider<TGridFunction>::template traits<ug::Vertex>::const_iterator const_iterator;
68 : const_iterator iterBegin = IteratorProvider<TGridFunction>::template begin<ug::Vertex>(*m_u, -1);
69 : const_iterator iterEnd = IteratorProvider<TGridFunction>::template end<ug::Vertex>(*m_u, -1);
70 0 : for( ; iterBegin != iterEnd; ++iterBegin, ++k)
71 : {
72 : ug::Vertex* v = *iterBegin;
73 0 : index[k] = m_aaVrtIndex[v];
74 : }
75 :
76 : k = 0;
77 0 : std::vector<DoFIndex> ind(1);
78 :
79 0 : iterBegin = IteratorProvider<TGridFunction>::template begin<ug::Vertex>(*m_u, -1);
80 0 : iterEnd = IteratorProvider<TGridFunction>::template end<ug::Vertex>(*m_u, -1);
81 0 : for( ; iterBegin != iterEnd; ++iterBegin, ++k)
82 : {
83 : ug::Vertex* v = *iterBegin;
84 :
85 : // get vector holding all indices on the vertex
86 : m_u->inner_dof_indices(v, 0, ind);
87 0 : DoFRef(*m_u, ind[0]) = index[k];
88 : }
89 0 : }
90 :
91 :
92 : template <typename TElem>
93 0 : void count_sizes()
94 : {
95 : // get the grid associated to the solution
96 0 : Grid& grid = *m_u->domain()->grid();
97 :
98 : // get reference element
99 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
100 :
101 : // number of corners of element
102 : static const int numCo = ref_elem_type::numCorners;
103 :
104 : // get iterators
105 : typedef typename IteratorProvider<TGridFunction>::template traits<TElem>::const_iterator const_iterator;
106 : const_iterator iterBegin = IteratorProvider<TGridFunction>::template begin<TElem>(*m_u, -1);
107 : const_iterator iterEnd = IteratorProvider<TGridFunction>::template end<TElem>(*m_u, -1);
108 :
109 : // loop elements
110 0 : for( ; iterBegin != iterEnd; ++iterBegin)
111 : {
112 : // get the element
113 : TElem *elem = *iterBegin;
114 :
115 : // count number of elements and number of connections;
116 : // handle octahedrons separately by splitting into a top and bottom pyramid
117 : if(ref_elem_type::REFERENCE_OBJECT_ID != ROID_OCTAHEDRON)
118 : {
119 0 : ++m_numElem;
120 0 : m_numConn += numCo;
121 : }
122 : else
123 : {
124 : // counting top and bottom pyramid
125 0 : m_numElem += 2;
126 0 : m_numConn += 10;
127 : }
128 :
129 : // loop vertices of the element
130 0 : for(int i = 0; i < numCo; ++i)
131 : {
132 : // get vertex of the element
133 0 : Vertex* v = GetVertex(elem, i);
134 :
135 : // if this vertex has already been counted, skip it
136 0 : if(grid.is_marked(v)) continue;
137 :
138 : // count vertex and mark it
139 0 : ++m_numVert;
140 : grid.mark(v);
141 : }
142 : }
143 0 : }
144 :
145 :
146 : template <typename TElem>
147 0 : void number_points_elementwise()
148 : {
149 : // get the grid associated to the solution
150 0 : Grid& grid = *m_u->domain()->grid();
151 :
152 : // get reference element
153 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
154 : // get iterators
155 : typedef typename IteratorProvider<TGridFunction>::template traits<TElem>::const_iterator const_iterator;
156 : const_iterator iterBegin = IteratorProvider<TGridFunction>::template begin<TElem>(*m_u, -1);
157 : const_iterator iterEnd = IteratorProvider<TGridFunction>::template end<TElem>(*m_u, -1);
158 :
159 : // loop all elements of the subset
160 0 : for( ; iterBegin != iterEnd; ++iterBegin)
161 : {
162 : // get the element
163 : TElem *elem = *iterBegin;
164 :
165 : // loop vertices of the element
166 0 : for(size_t i = 0; i < (size_t) ref_elem_type::numCorners; ++i)
167 : {
168 : // get vertex of element
169 : Vertex* v = GetVertex(elem, i);
170 :
171 : // if vertex has already be handled, skip it
172 0 : if(grid.is_marked(v)) continue;
173 :
174 : // mark the vertex as processed
175 : grid.mark(v);
176 :
177 : // number vertex
178 0 : m_aaVrtIndex[v] = m_n++;
179 : }
180 : }
181 0 : }
182 :
183 0 : void create_vtkoutput_ordering(){
184 : // check functions
185 0 : for(size_t fct = 0; fct < m_u->num_fct(); ++fct)
186 : {
187 : // check if function is defined everywhere
188 0 : if(!m_u->is_def_everywhere(fct))
189 0 : UG_THROW("only serial case implemented!");
190 : }
191 :
192 : // get the grid associated to the solution
193 0 : Grid& grid = *m_u->domain()->grid();
194 :
195 : // attach help indices
196 : AVrtIndex aVrtIndex;
197 : grid.attach_to_vertices(aVrtIndex);
198 : m_aaVrtIndex.access(grid, aVrtIndex);
199 :
200 0 : int dim = DimensionOfSubsets(*m_u->domain()->subset_handler());
201 :
202 : // Count needed sizes for vertices, elements and connections
203 : try{
204 : // reset all marks
205 0 : grid.begin_marking();
206 :
207 : // switch dimension
208 0 : switch(dim)
209 : {
210 0 : case 0: count_sizes<Vertex>(); break;
211 0 : case 1: count_sizes<RegularEdge>();
212 0 : count_sizes<ConstrainingEdge>(); break;
213 0 : case 2: count_sizes<Triangle>();
214 0 : count_sizes<Quadrilateral>();
215 0 : count_sizes<ConstrainingTriangle>();
216 0 : count_sizes<ConstrainingQuadrilateral>(); break;
217 0 : case 3: count_sizes<Tetrahedron>();
218 0 : count_sizes<Pyramid>();
219 0 : count_sizes<Prism>();
220 0 : count_sizes<Octahedron>();
221 0 : count_sizes<Hexahedron>(); break;
222 0 : default: UG_THROW(name() << "::create_vtkoutput_ordering: Dimension " << dim << " is not supported.");
223 : }
224 :
225 : // signal end of marking
226 0 : grid.end_marking();
227 : }
228 0 : UG_CATCH_THROW(name() << "::create_vtkoutput_ordering: Can not count piece sizes.");
229 :
230 : // start marking of vertices
231 0 : grid.begin_marking();
232 :
233 : // switch dimension
234 0 : if(m_numVert > 0){
235 0 : switch(dim){
236 0 : case 0: number_points_elementwise<Vertex>(); break;
237 0 : case 1: number_points_elementwise<RegularEdge>();
238 0 : number_points_elementwise<ConstrainingEdge>(); break;
239 0 : case 2: number_points_elementwise<Triangle>();
240 0 : number_points_elementwise<Quadrilateral>();
241 0 : number_points_elementwise<ConstrainingTriangle>();
242 0 : number_points_elementwise<ConstrainingQuadrilateral>(); break;
243 0 : case 3: number_points_elementwise<Tetrahedron>();
244 0 : number_points_elementwise<Pyramid>();
245 0 : number_points_elementwise<Prism>();
246 0 : number_points_elementwise<Octahedron>();
247 0 : number_points_elementwise<Hexahedron>(); break;
248 : default: UG_THROW(name() << "::create_vtkoutput_ordering: Dimension " << dim << " is not supported.");
249 : }
250 : }
251 : // signal end of marking the grid
252 0 : grid.end_marking();
253 0 : }
254 :
255 0 : SmartPtr<TGridFunctionNumberData> get(){
256 0 : return SmartPtr<TGridFunctionNumberData>(new TGridFunctionNumberData(m_u, m_name));
257 : }
258 :
259 : const char* name() const {return "GridPointsOrdering";}
260 :
261 : private:
262 : SmartPtr<TGridFunction> m_u;
263 : const char* m_name;
264 : size_t m_n;
265 : int m_numVert;
266 : int m_numElem;
267 : int m_numConn;
268 :
269 : Grid::VertexAttachmentAccessor<AVrtIndex> m_aaVrtIndex;
270 : };
271 :
272 : } //namespace
|
