Line data Source code
1 : /*
2 : * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
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__USER_DATA__USER_DATA_IMPL__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__USER_DATA__USER_DATA_IMPL__
35 :
36 : #include "user_data.h"
37 : #include "lib_disc/common/groups_util.h"
38 :
39 : namespace ug{
40 :
41 : ////////////////////////////////////////////////////////////////////////////////
42 : // ICplUserData
43 : ////////////////////////////////////////////////////////////////////////////////
44 :
45 : template <int dim>
46 0 : ICplUserData<dim>::ICplUserData()
47 0 : : m_locPosDim(-1), m_timePoint(0), m_defaultTimePoint(-1), m_si(-1)
48 : {
49 : m_vNumIP.clear();
50 : m_vMayChange.clear();
51 : m_locPosDim = -1;
52 : m_pvLocIP1d.clear(); m_pvLocIP2d.clear(); m_pvLocIP3d.clear();
53 0 : m_vTime.clear(); m_vTime.push_back(0.0);
54 0 : }
55 :
56 : template <int dim>
57 0 : void ICplUserData<dim>::clear()
58 : {
59 0 : local_ip_series_to_be_cleared();
60 : m_vNumIP.clear();
61 : m_vMayChange.clear();
62 0 : m_locPosDim = -1;
63 : m_pvLocIP1d.clear(); m_pvLocIP2d.clear(); m_pvLocIP3d.clear();
64 0 : m_timePoint = 0;
65 0 : m_vTime.clear(); m_vTime.push_back(0.0);
66 0 : m_si = -1;
67 0 : }
68 :
69 : template <int dim>
70 : template <int ldim>
71 0 : size_t ICplUserData<dim>::register_local_ip_series(const MathVector<ldim>* vPos,
72 : const size_t numIP,
73 : const int timePointSpec,
74 : bool bMayChange)
75 : {
76 : // check, that dimension is ok.
77 0 : if(m_locPosDim == -1) m_locPosDim = ldim;
78 0 : else if(m_locPosDim != ldim)
79 0 : UG_THROW("Local IP dimension conflict");
80 :
81 : // get the "right" time point specification
82 0 : int theTimePoint = (m_defaultTimePoint >= 0)? m_defaultTimePoint : timePointSpec;
83 :
84 : // get local positions
85 : std::vector<const MathVector<ldim>*>& vvIP = get_local_ips(Int2Type<ldim>());
86 :
87 : // search for ips
88 : // we only identify ip series if the local ip positions will not change
89 0 : if(!bMayChange && numIP != 0)
90 0 : for(size_t s = 0; s < vvIP.size(); ++s)
91 : {
92 : // return series number iff exists and local ips remain constant
93 0 : if(!m_vMayChange[s])
94 0 : if(vvIP[s] == vPos && m_vNumIP[s] == numIP && m_vTimePoint[s] == theTimePoint)
95 0 : return s;
96 : }
97 :
98 : // if series not yet registered, add it
99 0 : vvIP.push_back(vPos);
100 0 : m_vNumIP.push_back(numIP);
101 0 : m_vTimePoint.push_back(theTimePoint);
102 0 : m_vMayChange.push_back(bMayChange);
103 :
104 : // invoke callback:
105 : // This callback is called, whenever the local_ip_series have changed. It
106 : // allows derived classes to react on this changes. For example, the data
107 : // linker must himself request local_ip_series from the data inputs of
108 : // the linker. In addition value fields and derivative fields must be adjusted
109 : // in UserData<TData, dim> etc.
110 0 : local_ip_series_added(m_vNumIP.size() - 1);
111 :
112 : // return new series id
113 0 : return m_vNumIP.size() - 1;
114 : }
115 :
116 :
117 : template <int dim>
118 : template <int ldim>
119 0 : void ICplUserData<dim>::set_local_ips(const size_t seriesID,
120 : const MathVector<ldim>* vPos,
121 : const size_t numIP)
122 : {
123 : // check series id
124 0 : if(seriesID >= num_series())
125 0 : UG_THROW("Trying to set new ips for invalid seriesID "<<seriesID);
126 :
127 : // check that series is changeable
128 0 : if(!m_vMayChange[seriesID])
129 0 : UG_THROW("Local IP is not changable, but trying to set new ips.");
130 :
131 : // check, that dimension is ok.
132 0 : if(m_locPosDim == -1) m_locPosDim = ldim;
133 0 : else if(m_locPosDim != ldim)
134 0 : UG_THROW("Local IP dimension conflict");
135 :
136 : // get local positions
137 : std::vector<const MathVector<ldim>*>& vvIP = get_local_ips(Int2Type<ldim>());
138 :
139 : // check if still at same position and with same numIP. In that case the
140 : // positions have not changed. We have nothing to do
141 0 : if(vvIP[seriesID] == vPos && m_vNumIP[seriesID] == numIP) return;
142 :
143 : // remember new positions and numIP
144 0 : vvIP[seriesID] = vPos;
145 0 : m_vNumIP[seriesID] = numIP;
146 :
147 : // invoke callback:
148 : // This callback is called, whenever the local_ip_series have changed. It
149 : // allows derived classes to react on this changes. For example, the data
150 : // linker must himself request local_ip_series from the data inputs of
151 : // the linker. In addition value fields and derivative fields must be adjusted
152 : // in UserData<TData, dim> etc.
153 0 : local_ips_changed(seriesID, numIP);
154 : }
155 :
156 : template <int dim>
157 0 : void ICplUserData<dim>::set_time_point(const size_t seriesID,
158 : const int timePointSpec)
159 : {
160 : // check series id
161 0 : if(seriesID >= num_series())
162 0 : UG_THROW("Trying to set new ips for invalid seriesID "<<seriesID);
163 :
164 : // check that series is changeable
165 0 : if(!m_vMayChange[seriesID])
166 0 : UG_THROW("Time point specification is not changable, but trying to set a new one.");
167 :
168 : // set the new time point specification (if it is not prescribed by the object)
169 0 : m_vTimePoint[seriesID] = (m_defaultTimePoint >= 0)? m_defaultTimePoint : timePointSpec;
170 :
171 : //TODO: Should we call the callback here? (No data sizes are changed!)
172 0 : }
173 :
174 : template <int dim>
175 : template <int ldim>
176 0 : const MathVector<ldim>* ICplUserData<dim>::local_ips(size_t s) const
177 : {
178 : // check, that dimension is ok.
179 0 : if(m_locPosDim != ldim) UG_THROW("Local IP dimension conflict");
180 :
181 : UG_ASSERT(s < num_series(), "Wrong series id");
182 :
183 : // NOTE: local ips may be NULL, if no ip position given, i.e. num_ip(s) == 0
184 0 : return get_local_ips(Int2Type<ldim>())[s];
185 : }
186 :
187 : template <int dim>
188 : template <int ldim>
189 : const MathVector<ldim>& ICplUserData<dim>::local_ip(size_t s, size_t ip) const
190 : {
191 : // check, that dimension is ok.
192 : if(m_locPosDim != ldim) UG_THROW("Local IP dimension conflict");
193 :
194 : UG_ASSERT(s < num_series(), "Wrong series id");
195 : UG_ASSERT(ip < num_ip(s), "Invalid index.");
196 :
197 : return get_local_ips(Int2Type<ldim>())[s][ip];
198 : }
199 :
200 : template <int dim>
201 : int ICplUserData<dim>::time_point_specification(size_t s) const
202 : {
203 : UG_ASSERT(s < num_series(), "Wrong series id");
204 :
205 : return m_vTimePoint[s];
206 : }
207 :
208 : template <int dim>
209 : size_t ICplUserData<dim>::time_point(size_t s) const
210 : {
211 : UG_ASSERT(s < num_series(), "Wrong series id:" << s << ">=" << num_series());
212 :
213 : // size_t time_spec;
214 : // if ((time_spec = m_vTimePoint[s]) >= 0)
215 0 : if (m_vTimePoint[s] >= 0)
216 0 : return m_vTimePoint[s];
217 0 : return m_timePoint;
218 : }
219 :
220 : template <int dim>
221 : bool ICplUserData<dim>::at_current_time(size_t s) const
222 : {
223 : UG_ASSERT(s < num_series(), "Wrong series id:" << s << ">=" << num_series());
224 :
225 : int time_spec;
226 0 : if ((time_spec = m_vTimePoint[s]) >= 0)
227 0 : return ((size_t) time_spec) == m_timePoint;
228 : return true;
229 : }
230 :
231 : template <int dim>
232 0 : void ICplUserData<dim>::set_global_ips(size_t s, const MathVector<dim>* vPos, size_t numIP)
233 : {
234 : UG_ASSERT(s < num_series(), "Wrong series id: "<<s<<" (numSeries: "<<num_series()<<")");
235 :
236 : // check number of ips (must match local ip number)
237 0 : if(numIP != num_ip(s))
238 0 : UG_THROW("UserData::set_global_ips: Num Local IPs is " << num_ip(s)
239 : << ", but trying to set Num Global IPs: " << numIP <<
240 : " for series "<< s);
241 :
242 : // remember global positions
243 0 : m_vvGlobPos[s] = vPos;
244 :
245 : // invoke callback:
246 : // this callback is called every time the global position changes. It gives
247 : // derived classes the possibility to react on this fact. E.g. the data
248 : // linker must forward the global positions to its own imports.
249 0 : global_ips_changed(s, vPos, numIP);
250 0 : }
251 :
252 : template <int dim>
253 : inline void ICplUserData<dim>::check_s(size_t s) const
254 : {
255 : UG_ASSERT(s < num_series(), "Wrong series id");
256 : UG_ASSERT(s < m_vvGlobPos.size(), "Invalid index.");
257 : }
258 :
259 : template <int dim>
260 : inline void ICplUserData<dim>::check_s_ip(size_t s, size_t ip) const
261 : {
262 : check_s(s);
263 : UG_ASSERT(ip < num_ip(s), "Invalid index.");
264 : UG_ASSERT(m_vvGlobPos[s] != NULL, "Global IP not set.");
265 : }
266 :
267 : ////////////////////////////////////////////////////////////////////////////////
268 : // UserData
269 : ////////////////////////////////////////////////////////////////////////////////
270 :
271 : template <typename TData, int dim, typename TRet>
272 0 : void CplUserData<TData,dim,TRet>::
273 : register_storage_callback(DataImport<TData,dim>* obj, void (DataImport<TData,dim>::*func)())
274 : {
275 : typedef std::pair<DataImport<TData,dim>*, CallbackFct> Pair;
276 : // m_vCallback.push_back(Pair(obj,func));
277 0 : m_vCallback.push_back(Pair(obj, boost::bind(func, obj)));
278 0 : }
279 :
280 : template <typename TData, int dim, typename TRet>
281 0 : void CplUserData<TData,dim,TRet>::
282 : unregister_storage_callback(DataImport<TData,dim>* obj)
283 : {
284 : typedef typename std::vector<std::pair<DataImport<TData,dim>*, CallbackFct> > VecType;
285 : typedef typename VecType::iterator iterator;
286 : iterator iter = m_vCallback.begin();
287 0 : while(iter != m_vCallback.end())
288 : {
289 0 : if((*iter).first == obj) iter = m_vCallback.erase(iter);
290 : else ++iter;
291 : }
292 0 : }
293 :
294 : template <typename TData, int dim, typename TRet>
295 : void CplUserData<TData,dim,TRet>::
296 : call_storage_callback() const
297 : {
298 : typedef typename std::vector<std::pair<DataImport<TData,dim>*, CallbackFct> > VecType;
299 : typedef typename VecType::const_iterator iterator;
300 0 : for(iterator iter = m_vCallback.begin(); iter != m_vCallback.end(); ++iter)
301 : {
302 : // (((*iter).first)->*((*iter).second))();
303 0 : ((*iter).second)();
304 : }
305 : }
306 :
307 : template <typename TData, int dim, typename TRet>
308 : inline void CplUserData<TData,dim,TRet>::check_series(size_t s) const
309 : {
310 : UG_ASSERT(s < num_series(), "Wrong series id"<<s);
311 : UG_ASSERT(s < m_vvValue.size(), "Invalid index "<<s);
312 : }
313 :
314 : template <typename TData, int dim, typename TRet>
315 : inline void CplUserData<TData,dim,TRet>::check_series_ip(size_t s, size_t ip) const
316 : {
317 : check_series(s);
318 : UG_ASSERT(ip < num_ip(s), "Invalid index "<<ip);
319 : UG_ASSERT(ip < m_vvValue[s].size(), "Invalid index "<<ip);
320 : }
321 :
322 : template <typename TData, int dim, typename TRet>
323 0 : void CplUserData<TData,dim,TRet>::local_ip_series_added(const size_t seriesID)
324 : {
325 : const size_t s = seriesID;
326 :
327 : // check, that only increasing the data, this is important to guarantee,
328 : // that the allocated memory pointer remain valid. They are used outside of
329 : // the class as well to allow fast access to the data.
330 0 : if(s < m_vvValue.size())
331 0 : UG_THROW("Decrease is not implemented. Series: "<<s<<
332 : ", currNumSeries: "<<m_vvValue.size());
333 :
334 : // increase number of series if needed
335 0 : m_vvValue.resize(s+1);
336 0 : m_vvBoolFlag.resize(s+1);
337 :
338 : // allocate new storage
339 0 : m_vvValue[s].resize(num_ip(s));
340 0 : m_vvBoolFlag[s].resize(num_ip(s), true);
341 0 : value_storage_changed(s);
342 : call_storage_callback();
343 :
344 : // call base class callback
345 : base_type::local_ip_series_added(seriesID);
346 0 : }
347 :
348 : template <typename TData, int dim, typename TRet>
349 0 : void CplUserData<TData,dim,TRet>::local_ip_series_to_be_cleared()
350 : {
351 : // free the memory
352 : // clear all series
353 : m_vvValue.clear();
354 : m_vvBoolFlag.clear();
355 :
356 : // call base class callback (if implementation given)
357 : // base_type::local_ip_series_to_be_cleared();
358 0 : }
359 :
360 : template <typename TData, int dim, typename TRet>
361 0 : void CplUserData<TData,dim,TRet>::local_ips_changed(const size_t seriesID, const size_t newNumIP)
362 : {
363 : // resize only when more data is needed than actually allocated
364 0 : if(newNumIP >= m_vvValue[seriesID].size())
365 : {
366 : // resize
367 0 : m_vvValue[seriesID].resize(newNumIP);
368 0 : m_vvBoolFlag[seriesID].resize(newNumIP, true);
369 :
370 : // invoke callback
371 0 : value_storage_changed(seriesID);
372 : call_storage_callback();
373 : }
374 :
375 : // call base class callback (if implementation given)
376 : // base_type::local_ips_changed(seriesID);
377 0 : }
378 :
379 : ////////////////////////////////////////////////////////////////////////////////
380 : // DependentUserData
381 : ////////////////////////////////////////////////////////////////////////////////
382 :
383 : template <typename TData, int dim>
384 0 : void DependentUserData<TData,dim>::set_function_pattern(ConstSmartPtr<FunctionPattern> fctPatt)
385 : {
386 0 : this->m_fctGrp.set_function_pattern(fctPatt);
387 0 : extract_fct_grp();
388 0 : }
389 :
390 : template <typename TData, int dim>
391 0 : void DependentUserData<TData,dim>::set_functions(const char* symbFct)
392 : {
393 0 : set_functions(std::string(symbFct));
394 0 : }
395 :
396 : template <typename TData, int dim>
397 0 : void DependentUserData<TData,dim>::set_functions(const std::string& symbFct)
398 : {
399 0 : set_functions(TokenizeTrimString(symbFct));
400 0 : }
401 :
402 : template <typename TData, int dim>
403 : void DependentUserData<TData,dim>::set_functions(const std::vector<std::string>& symbFct)
404 : {
405 0 : m_SymbFct = symbFct;
406 0 : extract_fct_grp();
407 0 : }
408 :
409 : template <typename TData, int dim>
410 0 : void DependentUserData<TData,dim>::extract_fct_grp()
411 : {
412 : // if associated infos missing return
413 0 : ConstSmartPtr<FunctionPattern> spFctPatt = this->m_fctGrp.function_pattern();
414 0 : if(spFctPatt.invalid()) return;
415 :
416 : // if no function passed, clear functions
417 0 : if(m_SymbFct.size() == 1 && m_SymbFct[0].empty()) m_SymbFct.clear();
418 :
419 : // if functions passed with separator, but not all tokens filled, throw error
420 0 : for(size_t i = 0; i < m_SymbFct.size(); ++i)
421 : {
422 0 : if(m_SymbFct.empty())
423 0 : UG_THROW("Error while setting functions in a DependentUserData: passed "
424 : "function string lacks a "
425 : "function specification at position "<<i<<"(of "
426 : <<m_SymbFct.size()-1<<")");
427 : }
428 :
429 0 : if(m_SymbFct.empty()){
430 0 : this->m_fctGrp.clear();
431 0 : return;
432 : }
433 :
434 : // create function group of this elem disc
435 : try{
436 0 : this->m_fctGrp.clear();
437 0 : this->m_fctGrp.add(m_SymbFct);
438 0 : }UG_CATCH_THROW("DependentUserData: Cannot find some symbolic function "
439 : "name.");
440 :
441 : // create a mapping between all functions and the function group of this
442 : // element disc.
443 : try{
444 0 : CreateFunctionIndexMapping(this->m_map, this->m_fctGrp, spFctPatt);
445 0 : }UG_CATCH_THROW("DependentUserData: Cannot create Function Index Mapping.");
446 :
447 0 : this->check_setup();
448 : }
449 :
450 : template <typename TData, int dim>
451 0 : void DependentUserData<TData,dim>::update_dof_sizes(const LocalIndices& ind)
452 : {
453 : // check size
454 0 : const FunctionIndexMapping& map = this->map();
455 : UG_ASSERT(map.num_fct() == this->num_fct(), "Number function mismatch.");
456 :
457 : // cache numFct and their numDoFs
458 0 : m_vvNumDoFPerFct.resize(map.num_fct());
459 0 : for(size_t fct = 0; fct < m_vvNumDoFPerFct.size(); ++fct)
460 0 : m_vvNumDoFPerFct[fct] = ind.num_dof(map[fct]);
461 :
462 : resize_deriv_array();
463 0 : }
464 :
465 : template <typename TData, int dim>
466 : void DependentUserData<TData,dim>::resize_deriv_array()
467 : {
468 : // resize num fct
469 0 : for(size_t s = 0; s < m_vvvvDeriv.size(); ++s)
470 0 : resize_deriv_array(s);
471 : }
472 :
473 : template <typename TData, int dim>
474 0 : void DependentUserData<TData,dim>::resize_deriv_array(const size_t s)
475 : {
476 : // resize ips
477 0 : m_vvvvDeriv[s].resize(num_ip(s));
478 :
479 0 : for(size_t ip = 0; ip < m_vvvvDeriv[s].size(); ++ip)
480 : {
481 : // resize num fct
482 0 : m_vvvvDeriv[s][ip].resize(m_vvNumDoFPerFct.size());
483 :
484 : // resize dofs
485 0 : for(size_t fct = 0; fct < m_vvNumDoFPerFct.size(); ++fct)
486 0 : m_vvvvDeriv[s][ip][fct].resize(m_vvNumDoFPerFct[fct]);
487 : }
488 0 : }
489 :
490 : template <typename TData, int dim>
491 0 : void DependentUserData<TData,dim>::set_zero(std::vector<std::vector<TData> > vvvDeriv[], const size_t nip)
492 : {
493 0 : for(size_t ip = 0; ip < nip; ++ip)
494 0 : for(size_t fct = 0; fct < vvvDeriv[ip].size(); ++fct)
495 0 : for(size_t sh = 0; sh < vvvDeriv[ip][fct].size(); ++sh)
496 : {
497 0 : vvvDeriv[ip][fct][sh] = 0.0;
498 : }
499 0 : }
500 :
501 : template <typename TData, int dim>
502 : inline void DependentUserData<TData,dim>::check_s_ip(size_t s, size_t ip) const
503 : {
504 : UG_ASSERT(s < this->num_series(), "Wrong series id"<<s);
505 : UG_ASSERT(s < m_vvvvDeriv.size(), "Invalid index "<<s);
506 : UG_ASSERT(ip < this->num_ip(s), "Invalid index "<<ip);
507 : UG_ASSERT(ip < m_vvvvDeriv[s].size(), "Invalid index "<<ip);
508 : }
509 :
510 : template <typename TData, int dim>
511 : inline void DependentUserData<TData,dim>::check_s_ip_fct(size_t s, size_t ip, size_t fct) const
512 : {
513 : check_s_ip(s,ip);
514 : UG_ASSERT(fct < m_vvvvDeriv[s][ip].size(), "Invalid index.");
515 : }
516 :
517 : template <typename TData, int dim>
518 : inline void DependentUserData<TData,dim>::check_s_ip_fct_dof(size_t s, size_t ip, size_t fct, size_t dof) const
519 : {
520 : check_s_ip_fct(s,ip,fct);
521 : UG_ASSERT(dof < m_vvvvDeriv[s][ip][fct].size(), "Invalid index.");
522 : }
523 :
524 : template <typename TData, int dim>
525 0 : void DependentUserData<TData,dim>::local_ip_series_added(const size_t seriesID)
526 : {
527 : // adjust data arrays
528 0 : m_vvvvDeriv.resize(seriesID+1);
529 :
530 : // forward change signal to base class
531 0 : base_type::local_ip_series_added(seriesID);
532 0 : }
533 :
534 : template <typename TData, int dim>
535 0 : void DependentUserData<TData,dim>::local_ip_series_to_be_cleared()
536 : {
537 : // adjust data arrays
538 : m_vvvvDeriv.clear();
539 :
540 : // forward change signal to base class
541 0 : base_type::local_ip_series_to_be_cleared();
542 0 : }
543 :
544 : template <typename TData, int dim>
545 0 : void DependentUserData<TData,dim>::local_ips_changed(const size_t seriesID, const size_t newNumIP)
546 : {
547 : UG_ASSERT(seriesID < m_vvvvDeriv.size(), "wrong series id.");
548 :
549 : // resize only when more data is needed than actually allocated
550 0 : if(newNumIP >= m_vvvvDeriv[seriesID].size())
551 0 : resize_deriv_array(seriesID);
552 :
553 : // call base class callback (if implementation given)
554 0 : base_type::local_ips_changed(seriesID, newNumIP);
555 0 : }
556 :
557 : } // end namespace ug
558 :
559 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__USER_DATA__USER_DATA_IMPL__ */
|
