Line data Source code
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 :
36 : namespace ug{
37 :
38 : template <typename TDomain, typename TElemDisc>
39 0 : DataEvaluatorBase<TDomain, TElemDisc>::
40 : DataEvaluatorBase(int discPart,
41 : const std::vector<TElemDisc*>& vElemDisc,
42 : ConstSmartPtr<FunctionPattern> fctPat,
43 : const bool bNonRegularGrid,
44 : LocalVectorTimeSeries* pLocTimeSeries,
45 : const std::vector<number>* pvScaleMass,
46 : const std::vector<number>* pvScaleStiff)
47 0 : : m_spFctPattern(fctPat)
48 : {
49 : // remember infos
50 0 : m_discPart = discPart;
51 0 : m_pLocTimeSeries = pLocTimeSeries;
52 0 : m_bNeedLocTimeSeries = false; // initially
53 0 : m_bUseHanging = false; // initially
54 :
55 0 : m_vElemDisc[PT_ALL] = vElemDisc;
56 :
57 : // create FunctionIndexMapping for each Disc
58 0 : for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
59 : {
60 : // get elem disc
61 0 : TElemDisc* disc = m_vElemDisc[PT_ALL][i];
62 :
63 : // handle time dependency
64 0 : if(pLocTimeSeries != NULL && pvScaleMass != NULL && pvScaleStiff != NULL){
65 0 : disc->set_time_dependent(*pLocTimeSeries, *pvScaleMass, *pvScaleStiff);
66 : }
67 0 : else if(pLocTimeSeries != NULL){
68 0 : disc->set_time_dependent(*pLocTimeSeries, std::vector<number>(), std::vector<number>());
69 : }
70 : else{
71 0 : disc->set_time_independent();
72 : }
73 :
74 : // checks
75 0 : disc->check_setup(bNonRegularGrid);
76 :
77 : // cache time dependency
78 0 : m_bNeedLocTimeSeries |= disc->local_time_series_needed();
79 :
80 : // cache grid type required
81 0 : if(bNonRegularGrid)
82 0 : 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 0 : m_vElemDisc[process].push_back(disc);
89 : }
90 0 : }
91 :
92 :
93 : ///////////////////////////////////////////////////////////////////////////////
94 : // DataEvaluatorBase Setup
95 : ///////////////////////////////////////////////////////////////////////////////
96 :
97 :
98 :
99 : template <typename TDomain, typename TElemDisc>
100 0 : void DataEvaluatorBase<TDomain, TElemDisc>::
101 : add_data_to_eval_data(std::vector<SmartPtr<ICplUserData<dim> > >& vEvalData,
102 : std::vector<SmartPtr<ICplUserData<dim> > >& vTryingToAdd)
103 : {
104 : // if empty, we're done
105 0 : if(vTryingToAdd.empty()) return;
106 :
107 : // search for element in already scheduled data
108 : typename std::vector<SmartPtr<ICplUserData<dim> > >::iterator it, itEnd;
109 0 : it = find(vEvalData.begin(), vEvalData.end(), vTryingToAdd.back());
110 :
111 : // if found, skip this data
112 0 : if(it != vEvalData.end())
113 : {
114 : vTryingToAdd.pop_back();
115 0 : 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 0 : it = find(vTryingToAdd.begin(), itEnd, *itEnd);
123 :
124 : // if found, return error of circle dependency
125 0 : if(it != itEnd)
126 0 : 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 0 : for(size_t i = 0; i < data->num_needed_data(); ++i)
132 : {
133 : // add each data separately
134 0 : vTryingToAdd.push_back(data->needed_data(i));
135 0 : add_data_to_eval_data(vEvalData, vTryingToAdd);
136 : }
137 :
138 : // add this data to the evaluation list
139 0 : vEvalData.push_back(data);
140 :
141 : // pop last one, since now added to eval list
142 0 : if(!vTryingToAdd.empty())
143 : vTryingToAdd.pop_back();
144 : }
145 :
146 : template <typename TDomain, typename TElemDisc>
147 0 : void DataEvaluatorBase<TDomain, TElemDisc>::extract_imports_and_userdata(int subsetIndex, int discPart)
148 : {
149 0 : 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 0 : for(size_t d = 0; d < m_vElemDisc[PT_ALL].size(); ++d)
175 : {
176 0 : TElemDisc* disc = m_vElemDisc[PT_ALL][d];
177 :
178 : // loop imports
179 0 : for(size_t i = 0; i < disc->num_imports(); ++i)
180 : {
181 : // get import
182 0 : IDataImport<dim>* iimp = &(disc->get_import(i));
183 :
184 : // skip non-given data (no need for evaluation)
185 0 : if(!iimp->data_given()) continue;
186 :
187 : // check part
188 0 : if( !(iimp->part() & discPart) ) continue;
189 :
190 : // check correct process type
191 0 : if(iimp->part() == MASS)
192 0 : if(!disc->is_time_dependent()) continue;
193 :
194 : // push export on stack of needed data
195 0 : vTryingToAdd.push_back(iimp->data());
196 :
197 : // add data and all dependency to evaluation list
198 : try{
199 0 : add_data_to_eval_data(vEvalData, vTryingToAdd);
200 : }
201 0 : 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 0 : if(!vTryingToAdd.empty())
206 0 : UG_THROW("DataEvaluatorBase:"
207 : " Internal Error, UserData queue not empty after adding.");
208 :
209 : // done if and only if zero-derivative
210 0 : 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 0 : m_vImport[PT_ALL][iimp->part()].push_back(iimp);
218 0 : 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 0 : 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 0 : ipData->set_function_pattern(m_spFctPattern);
253 : }
254 0 : UG_CATCH_THROW("DataEvaluatorBase: Cannot set FunctionPattern to UserData.");
255 :
256 : // sort data into const and non-solution dependent
257 0 : if(ipData->constant()) {m_vConstData.push_back(ipData); continue;}
258 0 : if(ipData->zero_derivative()){m_vPosData.push_back(ipData); continue;}
259 :
260 : // save as dependent data
261 0 : m_vDependentData.push_back(ipData);
262 : }
263 :
264 : // Handle time dependency
265 : // NOTE: constant data is not processed.
266 0 : if(m_pLocTimeSeries != NULL){
267 0 : for(size_t i = 0; i < m_vPosData.size(); ++i)
268 0 : m_vPosData[i]->set_times(m_pLocTimeSeries->times());
269 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i)
270 0 : m_vDependentData[i]->set_times(m_pLocTimeSeries->times());
271 : }
272 :
273 : // NOTE: constant data is not processed, since constant == independent of si
274 0 : for(size_t i = 0; i < m_vPosData.size(); ++i)
275 : m_vPosData[i]->set_subset(subsetIndex);
276 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i)
277 : m_vDependentData[i]->set_subset(subsetIndex);
278 0 : }
279 :
280 : template <typename TDomain, typename TElemDisc>
281 0 : void DataEvaluatorBase<TDomain, TElemDisc>::clear_extracted_data_and_mappings()
282 : {
283 0 : 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 0 : }
293 :
294 :
295 :
296 : template <typename TDomain, typename TElemDisc>
297 0 : void DataEvaluatorBase<TDomain, TElemDisc>::clear_positions_in_user_data()
298 : {
299 : // remove ip series for all used UserData
300 0 : for(size_t i = 0; i < m_vConstData.size(); ++i) m_vConstData[i]->clear();
301 0 : for(size_t i = 0; i < m_vPosData.size(); ++i) m_vPosData[i]->clear();
302 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i) m_vDependentData[i]->clear();
303 0 : }
304 :
305 :
306 : ///////////////////////////////////////////////////////////////////////////////
307 : // prepare / finish
308 : ///////////////////////////////////////////////////////////////////////////////
309 :
310 : template <typename TDomain, typename TElemDisc>
311 0 : void DataEvaluatorBase<TDomain, TElemDisc>::set_time_point(const size_t timePoint)
312 : {
313 0 : for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
314 0 : m_vElemDisc[PT_ALL][i]->set_time_point(timePoint);
315 :
316 : // NOTE: constant data is not processed.
317 0 : for(size_t i = 0; i < m_vPosData.size(); ++i)
318 : m_vPosData[i]->set_time_point(timePoint);
319 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i)
320 : m_vDependentData[i]->set_time_point(timePoint);
321 0 : }
322 :
323 :
324 :
325 : ///////////////////////////////////////////////////////////////////////////////
326 : // Error estimator's routines
327 : ///////////////////////////////////////////////////////////////////////////////
328 :
329 : template <typename TDomain, typename TElemDisc>
330 0 : void DataEvaluatorBase<TDomain, TElemDisc>::
331 : prepare_err_est_elem_loop(const ReferenceObjectID id, int si)
332 : {
333 : // prepare loop (elem disc set local ip series here)
334 : try{
335 0 : for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
336 0 : m_vElemDisc[PT_ALL][i]->do_prep_err_est_elem_loop(id, si);
337 : }
338 0 : 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 0 : extract_imports_and_userdata(si, m_discPart);
344 : }
345 0 : UG_CATCH_THROW("DataEvaluatorBase::prepare_err_est_elem_loop: "
346 : "Cannot extract imports and userdata.");
347 :
348 : // check setup of imports
349 : try{
350 0 : for(size_t i = 0; i < m_vImport[PT_ALL][MASS].size(); ++i)
351 0 : m_vImport[PT_ALL][MASS][i]->check_setup();
352 0 : for(size_t i = 0; i < m_vImport[PT_ALL][STIFF].size(); ++i)
353 0 : m_vImport[PT_ALL][STIFF][i]->check_setup();
354 0 : for(size_t i = 0; i < m_vImport[PT_ALL][RHS].size(); ++i)
355 0 : m_vImport[PT_ALL][RHS][i]->check_setup();
356 : }
357 0 : UG_CATCH_THROW("DataEvaluatorBase::prepare_err_est_elem_loop: Import not correctly implemented.");
358 :
359 : // prepare and check dependent data
360 : try{
361 0 : for(size_t i = 0; i < m_vDependentData.size(); ++i){
362 0 : m_vDependentData[i]->check_setup();
363 : }
364 : }
365 0 : 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 0 : for(size_t i = 0; i < m_vConstData.size(); ++i)
370 0 : m_vConstData[i]->compute((LocalVector*)NULL, NULL, NULL, false);
371 0 : }
372 :
373 : template <typename TDomain, typename TElemDisc>
374 0 : void DataEvaluatorBase<TDomain, TElemDisc>::finish_err_est_elem_loop()
375 : {
376 : // finish each elem error estimator disc
377 : try{
378 0 : for(size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
379 0 : m_vElemDisc[PT_ALL][i]->do_fsh_err_est_elem_loop();
380 : }
381 0 : UG_CATCH_THROW("DataEvaluatorBase::finish_err_est_elem_loop: Cannot finish error estimator element loop");
382 :
383 : // clear positions at user data
384 0 : clear_positions_in_user_data();
385 0 : }
386 :
387 : template <typename TDomain, typename TElemDisc>
388 0 : void DataEvaluatorBase<TDomain, TElemDisc>::prepare_err_est_elem
389 : (LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[],
390 : const LocalIndices& ind, bool bDeriv)
391 : {
392 : try
393 : {
394 0 : for (size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
395 0 : m_vElemDisc[PT_ALL][i]->do_prep_err_est_elem(u, elem, vCornerCoords);
396 : }
397 0 : UG_CATCH_THROW("DataEvaluatorBase::compute_elem_err_est: Cannot prepare element.");
398 :
399 : // evaluate position data
400 0 : for (size_t i = 0; i < m_vPosData.size(); ++i)
401 0 : 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 0 : for (size_t i = 0; i < m_vDependentData.size(); ++i)
414 : {
415 0 : u.access_by_map(m_vDependentData[i]->map());
416 0 : m_vDependentData[i]->compute(&u, elem, vCornerCoords, false);
417 : }
418 : }
419 0 : UG_CATCH_THROW("DataEvaluatorBase::compute_elem_err_est: Cannot compute data for Export or Linker.");
420 0 : }
421 :
422 :
423 : template <typename TDomain, typename TElemDisc>
424 0 : void DataEvaluatorBase<TDomain, TElemDisc>::
425 : compute_err_est_A_elem
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 0 : for (std::size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
437 0 : m_vElemDisc[PT_ALL][i]->do_compute_err_est_A_elem(u, elem, vCornerCoords, scaleStiff);
438 : }
439 0 : UG_CATCH_THROW("DataEvaluatorBase::compute_err_est_A_elem: Cannot assemble stiffness part of error estimator");
440 0 : }
441 :
442 : template <typename TDomain, typename TElemDisc>
443 0 : void DataEvaluatorBase<TDomain, TElemDisc>::
444 : compute_err_est_M_elem
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 0 : for (std::size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
456 0 : m_vElemDisc[PT_ALL][i]->do_compute_err_est_M_elem(u, elem, vCornerCoords, scaleMass);
457 : }
458 0 : UG_CATCH_THROW("DataEvaluatorBase::compute_err_est_A_elem: Cannot assemble stiffness part of error estimator");
459 0 : }
460 :
461 : template <typename TDomain, typename TElemDisc>
462 0 : void DataEvaluatorBase<TDomain, TElemDisc>::
463 : compute_err_est_rhs_elem
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 0 : for (std::size_t i = 0; i < m_vElemDisc[PT_ALL].size(); ++i)
474 0 : m_vElemDisc[PT_ALL][i]->do_compute_err_est_rhs_elem(elem, vCornerCoords, scaleStiff);
475 : }
476 0 : UG_CATCH_THROW("DataEvaluatorBase::compute_err_est_rhs_elem: Cannot assemble rhs part of error estimator");
477 0 : }
478 :
479 :
480 : } // namespace ug
481 :
482 :
483 : #endif
|
