Line data Source code
1 : /*
2 : * Copyright (c) 2015: G-CSC, Goethe University Frankfurt
3 : * Author: Markus Breit
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 UG__LIB_DISC__FUNCTION_SPACES__INTERPOLATE_INNER_H
34 : #define UG__LIB_DISC__FUNCTION_SPACES__INTERPOLATE_INNER_H
35 :
36 : #include "common/common.h"
37 : #include "common/util/smart_pointer.h"
38 :
39 : #include "lib_grid/tools/subset_group.h"
40 :
41 : #include "lib_disc/function_spaces/interpolate.h"
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 : #include "lib_disc/function_spaces/dof_position_util.h"
48 :
49 : #ifdef UG_FOR_LUA
50 : #include "bindings/lua/lua_user_data.h"
51 : #endif
52 :
53 : namespace ug{
54 :
55 :
56 : /**
57 : * This function interpolates a grid function on a specific subset and for
58 : * a specific element type. All inner dofs of all grid elements of the given
59 : * type in the given subset will be assigned an interpolated value.
60 : *
61 : * @param[in] spInterpolFunction data providing interpolation values
62 : * @param[out] spGridFct interpolated grid function
63 : * @param[in] fct symbolic name of function component
64 : * @param[in] si subset, where to interpolate
65 : * @param[in] time time point
66 : */
67 : template <typename TElem, typename TGridFunction>
68 0 : void InterpolateOnElementsInner(
69 : SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
70 : SmartPtr<TGridFunction> spGridFct, size_t fct, int si, number time, const MathVector<TGridFunction::dim> diff_pos)
71 : {
72 : // do nothing if no DoFs here
73 0 : if (!spGridFct->max_fct_dofs(fct, ROID_TETRAHEDRON, si))
74 0 : return;
75 :
76 : // domain type and position_type
77 : typedef typename TGridFunction::domain_type domain_type;
78 : typedef typename domain_type::position_type position_type;
79 :
80 : // get iterators
81 : typename TGridFunction::template traits<TElem>::const_iterator iterEnd, iter;
82 0 : iterEnd = spGridFct->template end<TElem>(si);
83 0 : iter = spGridFct->template begin<TElem>(si);
84 :
85 : // check if something to do:
86 0 : if (iter == iterEnd) return;
87 :
88 : // id of shape functions used
89 0 : LFEID id = spGridFct->local_finite_element_id(fct);
90 :
91 : // iterate over all elements
92 0 : for (; iter != iterEnd; ++iter)
93 : {
94 : // get element
95 : TElem* elem = *iter;
96 :
97 : // get multiindices of element
98 : std::vector<DoFIndex> ind;
99 : spGridFct->inner_dof_indices(elem, fct, ind);
100 :
101 : // global positions of DoFs
102 : std::vector<position_type> glob_pos;
103 0 : InnerDoFPosition<domain_type>(glob_pos, elem, *spGridFct->domain(), id);
104 :
105 : // check global positions size
106 : size_t ind_sz = ind.size();
107 : size_t gp_sz = glob_pos.size();
108 0 : if (ind_sz != gp_sz)
109 0 : UG_THROW("InterpolateOnElem: On subset " << si << ": Number of DoFs is "
110 : << ind_sz << ", but number of DoF positions is "
111 : << gp_sz << "." << std::endl);
112 :
113 : // loop all dofs
114 0 : for (size_t i = 0; i < ind_sz && i < gp_sz; ++i)
115 : {
116 : // value at position
117 : number val;
118 : position_type rel_pos=glob_pos[i];
119 : rel_pos-=diff_pos;
120 :
121 0 : (*spInterpolFunction)(val, rel_pos, time, si);
122 :
123 : // set value
124 0 : DoFRef(*spGridFct, ind[i]) = val;
125 : }
126 : }
127 : }
128 : /**
129 : * This function interpolates a grid function on a specific subset and for
130 : * a specific element type. All inner dofs of all grid elements of the given
131 : * type in the given subset will be assigned an interpolated value.
132 : *
133 : * @param[in] spInterpolFunction data providing interpolation values
134 : * @param[out] spGridFct interpolated grid function
135 : * @param[in] fct symbolic name of function component
136 : * @param[in] si subset, where to interpolate
137 : * @param[in] time time point
138 : */
139 : template <typename TElem, typename TGridFunction>
140 0 : void InterpolateOnElementsInner(
141 : SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
142 : SmartPtr<TGridFunction> spGridFct, size_t fct, int si, number time)
143 : {
144 : // do nothing if no DoFs here
145 0 : if (!spGridFct->max_fct_dofs(fct, ROID_TETRAHEDRON, si))
146 : return;
147 :
148 : // domain type and position_type
149 : typedef typename TGridFunction::domain_type domain_type;
150 : typedef typename domain_type::position_type position_type;
151 : typedef typename position_type::value_type value_type;
152 :
153 : // dimension of reference element
154 : const int dim = TGridFunction::dim;
155 :
156 : // get iterators
157 : /*typename TGridFunction::template traits<TElem>::const_iterator iterEnd, iter;
158 : iterEnd = spGridFct->template end<TElem>(si);
159 : iter = spGridFct->template begin<TElem>(si);
160 :
161 : // check if something to do:
162 : if (iter == iterEnd) return;
163 :
164 : // id of shape functions used
165 : LFEID id = spGridFct->local_finite_element_id(fct);
166 :
167 : // iterate over all elements
168 : for (; iter != iterEnd; ++iter)
169 : {
170 : // get element
171 : TElem* elem = *iter;
172 :
173 : // get multiindices of element
174 : std::vector<DoFIndex> ind;
175 : spGridFct->inner_dof_indices(elem, fct, ind);
176 :
177 : // global positions of DoFs
178 : std::vector<position_type> glob_pos;
179 : InnerDoFPosition<domain_type>(glob_pos, elem, *spGridFct->domain(), id);
180 :
181 : // check global positions size
182 : size_t ind_sz = ind.size();
183 : size_t gp_sz = glob_pos.size();
184 : if (ind_sz != gp_sz)
185 : UG_THROW("InterpolateOnElem: On subset " << si << ": Number of DoFs is "
186 : << ind_sz << ", but number of DoF positions is "
187 : << gp_sz << "." << std::endl);
188 :
189 : // loop all dofs
190 : for (size_t i = 0; i < ind_sz && i < gp_sz; ++i)
191 : {
192 : // value at position
193 : number val;
194 : (*spInterpolFunction)(val, glob_pos[i], time, si);
195 :
196 : // set value
197 : DoFRef(*spGridFct, ind[i]) = val;
198 : }
199 : }*/
200 :
201 0 : MathVector<dim>* diff_pos=new MathVector<dim, value_type>();
202 0 : InterpolateOnElementsInner<TElem,TGridFunction>(spInterpolFunction, spGridFct, fct, si,time, *diff_pos);
203 :
204 : }
205 :
206 :
207 : /**
208 : * This function interpolates a grid function on a subset group. All inner
209 : * dofs of all grid elements in the subset group will be assigned an inter-
210 : * polated value.
211 : *
212 : * @param[in] spInterpolFunction data providing interpolation values
213 : * @param[out] spGridFct interpolated grid function
214 : * @param[in] fct symbolic name of function component
215 : * @param[in] ssGrp subsets, where to interpolate
216 : * @param[in] time time point
217 : */
218 : template <typename TGridFunction>
219 0 : void InterpolateOnElementsInner
220 : (
221 : SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
222 : SmartPtr<TGridFunction> spGridFct,
223 : size_t fct,
224 : number time,
225 : const SubsetGroup& ssGrp
226 : )
227 : {
228 : // loop subsets
229 0 : for(size_t i = 0; i < ssGrp.size(); ++i)
230 : {
231 : // get subset index
232 : const int si = ssGrp[i];
233 :
234 : // skip if function is not defined in subset
235 0 : if (!spGridFct->is_def_in_subset(fct, si)) continue;
236 :
237 : // switch dimensions
238 : try
239 : {
240 0 : const int dim = ssGrp.dim(i);
241 :
242 0 : if (dim > TGridFunction::dim)
243 0 : UG_THROW("InterpolateOnElements: Dimension of subset is " << dim << ", but "
244 : " World Dimension is " << TGridFunction::dim << ". Cannot interpolate this.");
245 :
246 : // In a parallel scenario, the distribution CAN cause elements of of lower
247 : // dimension than the rest of their subset to be located disconnected from
248 : // the rest of the subset on a processor. For example, in 2D, think of a
249 : // (1D) boundary subset and a distribution where the boundary of a proc's
250 : // domain only touches the boundary subset in a vertex, but intersects with
251 : // the boundary subset in another place.
252 : // Therefore, we need to interpolate on all dimensions < dim and we will
253 : // only consider inner indices then, thus we simply switch-case backwards
254 : // and without breaks.
255 0 : switch (dim)
256 : {
257 : case DIM_SUBSET_EMPTY_GRID: break;
258 : case 3:
259 0 : if (spGridFct->max_fct_dofs(fct, VOLUME, si))
260 : {
261 : InterpolateOnElementsInner<Tetrahedron, TGridFunction>
262 0 : (spInterpolFunction, spGridFct, fct, si, time);
263 : InterpolateOnElementsInner<Hexahedron, TGridFunction>
264 0 : (spInterpolFunction, spGridFct, fct, si, time);
265 : InterpolateOnElementsInner<Prism, TGridFunction>
266 0 : (spInterpolFunction, spGridFct, fct, si, time);
267 : InterpolateOnElementsInner<Pyramid, TGridFunction>
268 0 : (spInterpolFunction, spGridFct, fct, si, time);
269 : InterpolateOnElementsInner<Octahedron, TGridFunction>
270 0 : (spInterpolFunction, spGridFct, fct, si, time);
271 : }
272 : case 2:
273 0 : if (spGridFct->max_fct_dofs(fct, FACE, si))
274 : {
275 : InterpolateOnElementsInner<Triangle, TGridFunction>
276 0 : (spInterpolFunction, spGridFct, fct, si, time);
277 : InterpolateOnElementsInner<Quadrilateral, TGridFunction>
278 0 : (spInterpolFunction, spGridFct, fct, si, time);
279 : InterpolateOnElementsInner<ConstrainingTriangle, TGridFunction>
280 0 : (spInterpolFunction, spGridFct, fct, si, time);
281 : InterpolateOnElementsInner<ConstrainingQuadrilateral, TGridFunction>
282 0 : (spInterpolFunction, spGridFct, fct, si, time);
283 : InterpolateOnElementsInner<ConstrainedTriangle, TGridFunction>
284 0 : (spInterpolFunction, spGridFct, fct, si, time);
285 : InterpolateOnElementsInner<ConstrainedQuadrilateral, TGridFunction>
286 0 : (spInterpolFunction, spGridFct, fct, si, time);
287 : }
288 : case 1:
289 0 : if (spGridFct->max_fct_dofs(fct, EDGE, si))
290 : {
291 : InterpolateOnElementsInner<RegularEdge, TGridFunction>
292 0 : (spInterpolFunction, spGridFct, fct, si, time);
293 : InterpolateOnElementsInner<ConstrainingEdge, TGridFunction>
294 0 : (spInterpolFunction, spGridFct, fct, si, time);
295 : InterpolateOnElementsInner<ConstrainedEdge, TGridFunction>
296 0 : (spInterpolFunction, spGridFct, fct, si, time);
297 : }
298 : case 0:
299 0 : if (spGridFct->max_fct_dofs(fct, VERTEX, si))
300 : {
301 : InterpolateOnElementsInner<RegularVertex, TGridFunction>
302 0 : (spInterpolFunction, spGridFct, fct, si, time);
303 : InterpolateOnElementsInner<ConstrainedVertex, TGridFunction>
304 0 : (spInterpolFunction, spGridFct, fct, si, time);
305 : }
306 : break;
307 0 : default: UG_THROW("InterpolateOnElements: Dimension " <<dim<<
308 : " not possible for world dim "<<3<<".");
309 : }
310 : }
311 0 : UG_CATCH_THROW("InterpolateOnElements: Cannot interpolate on elements.");
312 : }
313 0 : }
314 :
315 : ////////////////////////////////////////////////////////////////////////////////
316 : // Interpolate routine
317 : ////////////////////////////////////////////////////////////////////////////////
318 :
319 : /// interpolates a function on a subset
320 : /**
321 : * This function interpolates a Lagrange type GridFunction.
322 : * In contrast to its "big brother" Interpolate() in interpolate.h,
323 : * this function will only write to _inner_ DoFs, but to _all_ of them,
324 : * on any element located in any of the given subsets.
325 : * This is more intuitive and secure than writing to all (not necessarily
326 : * inner) DoFs of an element located in a subset, as this takes into
327 : * account the unique subsets of sub-elements.
328 : *
329 : * At the moment, this is only meant for Lagrange-type elements.
330 : * Please feel free to add functionality for other types if you like.
331 : *
332 : * @param[in] spInterpolFunction data providing interpolation values
333 : * @param[out] spGridFct interpolated grid function
334 : * @param[in] cmp symbolic name of function component
335 : * @param[in] subsets subsets, where to interpolate (NULL = everywhere)
336 : * @param[in] time time point
337 : */
338 : template <typename TGridFunction>
339 : void InterpolateInnerDiff(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
340 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
341 : const char* subsets, number time, const SmartPtr<CplUserData<MathVector<TGridFunction::dim>, TGridFunction::dim> > m_diff_pos)
342 : {
343 : // check, that values do not depend on a solution
344 : if (spInterpolFunction->requires_grid_fct())
345 : UG_THROW("Interpolate: The interpolation values depend on a grid function."
346 : " This is not allowed in the current implementation. Use constant,"
347 : " lua-callback or vrl-callback user data only (even within linkers).");
348 :
349 : // get function id of name
350 : const size_t fct = spGridFct->fct_id_by_name(cmp);
351 :
352 : // check that function found
353 : if (fct > spGridFct->num_fct())
354 : UG_THROW("Interpolate: Name of component '"<< cmp <<"' not found.");
355 :
356 : // check that type is Lagrange
357 : if (spGridFct->local_finite_element_id(fct).type() != LFEID::LAGRANGE)
358 : UG_THROW("This interpolation only allows Lagrange-type elements.\n"
359 : "Feel free to add functionality for other types as needed.");
360 :
361 : // check if fast P1 interpolation can be used
362 : const bool bUseP1Interpolation
363 : = spGridFct->local_finite_element_id(fct).order() == 1;
364 :
365 : // create subset group
366 : SubsetGroup ssGrp(spGridFct->domain()->subset_handler());
367 : if (subsets != NULL)
368 : ssGrp.add(TokenizeString(subsets));
369 : else
370 : ssGrp.add_all();
371 :
372 : // forward
373 : if (bUseP1Interpolation)
374 : InterpolateOnDiffVertices<TGridFunction>(spInterpolFunction, spGridFct, fct, time, ssGrp, m_diff_pos);
375 : else
376 : InterpolateOnElementsInner<TGridFunction>(spInterpolFunction, spGridFct, fct, time, ssGrp, m_diff_pos);
377 :
378 : // adjust parallel storage state
379 : #ifdef UG_PARALLEL
380 : spGridFct->set_storage_type(PST_CONSISTENT);
381 : #endif
382 : }
383 :
384 : /// interpolates a function on a subset
385 : /**
386 : * This function interpolates a Lagrange type GridFunction.
387 : * In contrast to its "big brother" Interpolate() in interpolate.h,
388 : * this function will only write to _inner_ DoFs, but to _all_ of them,
389 : * on any element located in any of the given subsets.
390 : * This is more intuitive and secure than writing to all (not necessarily
391 : * inner) DoFs of an element located in a subset, as this takes into
392 : * account the unique subsets of sub-elements.
393 : *
394 : * At the moment, this is only meant for Lagrange-type elements.
395 : * Please feel free to add functionality for other types if you like.
396 : *
397 : * @param[in] spInterpolFunction data providing interpolation values
398 : * @param[out] spGridFct interpolated grid function
399 : * @param[in] cmp symbolic name of function component
400 : * @param[in] subsets subsets, where to interpolate (NULL = everywhere)
401 : * @param[in] time time point
402 : */
403 : template <typename TGridFunction>
404 0 : void InterpolateInner(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
405 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
406 : const char* subsets, number time)
407 : {
408 : // check, that values do not depend on a solution
409 0 : if (spInterpolFunction->requires_grid_fct())
410 0 : UG_THROW("Interpolate: The interpolation values depend on a grid function."
411 : " This is not allowed in the current implementation. Use constant,"
412 : " lua-callback or vrl-callback user data only (even within linkers).");
413 :
414 : // get function id of name
415 : const size_t fct = spGridFct->fct_id_by_name(cmp);
416 :
417 : // check that function found
418 0 : if (fct > spGridFct->num_fct())
419 0 : UG_THROW("Interpolate: Name of component '"<< cmp <<"' not found.");
420 :
421 : // check that type is Lagrange
422 0 : if (spGridFct->local_finite_element_id(fct).type() != LFEID::LAGRANGE)
423 0 : UG_THROW("This interpolation only allows Lagrange-type elements.\n"
424 : "Feel free to add functionality for other types as needed.");
425 :
426 : // check if fast P1 interpolation can be used
427 : const bool bUseP1Interpolation
428 : = spGridFct->local_finite_element_id(fct).order() == 1;
429 :
430 : // create subset group
431 0 : SubsetGroup ssGrp(spGridFct->domain()->subset_handler());
432 0 : if (subsets != NULL)
433 0 : ssGrp.add(TokenizeString(subsets));
434 : else
435 0 : ssGrp.add_all();
436 :
437 : // forward
438 0 : if (bUseP1Interpolation)
439 0 : InterpolateOnVertices<TGridFunction>(spInterpolFunction, spGridFct, fct, time, ssGrp);
440 : else
441 0 : InterpolateOnElementsInner<TGridFunction>(spInterpolFunction, spGridFct, fct, time, ssGrp);
442 :
443 : // adjust parallel storage state
444 : #ifdef UG_PARALLEL
445 : spGridFct->set_storage_type(PST_CONSISTENT);
446 : #endif
447 0 : }
448 :
449 :
450 :
451 : template <typename TGridFunction>
452 0 : void InterpolateInner(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
453 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
454 : number time)
455 0 : {InterpolateInner(spInterpolFunction, spGridFct, cmp, NULL, time);}
456 : template <typename TGridFunction>
457 0 : void InterpolateInner(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
458 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
459 : const char* subsets)
460 0 : {InterpolateInner(spInterpolFunction, spGridFct, cmp, subsets, 0.0);}
461 : template <typename TGridFunction>
462 0 : void InterpolateInner(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
463 : SmartPtr<TGridFunction> spGridFct, const char* cmp)
464 0 : {InterpolateInner(spInterpolFunction, spGridFct, cmp, NULL, 0.0);}
465 :
466 : ///////////////
467 : // const data
468 : ///////////////
469 :
470 : template <typename TGridFunction>
471 0 : void InterpolateInner(number val,
472 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
473 : const char* subsets, number time)
474 : {
475 : static const int dim = TGridFunction::dim;
476 0 : SmartPtr<UserData<number, dim> > sp =
477 0 : make_sp(new ConstUserNumber<dim>(val));
478 0 : InterpolateInner(sp, spGridFct, cmp, subsets, time);
479 0 : }
480 : template <typename TGridFunction>
481 0 : void InterpolateInner(number val,
482 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
483 : number time)
484 0 : {InterpolateInner(val, spGridFct, cmp, NULL, time);}
485 : template <typename TGridFunction>
486 0 : void InterpolateInner(number val,
487 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
488 : const char* subsets)
489 0 : {InterpolateInner(val, spGridFct, cmp, subsets, 0.0);}
490 : template <typename TGridFunction>
491 0 : void InterpolateInner(number val,
492 : SmartPtr<TGridFunction> spGridFct, const char* cmp)
493 0 : {InterpolateInner(val, spGridFct, cmp, NULL, 0.0);}
494 :
495 : ///////////////
496 : // lua data
497 : ///////////////
498 :
499 : #ifdef UG_FOR_LUA
500 : template <typename TGridFunction>
501 0 : void InterpolateInner(const char* LuaFunction,
502 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
503 : const char* subsets, number time)
504 : {
505 : static const int dim = TGridFunction::dim;
506 0 : SmartPtr<UserData<number, dim> > sp =
507 : LuaUserDataFactory<number, dim>::create(LuaFunction);
508 0 : InterpolateInner(sp, spGridFct, cmp, subsets, time);
509 0 : }
510 : template <typename TGridFunction>
511 0 : void InterpolateInner(const char* LuaFunction,
512 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
513 : number time)
514 0 : {InterpolateInner(LuaFunction, spGridFct, cmp, NULL, time);}
515 : template <typename TGridFunction>
516 0 : void InterpolateInner(const char* LuaFunction,
517 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
518 : const char* subsets)
519 0 : {InterpolateInner(LuaFunction, spGridFct, cmp, subsets, 0.0);}
520 : template <typename TGridFunction>
521 0 : void InterpolateInner(const char* LuaFunction,
522 : SmartPtr<TGridFunction> spGridFct, const char* cmp)
523 0 : {InterpolateInner(LuaFunction, spGridFct, cmp, NULL, 0.0);}
524 : #endif
525 :
526 :
527 : } // namespace ug
528 :
529 : #endif // UG__LIB_DISC__FUNCTION_SPACES__INTERPOLATE_INNER_H
|
