Loading [MathJax]/extensions/tex2jax.js
ug4
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
interpolate_inner.h
Go to the documentation of this file.
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
53namespace ug{
54
55
67template <typename TElem, typename TGridFunction>
68void 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 if (!spGridFct->max_fct_dofs(fct, ROID_TETRAHEDRON, si))
74 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 iterEnd = spGridFct->template end<TElem>(si);
83 iter = spGridFct->template begin<TElem>(si);
84
85 // check if something to do:
86 if (iter == iterEnd) return;
87
88 // id of shape functions used
89 LFEID id = spGridFct->local_finite_element_id(fct);
90
91 // iterate over all elements
92 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 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 if (ind_sz != gp_sz)
109 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 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 (*spInterpolFunction)(val, rel_pos, time, si);
122
123 // set value
124 DoFRef(*spGridFct, ind[i]) = val;
125 }
126 }
127}
139template <typename TElem, typename TGridFunction>
140void 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 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 MathVector<dim>* diff_pos=new MathVector<dim, value_type>();
202 InterpolateOnElementsInner<TElem,TGridFunction>(spInterpolFunction, spGridFct, fct, si,time, *diff_pos);
203
204}
205
206
218template <typename TGridFunction>
219void 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 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 if (!spGridFct->is_def_in_subset(fct, si)) continue;
236
237 // switch dimensions
238 try
239 {
240 const int dim = ssGrp.dim(i);
241
242 if (dim > TGridFunction::dim)
243 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 switch (dim)
256 {
257 case DIM_SUBSET_EMPTY_GRID: break;
258 case 3:
259 if (spGridFct->max_fct_dofs(fct, VOLUME, si))
260 {
261 InterpolateOnElementsInner<Tetrahedron, TGridFunction>
262 (spInterpolFunction, spGridFct, fct, si, time);
263 InterpolateOnElementsInner<Hexahedron, TGridFunction>
264 (spInterpolFunction, spGridFct, fct, si, time);
265 InterpolateOnElementsInner<Prism, TGridFunction>
266 (spInterpolFunction, spGridFct, fct, si, time);
267 InterpolateOnElementsInner<Pyramid, TGridFunction>
268 (spInterpolFunction, spGridFct, fct, si, time);
269 InterpolateOnElementsInner<Octahedron, TGridFunction>
270 (spInterpolFunction, spGridFct, fct, si, time);
271 }
272 case 2:
273 if (spGridFct->max_fct_dofs(fct, FACE, si))
274 {
275 InterpolateOnElementsInner<Triangle, TGridFunction>
276 (spInterpolFunction, spGridFct, fct, si, time);
277 InterpolateOnElementsInner<Quadrilateral, TGridFunction>
278 (spInterpolFunction, spGridFct, fct, si, time);
279 InterpolateOnElementsInner<ConstrainingTriangle, TGridFunction>
280 (spInterpolFunction, spGridFct, fct, si, time);
281 InterpolateOnElementsInner<ConstrainingQuadrilateral, TGridFunction>
282 (spInterpolFunction, spGridFct, fct, si, time);
283 InterpolateOnElementsInner<ConstrainedTriangle, TGridFunction>
284 (spInterpolFunction, spGridFct, fct, si, time);
285 InterpolateOnElementsInner<ConstrainedQuadrilateral, TGridFunction>
286 (spInterpolFunction, spGridFct, fct, si, time);
287 }
288 case 1:
289 if (spGridFct->max_fct_dofs(fct, EDGE, si))
290 {
291 InterpolateOnElementsInner<RegularEdge, TGridFunction>
292 (spInterpolFunction, spGridFct, fct, si, time);
293 InterpolateOnElementsInner<ConstrainingEdge, TGridFunction>
294 (spInterpolFunction, spGridFct, fct, si, time);
295 InterpolateOnElementsInner<ConstrainedEdge, TGridFunction>
296 (spInterpolFunction, spGridFct, fct, si, time);
297 }
298 case 0:
299 if (spGridFct->max_fct_dofs(fct, VERTEX, si))
300 {
301 InterpolateOnElementsInner<RegularVertex, TGridFunction>
302 (spInterpolFunction, spGridFct, fct, si, time);
303 InterpolateOnElementsInner<ConstrainedVertex, TGridFunction>
304 (spInterpolFunction, spGridFct, fct, si, time);
305 }
306 break;
307 default: UG_THROW("InterpolateOnElements: Dimension " <<dim<<
308 " not possible for world dim "<<3<<".");
309 }
310 }
311 UG_CATCH_THROW("InterpolateOnElements: Cannot interpolate on elements.");
312 }
313}
314
316// Interpolate routine
318
320
338template <typename TGridFunction>
339void 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
385
403template <typename TGridFunction>
404void 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 if (spInterpolFunction->requires_grid_fct())
410 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 if (fct > spGridFct->num_fct())
419 UG_THROW("Interpolate: Name of component '"<< cmp <<"' not found.");
420
421 // check that type is Lagrange
422 if (spGridFct->local_finite_element_id(fct).type() != LFEID::LAGRANGE)
423 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 SubsetGroup ssGrp(spGridFct->domain()->subset_handler());
432 if (subsets != NULL)
433 ssGrp.add(TokenizeString(subsets));
434 else
435 ssGrp.add_all();
436
437 // forward
438 if (bUseP1Interpolation)
439 InterpolateOnVertices<TGridFunction>(spInterpolFunction, spGridFct, fct, time, ssGrp);
440 else
441 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}
448
449
450
451template <typename TGridFunction>
452void InterpolateInner(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
453 SmartPtr<TGridFunction> spGridFct, const char* cmp,
454 number time)
455{InterpolateInner(spInterpolFunction, spGridFct, cmp, NULL, time);}
456template <typename TGridFunction>
457void InterpolateInner(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
458 SmartPtr<TGridFunction> spGridFct, const char* cmp,
459 const char* subsets)
460{InterpolateInner(spInterpolFunction, spGridFct, cmp, subsets, 0.0);}
461template <typename TGridFunction>
462void InterpolateInner(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
463 SmartPtr<TGridFunction> spGridFct, const char* cmp)
464{InterpolateInner(spInterpolFunction, spGridFct, cmp, NULL, 0.0);}
465
467// const data
469
470template <typename TGridFunction>
471void 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 SmartPtr<UserData<number, dim> > sp =
477 make_sp(new ConstUserNumber<dim>(val));
478 InterpolateInner(sp, spGridFct, cmp, subsets, time);
479}
480template <typename TGridFunction>
481void InterpolateInner(number val,
482 SmartPtr<TGridFunction> spGridFct, const char* cmp,
483 number time)
484{InterpolateInner(val, spGridFct, cmp, NULL, time);}
485template <typename TGridFunction>
486void InterpolateInner(number val,
487 SmartPtr<TGridFunction> spGridFct, const char* cmp,
488 const char* subsets)
489{InterpolateInner(val, spGridFct, cmp, subsets, 0.0);}
490template <typename TGridFunction>
491void InterpolateInner(number val,
492 SmartPtr<TGridFunction> spGridFct, const char* cmp)
493{InterpolateInner(val, spGridFct, cmp, NULL, 0.0);}
494
496// lua data
498
499#ifdef UG_FOR_LUA
500template <typename TGridFunction>
501void 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 SmartPtr<UserData<number, dim> > sp =
507 LuaUserDataFactory<number, dim>::create(LuaFunction);
508 InterpolateInner(sp, spGridFct, cmp, subsets, time);
509}
510template <typename TGridFunction>
511void InterpolateInner(const char* LuaFunction,
512 SmartPtr<TGridFunction> spGridFct, const char* cmp,
513 number time)
514{InterpolateInner(LuaFunction, spGridFct, cmp, NULL, time);}
515template <typename TGridFunction>
516void InterpolateInner(const char* LuaFunction,
517 SmartPtr<TGridFunction> spGridFct, const char* cmp,
518 const char* subsets)
519{InterpolateInner(LuaFunction, spGridFct, cmp, subsets, 0.0);}
520template <typename TGridFunction>
521void InterpolateInner(const char* LuaFunction,
522 SmartPtr<TGridFunction> spGridFct, const char* cmp)
523{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
SmartPtr
Definition smart_pointer.h:108
ug::ConstUserNumber
constant scalar user data
Definition const_user_data.h:153
ug::CplUserData
Type based UserData.
Definition user_data.h:501
ug::LFEID
Identifier for Local Finite Elements.
Definition local_finite_element_id.h:98
ug::LFEID::LAGRANGE
@ LAGRANGE
Definition local_finite_element_id.h:104
ug::LuaUserDataFactory::create
static SmartPtr< LuaUserData< TData, dim, TRet > > create(const std::string &name)
Definition lua_user_data.h:223
ug::MathVector
a mathematical Vector with N entries.
Definition math_vector.h:97
ug::SubsetGroup
Group of subsets.
Definition subset_group.h:51
ug::SubsetGroup::size
size_t size() const
number of subsets in this group
Definition subset_group.h:122
ug::SubsetGroup::add_all
void add_all()
select all subsets of underlying subset handler
Definition subset_group.cpp:133
ug::SubsetGroup::dim
int dim(size_t i) const
dimension of subset
Definition subset_group.cpp:237
ug::SubsetGroup::add
void add(int si)
adds a subset by number to this group
Definition subset_group.cpp:64
ug::UserData
Type based UserData.
Definition user_data.h:143
const_user_data.h
dof_position_util.h
domain_util.h
ug::PST_CONSISTENT
@ PST_CONSISTENT
Definition parallel_storage_type.h:68
UG_CATCH_THROW
#define UG_CATCH_THROW(msg)
Definition error.h:64
UG_THROW
#define UG_THROW(msg)
Definition error.h:57
number
double number
Definition types.h:124
groups_util.h
local_finite_element_provider.h
lua_user_data.h
ug
the ug namespace
ug::InterpolateInnerDiff
void InterpolateInnerDiff(SmartPtr< UserData< number, TGridFunction::dim > > spInterpolFunction, SmartPtr< TGridFunction > spGridFct, const char *cmp, const char *subsets, number time, const SmartPtr< CplUserData< MathVector< TGridFunction::dim >, TGridFunction::dim > > m_diff_pos)
interpolates a function on a subset
Definition interpolate_inner.h:339
ug::ROID_TETRAHEDRON
@ ROID_TETRAHEDRON
Definition grid_base_objects.h:80
ug::InterpolateOnElementsInner
void InterpolateOnElementsInner(SmartPtr< UserData< number, TGridFunction::dim > > spInterpolFunction, SmartPtr< TGridFunction > spGridFct, size_t fct, int si, number time, const MathVector< TGridFunction::dim > diff_pos)
Definition interpolate_inner.h:68
ug::TokenizeString
void TokenizeString(const string &str, vector< string > &vToken, const char delimiter)
Definition string_util.cpp:56
ug::VOLUME
@ VOLUME
Definition grid_base_objects.h:63
ug::VERTEX
@ VERTEX
Definition grid_base_objects.h:60
ug::EDGE
@ EDGE
Definition grid_base_objects.h:61
ug::FACE
@ FACE
Definition grid_base_objects.h:62
ug::InterpolateInner
void InterpolateInner(SmartPtr< UserData< number, TGridFunction::dim > > spInterpolFunction, SmartPtr< TGridFunction > spGridFct, const char *cmp, const char *subsets, number time)
interpolates a function on a subset
Definition interpolate_inner.h:404
ug::DoFRef
number & DoFRef(TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
Definition multi_index.h:276
ug::DIM_SUBSET_EMPTY_GRID
@ DIM_SUBSET_EMPTY_GRID
Definition subset_dim_util.h:92
reference_mapping.h
smart_pointer.h
make_sp
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition smart_pointer.h:836
value_type
T value_type
Definition sparsematrix_interface.h:2
subset_group.h