Line data Source code
1 : /*
2 : * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Jan Friebertshäuser
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/ordering_strategies/algorithms/downwindorder.h"
35 : #include "common/common.h"
36 : #include <iostream>
37 :
38 : #include <utility> // for pair
39 : #include <map> // for graph structure
40 : #include "lib_disc/domain.h"
41 : #include "lib_disc/function_spaces/dof_position_util.h"
42 :
43 : namespace ug{
44 : /**
45 : * For Debugging:
46 : */
47 : DebugID LIB_DISC_ORDER("LIB_DISC_ORDER");
48 :
49 : // Numbers vertices based on a adjacency matrix.
50 0 : void NumeriereKnoten(const std::vector<std::vector<size_t> > &vvConnections,
51 : std::vector<bool> &vVisited, std::vector<size_t> & vAncestorsCount, std::vector<size_t> & vNewIndex, size_t & N, size_t v)
52 : {
53 : vVisited[v] = true;
54 0 : vNewIndex[v] = N;
55 0 : N++;
56 0 : const std::vector<size_t> connections = vvConnections[v];
57 : std::vector<size_t>::const_iterator AdjacencIter;
58 0 : for (AdjacencIter = connections.begin(); AdjacencIter != connections.end(); ++AdjacencIter)
59 : {
60 0 : vAncestorsCount[*AdjacencIter]--;
61 0 : if (vAncestorsCount[*AdjacencIter] == 0)
62 0 : NumeriereKnoten(vvConnections, vVisited, vAncestorsCount, vNewIndex, N, *AdjacencIter);
63 : }
64 0 : }
65 :
66 : // Calculates Downwind-Numbering for one DofDistribution only.
67 : template <typename TDomain>
68 0 : void OrderDownwindForDofDist(SmartPtr<DoFDistribution> dd, ConstSmartPtr<TDomain> domain,
69 : SmartPtr<UserData<MathVector<TDomain::dim>, TDomain::dim> > spVelocity,
70 : number time, int si, number threshold)
71 : {
72 : static const int dim = TDomain::dim;
73 : const size_t num_ind = dd->num_indices();
74 : typedef typename std::pair<MathVector<dim>, size_t> pos_type;
75 : typedef typename std::vector<std::vector<size_t> > adjacency_type;
76 :
77 : // get positions of indices
78 : typename std::vector<pos_type> vPositions;
79 0 : ExtractPositions(domain, dd, vPositions);
80 :
81 : // get adjacency vector of vectors
82 : adjacency_type vvConnections;
83 0 : dd->get_connections(vvConnections);
84 :
85 : // Check vector sizes match
86 0 : if (vvConnections.size() != num_ind)
87 0 : UG_THROW("OrderDownstreamForDofDist: "
88 : "Adjacency list of dimension " << num_ind << " expected, got "<< vvConnections.size());
89 :
90 0 : if (vPositions.size() != num_ind)
91 0 : UG_THROW("OrderDownstreamForDofDist: "
92 : "Position list of dimension " << num_ind << " expected, got "<< vPositions.size());
93 :
94 : // init helper structures
95 0 : std::vector<size_t> vNewIndex(num_ind, 0);
96 0 : std::vector<size_t> vAncestorsCount(num_ind, 0);
97 0 : std::vector<bool> vVisited(num_ind, false);
98 :
99 : // remove connections that are not in stream direction
100 : adjacency_type::iterator VertexIter;
101 : std::vector<size_t>::iterator AdjIter;
102 : std::vector<size_t> vAdjacency;
103 :
104 : MathVector<TDomain::dim> vVel1, vPos1, vPos2, vDir1_2;
105 : size_t i;
106 0 : for (VertexIter = vvConnections.begin(), i=0; VertexIter != vvConnections.end(); VertexIter++, i++)
107 : {
108 : UG_DLOG(LIB_DISC_ORDER, 2, "Filtering vertex " << i << " adjacency vector." <<std::endl);
109 : // count how many vertex were kept / removed per adjacency vector
110 : #ifdef UG_ENABLE_DEBUG_LOGS
111 : size_t initialcount = VertexIter->size();
112 : #endif
113 : size_t kept = 0, removed = 0;
114 : // get position and velocity of first trait
115 : vPos1 = vPositions.at(i).first;
116 0 : (*spVelocity)(vVel1, vPos1, time, si);
117 0 : if (VecLengthSq(vVel1) == 0 )
118 : {
119 : // if the velocity is zero at this trait it does not interfere with others
120 : // NOTE: otherwise this trait would be downwind-connected to all of it's neighbors
121 : // NOTE: VertexIter-> will access inner vector functions (*VertexIter) is the inner vector.
122 : removed = VertexIter->size();
123 : VertexIter->clear();
124 : }
125 : else {
126 : AdjIter = VertexIter->begin();
127 0 : while (AdjIter != VertexIter->end())
128 : {
129 : // get position of second trait
130 0 : vPos2 = vPositions.at(*AdjIter).first;
131 :
132 : // get difference vector as direction vector
133 : VecSubtract(vDir1_2, vPos2, vPos1);
134 :
135 : // compute angle between velocity and direction vector
136 0 : number anglex1_2 = VecAngle(vDir1_2, vVel1);
137 :
138 : // if angle is smaller then threshold continue else remove connection
139 0 : if (anglex1_2 <= threshold && i != *AdjIter)
140 : {
141 0 : vAncestorsCount.at(*AdjIter) += 1;
142 : ++AdjIter;
143 : kept++;
144 : } else {
145 : AdjIter = VertexIter->erase(AdjIter);
146 : removed++;
147 : }
148 : }
149 : }
150 : UG_DLOG(LIB_DISC_ORDER, 2, "Kept: " << kept << ", removed: " << removed << " of " << initialcount
151 : << " entries in adjacency matrix." << std::endl << std::endl);
152 : }
153 : // calculate downwindorder
154 : // Find vertexes without any ancestors and start NumeriereKnoten on them.
155 : size_t v,N;
156 0 : for (v=0, N=0; v < vvConnections.size(); v++)
157 : {
158 0 : if (vAncestorsCount[v] == 0 && !vVisited[v])
159 : {
160 0 : NumeriereKnoten(vvConnections, vVisited, vAncestorsCount, vNewIndex, N, v);
161 : }
162 : }
163 :
164 : // sanity check
165 0 : if (N < vvConnections.size()){
166 : size_t fails = 0;
167 0 : for (v=0; v < vvConnections.size(); v++) {
168 0 : if (!vVisited[v]) {
169 : UG_DLOG(LIB_DISC_ORDER, 2, v << "was not visited, has unresolved ancestors: " << vAncestorsCount[v] << std::endl);
170 0 : fails ++;
171 : }
172 : }
173 0 : UG_THROW("OrderDownwindForDist failed, " << fails << " traits unvisited." << std::endl);
174 : }
175 :
176 : // reorder traits
177 0 : dd->permute_indices(vNewIndex);
178 0 : }
179 :
180 : // Calculates Downwind Numbering for all DofDistributions of one Domain.
181 : template <typename TDomain>
182 0 : void OrderDownwind(ApproximationSpace<TDomain>& approxSpace,
183 : SmartPtr<UserData<MathVector<TDomain::dim>, TDomain::dim> > spVelocity, number threshold)
184 : {
185 : UG_LOG ("OrderDownwind: This function is obsolete and may cause problems. Avoid it! Alternatives: Ordering strategies in solvers etc.\n");
186 :
187 : // TODO: implement for variable time and subset
188 : number time = 0.0;
189 : int si = 0;
190 0 : std::vector<SmartPtr<DoFDistribution> > vDD = approxSpace.dof_distributions();
191 : UG_DLOG(LIB_DISC_ORDER, 2, "Starting DownwindOrdering." << std::endl);
192 0 : for(size_t i = 0; i < vDD.size(); ++i){
193 : UG_DLOG(LIB_DISC_ORDER, 2, "Ordering Domain Distribution " << i << "." << std::endl);
194 0 : OrderDownwindForDofDist<TDomain>(vDD[i], approxSpace.domain(), spVelocity, time, si, threshold);
195 : }
196 0 : }
197 :
198 : // Calculates Downwind Numbering for all DofDistributions of one Domain.
199 : template <typename TDomain>
200 0 : void OrderDownwind(ApproximationSpace<TDomain>& approxSpace,
201 : SmartPtr<UserData<MathVector<TDomain::dim>, TDomain::dim> > spVelocity)
202 : {
203 0 : OrderDownwind<TDomain>(approxSpace, spVelocity, PI/4.0);
204 0 : }
205 :
206 : // Calculates Downwind Numbering for all DofDistributions of one Domain.
207 : template <typename TDomain>
208 0 : void OrderDownwind(ApproximationSpace<TDomain>& approxSpace,
209 : const std::vector<number>& vVel)
210 : {
211 : static const int dim = TDomain::dim;
212 0 : if(vVel.size() != dim){
213 0 : UG_THROW("OrderDownstream: Velocity field of dimension " << dim << " expected, got "<< vVel.size());
214 : }
215 :
216 0 : OrderDownwind<TDomain>(approxSpace, SmartPtr<ConstUserVector<dim> >(new ConstUserVector<dim>(vVel)) );
217 0 : }
218 :
219 : // Calculates Downwind Numbering for all DofDistributions of one Domain.
220 : template <typename TDomain>
221 0 : void OrderDownwind(ApproximationSpace<TDomain>& approxSpace,
222 : const std::vector<number>& vVel, number threshold)
223 : {
224 : static const int dim = TDomain::dim;
225 0 : if(vVel.size() != dim){
226 0 : UG_THROW("OrderDownstream: Velocity field of dimension " << dim << " expected, got "<< vVel.size());
227 : }
228 :
229 0 : OrderDownwind<TDomain>(approxSpace, SmartPtr<ConstUserVector<dim> >(new ConstUserVector<dim>(vVel)), threshold );
230 0 : }
231 :
232 :
233 : #ifdef UG_FOR_LUA
234 :
235 : // Calculates Downwind Numbering for all DofDistributions of one Domain.
236 : template <typename TDomain>
237 0 : void OrderDownwind(ApproximationSpace<TDomain>& approxSpace, const char* strVelocity)
238 : {
239 : static const int dim = TDomain::dim;
240 :
241 0 : SmartPtr<UserData<MathVector<dim>, dim> > spVelocity
242 0 : = make_sp(new LuaUserData<MathVector<dim>, dim>(strVelocity));
243 :
244 0 : OrderDownwind<TDomain>(approxSpace, spVelocity);
245 0 : }
246 :
247 : // Calculates Downwind Numbering for all DofDistributions of one Domain.
248 : template <typename TDomain>
249 0 : void OrderDownwind(ApproximationSpace<TDomain>& approxSpace, const char* strVelocity, number threshold)
250 : {
251 : static const int dim = TDomain::dim;
252 :
253 0 : SmartPtr<UserData<MathVector<dim>, dim> > spVelocity
254 0 : = make_sp(new LuaUserData<MathVector<dim>, dim>(strVelocity));
255 :
256 0 : OrderDownwind<TDomain>(approxSpace, spVelocity, threshold);
257 0 : }
258 :
259 : #ifdef UG_DIM_1
260 : template void OrderDownwind<Domain1d>(ApproximationSpace<Domain1d>& approxSpace, const char* strVelocity);
261 : template void OrderDownwind<Domain1d>(ApproximationSpace<Domain1d>& approxSpace, const char* strVelocity, number threshold);
262 : #endif
263 : #ifdef UG_DIM_2
264 : template void OrderDownwind<Domain2d>(ApproximationSpace<Domain2d>& approxSpace, const char* strVelocity);
265 : template void OrderDownwind<Domain2d>(ApproximationSpace<Domain2d>& approxSpace, const char* strVelocity, number threshold);
266 : #endif
267 : #ifdef UG_DIM_3
268 : template void OrderDownwind<Domain3d>(ApproximationSpace<Domain3d>& approxSpace, const char* strVelocity);
269 : template void OrderDownwind<Domain3d>(ApproximationSpace<Domain3d>& approxSpace, const char* strVelocity, number threshold);
270 : #endif
271 :
272 : #endif
273 :
274 : #ifdef UG_DIM_1
275 : template void OrderDownwind<Domain1d>(ApproximationSpace<Domain1d>& approxSpace, SmartPtr<UserData<MathVector<Domain1d::dim>, Domain1d::dim> > spVelocity);
276 : template void OrderDownwind<Domain1d>(ApproximationSpace<Domain1d>& approxSpace, SmartPtr<UserData<MathVector<Domain1d::dim>, Domain1d::dim> > spVelocity, number threshold);
277 : template void OrderDownwind<Domain1d>(ApproximationSpace<Domain1d>& approxSpace, const std::vector<number>& vVel);
278 : template void OrderDownwind<Domain1d>(ApproximationSpace<Domain1d>& approxSpace, const std::vector<number>& vVel, number threshold);
279 : #endif
280 :
281 : #ifdef UG_DIM_2
282 : template void OrderDownwind<Domain2d>(ApproximationSpace<Domain2d>& approxSpace, SmartPtr<UserData<MathVector<Domain2d::dim>, Domain2d::dim> > spVelocity);
283 : template void OrderDownwind<Domain2d>(ApproximationSpace<Domain2d>& approxSpace, SmartPtr<UserData<MathVector<Domain2d::dim>, Domain2d::dim> > spVelocity, number threshold);
284 : template void OrderDownwind<Domain2d>(ApproximationSpace<Domain2d>& approxSpace, const std::vector<number>& vVel);
285 : template void OrderDownwind<Domain2d>(ApproximationSpace<Domain2d>& approxSpace, const std::vector<number>& vVel, number threshold);
286 : #endif
287 :
288 : #ifdef UG_DIM_3
289 : template void OrderDownwind<Domain3d>(ApproximationSpace<Domain3d>& approxSpace, SmartPtr<UserData<MathVector<Domain3d::dim>, Domain3d::dim> > spVelocity);
290 : template void OrderDownwind<Domain3d>(ApproximationSpace<Domain3d>& approxSpace, SmartPtr<UserData<MathVector<Domain3d::dim>, Domain3d::dim> > spVelocity, number threshold);
291 : template void OrderDownwind<Domain3d>(ApproximationSpace<Domain3d>& approxSpace, const std::vector<number>& vVel);
292 : template void OrderDownwind<Domain3d>(ApproximationSpace<Domain3d>& approxSpace, const std::vector<number>& vVel, number threshold);
293 : #endif
294 :
295 : } // end namespace ug
|
