Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Raphael Prohl
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_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__OBSTACLE_CONSTRAINT_INTERFACE_IMPL__
34 : #define __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__OBSTACLE_CONSTRAINT_INTERFACE_IMPL__
35 :
36 : #include "obstacle_constraint_interface.h"
37 : #include "lib_disc/function_spaces/dof_position_util.h"
38 :
39 : #ifdef UG_FOR_LUA
40 : #include "bindings/lua/lua_user_data.h"
41 : #endif
42 :
43 : namespace ug{
44 :
45 : template <typename TDomain, typename TAlgebra>
46 0 : void IObstacleConstraint<TDomain, TAlgebra>::
47 : init()
48 : {
49 0 : if(m_spDomain.invalid())
50 0 : UG_THROW("No domain set in 'IObstacleConstraint::init' \n");
51 :
52 0 : if(m_spDD.invalid())
53 0 : UG_THROW("DofDistribution not set in 'IObstacleConstraint'.");
54 :
55 : // build up a map of obstacle dofs and its corresponding obstacle values
56 0 : init_obstacle_dofs_with_values(1.0);
57 :
58 : UG_LOG("In IObstacleConstraint::init: "<<m_mObstacleValues.size()<< " obstacleDoFs tagged \n");
59 : UG_LOG("\n");
60 0 : }
61 :
62 : template <typename TDomain, typename TAlgebra>
63 0 : void IObstacleConstraint<TDomain, TAlgebra>::
64 : add(SmartPtr<UserData<number, dim, bool> > func, const char* function)
65 : {
66 0 : m_vCondNumberData.push_back(CondNumberData(func, function));
67 0 : }
68 :
69 : template <typename TDomain, typename TAlgebra>
70 0 : void IObstacleConstraint<TDomain, TAlgebra>::
71 : add(SmartPtr<UserData<number, dim, bool> > func, const char* function, const char* subsets)
72 : {
73 0 : m_vCondNumberData.push_back(CondNumberData(func, function, subsets));
74 0 : }
75 :
76 : template <typename TDomain, typename TAlgebra>
77 0 : void IObstacleConstraint<TDomain, TAlgebra>::
78 : add(SmartPtr<UserData<number, dim> > func, const char* function)
79 : {
80 0 : m_vNumberData.push_back(NumberData(func, function));
81 0 : }
82 :
83 : template <typename TDomain, typename TAlgebra>
84 0 : void IObstacleConstraint<TDomain, TAlgebra>::
85 : add(SmartPtr<UserData<number, dim> > func, const char* function, const char* subsets)
86 : {
87 0 : m_vNumberData.push_back(NumberData(func, function, subsets));
88 0 : }
89 :
90 : template <typename TDomain, typename TAlgebra>
91 0 : void IObstacleConstraint<TDomain, TAlgebra>::
92 : add(number value, const char* function)
93 : {
94 0 : m_vConstNumberData.push_back(ConstNumberData(value, function));
95 0 : }
96 :
97 : template <typename TDomain, typename TAlgebra>
98 0 : void IObstacleConstraint<TDomain, TAlgebra>::
99 : add(number value, const char* function, const char* subsets)
100 : {
101 0 : m_vConstNumberData.push_back(ConstNumberData(value, function, subsets));
102 0 : }
103 :
104 : template <typename TDomain, typename TAlgebra>
105 0 : void IObstacleConstraint<TDomain, TAlgebra>::
106 : add(SmartPtr<UserData<MathVector<dim>, dim> > func, const char* functions)
107 : {
108 0 : m_vVectorData.push_back(VectorData(func, functions));
109 0 : }
110 :
111 : template <typename TDomain, typename TAlgebra>
112 0 : void IObstacleConstraint<TDomain, TAlgebra>::
113 : add(SmartPtr<UserData<MathVector<dim>, dim> > func, const char* functions, const char* subsets)
114 : {
115 0 : m_vVectorData.push_back(VectorData(func, functions, subsets));
116 0 : }
117 :
118 : #ifdef UG_FOR_LUA
119 : template <typename TDomain, typename TAlgebra>
120 0 : void IObstacleConstraint<TDomain, TAlgebra>::
121 : add(const char* name, const char* function)
122 : {
123 0 : if(LuaUserData<number, dim>::check_callback_returns(name)){
124 0 : SmartPtr<UserData<number, dim> > sp =
125 : LuaUserDataFactory<number, dim>::create(name);
126 0 : add(sp, function);
127 : return;
128 : }
129 0 : if(LuaUserData<number, dim, bool>::check_callback_returns(name)){
130 0 : SmartPtr<UserData<number, dim, bool> > sp =
131 : LuaUserDataFactory<number, dim, bool>::create(name);
132 0 : add(sp, function);
133 : return;
134 : }
135 0 : if(LuaUserData<MathVector<dim>, dim>::check_callback_returns(name)){
136 0 : SmartPtr<UserData<MathVector<dim>, dim> > sp =
137 : LuaUserDataFactory<MathVector<dim>, dim>::create(name);
138 0 : add(sp, function);
139 : return;
140 : }
141 :
142 : // no match found
143 0 : if(!CheckLuaCallbackName(name))
144 0 : UG_THROW("IObstacleConstraint::add: Lua-Callback with name '"<<name<<
145 : "' does not exist.");
146 :
147 : // name exists but wrong signature
148 0 : UG_THROW("IObstacleConstraint::add: Cannot find matching callback "
149 : "signature. Use one of:\n"
150 : "a) Number - Callback\n"
151 : << (LuaUserData<number, dim>::signature()) << "\n" <<
152 : "b) Conditional Number - Callback\n"
153 : << (LuaUserData<number, dim, bool>::signature()) << "\n" <<
154 : "c) "<<dim<<"d Vector - Callback\n"
155 : << (LuaUserData<MathVector<dim>, dim>::signature()));
156 : }
157 :
158 : template <typename TDomain, typename TAlgebra>
159 0 : void IObstacleConstraint<TDomain, TAlgebra>::
160 : add(const char* name, const char* function, const char* subsets)
161 : {
162 0 : if(LuaUserData<number, dim>::check_callback_returns(name)){
163 0 : SmartPtr<UserData<number, dim> > sp =
164 : LuaUserDataFactory<number, dim>::create(name);
165 0 : add(sp, function, subsets);
166 : return;
167 : }
168 0 : if(LuaUserData<number, dim, bool>::check_callback_returns(name)){
169 0 : SmartPtr<UserData<number, dim, bool> > sp =
170 : LuaUserDataFactory<number, dim, bool>::create(name);
171 0 : add(sp, function, subsets);
172 : return;
173 : }
174 0 : if(LuaUserData<MathVector<dim>, dim>::check_callback_returns(name)){
175 0 : SmartPtr<UserData<MathVector<dim>, dim> > sp =
176 : LuaUserDataFactory<MathVector<dim>, dim>::create(name);
177 0 : add(sp, function, subsets);
178 : return;
179 : }
180 :
181 : // no match found
182 0 : if(!CheckLuaCallbackName(name))
183 0 : UG_THROW("IObstacleConstraint::add: Lua-Callback with name '"<<name<<
184 : "' does not exist.");
185 :
186 : // name exists but wrong signature
187 0 : UG_THROW("IObstacleConstraint::add: Cannot find matching callback "
188 : "signature. Use one of:\n"
189 : "a) Number - Callback\n"
190 : << (LuaUserData<number, dim>::signature()) << "\n" <<
191 : "b) Conditional Number - Callback\n"
192 : << (LuaUserData<number, dim, bool>::signature()) << "\n" <<
193 : "c) "<<dim<<"d Vector - Callback\n"
194 : << (LuaUserData<MathVector<dim>, dim>::signature()));
195 : }
196 : #endif
197 :
198 :
199 : template <typename TDomain, typename TAlgebra>
200 0 : void IObstacleConstraint<TDomain, TAlgebra>::
201 : clear()
202 : {
203 : m_vCondNumberData.clear();
204 : m_vNumberData.clear();
205 : m_vConstNumberData.clear();
206 : m_vVectorData.clear();
207 :
208 : m_mObstacleValues.clear();
209 : m_vActiveDofs.clear();
210 0 : }
211 :
212 : template <typename TDomain, typename TAlgebra>
213 0 : void IObstacleConstraint<TDomain, TAlgebra>::
214 : check_functions_and_subsets(FunctionGroup& functionGroup, SubsetGroup& subsetGroup,
215 : size_t numFct) const
216 : {
217 : // only number of functions allowed
218 0 : if(functionGroup.size() != numFct)
219 0 : UG_THROW("IObstacleConstraint:extract_data:"
220 : " Only "<<numFct<<" function(s) allowed in specification of a"
221 : " Obstacle Value, but the following functions given:"
222 : <<functionGroup);
223 :
224 : // get subsethandler
225 : ConstSmartPtr<ISubsetHandler> pSH = m_spDD->subset_handler();
226 :
227 : // loop subsets
228 0 : for(size_t si = 0; si < subsetGroup.size(); ++si)
229 : {
230 : // get subset index
231 : const int subsetIndex = subsetGroup[si];
232 :
233 : // check that subsetIndex is valid
234 0 : if(subsetIndex < 0 || subsetIndex >= pSH->num_subsets())
235 0 : UG_THROW("IObstacleConstraint:extract_data:"
236 : " Invalid Subset Index " << subsetIndex << ". (Valid is"
237 : " 0, .. , " << pSH->num_subsets() <<").");
238 :
239 : // check all functions
240 0 : for(size_t i=0; i < functionGroup.size(); ++i)
241 : {
242 : const size_t fct = functionGroup[i];
243 :
244 : // check if function exist
245 0 : if(fct >= m_spDD->num_fct())
246 0 : UG_THROW("IObstacleConstraint:extract_data:"
247 : " Function "<< fct << " does not exist in pattern.");
248 :
249 : // check that function is defined for segment
250 0 : if(!m_spDD->is_def_in_subset(fct, subsetIndex))
251 0 : UG_THROW("IObstacleConstraint:extract_data:"
252 : " Function "<<fct<<" not defined on subset "<<subsetIndex);
253 : }
254 : }
255 :
256 : // everything ok
257 0 : }
258 :
259 : template <typename TDomain, typename TAlgebra>
260 : template <typename TUserData, typename TScheduledUserData>
261 0 : void IObstacleConstraint<TDomain, TAlgebra>::
262 : extract_data(std::map<int, std::vector<TUserData*> >& mvUserDataObsSegment,
263 : std::vector<TScheduledUserData>& vUserData)
264 : {
265 : // clear the extracted data
266 : mvUserDataObsSegment.clear();
267 :
268 0 : for(size_t i = 0; i < vUserData.size(); ++i)
269 : {
270 : // create Function Group and Subset Group
271 0 : if (vUserData[i].bWholeDomain == false){
272 : try{
273 0 : vUserData[i].ssGrp = m_spDD->subset_grp_by_name(vUserData[i].ssName.c_str());
274 0 : }UG_CATCH_THROW(" Subsets '"<<vUserData[i].ssName<<"' not"
275 : " all contained in DoFDistribution.");
276 : }
277 : else{
278 0 : SubsetGroup ssGrp = SubsetGroup(m_spDD->subset_handler());
279 0 : ssGrp.add_all();
280 : vUserData[i].ssGrp = ssGrp;
281 : }
282 :
283 0 : FunctionGroup fctGrp;
284 : try{
285 0 : fctGrp = m_spDD->fct_grp_by_name(vUserData[i].fctName.c_str());
286 0 : }UG_CATCH_THROW(" Functions '"<<vUserData[i].fctName<<"' not"
287 : " all contained in DoFDistribution.");
288 :
289 : // check functions and subsets
290 0 : check_functions_and_subsets(fctGrp, vUserData[i].ssGrp, TUserData::numFct);
291 :
292 : // set functions
293 0 : if(fctGrp.size() != TUserData::numFct)
294 0 : UG_THROW("IObstacleConstraint: wrong number of fct");
295 :
296 0 : for(size_t fct = 0; fct < TUserData::numFct; ++fct)
297 : {
298 0 : vUserData[i].fct[fct] = fctGrp[fct];
299 : }
300 :
301 : // loop subsets
302 0 : for(size_t si = 0; si < vUserData[i].ssGrp.size(); ++si)
303 : {
304 : // get subset index and function
305 0 : const int subsetIndex = vUserData[i].ssGrp[si];
306 :
307 : // remember functor and function
308 0 : mvUserDataObsSegment[subsetIndex].push_back(&vUserData[i]);
309 : }
310 : }
311 0 : }
312 :
313 : template <typename TDomain, typename TAlgebra>
314 0 : void IObstacleConstraint<TDomain, TAlgebra>::
315 : extract_data()
316 : {
317 0 : extract_data(m_mNumberObsSegment, m_vNumberData);
318 0 : extract_data(m_mCondNumberObsSegment, m_vCondNumberData);
319 0 : extract_data(m_mConstNumberObsSegment, m_vConstNumberData);
320 0 : extract_data(m_mVectorObsSegment, m_vVectorData);
321 0 : }
322 :
323 : template <typename TDomain, typename TAlgebra>
324 0 : void IObstacleConstraint<TDomain, TAlgebra>::
325 : init_obstacle_dofs_with_values(number time)
326 : {
327 0 : extract_data();
328 :
329 : // reset map of obstacle values and vector of obstacle subset-indices
330 : m_mObstacleValues.clear();
331 0 : m_vObsSubsets.resize(0);
332 :
333 0 : init_obstacle_dofs_with_values<CondNumberData>(m_mCondNumberObsSegment, time);
334 0 : init_obstacle_dofs_with_values<NumberData>(m_mNumberObsSegment, time);
335 0 : init_obstacle_dofs_with_values<ConstNumberData>(m_mConstNumberObsSegment, time);
336 0 : init_obstacle_dofs_with_values<VectorData>(m_mVectorObsSegment, time);
337 0 : }
338 :
339 : template <typename TDomain, typename TAlgebra>
340 : template <typename TUserData>
341 0 : void IObstacleConstraint<TDomain, TAlgebra>::
342 : init_obstacle_dofs_with_values(const std::map<int, std::vector<TUserData*> >& mvUserData, number time)
343 : {
344 : typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
345 0 : for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
346 : {
347 : // get subset index
348 0 : const int si = (*iter).first;
349 :
350 : // store obstacle subsets
351 0 : m_vObsSubsets.push_back(si);
352 :
353 : // get vector of scheduled obstacle data on this subset
354 0 : const std::vector<TUserData*>& vUserData = (*iter).second;
355 :
356 : // gets obstacle values in each base element type
357 : try
358 : {
359 0 : if(m_spDD->max_dofs(VERTEX))
360 0 : init_obstacle_dofs_with_values<RegularVertex, TUserData>(vUserData, si, time);
361 0 : if(m_spDD->max_dofs(EDGE))
362 0 : init_obstacle_dofs_with_values<Edge, TUserData>(vUserData, si, time);
363 0 : if(m_spDD->max_dofs(FACE))
364 0 : init_obstacle_dofs_with_values<Face, TUserData>(vUserData, si, time);
365 0 : if(m_spDD->max_dofs(VOLUME))
366 0 : init_obstacle_dofs_with_values<Volume, TUserData>(vUserData, si, time);
367 : }
368 0 : UG_CATCH_THROW("IObstacleConstraint::init_obstacle_dofs_with_values:"
369 : " While calling 'obstacle_value' for TUserData, aborting.");
370 : }
371 0 : }
372 :
373 : template <typename TDomain, typename TAlgebra>
374 : template <typename TBaseElem, typename TUserData>
375 0 : void IObstacleConstraint<TDomain, TAlgebra>::
376 : init_obstacle_dofs_with_values(const std::vector<TUserData*>& vUserData, int si, number time)
377 : {
378 : // create Multiindex
379 : std::vector<DoFIndex> multInd;
380 :
381 : // readin value
382 : typename TUserData::value_type val;
383 :
384 : // position of dofs
385 : std::vector<position_type> vPos;
386 :
387 : // iterators
388 : typedef typename DoFDistribution::traits<TBaseElem>::const_iterator iter_type;
389 0 : iter_type iter = m_spDD->begin<TBaseElem>(si);
390 0 : iter_type iterEnd = m_spDD->end<TBaseElem>(si);
391 :
392 : // loop elements
393 0 : for( ; iter != iterEnd; iter++)
394 : {
395 : // get baseElem
396 : TBaseElem* elem = *iter;
397 :
398 : // loop obstacle functions on this segment
399 0 : for(size_t i = 0; i < vUserData.size(); ++i)
400 : {
401 0 : for(size_t f = 0; f < TUserData::numFct; ++f)
402 : {
403 : // get function index
404 0 : const size_t fct = vUserData[i]->fct[f];
405 :
406 : // get local finite element id
407 : const LFEID& lfeID = m_spDD->local_finite_element_id(fct);
408 :
409 : // get dof position
410 0 : InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
411 :
412 : // get multi indices
413 0 : m_spDD->inner_dof_indices(elem, fct, multInd);
414 :
415 : UG_ASSERT(multInd.size() == vPos.size(),
416 : "Mismatch: numInd="<<multInd.size()<<", numPos="
417 : <<vPos.size()<<" on "<<elem->reference_object_id());
418 :
419 : // loop dofs on element
420 0 : for(size_t j = 0; j < multInd.size(); ++j)
421 : {
422 : // check if function is an obstacle fct and read value
423 0 : if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
424 :
425 : // deposit obstacle values in a map
426 0 : m_mObstacleValues[multInd[j]] = val[f];
427 : }
428 : }
429 : }
430 : }
431 0 : }
432 :
433 : template <typename TDomain, typename TAlgebra>
434 : bool
435 : IObstacleConstraint<TDomain,TAlgebra>::
436 : is_obs_dof(const DoFIndex& dof)
437 : {
438 0 : if (m_mObstacleValues.find(dof) == m_mObstacleValues.end()){return false;}
439 : else {return true;}
440 : }
441 :
442 : template <typename TDomain, typename TAlgebra>
443 : void
444 0 : IObstacleConstraint<TDomain,TAlgebra>::
445 : adjust_restriction(matrix_type& R, ConstSmartPtr<DoFDistribution> ddCoarse,
446 : ConstSmartPtr<DoFDistribution> ddFine, int type, number time)
447 : {
448 : UG_LOG("IObstacleConstraint<TDomain,TAlgebra>::adjust_restrictionR \n");
449 :
450 0 : R.print();
451 :
452 : typedef typename vector<DoFIndex>::iterator iter_type;
453 : iter_type dofIter = m_vActiveDofs.begin();
454 : iter_type dofIterEnd = m_vActiveDofs.end();
455 0 : for( ; dofIter != dofIterEnd; dofIter++)
456 : {
457 0 : UG_LOG("IObstacleConstraint<TDomain,TAlgebra>::"
458 : "adjust_restrictionR::activeDof : " <<*dofIter<< "\n");
459 0 : SetCol(R, (*dofIter)[0], (*dofIter)[1]);
460 : }
461 :
462 0 : if (m_vActiveDofs.size() > 0)
463 : {
464 : UG_LOG("#OfActiveDofs: " <<m_vActiveDofs.size()<< "\n");
465 0 : R.print();
466 : }
467 : UG_LOG("IObstacleConstraint::adjust_restrictionR() \n");
468 0 : };
469 :
470 : } // end namespace ug
471 :
472 : #endif /* __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__OBSTACLE_CONSTRAINT_INTERFACE_IMPL__ */
|
