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__FUNCTION_SPACE__LEVEL_TRANSFER__
34 : #define __H__UG__LIB_DISC__FUNCTION_SPACE__LEVEL_TRANSFER__
35 :
36 : #include "lib_disc/function_spaces/grid_function.h"
37 : #include "lib_disc/reference_element/reference_mapping_provider.h"
38 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
39 : #include "lib_disc/function_spaces/dof_position_util.h"
40 :
41 : namespace ug{
42 :
43 : ////////////////////////////////////////////////////////////////////////////////
44 : // Prolongate
45 : ////////////////////////////////////////////////////////////////////////////////
46 :
47 : template <typename TDomain, typename TAlgebra>
48 0 : void ProlongateP1(GridFunction<TDomain, TAlgebra>& uFine,
49 : const GridFunction<TDomain, TAlgebra>& uCoarse)
50 : {
51 : typedef GridFunction<TDomain, TAlgebra> TGridFunction;
52 : typedef typename TGridFunction::template traits<Vertex>::const_iterator const_iterator;
53 :
54 : // get subsethandler and grid
55 0 : SmartPtr<MultiGrid> mg = uFine.domain()->grid();
56 :
57 : // get top level of gridfunctions
58 0 : const int fineTopLevel = uFine.dof_distribution()->grid_level().level();
59 0 : const int coarseTopLevel = uCoarse.dof_distribution()->grid_level().level();
60 :
61 : // check
62 0 : if(fineTopLevel == GridLevel::TOP || coarseTopLevel == GridLevel::TOP)
63 0 : UG_THROW("ProlongateP1: Top Level not supported.")
64 0 : if(fineTopLevel < coarseTopLevel)
65 0 : UG_THROW("ProlongateP1: fine level must be >= coarse level.");
66 :
67 : // storage
68 : std::vector<size_t> vFineIndex, vCoarseIndex;
69 :
70 : // loop elements
71 : const_iterator iterEnd = uFine.template end<Vertex>();
72 : const_iterator iter = uFine.template begin<Vertex>();
73 0 : for(; iter != iterEnd; ++iter)
74 : {
75 : // get vertex
76 : Vertex* vrt = *iter;
77 : const int vertexLevel = mg->get_level(vrt);
78 :
79 : // a) if not on the same level as the top level of the fine grid function
80 : // and the coarse grid function is already defined on the level
81 : // we can simply copy the values, since the coarse grid function and
82 : // the fine grid function are covering the identical part here
83 0 : if(vertexLevel != fineTopLevel && vertexLevel <= coarseTopLevel)
84 : {
85 : uFine.inner_algebra_indices(vrt, vFineIndex);
86 : uCoarse.inner_algebra_indices(vrt, vCoarseIndex);
87 :
88 0 : for(size_t i = 0; i < vFineIndex.size(); ++i)
89 0 : uFine[ vFineIndex[i] ] = uCoarse[ vCoarseIndex[i] ];
90 :
91 0 : continue;
92 0 : }
93 :
94 : // get parent and level where coarse grid function is defined
95 : GridObject* parent = mg->get_parent(vrt);
96 0 : const ReferenceObjectID parentBaseObjectID = parent->reference_object_id();
97 : int parentLevel = mg->get_level(parent);
98 0 : while(parentLevel > coarseTopLevel){
99 0 : parent = mg->get_parent(parent);
100 : parentLevel = mg->get_level(parent);
101 : }
102 :
103 : // b) if the parent, where the coarse grid function is defined and the
104 : // fine element are only one level separated, we can use an optimized
105 : // interpolation. This case will always apply if the two grid functions
106 : // are only one surface level separated.
107 0 : if(parentLevel == vertexLevel - 1)
108 : {
109 : // distinguish type of parent
110 0 : switch(parentBaseObjectID)
111 : {
112 0 : case ROID_VERTEX:
113 : {
114 : Vertex* pParent = static_cast<Vertex*>(parent);
115 : uFine.inner_algebra_indices(vrt, vFineIndex);
116 : uCoarse.inner_algebra_indices(pParent, vCoarseIndex);
117 :
118 0 : for(size_t i = 0; i < vFineIndex.size(); ++i)
119 0 : uFine[ vFineIndex[i] ] = uCoarse[ vCoarseIndex[i] ];
120 : }
121 : break;
122 : case ROID_EDGE:
123 : {
124 : uFine.inner_algebra_indices(vrt, vFineIndex);
125 0 : for(size_t i = 0; i < vFineIndex.size(); ++i)
126 0 : uFine[ vFineIndex[i] ] = 0.0;
127 :
128 : Edge* pParent = static_cast<Edge*>(parent);
129 0 : for(size_t i = 0; i < pParent->num_vertices(); ++i)
130 : {
131 0 : Vertex* edgeVrt = pParent->vertex(i);
132 : uCoarse.inner_algebra_indices(edgeVrt, vCoarseIndex);
133 :
134 0 : for(size_t i = 0; i < vFineIndex.size(); ++i)
135 0 : VecScaleAdd(uFine[ vFineIndex[i] ],
136 0 : 1.0, uFine[ vFineIndex[i] ],
137 : 0.5, uCoarse[ vCoarseIndex[i] ]);
138 : }
139 : }
140 : break;
141 : case ROID_QUADRILATERAL:
142 : {
143 : uFine.inner_algebra_indices(vrt, vFineIndex);
144 0 : for(size_t i = 0; i < vFineIndex.size(); ++i)
145 0 : uFine[ vFineIndex[i] ] = 0.0;
146 :
147 : Face* pParent = static_cast<Face*>(parent);
148 0 : for(size_t i = 0; i < pParent->num_vertices(); ++i)
149 : {
150 0 : Vertex* faceVrt = pParent->vertex(i);
151 : uCoarse.inner_algebra_indices(faceVrt, vCoarseIndex);
152 :
153 0 : for(size_t i = 0; i < vFineIndex.size(); ++i)
154 0 : VecScaleAdd(uFine[ vFineIndex[i] ],
155 0 : 1.0, uFine[ vFineIndex[i] ],
156 : 0.25, uCoarse[ vCoarseIndex[i] ]);
157 : }
158 : }
159 : break;
160 : case ROID_HEXAHEDRON:
161 : {
162 : uFine.inner_algebra_indices(vrt, vFineIndex);
163 0 : for(size_t i = 0; i < vFineIndex.size(); ++i)
164 0 : uFine[ vFineIndex[i] ] = 0.0;
165 :
166 : Volume* pParent = static_cast<Volume*>(parent);
167 0 : for(size_t i = 0; i < pParent->num_vertices(); ++i)
168 : {
169 0 : Vertex* hexVrt = pParent->vertex(i);
170 : uCoarse.inner_algebra_indices(hexVrt, vCoarseIndex);
171 :
172 0 : for(size_t i = 0; i < vFineIndex.size(); ++i)
173 0 : VecScaleAdd(uFine[ vFineIndex[i] ],
174 0 : 1.0, uFine[ vFineIndex[i] ],
175 : 0.125, uCoarse[ vCoarseIndex[i] ]);
176 : }
177 : }
178 : break;
179 : case ROID_TRIANGLE:
180 : case ROID_TETRAHEDRON:
181 : case ROID_PRISM:
182 : case ROID_PYRAMID:
183 : case ROID_OCTAHEDRON: /*nothing to do in those cases */ break;
184 0 : default: UG_THROW("Unexpected case appeared.");
185 : }
186 :
187 0 : continue;
188 0 : }
189 :
190 : // c) we must interpolate the values based on the trial space
191 0 : UG_THROW("This case not implemented.");
192 : }
193 0 : }
194 :
195 :
196 :
197 : template <typename TDomain, typename TAlgebra>
198 0 : void ProlongateElemwise(GridFunction<TDomain, TAlgebra>& uFine,
199 : const GridFunction<TDomain, TAlgebra>& uCoarse)
200 : {
201 : // dimension
202 : const int dim = TDomain::dim;
203 :
204 : // get subsethandler and grid
205 0 : SmartPtr<MultiGrid> mg = uFine.domain()->grid();
206 :
207 : // get top level of gridfunctions
208 0 : const int fineTopLevel = uFine.dof_distribution()->grid_level().level();
209 0 : const int coarseTopLevel = uCoarse.dof_distribution()->grid_level().level();
210 :
211 : // check
212 0 : if(fineTopLevel == GridLevel::TOP || coarseTopLevel == GridLevel::TOP)
213 0 : UG_THROW("ProlongateElemwise: Top Level not supported.")
214 0 : if(fineTopLevel < coarseTopLevel)
215 0 : UG_THROW("ProlongateElemwise: fine level must be >= coarse level.");
216 :
217 : // storage
218 : std::vector<DoFIndex> vCoarseMI, vFineMI;
219 :
220 : // vector of local finite element ids
221 : SmartPtr<DoFDistribution> fineDD = uFine.dof_distribution();
222 0 : std::vector<LFEID> vFineLFEID(fineDD->num_fct());
223 0 : for(size_t fct = 0; fct < fineDD->num_fct(); ++fct)
224 0 : vFineLFEID[fct] = fineDD->local_finite_element_id(fct);
225 : ConstSmartPtr<DoFDistribution> coarseDD = uCoarse.dof_distribution();
226 0 : std::vector<LFEID> vCoarseLFEID(coarseDD->num_fct());
227 0 : for(size_t fct = 0; fct < coarseDD->num_fct(); ++fct)
228 0 : vCoarseLFEID[fct] = coarseDD->local_finite_element_id(fct);
229 :
230 : // check fct
231 0 : if(vFineLFEID.size() != vCoarseLFEID.size())
232 0 : UG_THROW("ProlongateElemwise: Spaces must contain same number of functions.")
233 :
234 : // get flag if all trial spaces are equal
235 : // bool bSameLFEID = true;
236 0 : for(size_t fct = 0; fct < vFineLFEID.size(); ++fct){
237 : // if(vFineLFEID[fct] != vCoarseLFEID[fct])
238 : // bSameLFEID = false;
239 :
240 0 : if (vCoarseLFEID[fct].type() == LFEID::PIECEWISE_CONSTANT ||
241 : vFineLFEID[fct].type() == LFEID::PIECEWISE_CONSTANT)
242 0 : UG_THROW("Not implemented.")
243 : }
244 :
245 : // iterators
246 : typedef typename DoFDistribution::dim_traits<dim>::const_iterator const_iterator;
247 : typedef typename DoFDistribution::dim_traits<dim>::grid_base_object Element;
248 : const_iterator iter, iterBegin, iterEnd;
249 :
250 : // loop subsets on coarse level
251 0 : for(int si = 0; si < coarseDD->num_subsets(); ++si)
252 : {
253 0 : iterBegin = coarseDD->template begin<Element>(si);
254 0 : iterEnd = coarseDD->template end<Element>(si);
255 :
256 : // loop elem for coarse level subset
257 0 : for(iter = iterBegin; iter != iterEnd; ++iter)
258 : {
259 : // get element
260 0 : Element* coarseElem = *iter;
261 :
262 : // get children where fine grid function is defined
263 : std::vector<Element*> vChild;
264 : std::queue<Element*> qElem;
265 : qElem.push(coarseElem);
266 0 : while(!qElem.empty()){
267 0 : Element* elem = qElem.front(); qElem.pop();
268 0 : if(mg->get_level(elem) == fineTopLevel || !mg->has_children(elem)){
269 0 : vChild.push_back(elem);
270 : } else {
271 0 : for(size_t c = 0; c < mg->num_children<Element,Element>(elem); ++c){
272 0 : qElem.push(mg->get_child<Element,Element>(elem, c));
273 : }
274 : }
275 : }
276 :
277 : // type of father
278 0 : const ReferenceObjectID coarseROID = coarseElem->reference_object_id();
279 :
280 : // loop all components
281 0 : for(size_t fct = 0; fct < coarseDD->num_fct(); fct++)
282 : {
283 : // check that fct defined on subset
284 0 : if(!coarseDD->is_def_in_subset(fct, si)) continue;
285 :
286 : // get global indices
287 0 : coarseDD->dof_indices(coarseElem, fct, vCoarseMI);
288 :
289 : // get local finite element trial spaces
290 : const LocalShapeFunctionSet<dim>& lsfs
291 : = LocalFiniteElementProvider::get<dim>(coarseROID, vCoarseLFEID[fct]);
292 :
293 : // get corner coordinates
294 : std::vector<MathVector<dim> > vCornerCoarse;
295 0 : CollectCornerCoordinates(vCornerCoarse, *coarseElem, *uFine.domain());
296 :
297 : // loop children
298 0 : for(size_t c = 0; c < vChild.size(); ++c)
299 : {
300 0 : Element* child = vChild[c];
301 :
302 : // fine dof indices
303 0 : fineDD->dof_indices(child, fct, vFineMI);
304 :
305 : // global positions of fine dofs
306 : std::vector<MathVector<dim> > vDoFPos, vLocPos;
307 0 : DoFPosition(vDoFPos, child, *uFine.domain(), vFineLFEID[fct]);
308 :
309 : UG_ASSERT(vDoFPos.size() == vFineMI.size(), "numDoFPos ("
310 : <<vDoFPos.size()<<") != numDoFs ("<<vFineMI.size()<<").");
311 :
312 : // get Reference Mapping
313 : DimReferenceMapping<dim, dim>& map
314 : = ReferenceMappingProvider::get<dim, dim>(coarseROID, vCornerCoarse);
315 :
316 :
317 : // get local position of DoF
318 0 : vLocPos.resize(vDoFPos.size());
319 0 : for(size_t ip = 0; ip < vLocPos.size(); ++ip) VecSet(vLocPos[ip], 0.0);
320 0 : map.global_to_local(vLocPos, vDoFPos);
321 :
322 : // get all shape functions
323 : std::vector<std::vector<number> > vvShape;
324 :
325 : // evaluate coarse shape fct at fine local point
326 0 : lsfs.shapes(vvShape, vLocPos);
327 :
328 0 : for(size_t ip = 0; ip < vvShape.size(); ++ip){
329 0 : DoFRef(uFine, vFineMI[ip]) = 0.0;
330 0 : for(size_t sh = 0; sh < vvShape[ip].size(); ++sh){
331 0 : DoFRef(uFine, vFineMI[ip]) +=
332 0 : vvShape[ip][sh] * DoFRef(uCoarse, vCoarseMI[sh]);
333 : }
334 : }
335 : }
336 : }
337 : }
338 : }
339 0 : }
340 :
341 :
342 : template <typename TDomain, typename TAlgebra>
343 0 : void Prolongate(GridFunction<TDomain, TAlgebra>& uFine,
344 : const GridFunction<TDomain, TAlgebra>& uCoarse)
345 : {
346 : // grid functions must be from same Domain
347 0 : if(uFine.domain() != uCoarse.domain())
348 0 : UG_THROW("Prolongate: GridFunctions must have same Domain.");
349 :
350 : // grid functions must have same function pattern
351 0 : if(uFine.function_pattern().get() != uCoarse.function_pattern().get())
352 0 : UG_THROW("Prolongate: GridFunctions must have same Function Pattern.");
353 :
354 : // get grid levels
355 0 : const int fineTopLevel = uFine.dof_distribution()->grid_level().level();
356 0 : const int coarseTopLevel = uCoarse.dof_distribution()->grid_level().level();
357 0 : if(fineTopLevel == GridLevel::TOP || coarseTopLevel == GridLevel::TOP)
358 0 : UG_THROW("Prolongate: Top Level not supported.")
359 0 : if(fineTopLevel < coarseTopLevel)
360 0 : UG_THROW("Prolongate: fine level must be >= coarse level.");
361 :
362 : // loop functions
363 : bool bOnlyP1Fct = true;
364 0 : for(size_t fct = 0; fct < uFine.num_fct(); ++fct)
365 0 : if(uFine.local_finite_element_id(fct).type() != LFEID::LAGRANGE ||
366 : uFine.local_finite_element_id(fct).order() != 1)
367 : {
368 : bOnlyP1Fct = false; break;
369 : }
370 :
371 0 : if(bOnlyP1Fct &&
372 0 : (fineTopLevel == coarseTopLevel+1 || fineTopLevel == coarseTopLevel)){
373 0 : ProlongateP1(uFine, uCoarse);
374 : }
375 : else{
376 0 : ProlongateElemwise(uFine, uCoarse);
377 : }
378 :
379 : #ifdef UG_PARALLEL
380 : uFine.set_storage_type(uCoarse.get_storage_mask());
381 : #endif
382 0 : }
383 :
384 : ////////////////////////////////////////////////////////////////////////////////
385 : // Restrict
386 : ////////////////////////////////////////////////////////////////////////////////
387 :
388 :
389 : template <typename TDomain, typename TAlgebra>
390 0 : void RestrictP1(GridFunction<TDomain, TAlgebra>& uCoarse,
391 : const GridFunction<TDomain, TAlgebra>& uFine)
392 : {
393 : typedef GridFunction<TDomain, TAlgebra> TGridFunction;
394 : typedef typename TGridFunction::template traits<Vertex>::const_iterator const_iterator;
395 :
396 : // get subsethandler and grid
397 0 : SmartPtr<MultiGrid> mg = uCoarse.domain()->grid();
398 :
399 : // get top level of gridfunctions
400 0 : const int fineTopLevel = uFine.dof_distribution()->grid_level().level();
401 0 : const int coarseTopLevel = uCoarse.dof_distribution()->grid_level().level();
402 :
403 : // check
404 0 : if(fineTopLevel == GridLevel::TOP || coarseTopLevel == GridLevel::TOP)
405 0 : UG_THROW("RestrictP1: Top Level not supported.")
406 0 : if(fineTopLevel < coarseTopLevel)
407 0 : UG_THROW("RestrictP1: fine level must be >= coarse level.");
408 :
409 : // storage
410 : std::vector<size_t> vFineIndex, vCoarseIndex;
411 :
412 : // loop elements
413 : const_iterator iterEnd = uCoarse.template end<Vertex>();
414 : const_iterator iter = uCoarse.template begin<Vertex>();
415 0 : for(; iter != iterEnd; ++iter)
416 : {
417 : // get vertex
418 : Vertex* coarseVrt = *iter;
419 :
420 : // get children where fine grid function is defined
421 : Vertex* fineVrt = coarseVrt;
422 0 : while(mg->get_level(fineVrt) != fineTopLevel &&
423 : mg->has_children(fineVrt)){
424 : fineVrt = mg->get_child<Vertex,Vertex>(fineVrt, 0);
425 : }
426 :
427 : // copy values
428 : uFine.inner_algebra_indices(fineVrt, vFineIndex);
429 : uCoarse.inner_algebra_indices(coarseVrt, vCoarseIndex);
430 :
431 0 : for(size_t i = 0; i < vFineIndex.size(); ++i)
432 0 : uCoarse[ vCoarseIndex[i] ] = uFine[ vFineIndex[i] ];
433 : }
434 0 : }
435 :
436 :
437 :
438 : template <typename TDomain, typename TAlgebra>
439 : void RestrictElemwise(GridFunction<TDomain, TAlgebra>& uCoarse,
440 : const GridFunction<TDomain, TAlgebra>& uFine)
441 : {
442 : // dimension
443 : const int dim = TDomain::dim;
444 : const int locDim = TDomain::dim;
445 :
446 : // get subsethandler and grid
447 : SmartPtr<MultiGrid> mg = uCoarse.domain()->grid();
448 :
449 : // get top level of gridfunctions
450 : const int fineTopLevel = uFine.dof_distribution()->grid_level().level();
451 : const int coarseTopLevel = uCoarse.dof_distribution()->grid_level().level();
452 :
453 : // check
454 : if(fineTopLevel == GridLevel::TOP || coarseTopLevel == GridLevel::TOP)
455 : UG_THROW("RestrictElemwise: Top Level not supported.")
456 : if(fineTopLevel < coarseTopLevel)
457 : UG_THROW("RestrictElemwise: fine level must be >= coarse level.");
458 :
459 : // storage
460 : std::vector<DoFIndex> vCoarseMI, vFineMI;
461 :
462 : // vector of local finite element ids
463 : ConstSmartPtr<DoFDistribution> fineDD = uFine.dof_distribution();
464 : std::vector<LFEID> vFineLFEID(fineDD->num_fct());
465 : for(size_t fct = 0; fct < fineDD->num_fct(); ++fct)
466 : vFineLFEID[fct] = fineDD->local_finite_element_id(fct);
467 : SmartPtr<DoFDistribution> coarseDD = uCoarse.dof_distribution();
468 : std::vector<LFEID> vCoarseLFEID(coarseDD->num_fct());
469 : for(size_t fct = 0; fct < coarseDD->num_fct(); ++fct)
470 : vCoarseLFEID[fct] = coarseDD->local_finite_element_id(fct);
471 :
472 : // check fct
473 : if(vFineLFEID.size() != vCoarseLFEID.size())
474 : UG_THROW("RestrictElemwise: Spaces must contain same number of functions.")
475 :
476 : // get flag if all trial spaces are equal
477 : // bool bSameLFEID = true;
478 : for(size_t fct = 0; fct < vFineLFEID.size(); ++fct){
479 : // if(vFineLFEID[fct] != vCoarseLFEID[fct])
480 : // bSameLFEID = false;
481 :
482 : if (vCoarseLFEID[fct].type() == LFEID::PIECEWISE_CONSTANT ||
483 : vFineLFEID[fct].type() == LFEID::PIECEWISE_CONSTANT)
484 : UG_THROW("Not implemented.")
485 : }
486 :
487 : // iterators
488 : typedef typename DoFDistribution::dim_traits<locDim>::const_iterator const_iterator;
489 : typedef typename DoFDistribution::dim_traits<locDim>::grid_base_object Element;
490 : const_iterator iter, iterBegin, iterEnd;
491 :
492 : // loop subsets on coarse level
493 : for(int si = 0; si < coarseDD->num_subsets(); ++si)
494 : {
495 : iterBegin = coarseDD->template begin<Element>(si);
496 : iterEnd = coarseDD->template end<Element>(si);
497 :
498 : // loop elem for coarse level subset
499 : for(iter = iterBegin; iter != iterEnd; ++iter)
500 : {
501 : // get element
502 : Element* coarseElem = *iter;
503 :
504 : // get children where fine grid function is defined
505 : std::vector<Element*> vFineElem;
506 : std::queue<Element*> qElem;
507 : qElem.push(coarseElem);
508 : while(!qElem.empty()){
509 : Element* elem = qElem.front(); qElem.pop();
510 : if(mg->get_level(elem) == fineTopLevel || !mg->has_children(elem)){
511 : vFineElem.push_back(elem);
512 : } else {
513 : for(size_t c = 0; c < mg->num_children<Element,Element>(elem); ++c){
514 : qElem.push(mg->get_child<Element,Element>(elem, c));
515 : }
516 : }
517 : }
518 :
519 : // loop all components
520 : for(size_t fct = 0; fct < coarseDD->num_fct(); fct++)
521 : {
522 : // check that fct defined on subset
523 : if(!coarseDD->is_def_in_subset(fct, si)) continue;
524 :
525 : // get global indices
526 : coarseDD->inner_dof_indices(coarseElem, fct, vCoarseMI);
527 :
528 : // global positions of fine dofs
529 : std::vector<MathVector<dim> > vDoFPos;
530 : InnerDoFPosition(vDoFPos, coarseElem, *uCoarse.domain(), vCoarseLFEID[fct]);
531 :
532 : // loop dof points
533 : for(size_t ip = 0; ip < vDoFPos.size(); ++ip)
534 : {
535 : // loop children
536 : for(size_t c = 0; c < vFineElem.size(); ++c)
537 : {
538 : Element* fineElem = vFineElem[c];
539 :
540 : UG_THROW("This part does not work.");
541 : /* if(!ContainsPoint(fineElem, vDoFPos[ip], uFine.domain()->position_accessor())){
542 : if(c == vFineElem.size()-1)
543 : UG_THROW("Restrict: Cannot find child containing dof.");
544 : continue;
545 : }
546 : */
547 : // get corner coordinates
548 : std::vector<MathVector<dim> > vCornerFine;
549 : CollectCornerCoordinates(vCornerFine, *fineElem, *uFine.domain());
550 :
551 : // type of child
552 : const ReferenceObjectID fineROID = fineElem->reference_object_id();
553 :
554 : // get local finite element trial spaces
555 : const LocalShapeFunctionSet<locDim>& lsfs
556 : = LocalFiniteElementProvider::get<locDim>(fineROID, vFineLFEID[fct]);
557 :
558 : // get Reference Mapping
559 : DimReferenceMapping<locDim, dim>& map
560 : = ReferenceMappingProvider::get<locDim, dim>(fineROID, vCornerFine);
561 :
562 : // get local position of DoF
563 : MathVector<locDim> vLocPos;
564 : VecSet(vLocPos, 0.0);
565 : map.global_to_local(vLocPos, vDoFPos[ip]);
566 :
567 : // fine dof indices
568 : fineDD->dof_indices(fineElem, fct, vFineMI);
569 :
570 : // get all shape functions
571 : std::vector<number> vShape;
572 :
573 : // evaluate coarse shape fct at fine local point
574 : lsfs.shapes(vShape, vLocPos);
575 :
576 : // interpolate
577 : DoFRef(uCoarse, vCoarseMI[ip]) = 0.0;
578 : for(size_t sh = 0; sh < vShape.size(); ++sh)
579 : {
580 : DoFRef(uCoarse, vCoarseMI[ip]) +=
581 : vShape[sh] * DoFRef(uFine, vFineMI[sh]);
582 : }
583 : }
584 : }
585 : }
586 : }
587 : }
588 : }
589 :
590 :
591 : template <typename TDomain, typename TAlgebra>
592 0 : void Restrict(GridFunction<TDomain, TAlgebra>& uCoarse,
593 : const GridFunction<TDomain, TAlgebra>& uFine)
594 : {
595 : // grid functions must be from same Domain
596 0 : if(uCoarse.domain() != uFine.domain())
597 0 : UG_THROW("Restrict: GridFunctions must have same Domain.");
598 :
599 : // grid functions must have same function pattern
600 0 : if(uCoarse.function_pattern().get() != uFine.function_pattern().get())
601 0 : UG_THROW("Restrict: GridFunctions must have same Function Pattern.");
602 :
603 : // get grid levels
604 0 : const int coarseTopLevel = uCoarse.dof_distribution()->grid_level().level();
605 0 : const int fineTopLevel = uFine.dof_distribution()->grid_level().level();
606 0 : if(coarseTopLevel == GridLevel::TOP || fineTopLevel == GridLevel::TOP)
607 0 : UG_THROW("Restrict: Top Level not supported.")
608 0 : if(coarseTopLevel > fineTopLevel)
609 0 : UG_THROW("Restrict: fine level must be >= coarse level.");
610 :
611 : // loop functions
612 : bool bOnlyP1Fct = true;
613 0 : for(size_t fct = 0; fct < uCoarse.num_fct(); ++fct)
614 0 : if(uCoarse.local_finite_element_id(fct).type() != LFEID::LAGRANGE ||
615 : uCoarse.local_finite_element_id(fct).order() != 1)
616 : {
617 : bOnlyP1Fct = false; break;
618 : }
619 :
620 0 : if(bOnlyP1Fct &&
621 0 : (coarseTopLevel+1 == fineTopLevel || coarseTopLevel == fineTopLevel)){
622 0 : RestrictP1(uCoarse, uFine);
623 : }
624 : else{
625 0 : UG_THROW("Restrict: Only P1 implemented.")
626 : RestrictElemwise(uCoarse, uFine);
627 : }
628 :
629 : #ifdef UG_PARALLEL
630 : uCoarse.set_storage_type(uFine.get_storage_mask());
631 : #endif
632 0 : }
633 :
634 : } // end namespace ug
635 :
636 : #endif /* __H__UG__LIB_DISC__FUNCTION_SPACE__LEVEL_TRANSFER__ */
|
