Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
data_evaluator_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2011-2017: G-CSC, Goethe University Frankfurt
3 * Authors: Andreas Vogel, Arne Naegel
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__SPATIAL_DISC__DATA_EVALUATOR_IMPL__
34#define __H__UG__LIB_DISC__SPATIAL_DISC__DATA_EVALUATOR_IMPL__
35
36namespace ug{
37
38template <typename TDomain, typename TElemDisc>
40DataEvaluatorBase(int discPart,
41 const std::vector<TElemDisc*>& vElemDisc,
43 const bool bNonRegularGrid,
44 LocalVectorTimeSeries* pLocTimeSeries,
45 const std::vector<number>* pvScaleMass,
46 const std::vector<number>* pvScaleStiff)
47 : m_spFctPattern(fctPat)
48{
49// remember infos
50 m_discPart = discPart;
51 m_pLocTimeSeries = pLocTimeSeries;
52 m_bNeedLocTimeSeries = false; // initially
53 m_bUseHanging = false; // initially
54
55 m_vElemDisc[PT_ALL] = vElemDisc;
56
57// create FunctionIndexMapping for each Disc
58 for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
59 {
60 // get elem disc
61 TElemDisc* disc = m_vElemDisc[PT_ALL][i];
63 // handle time dependency
64 if(pLocTimeSeries != NULL && pvScaleMass != NULL && pvScaleStiff != NULL){
65 disc->set_time_dependent(*pLocTimeSeries, *pvScaleMass, *pvScaleStiff);
66 }
67 else if(pLocTimeSeries != NULL){
68 disc->set_time_dependent(*pLocTimeSeries, std::vector<number>(), std::vector<number>());
69 }
70 else{
71 disc->set_time_independent();
72 }
73
74 // checks
75 disc->check_setup(bNonRegularGrid);
76
77 // cache time dependency
78 m_bNeedLocTimeSeries |= disc->local_time_series_needed();
79
80 // cache grid type required
81 if(bNonRegularGrid)
82 m_bUseHanging |= disc->use_hanging();
83
84 // sort by process type
85 ProcessType process;
86 if(!disc->is_time_dependent()) process = PT_STATIONARY;
87 else process = PT_INSTATIONARY;
88 m_vElemDisc[process].push_back(disc);
89 }
90}
91
92
94// DataEvaluatorBase Setup
96
98
99template <typename TDomain, typename TElemDisc>
102 std::vector<SmartPtr<ICplUserData<dim> > >& vTryingToAdd)
103{
104// if empty, we're done
105 if(vTryingToAdd.empty()) return;
106
107// search for element in already scheduled data
108 typename std::vector<SmartPtr<ICplUserData<dim> > >::iterator it, itEnd;
109 it = find(vEvalData.begin(), vEvalData.end(), vTryingToAdd.back());
110
111// if found, skip this data
112 if(it != vEvalData.end())
113 {
114 vTryingToAdd.pop_back();
115 return;
116 }
117
118// search if element already contained in list. Then, the element
119// did start the adding procedure before and a circle dependency
120// is found
121 itEnd = vTryingToAdd.end(); itEnd--;
122 it = find(vTryingToAdd.begin(), itEnd, *itEnd);
123
124// if found, return error of circle dependency
125 if(it != itEnd)
126 UG_THROW("DataEvaluatorBase::add_data_to_eval_data:"
127 " Circle dependency of data detected for UserData.");
128
129// add all dependent datas
130 SmartPtr<ICplUserData<dim> > data = vTryingToAdd.back();
131 for(size_t i = 0; i < data->num_needed_data(); ++i)
132 {
133 // add each data separately
134 vTryingToAdd.push_back(data->needed_data(i));
135 add_data_to_eval_data(vEvalData, vTryingToAdd);
136 }
137
138// add this data to the evaluation list
139 vEvalData.push_back(data);
140
141// pop last one, since now added to eval list
142 if(!vTryingToAdd.empty())
143 vTryingToAdd.pop_back();
144}
145
146template <typename TDomain, typename TElemDisc>
148{
149 clear_extracted_data_and_mappings();
150
151// queue for all user data needed
152 std::vector<SmartPtr<ICplUserData<dim> > > vEvalData;
153 std::vector<SmartPtr<ICplUserData<dim> > > vTryingToAdd;
154
155// In the next loop we extract all needed UserData:
156// We only process the DataImport if there has been set data to the import
157// since otherwise no evaluation is needed.
158// If there is data given, we get the connected UserData and add it to the vector
159// of EvaluationData. This simply adds the UserData to the queue for UserData, if
160// the data does not depend on other Data. But if the UserData itself has
161// dependencies to other UserData, this data is added first (in a recursive
162// process). Of course, no circle dependency between UserData is allowed.
163
164// In the same loop over the data imports, we schedule the DataImports for
165// evaluation and compute the correct FunctionMapping for the linearization
166// of the defect and the Data, the Import is connected to:
167// If the UserData does not depend on the primary unknowns, we're done. Else
168// we have to setup the Function mappings between the common function group
169// and the DataImport-FunctionGroup. This is simply the same function map as
170// for the element discretization, since the DataImport depends by definition
171// from and only from the primary variables of its associated IElemDisc.
172
173// loop elem discs
174 for(size_t d = 0; d < m_vElemDisc[PT_ALL].size(); ++d)
175 {
176 TElemDisc* disc = m_vElemDisc[PT_ALL][d];
177
178 // loop imports
179 for(size_t i = 0; i < disc->num_imports(); ++i)
180 {
181 // get import
182 IDataImport<dim>* iimp = &(disc->get_import(i));
183
184 // skip non-given data (no need for evaluation)
185 if(!iimp->data_given()) continue;
186
187 // check part
188 if( !(iimp->part() & discPart) ) continue;
189
190 // check correct process type
191 if(iimp->part() == MASS)
192 if(!disc->is_time_dependent()) continue;
193
194 // push export on stack of needed data
195 vTryingToAdd.push_back(iimp->data());
196
197 // add data and all dependency to evaluation list
198 try{
199 add_data_to_eval_data(vEvalData, vTryingToAdd);
200 }
201 UG_CATCH_THROW("DataEvaluatorBase:"
202 " Circle dependency of data detected for UserData.");
203
204 // check that queue is empty now, else some internal error occured
205 if(!vTryingToAdd.empty())
206 UG_THROW("DataEvaluatorBase:"
207 " Internal Error, UserData queue not empty after adding.");
208
209 // done if and only if zero-derivative
210 if(iimp->zero_derivative()) continue;
211
212 // remember Import
213 ProcessType process;
214 if(!disc->is_time_dependent()) process = PT_STATIONARY;
215 else process = PT_INSTATIONARY;
216
217 m_vImport[PT_ALL][iimp->part()].push_back(iimp);
218 m_vImport[process][iimp->part()].push_back(iimp);
219 }
220 }
221
222// Now, we have processed all imports, that must be evaluated and have a long
223// vector of UserData that is connected to those imports. The UserData is already
224// sorted in this way: Data that depends on other data appears after the data
225// it depends on. This is important since we will schedule now the data for
226// evaluation and the data, that is needed by other data, will be computed
227// first. In addition, the data linker have to update their FunctionGroup and
228// must be sure that the data they depend on has already a correct FunctionGroup
229// set. This all is ensured by the (already produced) correct ordering.
230//
231// In the next loop we process all UserData, that will be evaluated during
232// assembling (i.e. is connected to an Import). First, we check if the data
233// is constant. If so simply add it the the Constant Data vector; nothing more
234// has to be done here. Else we check if the data depends on the primary
235// unknowns. If this is not the case, the UserData must be a Position-dependent
236// data, but not constant. Thus, schedule it at the Position Data vector.
237// If the data depends on the primary unknowns we must proceed as follows:
238// First, we update the FunctionGroup of the Data, since it could be a linker
239// and having an incorrect FunctionGroup (iff the FunctionGroup of the data
240// the linker depends on has been changed). Then we create the function
241// mapping between the functions the linker depends on and the common Function
242// Group.
243
244// loop all needed user data and group it
245 for(size_t i = 0; i < vEvalData.size(); ++i)
246 {
247 // get the user data
248 SmartPtr<ICplUserData<dim> > ipData = vEvalData[i];
249
250 // update function pattern (this will update functionGroups and Map of Data)
251 try{
252 ipData->set_function_pattern(m_spFctPattern);
253 }
254 UG_CATCH_THROW("DataEvaluatorBase: Cannot set FunctionPattern to UserData.");
255
256 // sort data into const and non-solution dependent
257 if(ipData->constant()) {m_vConstData.push_back(ipData); continue;}
258 if(ipData->zero_derivative()){m_vPosData.push_back(ipData); continue;}
259
260 // save as dependent data
261 m_vDependentData.push_back(ipData);
262 }
263
264// Handle time dependency
265// NOTE: constant data is not processed.
266 if(m_pLocTimeSeries != NULL){
267 for(size_t i = 0; i < m_vPosData.size(); ++i)
268 m_vPosData[i]->set_times(m_pLocTimeSeries->times());
269 for(size_t i = 0; i < m_vDependentData.size(); ++i)
270 m_vDependentData[i]->set_times(m_pLocTimeSeries->times());
271 }
272
273// NOTE: constant data is not processed, since constant == independent of si
274 for(size_t i = 0; i < m_vPosData.size(); ++i)
275 m_vPosData[i]->set_subset(subsetIndex);
276 for(size_t i = 0; i < m_vDependentData.size(); ++i)
277 m_vDependentData[i]->set_subset(subsetIndex);
278}
279
280template <typename TDomain, typename TElemDisc>
282{
283 for(int type = 0; type < MAX_PROCESS; ++type){
284 m_vImport[type][MASS].clear();
285 m_vImport[type][STIFF].clear();
286 m_vImport[type][RHS].clear();
287 }
288
289 m_vConstData.clear();
290 m_vPosData.clear();
291 m_vDependentData.clear();
292}
293
294
295
296template <typename TDomain, typename TElemDisc>
298{
299// remove ip series for all used UserData
300 for(size_t i = 0; i < m_vConstData.size(); ++i) m_vConstData[i]->clear();
301 for(size_t i = 0; i < m_vPosData.size(); ++i) m_vPosData[i]->clear();
302 for(size_t i = 0; i < m_vDependentData.size(); ++i) m_vDependentData[i]->clear();
303}
304
305
307// prepare / finish
309
310template <typename TDomain, typename TElemDisc>
312{
313 for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
314 m_vElemDisc[PT_ALL][i]->set_time_point(timePoint);
315
316 // NOTE: constant data is not processed.
317 for(size_t i = 0; i < m_vPosData.size(); ++i)
318 m_vPosData[i]->set_time_point(timePoint);
319 for(size_t i = 0; i < m_vDependentData.size(); ++i)
320 m_vDependentData[i]->set_time_point(timePoint);
321}
322
323
324
326// Error estimator's routines
328
329template <typename TDomain, typename TElemDisc>
332{
333// prepare loop (elem disc set local ip series here)
334 try{
335 for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
336 m_vElemDisc[PT_ALL][i]->do_prep_err_est_elem_loop(id, si);
337 }
338 UG_CATCH_THROW("DataEvaluatorBase::prepare_err_est_elem_loop: "
339 "Cannot prepare element loop.");
340
341// extract data imports and user data
342 try{
343 extract_imports_and_userdata(si, m_discPart);
344 }
345 UG_CATCH_THROW("DataEvaluatorBase::prepare_err_est_elem_loop: "
346 "Cannot extract imports and userdata.");
347
348// check setup of imports
349 try{
350 for(size_t i = 0; i < m_vImport[PT_ALL][MASS].size(); ++i)
351 m_vImport[PT_ALL][MASS][i]->check_setup();
352 for(size_t i = 0; i < m_vImport[PT_ALL][STIFF].size(); ++i)
353 m_vImport[PT_ALL][STIFF][i]->check_setup();
354 for(size_t i = 0; i < m_vImport[PT_ALL][RHS].size(); ++i)
355 m_vImport[PT_ALL][RHS][i]->check_setup();
356 }
357 UG_CATCH_THROW("DataEvaluatorBase::prepare_err_est_elem_loop: Import not correctly implemented.");
358
359// prepare and check dependent data
360 try{
361 for(size_t i = 0; i < m_vDependentData.size(); ++i){
362 m_vDependentData[i]->check_setup();
363 }
364 }
365 UG_CATCH_THROW("DataEvaluatorBase::prepare_err_est_elem_loop: Dependent UserData "
366 " (e.g. Linker or Export) is not ready for evaluation.");
367
368// evaluate constant data
369 for(size_t i = 0; i < m_vConstData.size(); ++i)
370 m_vConstData[i]->compute((LocalVector*)NULL, NULL, NULL, false);
371}
372
373template <typename TDomain, typename TElemDisc>
375{
376// finish each elem error estimator disc
377 try{
378 for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
379 m_vElemDisc[PT_ALL][i]->do_fsh_err_est_elem_loop();
380 }
381 UG_CATCH_THROW("DataEvaluatorBase::finish_err_est_elem_loop: Cannot finish error estimator element loop");
382
383// clear positions at user data
384 clear_positions_in_user_data();
385}
386
387template <typename TDomain, typename TElemDisc>
389(LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[],
390const LocalIndices& ind, bool bDeriv)
391{
392 try
393 {
394 for (size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
395 m_vElemDisc[PT_ALL][i]->do_prep_err_est_elem(u, elem, vCornerCoords);
396 }
397 UG_CATCH_THROW("DataEvaluatorBase::compute_elem_err_est: Cannot prepare element.");
398
399// evaluate position data
400 for (size_t i = 0; i < m_vPosData.size(); ++i)
401 m_vPosData[i]->compute(&u, elem, vCornerCoords, false);
402
403// process dependent data:
404// We can not simply compute exports first, then Linker, because an export
405// itself could depend on other data if implemented somehow in the IElemDisc
406// (e.g. using data from some DataImport). Thus, we have to loop the sorted
407// vector of all dependent data (that is correctly sorted the way that always
408// needed data has previously computed).
409
410// compute the data
411 try
412 {
413 for (size_t i = 0; i < m_vDependentData.size(); ++i)
414 {
415 u.access_by_map(m_vDependentData[i]->map());
416 m_vDependentData[i]->compute(&u, elem, vCornerCoords, false);
417 }
418 }
419 UG_CATCH_THROW("DataEvaluatorBase::compute_elem_err_est: Cannot compute data for Export or Linker.");
420}
421
422
423template <typename TDomain, typename TElemDisc>
426( LocalVector& u,
427 GridObject* elem,
428 const MathVector<dim> vCornerCoords[],
429 const LocalIndices& ind,
430 const number scaleMass,
431 const number scaleStiff)
432{
433 UG_ASSERT(m_discPart & STIFF, "Using compute_err_est_A_elem, but not STIFF requested.");
434 try
435 {
436 for (std::size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
437 m_vElemDisc[PT_ALL][i]->do_compute_err_est_A_elem(u, elem, vCornerCoords, scaleStiff);
438 }
439 UG_CATCH_THROW("DataEvaluatorBase::compute_err_est_A_elem: Cannot assemble stiffness part of error estimator");
440}
441
442template <typename TDomain, typename TElemDisc>
445( LocalVector& u,
446 GridObject* elem,
447 const MathVector<dim> vCornerCoords[],
448 const LocalIndices& ind,
449 const number scaleMass,
450 const number scaleStiff)
451{
452 UG_ASSERT(m_discPart & MASS, "Using compute_err_est_M_elem, but not MASS requested.");
453 try
454 {
455 for (std::size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
456 m_vElemDisc[PT_ALL][i]->do_compute_err_est_M_elem(u, elem, vCornerCoords, scaleMass);
457 }
458 UG_CATCH_THROW("DataEvaluatorBase::compute_err_est_A_elem: Cannot assemble stiffness part of error estimator");
459}
460
461template <typename TDomain, typename TElemDisc>
464( GridObject* elem,
465 const MathVector<dim> vCornerCoords[],
466 const LocalIndices& ind,
467 const number scaleMass,
468 const number scaleStiff)
469{
470 UG_ASSERT(m_discPart & RHS, "Using compute_err_est_rhs_elem, but not RHS requested.");
471 try
472 {
473 for (std::size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
474 m_vElemDisc[PT_ALL][i]->do_compute_err_est_rhs_elem(elem, vCornerCoords, scaleStiff);
475 }
476 UG_CATCH_THROW("DataEvaluatorBase::compute_err_est_rhs_elem: Cannot assemble rhs part of error estimator");
477}
478
479
480} // namespace ug
481
482
483#endif
Definition smart_pointer.h:296
Definition smart_pointer.h:108
bool m_bUseHanging
flag if hanging nodes are used
Definition data_evaluator.h:128
void prepare_err_est_elem_loop(const ReferenceObjectID id, int si)
prepares the element loop for all IElemDiscs for the computation of the error estimator
Definition data_evaluator_impl.h:331
bool m_bNeedLocTimeSeries
flag indicating if any elem disc needs local time series
Definition data_evaluator.h:131
LocalVectorTimeSeries * m_pLocTimeSeries
local time series (non-const since mapping may change)
Definition data_evaluator.h:134
std::vector< TElemDisc * > m_vElemDisc[MAX_PROCESS]
elem disc data
Definition data_evaluator.h:122
void clear_positions_in_user_data()
clears all requested positions in user data
Definition data_evaluator_impl.h:297
void clear_extracted_data_and_mappings()
clears imports and user data and mappings betweem commonFctGrp and local
Definition data_evaluator_impl.h:281
DataEvaluatorBase(int discPart, const std::vector< TElemDisc * > &vElemDisc, ConstSmartPtr< FunctionPattern > fctPat, const bool bNonRegularGrid, LocalVectorTimeSeries *locTimeSeries=NULL, const std::vector< number > *vScaleMass=NULL, const std::vector< number > *vScaleStiff=NULL)
sets the elem discs to evaluate
Definition data_evaluator_impl.h:40
int m_discPart
disc part needed (MASS and/or STIFF and/or RHS)
Definition data_evaluator.h:119
void add_data_to_eval_data(std::vector< SmartPtr< ICplUserData< dim > > > &vEvalData, std::vector< SmartPtr< ICplUserData< dim > > > &vTryingToAdd)
tries to add the last entry of vTryingToAdd to the eval data
Definition data_evaluator_impl.h:101
void prepare_err_est_elem(LocalVector &u, GridObject *elem, const MathVector< dim > vCornerCoords[], const LocalIndices &ind, bool bDeriv=false)
prepares the element for all IElemDiscs
Definition data_evaluator_impl.h:389
void set_time_point(const size_t timePoint)
sets the time point for data evaluation
Definition data_evaluator_impl.h:311
void extract_imports_and_userdata(int subsetIndex, int discPart)
extracts imports and userdata from IElemDiscs
Definition data_evaluator_impl.h:147
void compute_err_est_rhs_elem(GridObject *elem, const MathVector< dim > vCornerCoords[], const LocalIndices &ind, const number scaleMass=(number) 1.0, const number scaleStiff=(number) 1.0)
compute contributions of the local error indicators in one element for all IElemDiscs
Definition data_evaluator_impl.h:464
void compute_err_est_M_elem(LocalVector &u, GridObject *elem, const MathVector< dim > vCornerCoords[], const LocalIndices &ind, const number scaleMass=(number) 1.0, const number scaleStiff=(number) 1.0)
compute contributions of the local error indicators in one element for all IElemDiscs
Definition data_evaluator_impl.h:445
void finish_err_est_elem_loop()
finishes the error estimator element loop for all IElemDiscs
Definition data_evaluator_impl.h:374
void compute_err_est_A_elem(LocalVector &u, GridObject *elem, const MathVector< dim > vCornerCoords[], const LocalIndices &ind, const number scaleMass=(number) 1.0, const number scaleStiff=(number) 1.0)
compute contributions of the local error indicators in one element for all IElemDiscs
Definition data_evaluator_impl.h:426
The base class for all geometric objects, such as vertices, edges, faces, volumes,...
Definition grid_base_objects.h:157
Base class for UserData.
Definition user_data.h:260
Base class for data import.
Definition data_import.h:72
DiscPart part() const
returns if import is located in mass part (for time dependent problems)
Definition data_import.h:104
bool zero_derivative() const
returns if data depends on unknown functions
Definition data_import.h:119
virtual bool data_given() const =0
returns if data is set
virtual SmartPtr< ICplUserData< dim > > data()=0
returns the connected user data
Definition local_algebra.h:50
Definition local_algebra.h:198
void access_by_map(const FunctionIndexMapping &funcMap)
access only part of the functions using mapping (restrict functions)
Definition local_algebra.h:306
time series of local vectors
Definition solution_time_series.h:167
a mathematical Vector with N entries.
Definition math_vector.h:97
function util rates kinetic compute(ConvRateSetup)
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#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
IndexLayout::Interface::iterator find(IndexLayout::Interface &interface, size_t i)
Definition parallel_index_layout.h:77
ReferenceObjectID
these ids are used to identify the shape of a geometric object.
Definition grid_base_objects.h:74
ProcessType
Definition data_evaluator.h:42
@ MAX_PROCESS
Definition data_evaluator.h:42
@ PT_INSTATIONARY
Definition data_evaluator.h:42
@ PT_ALL
Definition data_evaluator.h:42
@ PT_STATIONARY
Definition data_evaluator.h:42
@ STIFF
Definition data_import.h:46
@ RHS
Definition data_import.h:47
@ MASS
Definition data_import.h:45