Line data Source code
1 : /*
2 : * Copyright (c) 2010-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_SPACES__INTERPOLATE__
34 : #define __H__UG__LIB_DISC__FUNCTION_SPACES__INTERPOLATE__
35 :
36 : #include "grid_function_global_user_data.h"
37 : #include "common/common.h"
38 : #include "common/util/smart_pointer.h"
39 :
40 : #include "lib_grid/tools/subset_group.h"
41 :
42 : #include "lib_disc/domain_util.h"
43 : #include "lib_disc/common/groups_util.h"
44 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
45 : #include "lib_disc/spatial_disc/user_data/const_user_data.h"
46 : #include "lib_disc/reference_element/reference_mapping.h"
47 :
48 : #ifdef UG_FOR_LUA
49 : #include "bindings/lua/lua_user_data.h"
50 : #endif
51 :
52 : namespace ug{
53 :
54 : ////////////////////////////////////////////////////////////////////////////////
55 : // Interpolate on Vertices only
56 : ////////////////////////////////////////////////////////////////////////////////
57 :
58 : /**
59 : * This function interpolates a grid function on a vertex loop. Thus, it can only
60 : * be used if all degrees of freedom are located in the vertices only (e.g. P1
61 : * finite elements). In those cases it is faster than the element by element
62 : * loop.
63 : *
64 : * @param[in] spInterpolFunction data providing interpolation values
65 : * @param[out] spGridFct interpolated grid function
66 : * @param[in] fct symbolic name of function component
67 : * @param[in] ssGrp subsets, where to interpolate
68 : * @param[in] time time point
69 : * @param[in] diff_pos different vector between Interpolatin values and spGridFct.
70 : */
71 :
72 : template <typename TGridFunction>
73 0 : void InterpolateOnDiffVertices(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
74 : SmartPtr<TGridFunction> spGridFct,
75 : size_t fct,
76 : number time,
77 : const SubsetGroup& ssGrp,
78 : const MathVector<TGridFunction::dim> diff_pos)
79 : {
80 :
81 : // domain type and position_type
82 : typedef typename TGridFunction::domain_type domain_type;
83 : typedef typename domain_type::position_type position_type;
84 :
85 : //std::cout<<"Interpolate diff_vector: "<<diff_pos;
86 :
87 : // get position accessor (of interpolated grid function)
88 : const typename domain_type::position_accessor_type& aaPos
89 0 : = spGridFct->domain()->position_accessor();
90 :
91 :
92 : std::vector<DoFIndex> ind;
93 : typename TGridFunction::template dim_traits<0>::const_iterator iterEnd, iter;
94 :
95 0 : for(size_t i = 0; i < ssGrp.size(); ++i)
96 : {
97 : // get subset index
98 : const int si = ssGrp[i];
99 :
100 : // skip if function is not defined in subset
101 0 : if(!spGridFct->is_def_in_subset(fct, si)) continue;
102 :
103 : // iterate over all elements
104 0 : iterEnd = spGridFct->template end<Vertex>(si);
105 0 : iter = spGridFct->template begin<Vertex>(si);
106 0 : for(; iter != iterEnd; ++iter)
107 : {
108 : // get element
109 : Vertex* vrt = *iter;
110 :
111 : // global position
112 : position_type glob_pos = aaPos[vrt]; // position (of interpolated grid function)
113 : position_type rel_pos=glob_pos;
114 : rel_pos -=diff_pos;
115 :
116 : // value at position
117 : number val;
118 0 : (*spInterpolFunction)(val, rel_pos, time, si);
119 :
120 : // get multiindices of element
121 : spGridFct->dof_indices(vrt, fct, ind);
122 :
123 : // loop all dofs
124 0 : for(size_t i = 0; i < ind.size(); ++i)
125 : {
126 : // set value
127 0 : DoFRef(*spGridFct, ind[i]) = val;
128 : }
129 : }
130 : }
131 0 : }
132 : //getting value of spInterpolFunction at position
133 :
134 :
135 : template <typename TGridFunction>
136 : number get_number_on_coords(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
137 : typename TGridFunction::domain_type::position_type pos,
138 : number time,
139 : const int si
140 : ){
141 : number val;
142 : (*spInterpolFunction)(val, pos, time, si);
143 :
144 : return val;
145 : }
146 :
147 : /**
148 : * This function interpolates a grid function on a vertex loop. Thus, it can only
149 : * be used if all degrees of freedom are located in the vertices only (e.g. P1
150 : * finite elements). In those cases it is faster than the element by element
151 : * loop.
152 : *
153 : * @param[in] spInterpolFunction data providing interpolation values
154 : * @param[out] spGridFct interpolated grid function
155 : * @param[in] fct symbolic name of function component
156 : * @param[in] ssGrp subsets, where to interpolate
157 : * @param[in] time time point
158 : */
159 : template <typename TGridFunction>
160 0 : void InterpolateOnVertices(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
161 : SmartPtr<TGridFunction> spGridFct,
162 : size_t fct,
163 : number time,
164 : const SubsetGroup& ssGrp)
165 : {
166 : // dimension of reference element
167 : const int dim = TGridFunction::dim;
168 :
169 : MathVector<dim> diff_pos(0.0);
170 0 : InterpolateOnDiffVertices<TGridFunction>(spInterpolFunction, spGridFct, fct, time, ssGrp, diff_pos);
171 0 : }
172 :
173 :
174 :
175 : ////////////////////////////////////////////////////////////////////////////////
176 : // Interpolate on Elements
177 : ////////////////////////////////////////////////////////////////////////////////
178 :
179 : /**
180 : * This function interpolates a grid function on an element by element loop. On
181 : * each element the all associated (up to the boundary of the element) are
182 : * interpolated and the values are stored in the grid function.
183 : *
184 : * @param[in] spInterpolFunction data providing interpolation values
185 : * @param[out] spGridFct interpolated grid function
186 : * @param[in] fct symbolic name of function component
187 : * @param[in] si subset, where to interpolate
188 : * @param[in] time time point
189 : */
190 : template <typename TElem, typename TGridFunction>
191 0 : void InterpolateOnDiffElements(
192 : SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
193 : SmartPtr<TGridFunction> spGridFct, size_t fct, int si, number time,
194 : const MathVector<TGridFunction::dim> diff_pos)
195 : {
196 : // get reference element type
197 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
198 : const ReferenceObjectID roid = ref_elem_type::REFERENCE_OBJECT_ID;
199 :
200 : // dimension of reference element
201 : const int dim = ref_elem_type::dim;
202 :
203 :
204 : // domain type and position_type
205 : typedef typename TGridFunction::domain_type domain_type;
206 : typedef typename domain_type::position_type position_type;
207 :
208 : // get iterators
209 : typename TGridFunction::template traits<TElem>::const_iterator iterEnd, iter;
210 0 : iterEnd = spGridFct->template end<TElem>(si);
211 0 : iter = spGridFct->template begin<TElem>(si);
212 :
213 : // check if something to do:
214 0 : if(iter == iterEnd) return;
215 :
216 : // id of shape functions used
217 0 : LFEID id = spGridFct->local_finite_element_id(fct);
218 :
219 : // get trial space
220 : const LocalShapeFunctionSet<dim>& trialSpace =
221 : LocalFiniteElementProvider::get<dim>(roid, id);
222 :
223 : // number of dofs on element
224 0 : const size_t nsh = trialSpace.num_sh();
225 :
226 : // load local positions of dofs for the trial space on element
227 0 : std::vector<MathVector<dim> > loc_pos(nsh);
228 0 : for(size_t i = 0; i < nsh; ++i)
229 0 : if(!trialSpace.position(i, loc_pos[i]))
230 0 : UG_THROW("InterpolateOnElem: Cannot find meaningful"
231 : " local positions of dofs.");
232 :
233 : // create a reference mapping
234 : ReferenceMapping<ref_elem_type, domain_type::dim> mapping;
235 :
236 : // iterate over all elements
237 0 : for( ; iter != iterEnd; ++iter)
238 : {
239 : // get element
240 : TElem* elem = *iter;
241 :
242 : // get all corner coordinates
243 : std::vector<position_type> vCorner;
244 0 : CollectCornerCoordinates(vCorner, *elem, *spGridFct->domain());
245 :
246 : // update the reference mapping for the corners
247 0 : mapping.update(&vCorner[0]);
248 :
249 : // get multiindices of element
250 : std::vector<DoFIndex> ind;
251 : spGridFct->dof_indices(elem, fct, ind);
252 :
253 : // check multi indices
254 0 : if(ind.size() != nsh)
255 0 : UG_THROW("InterpolateOnElem: On subset "<<si<<": Number of shapes is "
256 : <<nsh<<", but got "<<ind.size()<<" multi indices.");
257 :
258 : // loop all dofs
259 0 : for(size_t i = 0; i < nsh; ++i)
260 : {
261 : // global position
262 : position_type glob_pos;
263 :
264 :
265 : // map local dof position to global position
266 0 : mapping.local_to_global(glob_pos, loc_pos[i]);
267 :
268 : position_type rel_pos=glob_pos;
269 : rel_pos -=diff_pos;
270 :
271 : // value at position
272 : number val;
273 0 : (*spInterpolFunction)(val, rel_pos, time, si);
274 :
275 : // set value
276 0 : DoFRef(*spGridFct, ind[i]) = val;
277 : }
278 : }
279 0 : }
280 :
281 : /**
282 : * This function interpolates a grid function on an element by element loop. On
283 : * each element the all associated (up to the boundary of the element) are
284 : * interpolated and the values are stored in the grid function.
285 : *
286 : * @param[in] spInterpolFunction data providing interpolation values
287 : * @param[out] spGridFct interpolated grid function
288 : * @param[in] fct symbolic name of function component
289 : * @param[in] ssGrp subsets, where to interpolate
290 : * @param[in] time time point
291 : */
292 : template <typename TGridFunction>
293 0 : void InterpolateOnDiffElements(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
294 : SmartPtr<TGridFunction> spGridFct,
295 : size_t fct,
296 : number time,
297 : const SubsetGroup& ssGrp,const MathVector<TGridFunction::dim> diff_pos)
298 : {
299 : // loop subsets
300 0 : for(size_t i = 0; i < ssGrp.size(); ++i)
301 : {
302 : // get subset index
303 : const int si = ssGrp[i];
304 :
305 : // skip if function is not defined in subset
306 0 : if(!spGridFct->is_def_in_subset(fct, si)) continue;
307 :
308 : // switch dimensions
309 : try
310 : {
311 0 : const int dim = ssGrp.dim(i);
312 :
313 0 : if(dim > TGridFunction::dim)
314 0 : UG_THROW("InterpolateOnElements: Dimension of subset is "<<dim<<", but "
315 : " World Dimension is "<<TGridFunction::dim<<". Cannot interpolate this.");
316 :
317 : // FIXME (at least for Lagrange, order > 1, parallel)
318 : // In a parallel scenario, the distribution CAN cause elements of of lower
319 : // dimension than the rest of their subset to be located disconnected from
320 : // the rest of the subset on a processor. For example, in 2D, think of a
321 : // (1D) boundary subset and a distribution where the boundary of a proc's
322 : // domain only touches the boundary subset in a vertex, but intersects with
323 : // the boundary subset in another place.
324 : // This vertex will not be considered during interpolation even though it
325 : // should be!
326 0 : switch(dim)
327 : {
328 : case DIM_SUBSET_EMPTY_GRID: break;
329 : case 0: /* \TODO: do nothing may be wrong */ break;
330 : case 1:
331 0 : InterpolateOnDiffElements<RegularEdge, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time,diff_pos);
332 0 : break;
333 : case 2:
334 0 : InterpolateOnDiffElements<Triangle, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time, diff_pos);
335 0 : InterpolateOnDiffElements<Quadrilateral, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time, diff_pos);
336 0 : break;
337 : case 3:
338 0 : InterpolateOnDiffElements<Tetrahedron, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time, diff_pos);
339 0 : InterpolateOnDiffElements<Hexahedron, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time, diff_pos);
340 0 : InterpolateOnDiffElements<Prism, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time, diff_pos);
341 0 : InterpolateOnDiffElements<Pyramid, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time, diff_pos);
342 0 : InterpolateOnDiffElements<Octahedron, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time, diff_pos);
343 0 : break;
344 0 : default: UG_THROW("InterpolateOnElements: Dimension " <<dim<<
345 : " not possible for world dim "<<3<<".");
346 : }
347 : }
348 0 : UG_CATCH_THROW("InterpolateOnElements: Cannot interpolate on elements.");
349 : }
350 0 : }
351 :
352 : /**
353 : * This function interpolates a grid function on an element by element loop. On
354 : * each element the all associated (up to the boundary of the element) are
355 : * interpolated and the values are stored in the grid function.
356 : *
357 : * @param[in] spInterpolFunction data providing interpolation values
358 : * @param[out] spGridFct interpolated grid function
359 : * @param[in] fct symbolic name of function component
360 : * @param[in] si subset, where to interpolate
361 : * @param[in] time time point
362 : */
363 : template <typename TElem, typename TGridFunction>
364 : void InterpolateOnElements(
365 : SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
366 : SmartPtr<TGridFunction> spGridFct, size_t fct, int si, number time)
367 : {
368 : // dimension of reference element
369 : const int dim = TGridFunction::dim;
370 :
371 : MathVector<dim> diff_pos(0.0);
372 : InterpolateOnDiffElements<TElem,TGridFunction>(spInterpolFunction, spGridFct, fct, si,time, diff_pos);
373 : }
374 :
375 : /**
376 : * This function interpolates a grid function on an element by element loop. On
377 : * each element the all associated (up to the boundary of the element) are
378 : * interpolated and the values are stored in the grid function.
379 : *
380 : * @param[in] spInterpolFunction data providing interpolation values
381 : * @param[out] spGridFct interpolated grid function
382 : * @param[in] fct symbolic name of function component
383 : * @param[in] ssGrp subsets, where to interpolate
384 : * @param[in] time time point
385 : */
386 : template <typename TGridFunction>
387 : void InterpolateOnElements(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
388 : SmartPtr<TGridFunction> spGridFct,
389 : size_t fct,
390 : number time,
391 : const SubsetGroup& ssGrp)
392 : {
393 : // loop subsets
394 : for(size_t i = 0; i < ssGrp.size(); ++i)
395 : {
396 : // get subset index
397 : const int si = ssGrp[i];
398 :
399 : // skip if function is not defined in subset
400 : if(!spGridFct->is_def_in_subset(fct, si)) continue;
401 :
402 : // switch dimensions
403 : try
404 : {
405 : const int dim = ssGrp.dim(i);
406 :
407 : if(dim > TGridFunction::dim)
408 : UG_THROW("InterpolateOnElements: Dimension of subset is "<<dim<<", but "
409 : " World Dimension is "<<TGridFunction::dim<<". Cannot interpolate this.");
410 :
411 : // FIXME (at least for Lagrange, order > 1, parallel)
412 : // In a parallel scenario, the distribution CAN cause elements of of lower
413 : // dimension than the rest of their subset to be located disconnected from
414 : // the rest of the subset on a processor. For example, in 2D, think of a
415 : // (1D) boundary subset and a distribution where the boundary of a proc's
416 : // domain only touches the boundary subset in a vertex, but intersects with
417 : // the boundary subset in another place.
418 : // This vertex will not be considered during interpolation even though it
419 : // should be!
420 : switch(dim)
421 : {
422 : case DIM_SUBSET_EMPTY_GRID: break;
423 : case 0: /* \TODO: do nothing may be wrong */ break;
424 : case 1:
425 : InterpolateOnElements<RegularEdge, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time);
426 : break;
427 : case 2:
428 : InterpolateOnElements<Triangle, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time);
429 : InterpolateOnElements<Quadrilateral, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time);
430 : break;
431 : case 3:
432 : InterpolateOnElements<Tetrahedron, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time);
433 : InterpolateOnElements<Hexahedron, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time);
434 : InterpolateOnElements<Prism, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time);
435 : InterpolateOnElements<Pyramid, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time);
436 : InterpolateOnElements<Octahedron, TGridFunction>(spInterpolFunction, spGridFct, fct, si, time);
437 : break;
438 : default: UG_THROW("InterpolateOnElements: Dimension " <<dim<<
439 : " not possible for world dim "<<3<<".");
440 : }
441 : }
442 : UG_CATCH_THROW("InterpolateOnElements: Cannot interpolate on elements.");
443 : }
444 : }
445 :
446 : ////////////////////////////////////////////////////////////////////////////////
447 : // Interpolate routine
448 : ////////////////////////////////////////////////////////////////////////////////
449 :
450 : /// interpolates a function on a subset
451 : /**
452 : * This function interpolates a GridFunction. To evaluate the function on every
453 : * point a functor must be passed.
454 : *
455 : * @param[in] spInterpolFunction data providing interpolation values
456 : * @param[out] spGridFct interpolated grid function
457 : * @param[in] cmp symbolic name of function component
458 : * @param[in] subsets subsets, where to interpolate (NULL = everywhere)
459 : * @param[in] time time point
460 : */
461 : template <typename TGridFunction>
462 0 : void Interpolate(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
463 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
464 : const char* subsets, number time)
465 : {
466 : // dimension of reference element
467 : const int dim = TGridFunction::dim;
468 : MathVector<dim> diff_pos(0.0);
469 0 : Interpolate(spInterpolFunction, spGridFct, cmp, subsets, time, diff_pos);
470 0 : }
471 :
472 : /// interpolates a function on a subset
473 : /**
474 : * This function interpolates a GridFunction. To evaluate the function on every
475 : * point a functor must be passed.
476 : *
477 : * @param[in] spInterpolFunction data providing interpolation values
478 : * @param[out] spGridFct interpolated grid function
479 : * @param[in] fct id of function component
480 : * @param[in] subsets subsets, where to interpolate (NULL = everywhere)
481 : * @param[in] time time point
482 : */
483 : template <typename TGridFunction>
484 0 : void Interpolate(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
485 : SmartPtr<TGridFunction> spGridFct, const size_t fct,
486 : const SubsetGroup& ssGrp, number time, const MathVector<TGridFunction::dim> diff_pos)
487 : {
488 :
489 : // check, that values do not depend on a solution
490 0 : if(spInterpolFunction->requires_grid_fct())
491 0 : UG_THROW("Interpolate: The interpolation values depend on a grid function."
492 : " This is not allowed in the current implementation. Use constant,"
493 : " lua-callback or vrl-callback user data only (even within linkers).");
494 :
495 : // check if fast P1 interpolation can be used
496 : // \TODO: This should be improved. Manifold admissible if space continuous
497 : bool bUseP1Interpolation = false;
498 0 : if(spGridFct->local_finite_element_id(fct).type() == LFEID::LAGRANGE &&
499 : spGridFct->local_finite_element_id(fct).order() == 1)
500 : bUseP1Interpolation = true;
501 :
502 : //forward
503 : if(bUseP1Interpolation){
504 0 : InterpolateOnDiffVertices<TGridFunction>(spInterpolFunction, spGridFct, fct, time, ssGrp, diff_pos);
505 : }else{
506 0 : InterpolateOnDiffElements<TGridFunction>(spInterpolFunction, spGridFct, fct, time, ssGrp, diff_pos);
507 : }
508 : // adjust parallel storage state
509 : #ifdef UG_PARALLEL
510 : spGridFct->set_storage_type(PST_CONSISTENT);
511 : #endif
512 0 : }
513 :
514 : /// interpolates a function on a subset
515 : /**
516 : * This function interpolates a GridFunction. To evaluate the function on every
517 : * point a functor must be passed.
518 : *
519 : * @param[in] spInterpolFunction data providing interpolation values
520 : * @param[out] spGridFct interpolated grid function
521 : * @param[in] cmp symbolic name of function component
522 : * @param[in] subsets subsets, where to interpolate (NULL = everywhere)
523 : * @param[in] time time point
524 : */
525 : template <typename TGridFunction>
526 0 : void Interpolate(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
527 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
528 : const char* subsets, number time, const MathVector<TGridFunction::dim> diff_pos)
529 : {
530 :
531 :
532 : // get function id of name
533 : const size_t fct = spGridFct->fct_id_by_name(cmp);
534 :
535 : // check that function found
536 0 : if(fct > spGridFct->num_fct())
537 0 : UG_THROW("Interpolate: Name of component '"<<cmp<<"' not found.");
538 :
539 0 : Interpolate(spInterpolFunction, spGridFct, fct, subsets, time, diff_pos);
540 0 : }
541 :
542 :
543 : /// interpolates a function on a subset
544 : /**
545 : * This function interpolates a GridFunction. To evaluate the function on every
546 : * point a functor must be passed.
547 : *
548 : * @param[in] spInterpolFunction data providing interpolation values
549 : * @param[out] spGridFct interpolated grid function
550 : * @param[in] cmp symbolic name of function component
551 : * @param[in] subsets subsets, where to interpolate (NULL = everywhere)
552 : * @param[in] time time point
553 : */
554 : template <typename TGridFunction>
555 0 : void Interpolate(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
556 : SmartPtr<TGridFunction> spGridFct, const size_t fct,
557 : const char* subsets, number time, const MathVector<TGridFunction::dim> diff_pos)
558 : {
559 : const bool bAllowManyfoldInterpolation =
560 : (spGridFct->local_finite_element_id(fct).type() == LFEID::LAGRANGE);
561 :
562 : // create subset group
563 0 : SubsetGroup ssGrp(spGridFct->domain()->subset_handler());
564 0 : if(subsets != NULL)
565 : {
566 0 : ssGrp.add(TokenizeString(subsets));
567 0 : if(!bAllowManyfoldInterpolation)
568 0 : if(!SameDimensionsInAllSubsets(ssGrp))
569 0 : UG_THROW("Interpolate: Subsets '"<<subsets<<"' do not have same dimension."
570 : "Can not integrate on subsets of different dimensions.");
571 : }
572 : else
573 : {
574 : // add all subsets and remove lower dim subsets afterwards
575 0 : ssGrp.add_all();
576 0 : if(!bAllowManyfoldInterpolation)
577 0 : RemoveLowerDimSubsets(ssGrp);
578 : }
579 :
580 0 : Interpolate(spInterpolFunction, spGridFct, fct, ssGrp, time, diff_pos);
581 0 : }
582 :
583 :
584 : template <typename TGridFunction>
585 0 : void Interpolate(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
586 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
587 : number time)
588 : {
589 0 : Interpolate(spInterpolFunction, spGridFct, cmp, NULL, time);
590 0 : }
591 : template <typename TGridFunction>
592 0 : void Interpolate(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
593 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
594 : const char* subsets)
595 : {
596 0 : Interpolate(spInterpolFunction, spGridFct, cmp, subsets, 0.0);
597 0 : }
598 : template <typename TGridFunction>
599 0 : void Interpolate(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
600 : SmartPtr<TGridFunction> spGridFct, const char* cmp)
601 : {
602 0 : Interpolate(spInterpolFunction, spGridFct, cmp, NULL, 0.0);
603 0 : }
604 : //interpolate with diff_vector
605 : template <typename TGridFunction>
606 0 : void Interpolate(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
607 : SmartPtr<TGridFunction> spGridFct, const char* cmp, const MathVector<TGridFunction::dim>& m_diff_pos)
608 : {
609 0 : Interpolate(spInterpolFunction, spGridFct, cmp, NULL, 0.0, m_diff_pos);
610 0 : }
611 :
612 : ///////////////
613 : // const data
614 : ///////////////
615 :
616 : template <typename TGridFunction>
617 0 : void Interpolate(number val,
618 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
619 : const char* subsets, number time)
620 : {
621 : static const int dim = TGridFunction::dim;
622 0 : SmartPtr<UserData<number, dim> > sp =
623 0 : make_sp(new ConstUserNumber<dim>(val));
624 0 : Interpolate(sp, spGridFct, cmp, subsets, time);
625 0 : }
626 : template <typename TGridFunction>
627 0 : void Interpolate(number val,
628 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
629 : number time)
630 0 : {Interpolate(val, spGridFct, cmp, NULL, time);}
631 : template <typename TGridFunction>
632 0 : void Interpolate(number val,
633 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
634 : const char* subsets)
635 0 : {Interpolate(val, spGridFct, cmp, subsets, 0.0);}
636 : template <typename TGridFunction>
637 0 : void Interpolate(number val,
638 : SmartPtr<TGridFunction> spGridFct, const char* cmp)
639 0 : {Interpolate(val, spGridFct, cmp, NULL, 0.0);}
640 :
641 : //interpolate with diff_vector
642 : template <typename TGridFunction>
643 : void Interpolate(number val,
644 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
645 : const char* subsets, number time,const SmartPtr<CplUserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > m_diff_pos)
646 : {
647 : static const int dim = TGridFunction::dim;
648 : SmartPtr<UserData<number, dim> > sp =
649 : make_sp(new ConstUserNumber<dim>(val));
650 :
651 : InterpolateDiff(sp, spGridFct, cmp, subsets, time,m_diff_pos);
652 : }
653 :
654 :
655 :
656 : ///////////////
657 : // lua data
658 : ///////////////
659 :
660 : #ifdef UG_FOR_LUA
661 : // function-name as string
662 : template <typename TGridFunction>
663 0 : void Interpolate(const char* LuaFunction,
664 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
665 : const char* subsets, number time)
666 : {
667 : static const int dim = TGridFunction::dim;
668 0 : SmartPtr<UserData<number, dim> > sp =
669 : LuaUserDataFactory<number, dim>::create(LuaFunction);
670 0 : Interpolate(sp, spGridFct, cmp, subsets, time);
671 0 : }
672 : template <typename TGridFunction>
673 0 : void Interpolate(const char* LuaFunction,
674 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
675 : number time)
676 0 : {Interpolate(LuaFunction, spGridFct, cmp, NULL, time);}
677 : template <typename TGridFunction>
678 0 : void Interpolate(const char* LuaFunction,
679 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
680 : const char* subsets)
681 0 : {Interpolate(LuaFunction, spGridFct, cmp, subsets, 0.0);}
682 : template <typename TGridFunction>
683 0 : void Interpolate(const char* LuaFunction,
684 : SmartPtr<TGridFunction> spGridFct, const char* cmp)
685 0 : {Interpolate(LuaFunction, spGridFct, cmp, NULL, 0.0);}
686 :
687 : // function as function handle
688 : template <typename TGridFunction>
689 0 : void Interpolate(LuaFunctionHandle LuaFunction,
690 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
691 : const char* subsets, number time)
692 : {
693 : static const int dim = TGridFunction::dim;
694 0 : SmartPtr<UserData<number, dim> > sp =
695 0 : make_sp(new LuaUserData<number, dim>(LuaFunction));
696 0 : Interpolate(sp, spGridFct, cmp, subsets, time);
697 0 : }
698 : template <typename TGridFunction>
699 0 : void Interpolate(LuaFunctionHandle LuaFunction,
700 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
701 : number time)
702 0 : {Interpolate(LuaFunction, spGridFct, cmp, NULL, time);}
703 : template <typename TGridFunction>
704 0 : void Interpolate(LuaFunctionHandle LuaFunction,
705 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
706 : const char* subsets)
707 0 : {Interpolate(LuaFunction, spGridFct, cmp, subsets, 0.0);}
708 : template <typename TGridFunction>
709 0 : void Interpolate(LuaFunctionHandle LuaFunction,
710 : SmartPtr<TGridFunction> spGridFct, const char* cmp)
711 0 : {Interpolate(LuaFunction, spGridFct, cmp, NULL, 0.0);}
712 :
713 : //interpolate with Diff-vector:
714 : template <typename TGridFunction>
715 : void InterpolateDiff(const char* LuaFunction,
716 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
717 : const char* subsets, number time, SmartPtr<CplUserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > m_diff_pos )
718 : {
719 : static const int dim = TGridFunction::dim;
720 : SmartPtr<UserData<number, dim> > sp =
721 : LuaUserDataFactory<number, dim>::create(LuaFunction);
722 : InterpolateDiff(sp, spGridFct, cmp,subsets, time, m_diff_pos);
723 : }
724 :
725 : template <typename TGridFunction>
726 : void InterpolateDiff(LuaFunctionHandle LuaFunction,
727 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
728 : const char* subsets, number time,SmartPtr<CplUserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > m_diff_pos)
729 : {
730 : static const int dim = TGridFunction::dim;
731 : SmartPtr<UserData<number, dim> > sp =
732 : make_sp(new LuaUserData<number, dim>(LuaFunction));
733 : InterpolateDiff(sp, spGridFct, cmp, subsets, time, m_diff_pos);
734 : }
735 : template <typename TGridFunction>
736 : void Interpolate(LuaFunctionHandle LuaFunction,
737 : SmartPtr<TGridFunction> spGridFct, const char* cmp,const SmartPtr<CplUserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > m_diff_pos)
738 : {InterpolateDiff(LuaFunction, spGridFct, cmp, NULL, 0.0, m_diff_pos);}
739 :
740 :
741 :
742 : #endif
743 : #ifdef UG_PARALLEL
744 : /// interpolates a gridFunction across similiar geometries with different meshing, where distribution across processes is not
745 : /// consistent. This is inefficient and should not be used in hot loops.
746 : /**
747 : * Creates GlobalGridFunctionNumber and transfers unknown values across different vertex configurations on different processes
748 : *
749 : * @param[in] spGfSource data providing interpolation values
750 : * @param[out] spGfTarget interpolated grid function
751 : * @param[in] fct symbolic name of function component to be interpolated
752 : * @param[in] time time point
753 : */
754 : template <typename TGridFunction>
755 : void InterpolateGlobalGridFunctionAcrossProcesses(SmartPtr<TGridFunction> spGfSource,
756 : SmartPtr<TGridFunction> spGfTarget,
757 : const char* cmp,
758 : number time)
759 : {
760 :
761 : static const int dim = TGridFunction::dim;
762 : SmartPtr<GlobalGridFunctionNumberData<TGridFunction> > data = make_sp(new GlobalGridFunctionNumberData<TGridFunction>(spGfSource, cmp)) ;
763 : // domain type and position_type
764 :
765 : // check if we have multiple processes - if no we can just forward globalGfUserData to standard interpolate
766 : if (pcl::NumProcs() == 1) {
767 : Interpolate(data, spGfTarget, cmp, time);
768 : return;
769 : }
770 : typedef typename TGridFunction::domain_type domain_type;
771 : typedef typename domain_type::position_type position_type;
772 :
773 :
774 : // get position accessor (of interpolated grid function)
775 : const typename domain_type::position_accessor_type& aaPos
776 : = spGfTarget->domain()->position_accessor();
777 :
778 :
779 : std::vector<DoFIndex> ind;
780 : typename TGridFunction::template dim_traits<0>::const_iterator iterEnd, iter;
781 : const size_t fct = spGfSource->fct_id_by_name(cmp);
782 :
783 :
784 : pcl::ProcessCommunicator com;
785 : // each process iterates over its own local elements
786 : for (int i = 0; i < pcl::NumProcs(); ++i){
787 : // iterate over all elements
788 : position_type glob_pos;
789 : bool finished = false;
790 : // check for active Proc
791 : if (pcl::ProcRank() == i) {
792 : iterEnd = spGfTarget->template end<Vertex>();
793 : iter = spGfTarget->template begin<Vertex>();
794 : for(; iter != iterEnd; )
795 : {
796 : // get vertex
797 : Vertex* vrt = *iter;
798 :
799 : // global position (in case this position is not contained on the same proc for both distributions, we have to
800 : // search all procs)
801 : glob_pos = aaPos[vrt]; // position (of interpolated grid function)
802 : // current active proc sends current evaluation position
803 : com.broadcast(glob_pos, i);
804 :
805 :
806 : // value at position
807 : number val;
808 : // evaluate source grid function value for cmp and given position
809 : // evaluate_global() is already mpi - parallelized and thus needs to be called with the same
810 : // poosition from all procs - see logic for non active proc below
811 : data->evaluate_global(val, glob_pos);
812 :
813 : // set local function value
814 : // get multiindices of element
815 : spGfTarget->dof_indices(vrt, fct, ind);
816 :
817 : // loop all dofs
818 : for(size_t i = 0; i < ind.size(); ++i)
819 : {
820 : // set value
821 : DoFRef(*spGfTarget, ind[i]) = val;
822 : }
823 : ++iter;
824 : // check if active proc is finished with element loop and broadcast result
825 : finished = (iter == iterEnd);
826 : com.broadcast(finished, i);
827 : }
828 :
829 :
830 : }
831 : // if we are not the active proc, we wait for the active proc to finish iterating over its local elements
832 : // and answer the evaluate_global function call with the corresponding broadcasted position
833 : else {
834 : while (!finished){
835 : // while active proc has elements, we recieve current position and evaluate
836 : com.broadcast(glob_pos, i);
837 :
838 : number val;
839 : // matching evaluate_global call for active Proc
840 : data->evaluate_global(val, glob_pos);
841 : com.broadcast(finished, i);
842 : // finish loop when active proc has no elements left
843 : }
844 : }
845 :
846 :
847 : }
848 : // resulting vector is consistent
849 : spGfTarget->set_storage_type(PST_CONSISTENT);
850 :
851 : }
852 : // implement as overload for Interpolate
853 : template <typename TGridFunction>
854 : void Interpolate(SmartPtr<TGridFunction> spGfSource,
855 : SmartPtr<TGridFunction> spGfTarget,
856 : const char* cmp,
857 : number time)
858 : {InterpolateGlobalGridFunctionAcrossProcesses(spGfSource, spGfTarget,cmp, time);}
859 : #endif
860 :
861 : } // namespace ug
862 :
863 : #endif /*__H__UG__LIB_DISC__FUNCTION_SPACES__INTERPOLATE__*/
|
