Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Authors: Andreas Vogel, Christian Wehner
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 : /*
34 : * an adaption of interpolate.h to maximum error computation
35 : */
36 :
37 : #ifndef __H__UG__LIB_DISC__MAX__ERROR__
38 : #define __H__UG__LIB_DISC__MAX__ERROR__
39 :
40 : #include "common/common.h"
41 : #include "common/util/smart_pointer.h"
42 :
43 : #include "lib_grid/tools/subset_group.h"
44 :
45 : #include "lib_disc/common/groups_util.h"
46 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
47 : #include "lib_disc/spatial_disc/user_data/const_user_data.h"
48 : #include "lib_disc/reference_element/reference_mapping.h"
49 : #include "lib_disc/domain_util.h"
50 :
51 : #ifdef UG_FOR_LUA
52 : #include "bindings/lua/lua_user_data.h"
53 : #endif
54 :
55 : namespace ug{
56 :
57 : ////////////////////////////////////////////////////////////////////////////////
58 : // Compute max error on Vertices only
59 : ////////////////////////////////////////////////////////////////////////////////
60 :
61 : /**
62 : * This function computes max error of a grid function on a vertex loop. Thus, it can only
63 : * be used if all degrees of freedom are located in the vertices only (e.g. P1
64 : * finite elements). In those cases it is faster than the element by element
65 : * loop.
66 : *
67 : * @param[in] globalMaxError reference to global max error
68 : * @param[in] spInterpolFunction data providing exact solution
69 : * @param[out] spGridFct grid function
70 : * @param[in] fct symbolic name of function component
71 : * @param[in] ssGrp subsets, where to compute max error
72 : * @param[in] time time point
73 : */
74 : template <typename TGridFunction>
75 0 : void MaxErrorOnVertices(number& globalMaxError,SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
76 : SmartPtr<TGridFunction> spGridFct,
77 : size_t fct,
78 : number time,
79 : const SubsetGroup& ssGrp)
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 : // get position accessor
86 : const typename domain_type::position_accessor_type& aaPos
87 0 : = spGridFct->domain()->position_accessor();
88 :
89 : std::vector<DoFIndex> ind;
90 : typename TGridFunction::template dim_traits<0>::const_iterator iterEnd, iter;
91 :
92 0 : for(size_t i = 0; i < ssGrp.size(); ++i)
93 : {
94 : // get subset index
95 : const int si = ssGrp[i];
96 :
97 : // skip if function is not defined in subset
98 0 : if(!spGridFct->is_def_in_subset(fct, si)) continue;
99 :
100 : // iterate over all elements
101 0 : iterEnd = spGridFct->template end<Vertex>(si);
102 0 : iter = spGridFct->template begin<Vertex>(si);
103 0 : for(; iter != iterEnd; ++iter)
104 : {
105 : // get element
106 : Vertex* vrt = *iter;
107 :
108 : // global position
109 : position_type glob_pos = aaPos[vrt];
110 :
111 : // value at position
112 : number val;
113 0 : (*spInterpolFunction)(val, glob_pos, time, si);
114 :
115 : // get multiindices of element
116 : spGridFct->dof_indices(vrt, fct, ind);
117 :
118 : // loop all dofs
119 0 : for(size_t i = 0; i < ind.size(); ++i)
120 : {
121 : // compute error
122 0 : number localError = std::abs(DoFRef((*spGridFct), ind[i])-val);
123 0 : if (localError>globalMaxError) globalMaxError=localError;
124 : }
125 : }
126 : }
127 0 : }
128 :
129 : ////////////////////////////////////////////////////////////////////////////////
130 : // Maximum error on Elements
131 : ////////////////////////////////////////////////////////////////////////////////
132 :
133 : /**
134 : * This function computes maximum error of a grid function on an element by element loop.
135 : *
136 : * @param[in] globalMaxError reference to global max error
137 : * @param[in] spInterpolFunction data providing exact solution
138 : * @param[out] spGridFct interpolated grid function
139 : * @param[in] fct symbolic name of function component
140 : * @param[in] si subset
141 : * @param[in] time time point
142 : */
143 : template <typename TElem, typename TGridFunction>
144 0 : void MaxErrorOnElements(
145 : number& globalMaxError,
146 : SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
147 : SmartPtr<TGridFunction> spGridFct, size_t fct, int si, number time)
148 : {
149 : // get reference element type
150 : typedef typename reference_element_traits<TElem>::reference_element_type
151 : ref_elem_type;
152 : const ReferenceObjectID roid = ref_elem_type::REFERENCE_OBJECT_ID;
153 :
154 : // dimension of reference element
155 : const int dim = ref_elem_type::dim;
156 :
157 : // domain type and position_type
158 : typedef typename TGridFunction::domain_type domain_type;
159 : typedef typename domain_type::position_type position_type;
160 :
161 : // get iterators
162 : typename TGridFunction::template traits<TElem>::const_iterator iterEnd, iter;
163 0 : iterEnd = spGridFct->template end<TElem>(si);
164 0 : iter = spGridFct->template begin<TElem>(si);
165 :
166 : // check if something to do:
167 0 : if(iter == iterEnd) return;
168 :
169 : // id of shape functions used
170 0 : LFEID id = spGridFct->local_finite_element_id(fct);
171 :
172 : // get trial space
173 : const LocalShapeFunctionSet<dim>& trialSpace =
174 : LocalFiniteElementProvider::get<dim>(roid, id);
175 :
176 : // number of dofs on element
177 0 : const size_t nsh = trialSpace.num_sh();
178 :
179 : // load local positions of dofs for the trial space on element
180 0 : std::vector<MathVector<dim> > loc_pos(nsh);
181 0 : for(size_t i = 0; i < nsh; ++i)
182 0 : if(!trialSpace.position(i, loc_pos[i]))
183 0 : UG_THROW("MaxErrorOnElem: Cannot find meaningful"
184 : " local positions of dofs.");
185 :
186 : // create a reference mapping
187 : ReferenceMapping<ref_elem_type, domain_type::dim> mapping;
188 :
189 : // iterate over all elements
190 0 : for( ; iter != iterEnd; ++iter)
191 : {
192 : // get element
193 : TElem* elem = *iter;
194 :
195 : // get all corner coordinates
196 : std::vector<position_type> vCorner;
197 0 : CollectCornerCoordinates(vCorner, *elem, *spGridFct->domain());
198 :
199 : // update the reference mapping for the corners
200 0 : mapping.update(&vCorner[0]);
201 :
202 : // get multiindices of element
203 : std::vector<DoFIndex> ind;
204 : spGridFct->dof_indices(elem, fct, ind);
205 :
206 : // check multi indices
207 0 : if(ind.size() != nsh)
208 0 : UG_THROW("MaxErrorOnElem: On subset "<<si<<": Number of shapes is "
209 : <<nsh<<", but got "<<ind.size()<<" multi indices.");
210 :
211 : // loop all dofs
212 0 : for(size_t i = 0; i < nsh; ++i)
213 : {
214 : // global position
215 : position_type glob_pos;
216 :
217 : // map local dof position to global position
218 0 : mapping.local_to_global(glob_pos, loc_pos[i]);
219 :
220 : // value at position
221 : number val;
222 0 : (*spInterpolFunction)(val, glob_pos, time, si);
223 :
224 : // compute error
225 0 : number localError = std::abs(DoFRef((*spGridFct), ind[i])-val);
226 0 : if (localError>globalMaxError) globalMaxError=localError;
227 : }
228 : }
229 0 : }
230 :
231 : /**
232 : * This function computes maximum error of grid function on an element by element loop.
233 : *
234 : * @param[in] globalMaxError reference to global max error
235 : * @param[in] spInterpolFunction data providing exact solution
236 : * @param[out] spGridFct grid function
237 : * @param[in] fct symbolic name of function component
238 : * @param[in] ssGrp subsets
239 : * @param[in] time time point
240 : */
241 : template <typename TGridFunction>
242 0 : void MaxErrorOnElements(
243 : number& globalMaxError,
244 : SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
245 : SmartPtr<TGridFunction> spGridFct,
246 : size_t fct,
247 : number time,
248 : const SubsetGroup& ssGrp)
249 : {
250 : // loop subsets
251 0 : for(size_t i = 0; i < ssGrp.size(); ++i)
252 : {
253 : // get subset index
254 : const int si = ssGrp[i];
255 :
256 : // skip if function is not defined in subset
257 0 : if(!spGridFct->is_def_in_subset(fct, si)) continue;
258 :
259 : // switch dimensions
260 : try
261 : {
262 0 : const int dim = ssGrp.dim(i);
263 :
264 0 : if(dim > TGridFunction::dim)
265 0 : UG_THROW("MaxErrorOnElements: Dimension of subset is "<<dim<<", but "
266 : " World Dimension is "<<TGridFunction::dim<<". Cannot interpolate this.");
267 :
268 0 : switch(dim)
269 : {
270 : case DIM_SUBSET_EMPTY_GRID: break;
271 : case 0: /* \TODO: do nothing may be wrong */ break;
272 : case 1:
273 0 : MaxErrorOnElements<RegularEdge, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
274 0 : break;
275 : case 2:
276 0 : MaxErrorOnElements<Triangle, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
277 0 : MaxErrorOnElements<Quadrilateral, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
278 0 : break;
279 : case 3:
280 0 : MaxErrorOnElements<Tetrahedron, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
281 0 : MaxErrorOnElements<Hexahedron, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
282 0 : MaxErrorOnElements<Prism, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
283 0 : MaxErrorOnElements<Pyramid, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
284 0 : MaxErrorOnElements<Octahedron, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
285 0 : break;
286 0 : default: UG_THROW("MaxErrorOnElements: Dimension " <<dim<<
287 : " not possible for world dim "<<3<<".");
288 : }
289 : }
290 0 : UG_CATCH_THROW("MaxErrorOnElements: Cannot interpolate on elements.");
291 : }
292 0 : }
293 :
294 : ////////////////////////////////////////////////////////////////////////////////
295 : // MaxError routine
296 : ////////////////////////////////////////////////////////////////////////////////
297 :
298 : /// computes maximum error of a grid function on a subset
299 : /**
300 : * This function computes the maximum error of a GridFunction. To evaluate the exact solution on every
301 : * point a functor must be passed.
302 : *
303 : * @param[in] spInterpolFunction data providing exact solution
304 : * @param[out] spGridFct interpolated grid function
305 : * @param[in] cmp symbolic name of function component
306 : * @param[in] subsets subsets (NULL = everywhere)
307 : * @param[in] time time point
308 : * @returns globalMaxError maximum norm of difference
309 : */
310 : template <typename TGridFunction>
311 0 : number MaxError(
312 : SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
313 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
314 : const char* subsets, number time)
315 : {
316 0 : number globalMaxError=0;
317 : // check, that values do not depend on a solution
318 0 : if(spInterpolFunction->requires_grid_fct())
319 0 : UG_THROW("MaxError: The interpolation values depend on a grid function."
320 : " This is not allowed in the current implementation. Use constant,"
321 : " lua-callback or vrl-callback user data only (even within linkers).");
322 :
323 : // get function id of name
324 : const size_t fct = spGridFct->fct_id_by_name(cmp);
325 :
326 : // check that function found
327 0 : if(fct > spGridFct->num_fct())
328 0 : UG_THROW("MaxError: Name of component '"<<cmp<<"' not found.");
329 :
330 : // check if fast P1 interpolation can be used
331 : // \TODO: This should be improved. Manifold admissible if space continuous
332 : bool bUseP1Interpolation = false;
333 0 : if(spGridFct->local_finite_element_id(fct).type() == LFEID::LAGRANGE &&
334 : spGridFct->local_finite_element_id(fct).order() == 1)
335 : bUseP1Interpolation = true;
336 : const bool bAllowManyfoldInterpolation =
337 : (spGridFct->local_finite_element_id(fct).type() == LFEID::LAGRANGE);
338 :
339 : // create subset group
340 0 : SubsetGroup ssGrp(spGridFct->domain()->subset_handler());
341 0 : if(subsets != NULL)
342 : {
343 0 : ssGrp.add(TokenizeString(subsets));
344 0 : if(!bAllowManyfoldInterpolation)
345 0 : if(!SameDimensionsInAllSubsets(ssGrp))
346 0 : UG_THROW("MaxError: Subsets '"<<subsets<<"' do not have same dimension."
347 : "Can not integrate on subsets of different dimensions.");
348 : }
349 : else
350 : {
351 : // add all subsets and remove lower dim subsets afterwards
352 0 : ssGrp.add_all();
353 0 : if(!bAllowManyfoldInterpolation)
354 0 : RemoveLowerDimSubsets(ssGrp);
355 : }
356 :
357 : // forward
358 0 : if(bUseP1Interpolation)
359 0 : MaxErrorOnVertices<TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, time, ssGrp);
360 : else
361 0 : MaxErrorOnElements<TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, time, ssGrp);
362 :
363 : // adjust parallel storage state
364 : #ifdef UG_PARALLEL
365 : spGridFct->set_storage_type(PST_CONSISTENT);
366 : #endif
367 0 : return globalMaxError;
368 : }
369 :
370 : template <typename TGridFunction>
371 0 : number MaxError(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
372 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
373 : number time)
374 : {
375 0 : return MaxError(spInterpolFunction, spGridFct, cmp, NULL, time);
376 : }
377 : template <typename TGridFunction>
378 0 : number MaxError(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
379 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
380 : const char* subsets)
381 : {
382 0 : return MaxError(spInterpolFunction, spGridFct, cmp, subsets, 0.0);
383 : }
384 : template <typename TGridFunction>
385 0 : number MaxError(SmartPtr<UserData<number, TGridFunction::dim> > spInterpolFunction,
386 : SmartPtr<TGridFunction> spGridFct, const char* cmp)
387 : {
388 0 : return MaxError(spInterpolFunction, spGridFct, cmp, NULL, 0.0);
389 : }
390 :
391 : ///////////////
392 : // const data
393 : ///////////////
394 :
395 : template <typename TGridFunction>
396 0 : number MaxError(number val,
397 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
398 : const char* subsets, number time)
399 : {
400 : static const int dim = TGridFunction::dim;
401 0 : SmartPtr<UserData<number, dim> > sp =
402 0 : make_sp(new ConstUserNumber<dim>(val));
403 0 : return MaxError(sp, spGridFct, cmp, subsets, time);
404 : }
405 : template <typename TGridFunction>
406 0 : number MaxError(number val,
407 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
408 : number time)
409 : {
410 0 : return MaxError(val, spGridFct, cmp, NULL, time);
411 : }
412 : template <typename TGridFunction>
413 0 : number MaxError(number val,
414 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
415 : const char* subsets)
416 : {
417 0 : return MaxError(val, spGridFct, cmp, subsets, 0.0);
418 : }
419 : template <typename TGridFunction>
420 0 : number MaxError(number val,
421 : SmartPtr<TGridFunction> spGridFct, const char* cmp)
422 : {
423 0 : return MaxError(val, spGridFct, cmp, NULL, 0.0);
424 : }
425 :
426 : ///////////////
427 : // lua data
428 : ///////////////
429 :
430 : #ifdef UG_FOR_LUA
431 : template <typename TGridFunction>
432 0 : number MaxError(const char* LuaFunction,
433 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
434 : const char* subsets, number time)
435 : {
436 : static const int dim = TGridFunction::dim;
437 0 : SmartPtr<UserData<number, dim> > sp =
438 : LuaUserDataFactory<number, dim>::create(LuaFunction);
439 0 : return MaxError(sp, spGridFct, cmp, subsets, time);
440 : }
441 : template <typename TGridFunction>
442 0 : number MaxError(const char* LuaFunction,
443 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
444 : number time)
445 : {
446 0 : return MaxError(LuaFunction, spGridFct, cmp, NULL, time);
447 : }
448 :
449 : template <typename TGridFunction>
450 0 : number MaxError(const char* LuaFunction,
451 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
452 : const char* subsets)
453 : {
454 0 : return MaxError(LuaFunction, spGridFct, cmp, subsets, 0.0);
455 : }
456 : template <typename TGridFunction>
457 0 : number MaxError(const char* LuaFunction,
458 : SmartPtr<TGridFunction> spGridFct, const char* cmp)
459 : {
460 0 : return MaxError(LuaFunction, spGridFct, cmp, NULL, 0.0);
461 : }
462 : #endif
463 :
464 :
465 : } // namespace ug
466 :
467 : #endif /*__H__UG__LIB_DISC__MAX__ERROR__*/
|
