Line data Source code
1 : /*
2 : * Copyright (c) 2011-2022: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel, 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 : #include <algorithm>
34 : #include <vector>
35 : #include <queue>
36 : #include <utility>
37 :
38 : #include "common/common.h"
39 : #include "lib_disc/function_spaces/dof_position_util.h"
40 : #include "lib_disc/reference_element/reference_element_util.h"
41 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
42 : #include "lib_disc/domain.h"
43 :
44 : #include "lib_disc/ordering_strategies/algorithms/lexorder_comparators.cpp"
45 : #include "lib_disc/ordering_strategies/algorithms/lexorder.h"
46 :
47 : namespace ug{
48 :
49 : template<int dim>
50 0 : void ComputeLexicographicOrder(std::vector<size_t>& vNewIndex,
51 : std::vector<std::pair<MathVector<dim>, size_t> >& vPos,
52 : size_t orderDim, bool increasing)
53 : {
54 : // a) order all indices
55 0 : if(vNewIndex.size() == vPos.size()){
56 : // sort indices based on their position
57 0 : if(increasing){
58 0 : if (orderDim == 0)
59 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDim<dim, 0>);
60 0 : else if (orderDim == 1)
61 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDim<dim, 1>);
62 0 : else if (orderDim == 2)
63 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDim<dim, 2>);
64 0 : else UG_THROW("Invalid sorting direction.");
65 : }
66 : else{
67 0 : if (orderDim == 0)
68 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDimDec<dim, 0>);
69 0 : else if (orderDim == 1)
70 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDimDec<dim, 1>);
71 0 : else if (orderDim == 2)
72 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDimDec<dim, 2>);
73 0 : else UG_THROW("Invalid sorting direction.");
74 : }
75 :
76 : // write mapping
77 0 : for (size_t i=0; i < vPos.size(); ++i)
78 0 : vNewIndex[vPos[i].second] = i;
79 : }
80 : // b) only some indices to order
81 : else{
82 : // create copy of pos
83 0 : std::vector<std::pair<MathVector<dim>, size_t> > vPosOrig(vPos);
84 :
85 : // sort indices based on their position
86 0 : if(increasing){
87 0 : if (orderDim == 0){
88 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDim<dim, 0>);
89 : }
90 0 : else if (orderDim == 1)
91 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDim<dim, 1>);
92 0 : else if (orderDim == 2)
93 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDim<dim, 2>);
94 0 : else UG_THROW("Invalid sorting direction.");
95 : }
96 : else{
97 0 : if (orderDim == 0){
98 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDimDec<dim, 0>);
99 : }
100 0 : else if (orderDim == 1)
101 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDimDec<dim, 1>);
102 0 : else if (orderDim == 2)
103 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDimDec<dim, 2>);
104 0 : else UG_THROW("Invalid sorting direction.");
105 : }
106 :
107 : // write mapping
108 0 : for (size_t i=0; i < vNewIndex.size(); ++i)
109 0 : vNewIndex[i] = i;
110 0 : for (size_t i=0; i < vPos.size(); ++i)
111 0 : vNewIndex[vPos[i].second] = vPosOrig[i].second;
112 0 : }
113 0 : }
114 :
115 : /// orders the dof distribution using Cuthill-McKee
116 : template <typename TDomain>
117 0 : void OrderLexForDofDist(SmartPtr<DoFDistribution> dd, ConstSmartPtr<TDomain> domain, size_t orderDim, bool increasing)
118 : {
119 : // Lex Ordering is only possible in this cases:
120 : // b) Same number of DoFs on each geometric object (or no DoFs on object)
121 : // --> in this case we can order all dofs
122 : // a) different trial spaces, but DoFs for each trial spaces only on separate
123 : // geometric objects. (e.g. one space only vertices, one space only on edges)
124 : // --> in this case we can order all geometric objects separately
125 :
126 : // a) check for same number of DoFs on every geometric object
127 : bool bEqualNumDoFOnEachGeomObj = true;
128 : int numDoFOnGeomObj = -1;
129 0 : for(int si = 0; si < dd->num_subsets(); ++si){
130 0 : for(int roid = 0; roid < NUM_REFERENCE_OBJECTS; ++roid){
131 0 : const int numDoF = dd->num_dofs((ReferenceObjectID)roid, si);
132 :
133 0 : if(numDoF == 0) continue;
134 :
135 0 : if(numDoFOnGeomObj == -1){
136 : numDoFOnGeomObj = numDoF;
137 : }
138 : else{
139 0 : if(numDoFOnGeomObj != numDoF)
140 : bEqualNumDoFOnEachGeomObj = false;
141 : }
142 : }
143 : }
144 :
145 : // b) check for non-mixed spaces
146 0 : std::vector<bool> bSingleSpaceUsage(NUM_REFERENCE_OBJECTS, true);
147 0 : std::vector<bool> vHasDoFs(NUM_REFERENCE_OBJECTS, false);
148 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct){
149 :
150 0 : LFEID lfeID = dd->local_finite_element_id(fct);
151 0 : const CommonLocalDoFSet& locDoF = LocalFiniteElementProvider::get_dofs(lfeID);
152 :
153 0 : for(int roid = 0; roid < NUM_REFERENCE_OBJECTS; ++roid){
154 : const int numDoF = locDoF.num_dof((ReferenceObjectID)roid);
155 :
156 0 : if(numDoF <= 0) continue;
157 :
158 0 : if(vHasDoFs[roid] == false){
159 : vHasDoFs[roid] = true;
160 : }
161 : else{
162 : bSingleSpaceUsage[roid] = false;
163 : }
164 : }
165 : }
166 0 : std::vector<bool> bSortableComp(dd->num_fct(), true);
167 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct){
168 :
169 0 : LFEID lfeID = dd->local_finite_element_id(fct);
170 0 : const CommonLocalDoFSet& locDoF = LocalFiniteElementProvider::get_dofs(lfeID);
171 :
172 0 : for(int roid = 0; roid < NUM_REFERENCE_OBJECTS; ++roid){
173 0 : if(locDoF.num_dof((ReferenceObjectID)roid) != 0){
174 0 : if(bSingleSpaceUsage[roid] == false)
175 : bSortableComp[fct] = false;
176 : }
177 : }
178 : }
179 :
180 : // get position attachment
181 : typedef typename std::pair<MathVector<TDomain::dim>, size_t> pos_type;
182 :
183 : // get positions of indices
184 : std::vector<pos_type> vPositions;
185 :
186 : // a) we can order globally
187 0 : if(bEqualNumDoFOnEachGeomObj)
188 : {
189 0 : ExtractPositions(domain, dd, vPositions);
190 :
191 : // get mapping: old -> new index
192 0 : std::vector<size_t> vNewIndex(dd->num_indices());
193 0 : ComputeLexicographicOrder<TDomain::dim>(vNewIndex, vPositions, orderDim, increasing);
194 :
195 : /*
196 : std::vector<bool> vCheck(dd->num_indices(), false);
197 : for (size_t i = 0; i < vNewIndex.size(); ++i)
198 : {
199 : UG_COND_THROW(vCheck.at(vNewIndex[i]), "Double mapping to index " << vNewIndex[i] << ".");
200 : vCheck.at(vNewIndex[i]) = true;
201 : }
202 : for (size_t i = 0; i < vCheck.size(); ++i)
203 : UG_COND_THROW(!vCheck[i], "Nothing maps to index " << i << ".");
204 : */
205 :
206 : // reorder indices
207 0 : dd->permute_indices(vNewIndex);
208 0 : }
209 : // b) we can only order some spaces
210 : else
211 : {
212 : UG_LOG("OrderLex: Cannot order globally, trying to order some components:\n");
213 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct){
214 0 : if(bSortableComp[fct] == false){
215 0 : UG_LOG("OrderLex: '"<<dd->name(fct)<<" NOT SORTED.\n");
216 0 : continue;
217 : }
218 :
219 0 : ExtractPositions(domain, dd, fct, vPositions);
220 :
221 : // get mapping: old -> new index
222 0 : std::vector<size_t> vNewIndex(dd->num_indices());
223 0 : ComputeLexicographicOrder<TDomain::dim>(vNewIndex, vPositions, orderDim, increasing);
224 :
225 : // reorder indices
226 0 : dd->permute_indices(vNewIndex);
227 :
228 0 : UG_LOG("OrderLex: '"<<dd->name(fct)<<" SORTED.\n");
229 : }
230 : }
231 0 : }
232 :
233 :
234 : /// orders the all DofDistributions of the ApproximationSpace using lexicographic order
235 : template <typename TDomain>
236 0 : void OrderLex(ApproximationSpace<TDomain>& approxSpace, const char *order)
237 : {
238 0 : std::vector<SmartPtr<DoFDistribution> > vDD = approxSpace.dof_distributions();
239 :
240 0 : size_t len = strlen(order);
241 :
242 0 : if(len == 0){
243 0 : UG_THROW("OrderLex: Empty direction!");
244 : }
245 :
246 : size_t pos = 0;
247 : bool increasing = true;
248 : char sign;
249 0 : while(pos < len){
250 0 : if(increasing){
251 : sign = '+';
252 : }
253 : else{
254 : sign = '-';
255 : }
256 :
257 0 : switch(order[pos]){
258 0 : case '+':
259 0 : ++pos;
260 : increasing = true;
261 0 : break;
262 0 : case '-':
263 0 : ++pos;
264 : increasing = false;
265 0 : break;
266 : case 'x':
267 0 : UG_LOG("OrderLex: LexOrdering in " << sign << "x direction.\n")
268 0 : for (size_t i = 0; i < vDD.size(); ++i)
269 0 : OrderLexForDofDist<TDomain>(vDD[i], approxSpace.domain(), 0, increasing);
270 0 : ++pos;
271 : increasing = true;
272 0 : break;
273 : case 'y':
274 0 : UG_LOG("OrderLex: LexOrdering in " << sign << "y direction.\n")
275 0 : for (size_t i = 0; i < vDD.size(); ++i)
276 0 : OrderLexForDofDist<TDomain>(vDD[i], approxSpace.domain(), 1, increasing);
277 0 : ++pos;
278 : increasing = true;
279 0 : break;
280 : case 'z':
281 0 : UG_LOG("OrderLex: LexOrdering in " << sign << "z direction.\n")
282 0 : for (size_t i = 0; i < vDD.size(); ++i)
283 0 : OrderLexForDofDist<TDomain>(vDD[i], approxSpace.domain(), 2, increasing);
284 0 : ++pos;
285 : increasing = true;
286 0 : break;
287 0 : default:
288 0 : UG_THROW("OrderLex: Invalid token in direction string, valid tokens: +, -, x, y, z");
289 : }
290 : }
291 0 : }
292 :
293 : #ifdef UG_DIM_1
294 : template void ComputeLexicographicOrder<1>(std::vector<size_t>& vNewIndex, std::vector<std::pair<MathVector<1>, size_t> >& vPos, size_t, bool);
295 : template void OrderLexForDofDist<Domain1d>(SmartPtr<DoFDistribution> dd, ConstSmartPtr<Domain1d> domain, size_t, bool);
296 : template void OrderLex<Domain1d>(ApproximationSpace<Domain1d>& approxSpace, const char *order);
297 : #endif
298 : #ifdef UG_DIM_2
299 : template void ComputeLexicographicOrder<2>(std::vector<size_t>& vNewIndex, std::vector<std::pair<MathVector<2>, size_t> >& vPos, size_t, bool);
300 : template void OrderLexForDofDist<Domain2d>(SmartPtr<DoFDistribution> dd, ConstSmartPtr<Domain2d> domain, size_t, bool);
301 : template void OrderLex<Domain2d>(ApproximationSpace<Domain2d>& approxSpace, const char *order);
302 : #endif
303 : #ifdef UG_DIM_3
304 : template void ComputeLexicographicOrder<3>(std::vector<size_t>& vNewIndex, std::vector<std::pair<MathVector<3>, size_t> >& vPos, size_t, bool);
305 : template void OrderLexForDofDist<Domain3d>(SmartPtr<DoFDistribution> dd, ConstSmartPtr<Domain3d> domain, size_t, bool);
306 : template void OrderLex<Domain3d>(ApproximationSpace<Domain3d>& approxSpace, const char *order);
307 : #endif
308 :
309 : }
|
