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 : #ifndef __H__UG__LIB_DISC__ORDERING_STRATEGIES_ALGORITHMS__RIVERORDER__
34 : #define __H__UG__LIB_DISC__ORDERING_STRATEGIES_ALGORITHMS__RIVERORDER__
35 :
36 : #include <boost/graph/adjacency_list.hpp>
37 : #include <boost/graph/graph_traits.hpp>
38 : #include <boost/graph/properties.hpp>
39 :
40 : #include <vector>
41 : #include <utility> //for pair
42 : #include <climits> //INT_MAX
43 :
44 : #include "lib_disc/domain.h"
45 : #include "lib_disc/function_spaces/grid_function.h"
46 :
47 : #include "lib_algebra/ordering_strategies/algorithms/IOrderingAlgorithm.h"
48 : #include "lib_algebra/ordering_strategies/algorithms/util.h"
49 :
50 : #include "common/error.h"
51 :
52 : namespace ug{
53 :
54 :
55 :
56 : template <typename TAlgebra, typename TDomain, typename O_t>
57 : class RiverOrdering : public IOrderingAlgorithm<TAlgebra, O_t>
58 : {
59 : public:
60 : typedef typename TAlgebra::matrix_type M_t;
61 : typedef typename TAlgebra::vector_type V_t;
62 : typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::bidirectionalS> G_t;
63 : typedef IOrderingAlgorithm<TAlgebra, O_t> baseclass;
64 :
65 : /// Grid function type for the solution
66 : typedef GridFunction<TDomain, TAlgebra> GridFunc_t;
67 :
68 : typedef typename boost::graph_traits<G_t>::vertex_descriptor vd;
69 : typedef typename boost::graph_traits<G_t>::adjacency_iterator adj_iter;
70 :
71 0 : RiverOrdering() : m_ssIdx(-1){}
72 :
73 : /// clone constructor
74 : RiverOrdering( const LexOrdering<TAlgebra, TDomain, O_t> &parent )
75 : : baseclass(), m_ssIdx(parent.m_ssIdx){}
76 :
77 0 : SmartPtr<IOrderingAlgorithm<TAlgebra, O_t> > clone()
78 : {
79 0 : return make_sp(new RiverOrdering<TAlgebra, TDomain, O_t>(*this));
80 : }
81 :
82 :
83 0 : vd get_source_vertex(std::vector<BOOL>& visited, G_t& g){
84 0 : for(unsigned i = 0; i < boost::num_vertices(g); ++i){
85 0 : if(!visited[i] && m_sources[i]){
86 0 : if(boost::in_degree(i, g) == 1){
87 0 : visited[i] = true;
88 : adj_iter nIt, nEnd;
89 0 : boost::tie(nIt, nEnd) = boost::adjacent_vertices(i, g);
90 0 : m_sources[*nIt] = true;
91 0 : boost::clear_vertex(i, g);
92 : return i;
93 : }
94 0 : else if(boost::in_degree(i, g) == 0){
95 0 : visited[i] = true;
96 0 : return i;
97 : }
98 : else{}
99 : }
100 : }
101 :
102 : return -1u;
103 : }
104 :
105 0 : void topological_ordering(O_t& o, G_t& g){
106 : size_t n = boost::num_vertices(g);
107 0 : std::vector<BOOL> visited(n, false);
108 0 : for(unsigned i = 0; i < n; ++i){
109 0 : auto v = get_source_vertex(visited, g);
110 0 : o[v] = i;
111 : }
112 0 : }
113 :
114 0 : void compute(){
115 0 : topological_ordering(o, g);
116 :
117 0 : g = G_t(0);
118 :
119 : #ifdef UG_DEBUG
120 : check();
121 : #endif
122 0 : }
123 :
124 0 : void check(){
125 0 : if(!is_permutation(o)){
126 0 : UG_THROW(name() << "::check: Not a permutation!");
127 : }
128 0 : }
129 :
130 0 : O_t& ordering(){
131 0 : return o;
132 : }
133 :
134 0 : void init(M_t* A, const V_t& V){
135 : const GridFunc_t* pGridF;
136 : size_t numSources = 0;
137 : std::string ssName;
138 : size_t n = 0;
139 :
140 : //graph construction
141 : try{
142 0 : if((pGridF = dynamic_cast<const GridFunc_t*>(&V)) == 0){
143 0 : UG_THROW(name() << "::init: No DoFDistribution specified.");
144 : }
145 :
146 : SmartPtr<DoFDistribution> dd = ((GridFunc_t*) pGridF)->dof_distribution();
147 :
148 : n = dd->num_indices();
149 :
150 0 : if(n != A->num_rows ()){
151 0 : UG_THROW(name() << "::init: #indices != #rows");
152 : }
153 :
154 0 : m_ssIdx = pGridF->domain()->subset_handler()->get_subset_index(m_ssName);
155 :
156 0 : if(m_ssIdx < 0)
157 0 : UG_THROW(name() << "::init: Invalid subset for sources selected! Call 'select_sources(const char*)'.");
158 :
159 0 : std::vector<std::vector<size_t> > vvConnection(n);
160 : try{
161 0 : dd->get_connections<ug::Edge>(vvConnection);
162 : }
163 0 : UG_CATCH_THROW(name() << "::init: No adjacency graph available!");
164 :
165 0 : g = G_t(n);
166 :
167 0 : for(unsigned i = 0; i < n; ++i){
168 0 : for(unsigned j = 0; j < vvConnection[i].size(); ++j){
169 0 : if(i != vvConnection[i][j]){ //no self loops!
170 0 : boost::add_edge(i, vvConnection[i][j], g);
171 : }
172 : }
173 : }
174 0 : }
175 0 : catch(...){
176 0 : throw;
177 : }
178 :
179 0 : o.resize(n);
180 0 : m_sources = std::vector<BOOL>(n, false);
181 :
182 : //select source vertices according to m_ssIdx
183 : typedef typename GridFunc_t::template traits<ug::Vertex>::const_iterator ugVertIt_t;
184 :
185 : ugVertIt_t ugVertIt = pGridF->template begin<ug::Vertex>();
186 : ugVertIt_t ugVertEnd = pGridF->template end<ug::Vertex>();
187 : size_t k = 0;
188 0 : for(; ugVertIt != ugVertEnd; ++ugVertIt, ++k){
189 : ug::Vertex* v = *ugVertIt;
190 :
191 0 : int si = pGridF->domain()->subset_handler()->get_subset_index(v);
192 :
193 0 : if(si == m_ssIdx){
194 0 : m_sources[k] = true;
195 : ++numSources;
196 : }
197 : }
198 :
199 : #ifdef UG_ENABLE_DEBUG_LOGS
200 : UG_LOG("Using " << name() << " (subset " << m_ssIdx << ", " << m_ssName
201 : << ", n=" << boost::num_vertices(g) << ", m=2*" << boost::num_edges(g)/2
202 : << ", s=" << numSources << ")\n");
203 : #endif
204 0 : }
205 :
206 0 : void init(M_t*){
207 0 : UG_THROW(name() << "::init: Cannot initialize smoother without a geometry. Specify the 2nd argument for init!");
208 : }
209 :
210 0 : void init(M_t*, const V_t&, const O_t&){
211 0 : UG_THROW(name() << "::init: Algorithm does not support induced subgraph version!");
212 : }
213 :
214 0 : void init(M_t*, const O_t&){
215 0 : UG_THROW(name() << "::init: Algorithm does not support induced subgraph version!");
216 : }
217 :
218 0 : virtual const char* name() const {return "RiverOrdering";}
219 :
220 0 : void select_sources(const char* ssName){
221 : //m_ssIdx = ssIdx;
222 0 : m_ssName = ssName;
223 0 : }
224 :
225 : private:
226 : G_t g;
227 : O_t o;
228 :
229 : int m_ssIdx;
230 : const char* m_ssName;
231 : std::vector<BOOL> m_sources;
232 : };
233 :
234 : } // end namespace ug
235 :
236 : #endif
|
