Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
max_error.h
Go to the documentation of this file.
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"
42
44
50
51#ifdef UG_FOR_LUA
53#endif
54
55namespace ug{
56
58// Compute max error on Vertices only
60
74template <typename TGridFunction>
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 = spGridFct->domain()->position_accessor();
88
89 std::vector<DoFIndex> ind;
90 typename TGridFunction::template dim_traits<0>::const_iterator iterEnd, iter;
91
92 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 if(!spGridFct->is_def_in_subset(fct, si)) continue;
99
100 // iterate over all elements
101 iterEnd = spGridFct->template end<Vertex>(si);
102 iter = spGridFct->template begin<Vertex>(si);
103 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 (*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 for(size_t i = 0; i < ind.size(); ++i)
120 {
121 // compute error
122 number localError = std::abs(DoFRef((*spGridFct), ind[i])-val);
123 if (localError>globalMaxError) globalMaxError=localError;
124 }
125 }
126 }
127}
128
130// Maximum error on Elements
132
143template <typename TElem, typename TGridFunction>
145 number& globalMaxError,
147 SmartPtr<TGridFunction> spGridFct, size_t fct, int si, number time)
148{
149// get 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 iterEnd = spGridFct->template end<TElem>(si);
164 iter = spGridFct->template begin<TElem>(si);
165
166// check if something to do:
167 if(iter == iterEnd) return;
168
169// id of shape functions used
170 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 const size_t nsh = trialSpace.num_sh();
178
179// load local positions of dofs for the trial space on element
180 std::vector<MathVector<dim> > loc_pos(nsh);
181 for(size_t i = 0; i < nsh; ++i)
182 if(!trialSpace.position(i, loc_pos[i]))
183 UG_THROW("MaxErrorOnElem: Cannot find meaningful"
184 " local positions of dofs.");
185
186// create a reference mapping
188
189// iterate over all elements
190 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 CollectCornerCoordinates(vCorner, *elem, *spGridFct->domain());
198
199 // update the reference mapping for the corners
200 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 if(ind.size() != nsh)
208 UG_THROW("MaxErrorOnElem: On subset "<<si<<": Number of shapes is "
209 <<nsh<<", but got "<<ind.size()<<" multi indices.");
210
211 // loop all dofs
212 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 mapping.local_to_global(glob_pos, loc_pos[i]);
219
220 // value at position
221 number val;
222 (*spInterpolFunction)(val, glob_pos, time, si);
223
224 // compute error
225 number localError = std::abs(DoFRef((*spGridFct), ind[i])-val);
226 if (localError>globalMaxError) globalMaxError=localError;
227 }
228 }
229}
230
241template <typename TGridFunction>
243 number& globalMaxError,
245 SmartPtr<TGridFunction> spGridFct,
246 size_t fct,
247 number time,
248 const SubsetGroup& ssGrp)
249{
250// loop subsets
251 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 if(!spGridFct->is_def_in_subset(fct, si)) continue;
258
259 // switch dimensions
260 try
261 {
262 const int dim = ssGrp.dim(i);
263
264 if(dim > TGridFunction::dim)
265 UG_THROW("MaxErrorOnElements: Dimension of subset is "<<dim<<", but "
266 " World Dimension is "<<TGridFunction::dim<<". Cannot interpolate this.");
267
268 switch(dim)
269 {
270 case DIM_SUBSET_EMPTY_GRID: break;
271 case 0: /* \TODO: do nothing may be wrong */ break;
272 case 1:
273 MaxErrorOnElements<RegularEdge, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
274 break;
275 case 2:
276 MaxErrorOnElements<Triangle, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
277 MaxErrorOnElements<Quadrilateral, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
278 break;
279 case 3:
280 MaxErrorOnElements<Tetrahedron, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
281 MaxErrorOnElements<Hexahedron, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
282 MaxErrorOnElements<Prism, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
283 MaxErrorOnElements<Pyramid, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
284 MaxErrorOnElements<Octahedron, TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, si, time);
285 break;
286 default: UG_THROW("MaxErrorOnElements: Dimension " <<dim<<
287 " not possible for world dim "<<3<<".");
288 }
289 }
290 UG_CATCH_THROW("MaxErrorOnElements: Cannot interpolate on elements.");
291 }
292}
293
295// MaxError routine
297
299
310template <typename TGridFunction>
313 SmartPtr<TGridFunction> spGridFct, const char* cmp,
314 const char* subsets, number time)
315{
316 number globalMaxError=0;
317// check, that values do not depend on a solution
318 if(spInterpolFunction->requires_grid_fct())
319 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 if(fct > spGridFct->num_fct())
328 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 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 SubsetGroup ssGrp(spGridFct->domain()->subset_handler());
341 if(subsets != NULL)
342 {
343 ssGrp.add(TokenizeString(subsets));
344 if(!bAllowManyfoldInterpolation)
346 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 ssGrp.add_all();
353 if(!bAllowManyfoldInterpolation)
355 }
356
357// forward
358 if(bUseP1Interpolation)
359 MaxErrorOnVertices<TGridFunction>(globalMaxError,spInterpolFunction, spGridFct, fct, time, ssGrp);
360 else
361 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 return globalMaxError;
368}
369
370template <typename TGridFunction>
372 SmartPtr<TGridFunction> spGridFct, const char* cmp,
373 number time)
374{
375 return MaxError(spInterpolFunction, spGridFct, cmp, NULL, time);
376}
377template <typename TGridFunction>
379 SmartPtr<TGridFunction> spGridFct, const char* cmp,
380 const char* subsets)
381{
382 return MaxError(spInterpolFunction, spGridFct, cmp, subsets, 0.0);
383}
384template <typename TGridFunction>
386 SmartPtr<TGridFunction> spGridFct, const char* cmp)
387{
388 return MaxError(spInterpolFunction, spGridFct, cmp, NULL, 0.0);
389}
390
392// const data
394
395template <typename TGridFunction>
397 SmartPtr<TGridFunction> spGridFct, const char* cmp,
398 const char* subsets, number time)
399{
400 static const int dim = TGridFunction::dim;
403 return MaxError(sp, spGridFct, cmp, subsets, time);
404}
405template <typename TGridFunction>
407 SmartPtr<TGridFunction> spGridFct, const char* cmp,
408 number time)
409{
410 return MaxError(val, spGridFct, cmp, NULL, time);
411}
412template <typename TGridFunction>
414 SmartPtr<TGridFunction> spGridFct, const char* cmp,
415 const char* subsets)
416{
417 return MaxError(val, spGridFct, cmp, subsets, 0.0);
418}
419template <typename TGridFunction>
421 SmartPtr<TGridFunction> spGridFct, const char* cmp)
422{
423 return MaxError(val, spGridFct, cmp, NULL, 0.0);
424}
425
427// lua data
429
430#ifdef UG_FOR_LUA
431template <typename TGridFunction>
432number 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;
439 return MaxError(sp, spGridFct, cmp, subsets, time);
440}
441template <typename TGridFunction>
442number MaxError(const char* LuaFunction,
443 SmartPtr<TGridFunction> spGridFct, const char* cmp,
444 number time)
445{
446 return MaxError(LuaFunction, spGridFct, cmp, NULL, time);
447}
448
449template <typename TGridFunction>
450number MaxError(const char* LuaFunction,
451 SmartPtr<TGridFunction> spGridFct, const char* cmp,
452 const char* subsets)
453{
454 return MaxError(LuaFunction, spGridFct, cmp, subsets, 0.0);
455}
456template <typename TGridFunction>
457number MaxError(const char* LuaFunction,
458 SmartPtr<TGridFunction> spGridFct, const char* cmp)
459{
460 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__*/
Definition smart_pointer.h:108
constant scalar user data
Definition const_user_data.h:153
virtual bool position(size_t i, MathVector< TDim > &pos) const =0
local position of DoF i
Identifier for Local Finite Elements.
Definition local_finite_element_id.h:98
@ LAGRANGE
Definition local_finite_element_id.h:104
virtual size_t num_sh() const
Definition local_dof_set.cpp:46
virtual base class for local shape function sets
Definition local_shape_function_set.h:70
static SmartPtr< LuaUserData< TData, dim, TRet > > create(const std::string &name)
Definition lua_user_data.h:223
Definition reference_mapping.h:65
void local_to_global(MathVector< worldDim > &globPos, const MathVector< dim > &locPos) const
map local coordinate to global coordinate
void update(const MathVector< worldDim > *vCornerCoord)
refresh mapping for new set of corners
Group of subsets.
Definition subset_group.h:51
size_t size() const
number of subsets in this group
Definition subset_group.h:122
void add_all()
select all subsets of underlying subset handler
Definition subset_group.cpp:133
int dim(size_t i) const
dimension of subset
Definition subset_group.cpp:237
void add(int si)
adds a subset by number to this group
Definition subset_group.cpp:64
Type based UserData.
Definition user_data.h:143
Base-class for all vertex-types.
Definition grid_base_objects.h:231
@ PST_CONSISTENT
Definition parallel_storage_type.h:68
void CollectCornerCoordinates(std::vector< typename TAAPos::ValueType > &vCornerCoordsOut, const TElem &elem, const TAAPos &aaPos, bool clearContainer=true)
returns the corner coordinates of a geometric object
Definition domain_util_impl.h:75
#define UG_CATCH_THROW(msg)
Definition error.h:64
#define UG_THROW(msg)
Definition error.h:57
double number
Definition types.h:124
the ug namespace
void MaxErrorOnVertices(number &globalMaxError, SmartPtr< UserData< number, TGridFunction::dim > > spInterpolFunction, SmartPtr< TGridFunction > spGridFct, size_t fct, number time, const SubsetGroup &ssGrp)
Definition max_error.h:75
ReferenceObjectID
these ids are used to identify the shape of a geometric object.
Definition grid_base_objects.h:74
void TokenizeString(const string &str, vector< string > &vToken, const char delimiter)
Definition string_util.cpp:56
bool SameDimensionsInAllSubsets(const SubsetGroup &subsetGroup)
Definition subset_group.cpp:299
void RemoveLowerDimSubsets(SubsetGroup &subsetGroup)
Definition subset_group.cpp:315
void MaxErrorOnElements(number &globalMaxError, SmartPtr< UserData< number, TGridFunction::dim > > spInterpolFunction, SmartPtr< TGridFunction > spGridFct, size_t fct, int si, number time)
Definition max_error.h:144
number MaxError(SmartPtr< UserData< number, TGridFunction::dim > > spInterpolFunction, SmartPtr< TGridFunction > spGridFct, const char *cmp, const char *subsets, number time)
computes maximum error of a grid function on a subset
Definition max_error.h:311
number & DoFRef(TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
Definition multi_index.h:276
@ DIM_SUBSET_EMPTY_GRID
Definition subset_dim_util.h:92
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition smart_pointer.h:836
traits for reference elements
Definition reference_element_traits.h:48