Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
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__OPERATOR__LINEAR_OPERATOR__STD_INJECTION_IMPL__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__STD_INJECTION_IMPL__
35 :
36 : #include "std_injection.h"
37 :
38 : namespace ug{
39 :
40 : /**
41 : * This functions assembles the interpolation matrix between to
42 : * grid levels using only the RegularVertex degrees of freedom.
43 : *
44 : * \param[out] mat Assembled interpolation matrix that interpolates u -> v
45 : * \param[in] approxSpace Approximation Space
46 : * \param[in] coarseLevel Coarse Level index
47 : * \param[in] fineLevel Fine Level index
48 : */
49 : template <typename TAlgebra>
50 0 : void AssembleInjectionForP1Lagrange(typename TAlgebra::matrix_type& mat,
51 : const DoFDistribution& coarseDD,
52 : const DoFDistribution& fineDD)
53 : {
54 : PROFILE_FUNC_GROUP("gmg");
55 : // Allow only lagrange P1 functions
56 0 : for(size_t fct = 0; fct < fineDD.num_fct(); ++fct)
57 0 : if(fineDD.local_finite_element_id(fct).type() != LFEID::LAGRANGE ||
58 : fineDD.local_finite_element_id(fct).order() != 1)
59 0 : UG_THROW("AssembleInjectionForP1Lagrange: "
60 : "Interpolation only implemented for Lagrange P1 functions.");
61 :
62 : // get MultiGrid
63 0 : const MultiGrid& grid = *coarseDD.multi_grid();
64 :
65 : // get number of dofs on different levels
66 : const size_t numFineDoFs = fineDD.num_indices();
67 : const size_t numCoarseDoFs = coarseDD.num_indices();
68 :
69 : // resize matrix
70 0 : mat.resize_and_clear(numCoarseDoFs, numFineDoFs);
71 :
72 : std::vector<size_t> coarseInd, fineInd;
73 :
74 : // RegularVertex iterators
75 : typedef DoFDistribution::traits<RegularVertex>::const_iterator const_iterator;
76 : const_iterator iter, iterBegin, iterEnd;
77 :
78 0 : iterBegin = fineDD.template begin<RegularVertex>();
79 0 : iterEnd = fineDD.template end<RegularVertex>();
80 :
81 : // loop nodes of fine subset
82 0 : for(iter = iterBegin; iter != iterEnd; ++iter)
83 : {
84 : // get father
85 : GridObject* geomObj = grid.get_parent(*iter);
86 0 : Vertex* vert = dynamic_cast<Vertex*>(geomObj);
87 :
88 : // Check if father is RegularVertex
89 0 : if(vert != NULL)
90 : {
91 : // get global indices
92 0 : coarseDD.inner_algebra_indices(vert, coarseInd);
93 : }
94 0 : else continue;
95 :
96 : // get global indices
97 0 : fineDD.inner_algebra_indices(*iter, fineInd);
98 :
99 0 : for(size_t i = 0; i < coarseInd.size(); ++i)
100 0 : mat(coarseInd[i], fineInd[i]) = 1.0;
101 : }
102 0 : }
103 :
104 : template <int dim, typename TAlgebra>
105 0 : void AssembleInjectionByAverageOfChildren(typename TAlgebra::matrix_type& mat,
106 : const DoFDistribution& coarseDD,
107 : const DoFDistribution& fineDD)
108 : {
109 : PROFILE_FUNC_GROUP("gmg");
110 : // get MultiGrid
111 0 : const MultiGrid& grid = *coarseDD.multi_grid();
112 :
113 : std::vector<size_t> coarseInd, fineInd;
114 :
115 : // RegularVertex iterators
116 : typedef typename DoFDistribution::dim_traits<dim>::const_iterator const_iterator;
117 : typedef typename DoFDistribution::dim_traits<dim>::grid_base_object Element;
118 : const_iterator iter, iterBegin, iterEnd;
119 :
120 0 : iterBegin = coarseDD.template begin<Element>();
121 0 : iterEnd = coarseDD.template end<Element>();
122 :
123 : // loop elements of coarse subset
124 0 : for(iter = iterBegin; iter != iterEnd; ++iter)
125 : {
126 : // get element
127 : Element* coarseElem = *iter;
128 :
129 : // get children
130 : const size_t numChild = grid.num_children<Element, Element>(coarseElem);
131 0 : if(numChild == 0) continue;
132 :
133 : // get global indices
134 0 : coarseDD.inner_algebra_indices(coarseElem, coarseInd);
135 :
136 0 : for(size_t c = 0; c < numChild; ++c)
137 : {
138 : Element* child = grid.get_child<Element, Element>(coarseElem, c);
139 :
140 0 : fineDD.inner_algebra_indices(child, fineInd);
141 :
142 0 : for(size_t i = 0; i < coarseInd.size(); ++i)
143 0 : mat(coarseInd[i], fineInd[i]) = 1.0/(numChild);
144 : }
145 : }
146 0 : }
147 :
148 : template <typename TAlgebra>
149 0 : void AssembleInjectionByAverageOfChildren(typename TAlgebra::matrix_type& mat,
150 : const DoFDistribution& coarseDD,
151 : const DoFDistribution& fineDD)
152 : {
153 : // get number of dofs on different levels
154 : const size_t numFineDoFs = fineDD.num_indices();
155 : const size_t numCoarseDoFs = coarseDD.num_indices();
156 :
157 : // resize matrix
158 0 : mat.resize_and_clear(numCoarseDoFs, numFineDoFs);
159 :
160 0 : if(coarseDD.max_dofs(VERTEX)) AssembleInjectionByAverageOfChildren<0, TAlgebra>(mat, coarseDD, fineDD);
161 0 : if(coarseDD.max_dofs(EDGE)) AssembleInjectionByAverageOfChildren<1, TAlgebra>(mat, coarseDD, fineDD);
162 0 : if(coarseDD.max_dofs(FACE)) AssembleInjectionByAverageOfChildren<2, TAlgebra>(mat, coarseDD, fineDD);
163 0 : if(coarseDD.max_dofs(VOLUME)) AssembleInjectionByAverageOfChildren<3, TAlgebra>(mat, coarseDD, fineDD);
164 0 : }
165 :
166 : template <typename TDomain, typename TAlgebra>
167 : template <typename TElem>
168 0 : void StdInjection<TDomain, TAlgebra>::
169 : set_identity_on_pure_surface(matrix_type& mat,
170 : const DoFDistribution& coarseDD, const DoFDistribution& fineDD)
171 : {
172 : PROFILE_FUNC_GROUP("gmg");
173 :
174 : std::vector<size_t> vCoarseIndex, vFineIndex;
175 0 : const MultiGrid& mg = *coarseDD.multi_grid();
176 :
177 : // iterators
178 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
179 : const_iterator iter, iterBegin, iterEnd;
180 :
181 : // loop subsets on fine level
182 0 : for(int si = 0; si < coarseDD.num_subsets(); ++si)
183 : {
184 0 : iterBegin = coarseDD.template begin<TElem>(si);
185 0 : iterEnd = coarseDD.template end<TElem>(si);
186 :
187 : // loop vertices for fine level subset
188 0 : for(iter = iterBegin; iter != iterEnd; ++iter)
189 : {
190 : // get element
191 : TElem* coarseElem = *iter;
192 :
193 : const size_t numChild = mg.num_children<TElem, TElem>(coarseElem);
194 0 : if(numChild != 0) continue;
195 :
196 : // get indices
197 0 : coarseDD.inner_algebra_indices(coarseElem, vCoarseIndex);
198 0 : fineDD.inner_algebra_indices(coarseElem, vFineIndex);
199 : UG_ASSERT(vCoarseIndex.size() == vFineIndex.size(), "Size mismatch");
200 :
201 : // set identity
202 0 : for(size_t i = 0; i < vCoarseIndex.size(); ++i)
203 0 : mat(vCoarseIndex[i], vFineIndex[i]) = 1.0;
204 : }
205 : }
206 0 : }
207 :
208 : template <typename TDomain, typename TAlgebra>
209 0 : void StdInjection<TDomain, TAlgebra>::
210 : set_identity_on_pure_surface(matrix_type& mat,
211 : const DoFDistribution& coarseDD, const DoFDistribution& fineDD)
212 : {
213 0 : if(coarseDD.max_dofs(VERTEX)) set_identity_on_pure_surface<Vertex>(mat, coarseDD, fineDD);
214 0 : if(coarseDD.max_dofs(EDGE)) set_identity_on_pure_surface<Edge>(mat, coarseDD, fineDD);
215 0 : if(coarseDD.max_dofs(FACE)) set_identity_on_pure_surface<Face>(mat, coarseDD, fineDD);
216 0 : if(coarseDD.max_dofs(VOLUME)) set_identity_on_pure_surface<Volume>(mat, coarseDD, fineDD);
217 0 : }
218 :
219 : ////////////////////////////////////////////////////////////////////////////////
220 : ////////////////////////////////////////////////////////////////////////////////
221 : // StdInjection
222 : ////////////////////////////////////////////////////////////////////////////////
223 : ////////////////////////////////////////////////////////////////////////////////
224 :
225 :
226 : template <typename TDomain, typename TAlgebra>
227 : void StdInjection<TDomain, TAlgebra>::
228 : set_approximation_space(SmartPtr<ApproximationSpace<TDomain> > approxSpace)
229 : {
230 0 : m_spApproxSpace = approxSpace;
231 : }
232 :
233 : template <typename TDomain, typename TAlgebra>
234 0 : void StdInjection<TDomain, TAlgebra>::
235 : set_levels(GridLevel coarseLevel, GridLevel fineLevel)
236 : {
237 0 : m_fineLevel = fineLevel;
238 0 : m_coarseLevel = coarseLevel;
239 :
240 0 : if(m_fineLevel.level() - m_coarseLevel.level() != 1)
241 0 : UG_THROW("StdInjection::set_levels:"
242 : " Can only project between successive level.");
243 0 : }
244 :
245 : template <typename TDomain, typename TAlgebra>
246 0 : void StdInjection<TDomain, TAlgebra>::init()
247 : {
248 : PROFILE_FUNC_GROUP("gmg");
249 0 : if(!m_spApproxSpace.valid())
250 0 : UG_THROW("StdInjection::init: "
251 : "Approximation Space not set. Cannot init Projection.");
252 :
253 : // check only lagrange P1 functions
254 : bool P1LagrangeOnly = true;
255 0 : for(size_t fct = 0; fct < m_spApproxSpace->num_fct(); ++fct)
256 0 : if(m_spApproxSpace->local_finite_element_id(fct).type() != LFEID::LAGRANGE ||
257 : m_spApproxSpace->local_finite_element_id(fct).order() != 1)
258 : P1LagrangeOnly = false;
259 :
260 : try{
261 0 : if(P1LagrangeOnly)
262 : {
263 : AssembleInjectionForP1Lagrange<TAlgebra>
264 0 : (m_matrix,
265 0 : *m_spApproxSpace->dof_distribution(m_coarseLevel),
266 0 : *m_spApproxSpace->dof_distribution(m_fineLevel));
267 : }
268 : else
269 : {
270 : AssembleInjectionByAverageOfChildren<TAlgebra>
271 0 : (m_matrix,
272 0 : *m_spApproxSpace->dof_distribution(m_coarseLevel),
273 0 : *m_spApproxSpace->dof_distribution(m_fineLevel));
274 : }
275 :
276 0 : } UG_CATCH_THROW("StdInjection::init():"
277 : " Cannot assemble interpolation matrix.");
278 :
279 0 : if(m_coarseLevel.is_surface()){
280 0 : set_identity_on_pure_surface(m_matrix, *m_spApproxSpace->dof_distribution(m_coarseLevel), *m_spApproxSpace->dof_distribution(m_fineLevel));
281 : }
282 :
283 : #ifdef UG_PARALLEL
284 : m_matrix.set_storage_type(PST_CONSISTENT);
285 : #endif
286 :
287 0 : m_bInit = true;
288 0 : }
289 :
290 : template <typename TDomain, typename TAlgebra>
291 0 : void StdInjection<TDomain, TAlgebra>::
292 : prolongate(vector_type& uFine, const vector_type& uCoarse)
293 : {
294 : PROFILE_FUNC_GROUP("gmg");
295 : // Check, that operator is initiallized
296 0 : if(!m_bInit)
297 0 : UG_THROW("StdInjection::apply:"
298 : " Operator not initialized.");
299 :
300 : // Some Assertions
301 : UG_ASSERT(uFine.size() >= m_matrix.num_rows(),
302 : "Vector [size= " << uFine.size() << "] and Rows [size= "
303 : << m_matrix.num_rows() <<"] sizes have to match!");
304 : UG_ASSERT(uCoarse.size() >= m_matrix.num_cols(), "Vector [size= "
305 : << uCoarse.size() << "] and Cols [size= " <<
306 : m_matrix.num_cols() <<"] sizes have to match!");
307 :
308 : // Apply matrix
309 0 : if(!m_matrix.apply_transposed(uFine, uCoarse))
310 : UG_THROW("StdInjection::apply: Cannot apply matrix.");
311 0 : }
312 :
313 : template <typename TDomain, typename TAlgebra>
314 0 : void StdInjection<TDomain, TAlgebra>::
315 : do_restrict(vector_type& uCoarse, const vector_type& uFine)
316 : {
317 : PROFILE_FUNC_GROUP("gmg");
318 : // Check, that operator is initialized
319 0 : if(!m_bInit)
320 0 : UG_THROW("StdInjection::apply_transposed:"
321 : "Operator not initialized.");
322 :
323 : // Some Assertions
324 : UG_ASSERT(uFine.size() >= m_matrix.num_cols(),
325 : "Vector [size= " << uFine.size() << "] and Cols [size= "
326 : << m_matrix.num_cols() <<"] sizes have to match!");
327 : UG_ASSERT(uCoarse.size() >= m_matrix.num_rows(), "Vector [size= "
328 : << uCoarse.size() << "] and Rows [size= " <<
329 : m_matrix.num_rows() <<"] sizes have to match!");
330 :
331 : // Apply matrix
332 : try{
333 0 : m_matrix.apply_ignore_zero_rows(uCoarse, 1.0, uFine);
334 : }
335 : UG_CATCH_THROW("StdInjection::apply_transposed:"
336 : " Cannot apply transposed matrix.");
337 0 : }
338 :
339 :
340 : template <typename TDomain, typename TAlgebra>
341 : SmartPtr<ITransferOperator<TDomain, TAlgebra> >
342 0 : StdInjection<TDomain, TAlgebra>::clone()
343 : {
344 0 : SmartPtr<StdInjection> op(new StdInjection);
345 0 : op->set_approximation_space(m_spApproxSpace);
346 0 : return op;
347 : }
348 :
349 : } // end namespace ug
350 :
351 : #endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__STD_INJECTION_IMPL__ */
|
