Line data Source code
1 : /*
2 : * Copyright (c) 2012-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__CONSTRAINTS__LAGRANGE_DIRICHLET_BOUNDARY_IMPL__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__CONSTRAINTS__LAGRANGE_DIRICHLET_BOUNDARY_IMPL__
35 :
36 : #include "lagrange_dirichlet_boundary.h"
37 : #include "lib_disc/function_spaces/grid_function.h"
38 : #include "lib_disc/function_spaces/dof_position_util.h"
39 :
40 : #ifdef UG_FOR_LUA
41 : #include "bindings/lua/lua_user_data.h"
42 : #endif
43 :
44 : namespace ug{
45 :
46 :
47 : ////////////////////////////////////////////////////////////////////////////////
48 : // setup
49 : ////////////////////////////////////////////////////////////////////////////////
50 :
51 : template <typename TDomain, typename TAlgebra>
52 0 : void DirichletBoundary<TDomain, TAlgebra>::
53 : set_approximation_space(SmartPtr<ApproximationSpace<TDomain> > approxSpace)
54 : {
55 0 : base_type::set_approximation_space(approxSpace);
56 0 : m_spApproxSpace = approxSpace;
57 0 : m_spDomain = approxSpace->domain();
58 0 : m_aaPos = m_spDomain->position_accessor();
59 0 : }
60 :
61 : template <typename TDomain, typename TAlgebra>
62 0 : void DirichletBoundary<TDomain, TAlgebra>::
63 : clear()
64 : {
65 : m_vBNDNumberData.clear();
66 : m_vNumberData.clear();
67 : m_vConstNumberData.clear();
68 : m_vVectorData.clear();
69 0 : }
70 :
71 : template <typename TDomain, typename TAlgebra>
72 0 : void DirichletBoundary<TDomain, TAlgebra>::
73 : add(SmartPtr<UserData<number, dim, bool> > func, const char* function, const char* subsets)
74 : {
75 0 : m_vBNDNumberData.push_back(CondNumberData(func, function, subsets));
76 0 : }
77 :
78 : template <typename TDomain, typename TAlgebra>
79 0 : void DirichletBoundary<TDomain, TAlgebra>::
80 : add(SmartPtr<UserData<number, dim, bool> > func, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
81 : {
82 : std::string function;
83 0 : for(size_t i = 0; i < Fcts.size(); ++i){
84 0 : if(i > 0) function.append(",");
85 : function.append(Fcts[i]);
86 : }
87 : std::string subsets;
88 0 : for(size_t i = 0; i < Subsets.size(); ++i){
89 0 : if(i > 0) subsets.append(",");
90 : subsets.append(Subsets[i]);
91 : }
92 :
93 0 : add(func, function.c_str(), subsets.c_str());
94 0 : }
95 :
96 : template <typename TDomain, typename TAlgebra>
97 0 : void DirichletBoundary<TDomain, TAlgebra>::
98 : add(SmartPtr<UserData<number, dim> > func, const char* function, const char* subsets)
99 : {
100 0 : m_vNumberData.push_back(NumberData(func, function, subsets));
101 0 : }
102 :
103 : template <typename TDomain, typename TAlgebra>
104 0 : void DirichletBoundary<TDomain, TAlgebra>::
105 : add(SmartPtr<UserData<number, dim> > func, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
106 : {
107 : std::string function;
108 0 : for(size_t i = 0; i < Fcts.size(); ++i){
109 0 : if(i > 0) function.append(",");
110 : function.append(Fcts[i]);
111 : }
112 : std::string subsets;
113 0 : for(size_t i = 0; i < Subsets.size(); ++i){
114 0 : if(i > 0) subsets.append(",");
115 : subsets.append(Subsets[i]);
116 : }
117 :
118 0 : add(func, function.c_str(), subsets.c_str());
119 0 : }
120 :
121 : template <typename TDomain, typename TAlgebra>
122 0 : void DirichletBoundary<TDomain, TAlgebra>::
123 : add(number value, const char* function, const char* subsets)
124 : {
125 0 : m_vConstNumberData.push_back(ConstNumberData(value, function, subsets));
126 0 : }
127 :
128 : template <typename TDomain, typename TAlgebra>
129 0 : void DirichletBoundary<TDomain, TAlgebra>::
130 : add(number value, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
131 : {
132 : std::string function;
133 0 : for(size_t i = 0; i < Fcts.size(); ++i){
134 0 : if(i > 0) function.append(",");
135 : function.append(Fcts[i]);
136 : }
137 : std::string subsets;
138 0 : for(size_t i = 0; i < Subsets.size(); ++i){
139 0 : if(i > 0) subsets.append(",");
140 : subsets.append(Subsets[i]);
141 : }
142 :
143 0 : add(value, function.c_str(), subsets.c_str());
144 0 : }
145 :
146 : template <typename TDomain, typename TAlgebra>
147 0 : void DirichletBoundary<TDomain, TAlgebra>::
148 : add(SmartPtr<UserData<MathVector<dim>, dim> > func, const char* functions, const char* subsets)
149 : {
150 0 : m_vVectorData.push_back(VectorData(func, functions, subsets));
151 0 : }
152 :
153 : template <typename TDomain, typename TAlgebra>
154 0 : void DirichletBoundary<TDomain, TAlgebra>::
155 : add(SmartPtr<UserData<MathVector<dim>, dim> > func, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
156 : {
157 : std::string function;
158 0 : for(size_t i = 0; i < Fcts.size(); ++i){
159 0 : if(i > 0) function.append(",");
160 : function.append(Fcts[i]);
161 : }
162 : std::string subsets;
163 0 : for(size_t i = 0; i < Subsets.size(); ++i){
164 0 : if(i > 0) subsets.append(",");
165 : subsets.append(Subsets[i]);
166 : }
167 :
168 0 : add(func, function.c_str(), subsets.c_str());
169 0 : }
170 :
171 : #ifdef UG_FOR_LUA
172 : template <typename TDomain, typename TAlgebra>
173 0 : void DirichletBoundary<TDomain, TAlgebra>::
174 : add(const char* name, const char* function, const char* subsets)
175 : {
176 0 : if(LuaUserData<number, dim>::check_callback_returns(name)){
177 0 : SmartPtr<UserData<number, dim> > sp =
178 : LuaUserDataFactory<number, dim>::create(name);
179 0 : add(sp, function, subsets);
180 : return;
181 : }
182 0 : if(LuaUserData<number, dim, bool>::check_callback_returns(name)){
183 0 : SmartPtr<UserData<number, dim, bool> > sp =
184 : LuaUserDataFactory<number, dim, bool>::create(name);
185 0 : add(sp, function, subsets);
186 : return;
187 : }
188 0 : if(LuaUserData<MathVector<dim>, dim>::check_callback_returns(name)){
189 0 : SmartPtr<UserData<MathVector<dim>, dim> > sp =
190 : LuaUserDataFactory<MathVector<dim>, dim>::create(name);
191 0 : add(sp, function, subsets);
192 : return;
193 : }
194 :
195 : // no match found
196 0 : if(!CheckLuaCallbackName(name))
197 0 : UG_THROW("LagrangeDirichlet::add: Lua-Callback with name '"<<name<<
198 : "' does not exist.");
199 :
200 : // name exists but wrong signature
201 0 : UG_THROW("LagrangeDirichlet::add: Cannot find matching callback "
202 : "signature. Use one of:\n"
203 : "a) Number - Callback\n"
204 : << (LuaUserData<number, dim>::signature()) << "\n" <<
205 : "b) Conditional Number - Callback\n"
206 : << (LuaUserData<number, dim, bool>::signature()) << "\n" <<
207 : "c) "<<dim<<"d Vector - Callback\n"
208 : << (LuaUserData<MathVector<dim>, dim>::signature()));
209 : }
210 :
211 : template <typename TDomain, typename TAlgebra>
212 0 : void DirichletBoundary<TDomain, TAlgebra>::
213 : add(const char* name, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
214 : {
215 : std::string function;
216 0 : for(size_t i = 0; i < Fcts.size(); ++i){
217 0 : if(i > 0) function.append(",");
218 : function.append(Fcts[i]);
219 : }
220 : std::string subsets;
221 0 : for(size_t i = 0; i < Subsets.size(); ++i){
222 0 : if(i > 0) subsets.append(",");
223 : subsets.append(Subsets[i]);
224 : }
225 :
226 0 : add(name, function.c_str(), subsets.c_str());
227 0 : }
228 :
229 : template <typename TDomain, typename TAlgebra>
230 0 : void DirichletBoundary<TDomain, TAlgebra>::
231 : add(LuaFunctionHandle fct, const char* function, const char* subsets)
232 : {
233 0 : if(LuaUserData<number, dim>::check_callback_returns(fct)){
234 0 : SmartPtr<UserData<number, dim> > sp =
235 0 : make_sp(new LuaUserData<number, dim>(fct));
236 0 : add(sp, function, subsets);
237 : return;
238 : }
239 0 : if(LuaUserData<number, dim, bool>::check_callback_returns(fct)){
240 0 : SmartPtr<UserData<number, dim, bool> > sp =
241 0 : make_sp(new LuaUserData<number, dim, bool>(fct));
242 0 : add(sp, function, subsets);
243 : return;
244 : }
245 0 : if(LuaUserData<MathVector<dim>, dim>::check_callback_returns(fct)){
246 0 : SmartPtr<UserData<MathVector<dim>, dim> > sp =
247 0 : make_sp(new LuaUserData<MathVector<dim>, dim>(fct));
248 0 : add(sp, function, subsets);
249 : return;
250 : }
251 :
252 : // name exists but wrong signature
253 0 : UG_THROW("LagrangeDirichlet::add: Cannot find matching callback "
254 : "signature. Use one of:\n"
255 : "a) Number - Callback\n"
256 : << (LuaUserData<number, dim>::signature()) << "\n" <<
257 : "b) Conditional Number - Callback\n"
258 : << (LuaUserData<number, dim, bool>::signature()) << "\n" <<
259 : "c) "<<dim<<"d Vector - Callback\n"
260 : << (LuaUserData<MathVector<dim>, dim>::signature()));
261 : }
262 :
263 : template <typename TDomain, typename TAlgebra>
264 0 : void DirichletBoundary<TDomain, TAlgebra>::
265 : add(LuaFunctionHandle fct, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
266 : {
267 : std::string function;
268 0 : for(size_t i = 0; i < Fcts.size(); ++i){
269 0 : if(i > 0) function.append(",");
270 : function.append(Fcts[i]);
271 : }
272 : std::string subsets;
273 0 : for(size_t i = 0; i < Subsets.size(); ++i){
274 0 : if(i > 0) subsets.append(",");
275 : subsets.append(Subsets[i]);
276 : }
277 :
278 0 : add(fct, function.c_str(), subsets.c_str());
279 0 : }
280 : #endif
281 :
282 : template <typename TDomain, typename TAlgebra>
283 0 : void DirichletBoundary<TDomain, TAlgebra>::
284 : add(const char* functions, const char* subsets)
285 : {
286 0 : m_vOldNumberData.push_back(OldNumberData(functions, subsets));
287 0 : }
288 :
289 : template <typename TDomain, typename TAlgebra>
290 0 : void DirichletBoundary<TDomain, TAlgebra>::
291 : add(const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
292 : {
293 : std::string function;
294 0 : for(size_t i = 0; i < Fcts.size(); ++i){
295 0 : if(i > 0) function.append(",");
296 : function.append(Fcts[i]);
297 : }
298 : std::string subsets;
299 0 : for(size_t i = 0; i < Subsets.size(); ++i){
300 0 : if(i > 0) subsets.append(",");
301 : subsets.append(Subsets[i]);
302 : }
303 :
304 0 : add(function.c_str(), subsets.c_str());
305 0 : }
306 :
307 : template <typename TDomain, typename TAlgebra>
308 0 : void DirichletBoundary<TDomain, TAlgebra>::
309 : check_functions_and_subsets(FunctionGroup& functionGroup, SubsetGroup& subsetGroup, size_t numFct) const
310 : {
311 : // only number of functions allowed
312 0 : if(functionGroup.size() != numFct)
313 0 : UG_THROW("DirichletBoundary:extract_data:"
314 : " Only "<<numFct<<" function(s) allowed in specification of a"
315 : " Dirichlet Value, but the following functions given:"
316 : <<functionGroup);
317 :
318 : // get subsethandler
319 0 : ConstSmartPtr<ISubsetHandler> pSH = m_spApproxSpace->subset_handler();
320 :
321 : // loop subsets
322 0 : for(size_t si = 0; si < subsetGroup.size(); ++si)
323 : {
324 : // get subset index
325 : const int subsetIndex = subsetGroup[si];
326 :
327 : // check that subsetIndex is valid
328 0 : if(subsetIndex < 0 || subsetIndex >= pSH->num_subsets())
329 0 : UG_THROW("DirichletBoundary:extract_data:"
330 : " Invalid Subset Index " << subsetIndex << ". (Valid is"
331 : " 0, .. , " << pSH->num_subsets() <<").");
332 :
333 : // check all functions
334 0 : for(size_t i=0; i < functionGroup.size(); ++i)
335 : {
336 : const size_t fct = functionGroup[i];
337 :
338 : // check if function exist
339 0 : if(fct >= m_spApproxSpace->num_fct())
340 0 : UG_THROW("DirichletBoundary:extract_data:"
341 : " Function "<< fct << " does not exist in pattern.");
342 :
343 : // check that function is defined for segment
344 0 : if(!m_spApproxSpace->is_def_in_subset(fct, subsetIndex))
345 0 : UG_THROW("DirichletBoundary:extract_data:"
346 : " Function "<<fct<<" not defined on subset "<<subsetIndex);
347 : }
348 : }
349 :
350 : // everything ok
351 0 : }
352 :
353 : template <typename TDomain, typename TAlgebra>
354 : template <typename TUserData, typename TScheduledUserData>
355 0 : void DirichletBoundary<TDomain, TAlgebra>::
356 : extract_data(std::map<int, std::vector<TUserData*> >& mvUserDataBndSegment,
357 : std::vector<TScheduledUserData>& vUserData)
358 : {
359 : // clear the extracted data
360 : mvUserDataBndSegment.clear();
361 :
362 0 : for(size_t i = 0; i < vUserData.size(); ++i)
363 : {
364 : // create Function Group and Subset Group
365 : try
366 : {
367 0 : if (! m_bInvertSubsetSelection)
368 0 : vUserData[i].ssGrp = m_spApproxSpace->subset_grp_by_name
369 : (vUserData[i].ssName.c_str());
370 : else
371 0 : vUserData[i].ssGrp = m_spApproxSpace->all_subsets_grp_except_for
372 : (vUserData[i].ssName.c_str());
373 : }
374 0 : UG_CATCH_THROW(" Subsets '"<<vUserData[i].ssName<<"' not"
375 : " all contained in ApproximationSpace.");
376 :
377 0 : FunctionGroup fctGrp;
378 : try{
379 0 : fctGrp = m_spApproxSpace->fct_grp_by_name(vUserData[i].fctName.c_str());
380 0 : }UG_CATCH_THROW(" Functions '"<<vUserData[i].fctName<<"' not"
381 : " all contained in ApproximationSpace.");
382 :
383 : // check functions and subsets
384 0 : check_functions_and_subsets(fctGrp, vUserData[i].ssGrp, TUserData::numFct);
385 :
386 : // set functions
387 0 : if(fctGrp.size() != TUserData::numFct)
388 0 : UG_THROW("LagrangeDirichletBoundary: wrong number of fct");
389 :
390 0 : for(size_t fct = 0; fct < TUserData::numFct; ++fct)
391 : {
392 0 : vUserData[i].fct[fct] = fctGrp[fct];
393 : }
394 :
395 : // loop subsets
396 0 : for(size_t si = 0; si < vUserData[i].ssGrp.size(); ++si)
397 : {
398 : // get subset index and function
399 0 : const int subsetIndex = vUserData[i].ssGrp[si];
400 :
401 : // remember functor and function
402 0 : mvUserDataBndSegment[subsetIndex].push_back(&vUserData[i]);
403 : }
404 : }
405 0 : }
406 :
407 :
408 : template <typename TDomain, typename TAlgebra>
409 0 : void DirichletBoundary<TDomain, TAlgebra>::
410 : extract_data()
411 : {
412 : // check that function pattern exists
413 0 : if(!m_spApproxSpace.valid())
414 0 : UG_THROW("DirichletBoundary:extract_data: "
415 : " Approximation Space not set.");
416 :
417 0 : extract_data(m_mNumberBndSegment, m_vNumberData);
418 0 : extract_data(m_mBNDNumberBndSegment, m_vBNDNumberData);
419 0 : extract_data(m_mConstNumberBndSegment, m_vConstNumberData);
420 0 : extract_data(m_mVectorBndSegment, m_vVectorData);
421 0 : extract_data(m_mOldNumberBndSegment, m_vOldNumberData);
422 0 : }
423 :
424 : ////////////////////////////////////////////////////////////////////////////////
425 : // assemble_dirichlet_rows
426 : ////////////////////////////////////////////////////////////////////////////////
427 :
428 : template <typename TDomain, typename TAlgebra>
429 : void DirichletBoundary<TDomain, TAlgebra>::
430 : assemble_dirichlet_rows(matrix_type& mat, ConstSmartPtr<DoFDistribution> dd, number time)
431 : {
432 : extract_data();
433 :
434 : // loop boundary subsets
435 : typename std::map<int, std::vector<CondNumberData*> >::const_iterator iter;
436 : for(iter = m_mBNDNumberBndSegment.begin(); iter != m_mBNDNumberBndSegment.end(); ++iter)
437 : {
438 : int si = (*iter).first;
439 : const std::vector<CondNumberData*>& userData = (*iter).second;
440 :
441 : DoFDistribution::traits<Vertex>::const_iterator iterBegin = dd->begin<Vertex>(si);
442 : DoFDistribution::traits<Vertex>::const_iterator iterEnd = dd->end<Vertex>(si);
443 :
444 : // create Multiindex
445 : std::vector<DoFIndex> multInd;
446 :
447 : // for readin
448 : MathVector<1> val;
449 : position_type corner;
450 :
451 : // loop vertices
452 : for(DoFDistribution::traits<Vertex>::const_iterator iter = iterBegin; iter != iterEnd; iter++)
453 : {
454 : // get vertex
455 : Vertex* vertex = *iter;
456 :
457 : // get corner position
458 : corner = m_aaPos[vertex];
459 :
460 : // loop dirichlet functions on this segment
461 : for(size_t i = 0; i < userData.size(); ++i)
462 : {
463 : // check if function is dirichlet
464 : if(!(*userData[i])(val, corner, time, si)) continue;
465 :
466 : // get function index
467 : const size_t fct = userData[i]->fct[0];
468 :
469 : // get multi indices
470 : if(dd->inner_dof_indices(vertex, fct, multInd) != 1)
471 : return;
472 :
473 : this->m_spAssTuner->set_dirichlet_row(mat, multInd[0]);
474 : }
475 : }
476 : }
477 : }
478 :
479 : ////////////////////////////////////////////////////////////////////////////////
480 : // adjust TRANSFER
481 : ////////////////////////////////////////////////////////////////////////////////
482 :
483 :
484 : template <typename TDomain, typename TAlgebra>
485 0 : void DirichletBoundary<TDomain, TAlgebra>::
486 : adjust_prolongation(matrix_type& P,
487 : ConstSmartPtr<DoFDistribution> ddFine,
488 : ConstSmartPtr<DoFDistribution> ddCoarse,
489 : int type,
490 : number time)
491 : {
492 : #ifdef LAGRANGE_DIRICHLET_ADJ_TRANSFER_FIX
493 : if (!m_bAdjustTransfers)
494 : {
495 : std::cerr << "Avoiding adjust_prolongation" << std::endl;
496 : return;
497 : }
498 : #endif
499 0 : extract_data();
500 :
501 0 : adjust_prolongation<CondNumberData>(m_mBNDNumberBndSegment, P, ddFine, ddCoarse, time);
502 0 : adjust_prolongation<NumberData>(m_mNumberBndSegment, P, ddFine, ddCoarse, time);
503 0 : adjust_prolongation<ConstNumberData>(m_mConstNumberBndSegment, P, ddFine, ddCoarse, time);
504 :
505 0 : adjust_prolongation<VectorData>(m_mVectorBndSegment, P, ddFine, ddCoarse, time);
506 :
507 0 : adjust_prolongation<OldNumberData>(m_mOldNumberBndSegment, P, ddFine, ddCoarse, time);
508 0 : }
509 :
510 : template <typename TDomain, typename TAlgebra>
511 : template <typename TUserData>
512 0 : void DirichletBoundary<TDomain, TAlgebra>::
513 : adjust_prolongation(const std::map<int, std::vector<TUserData*> >& mvUserData,
514 : matrix_type& P,
515 : ConstSmartPtr<DoFDistribution> ddFine,
516 : ConstSmartPtr<DoFDistribution> ddCoarse,
517 : number time)
518 : {
519 : // loop boundary subsets
520 : typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
521 0 : for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
522 : {
523 : // get subset index
524 0 : const int si = (*iter).first;
525 :
526 : // get vector of scheduled dirichlet data on this subset
527 0 : const std::vector<TUserData*>& vUserData = (*iter).second;
528 :
529 : // adapt jacobian for dofs in each base element type
530 : try
531 : {
532 0 : if(ddFine->max_dofs(VERTEX)) adjust_prolongation<RegularVertex, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
533 0 : if(ddFine->max_dofs(EDGE)) adjust_prolongation<Edge, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
534 0 : if(ddFine->max_dofs(FACE)) adjust_prolongation<Face, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
535 0 : if(ddFine->max_dofs(VOLUME)) adjust_prolongation<Volume, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
536 : }
537 0 : UG_CATCH_THROW("DirichletBoundary::adjust_prolongation:"
538 : " While calling 'adapt_prolongation' for TUserData, aborting.");
539 : }
540 0 : }
541 :
542 : template <typename TDomain, typename TAlgebra>
543 : template <typename TBaseElem, typename TUserData>
544 0 : void DirichletBoundary<TDomain, TAlgebra>::
545 : adjust_prolongation(const std::vector<TUserData*>& vUserData, int si,
546 : matrix_type& P,
547 : ConstSmartPtr<DoFDistribution> ddFine,
548 : ConstSmartPtr<DoFDistribution> ddCoarse,
549 : number time)
550 : {
551 : // create Multiindex
552 : std::vector<DoFIndex> vFineDoF, vCoarseDoF;
553 :
554 : // dummy for readin
555 : typename TUserData::value_type val;
556 :
557 : // position of dofs
558 : std::vector<position_type> vPos;
559 :
560 : // iterators
561 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
562 0 : iter = ddFine->begin<TBaseElem>(si);
563 0 : iterEnd = ddFine->end<TBaseElem>(si);
564 :
565 : // loop elements
566 0 : for( ; iter != iterEnd; iter++)
567 : {
568 : // get vertex
569 : TBaseElem* elem = *iter;
570 0 : GridObject* parent = m_spDomain->grid()->get_parent(elem);
571 0 : if(!parent) continue;
572 0 : if(!ddCoarse->is_contained(parent)) continue;
573 :
574 : // loop dirichlet functions on this segment
575 0 : for(size_t i = 0; i < vUserData.size(); ++i)
576 : {
577 0 : for(size_t f = 0; f < TUserData::numFct; ++f)
578 : {
579 : // get function index
580 0 : const size_t fct = vUserData[i]->fct[f];
581 :
582 : // get local finite element id
583 : const LFEID& lfeID = ddFine->local_finite_element_id(fct);
584 :
585 : // get multi indices
586 0 : ddFine->inner_dof_indices(elem, fct, vFineDoF);
587 0 : ddCoarse->inner_dof_indices(parent, fct, vCoarseDoF);
588 :
589 : // get dof position
590 : if(TUserData::isConditional){
591 0 : InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
592 : UG_ASSERT(vFineDoF.size() == vPos.size(), "Size mismatch");
593 : }
594 :
595 : // loop dofs on element
596 0 : for(size_t j = 0; j < vFineDoF.size(); ++j)
597 : {
598 : // check if function is dirichlet
599 : if(TUserData::isConditional){
600 0 : if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
601 : }
602 :
603 : SetRow(P, vFineDoF[j], 0.0);
604 : }
605 :
606 0 : if(vFineDoF.size() > 0){
607 0 : for(size_t k = 0; k < vCoarseDoF.size(); ++k){
608 0 : DoFRef(P, vFineDoF[0], vCoarseDoF[k]) = 1.0;
609 : }
610 : }
611 : }
612 : }
613 : }
614 0 : }
615 :
616 : template <typename TDomain, typename TAlgebra>
617 0 : void DirichletBoundary<TDomain, TAlgebra>::
618 : adjust_restriction(matrix_type& R,
619 : ConstSmartPtr<DoFDistribution> ddCoarse,
620 : ConstSmartPtr<DoFDistribution> ddFine,
621 : int type,
622 : number time)
623 : {
624 : #ifdef LAGRANGE_DIRICHLET_ADJ_TRANSFER_FIX
625 : if (!m_bAdjustTransfers)
626 : {
627 : std::cerr << "Avoiding adjust_restriction" << std::endl;
628 : return;
629 : }
630 : #endif
631 0 : extract_data();
632 :
633 0 : adjust_restriction<CondNumberData>(m_mBNDNumberBndSegment, R, ddCoarse, ddFine, time);
634 0 : adjust_restriction<NumberData>(m_mNumberBndSegment, R, ddCoarse, ddFine, time);
635 0 : adjust_restriction<ConstNumberData>(m_mConstNumberBndSegment, R, ddCoarse, ddFine, time);
636 :
637 0 : adjust_restriction<VectorData>(m_mVectorBndSegment, R, ddCoarse, ddFine, time);
638 :
639 0 : adjust_restriction<OldNumberData>(m_mOldNumberBndSegment, R, ddCoarse, ddFine, time);
640 0 : }
641 :
642 : template <typename TDomain, typename TAlgebra>
643 : template <typename TUserData>
644 0 : void DirichletBoundary<TDomain, TAlgebra>::
645 : adjust_restriction(const std::map<int, std::vector<TUserData*> >& mvUserData,
646 : matrix_type& R,
647 : ConstSmartPtr<DoFDistribution> ddCoarse,
648 : ConstSmartPtr<DoFDistribution> ddFine,
649 : number time)
650 : {
651 : // loop boundary subsets
652 : typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
653 0 : for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
654 : {
655 : // get subset index
656 0 : const int si = (*iter).first;
657 :
658 : // get vector of scheduled dirichlet data on this subset
659 0 : const std::vector<TUserData*>& vUserData = (*iter).second;
660 :
661 : // adapt jacobian for dofs in each base element type
662 : try
663 : {
664 0 : if(ddFine->max_dofs(VERTEX)) adjust_restriction<RegularVertex, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
665 0 : if(ddFine->max_dofs(EDGE)) adjust_restriction<Edge, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
666 0 : if(ddFine->max_dofs(FACE)) adjust_restriction<Face, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
667 0 : if(ddFine->max_dofs(VOLUME)) adjust_restriction<Volume, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
668 : }
669 0 : UG_CATCH_THROW("DirichletBoundary::adjust_restriction:"
670 : " While calling 'adjust_restriction' for TUserData, aborting.");
671 : }
672 0 : }
673 :
674 : template <typename TDomain, typename TAlgebra>
675 : template <typename TBaseElem, typename TUserData>
676 0 : void DirichletBoundary<TDomain, TAlgebra>::
677 : adjust_restriction(const std::vector<TUserData*>& vUserData, int si,
678 : matrix_type& R,
679 : ConstSmartPtr<DoFDistribution> ddCoarse,
680 : ConstSmartPtr<DoFDistribution> ddFine,
681 : number time)
682 : {
683 : // create Multiindex
684 : std::vector<DoFIndex> vFineDoF, vCoarseDoF;
685 :
686 : // dummy for readin
687 : typename TUserData::value_type val;
688 :
689 : // position of dofs
690 : std::vector<position_type> vPos;
691 :
692 : // iterators
693 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
694 0 : iter = ddFine->begin<TBaseElem>(si);
695 0 : iterEnd = ddFine->end<TBaseElem>(si);
696 :
697 : // loop elements
698 0 : for( ; iter != iterEnd; iter++)
699 : {
700 : // get vertex
701 : TBaseElem* elem = *iter;
702 0 : GridObject* parent = m_spDomain->grid()->get_parent(elem);
703 0 : if(!parent) continue;
704 0 : if(!ddCoarse->is_contained(parent)) continue;
705 :
706 : // loop dirichlet functions on this segment
707 0 : for(size_t i = 0; i < vUserData.size(); ++i)
708 : {
709 0 : for(size_t f = 0; f < TUserData::numFct; ++f)
710 : {
711 : // get function index
712 0 : const size_t fct = vUserData[i]->fct[f];
713 :
714 : // get local finite element id
715 : const LFEID& lfeID = ddFine->local_finite_element_id(fct);
716 :
717 : // get multi indices
718 0 : ddFine->inner_dof_indices(elem, fct, vFineDoF);
719 0 : ddCoarse->inner_dof_indices(parent, fct, vCoarseDoF);
720 :
721 : // get dof position
722 : if(TUserData::isConditional){
723 0 : InnerDoFPosition<TDomain>(vPos, parent, *m_spDomain, lfeID);
724 : UG_ASSERT(vCoarseDoF.size() == vPos.size(), "Size mismatch");
725 : }
726 :
727 : // loop dofs on element
728 0 : for(size_t j = 0; j < vCoarseDoF.size(); ++j)
729 : {
730 : // check if function is dirichlet
731 : if(TUserData::isConditional){
732 0 : if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
733 : }
734 :
735 : SetRow(R, vCoarseDoF[j], 0.0);
736 : }
737 :
738 0 : if(vFineDoF.size() > 0){
739 0 : for(size_t k = 0; k < vCoarseDoF.size(); ++k){
740 0 : DoFRef(R, vCoarseDoF[k], vFineDoF[0]) = 1.0;
741 : }
742 : }
743 : }
744 : }
745 : }
746 0 : }
747 :
748 : ////////////////////////////////////////////////////////////////////////////////
749 : // adjust JACOBIAN
750 : ////////////////////////////////////////////////////////////////////////////////
751 :
752 : template <typename TDomain, typename TAlgebra>
753 0 : void DirichletBoundary<TDomain, TAlgebra>::
754 : adjust_jacobian(matrix_type& J, const vector_type& u,
755 : ConstSmartPtr<DoFDistribution> dd, int type, number time,
756 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
757 : const number s_a0)
758 : {
759 0 : extract_data();
760 :
761 0 : adjust_jacobian<CondNumberData>(m_mBNDNumberBndSegment, J, u, dd, time);
762 0 : adjust_jacobian<NumberData>(m_mNumberBndSegment, J, u, dd, time);
763 0 : adjust_jacobian<ConstNumberData>(m_mConstNumberBndSegment, J, u, dd, time);
764 :
765 0 : adjust_jacobian<VectorData>(m_mVectorBndSegment, J, u, dd, time);
766 :
767 0 : adjust_jacobian<OldNumberData>(m_mOldNumberBndSegment, J, u, dd, time);
768 0 : }
769 :
770 :
771 : template <typename TDomain, typename TAlgebra>
772 : template <typename TUserData>
773 0 : void DirichletBoundary<TDomain, TAlgebra>::
774 : adjust_jacobian(const std::map<int, std::vector<TUserData*> >& mvUserData,
775 : matrix_type& J, const vector_type& u,
776 : ConstSmartPtr<DoFDistribution> dd, number time)
777 : {
778 : // loop boundary subsets
779 : typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
780 0 : for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
781 : {
782 : // get subset index
783 0 : const int si = (*iter).first;
784 :
785 : // get vector of scheduled dirichlet data on this subset
786 0 : const std::vector<TUserData*>& vUserData = (*iter).second;
787 :
788 : // adapt jacobian for dofs in each base element type
789 : try
790 : {
791 0 : if(dd->max_dofs(VERTEX))
792 0 : adjust_jacobian<RegularVertex, TUserData>(vUserData, si, J, u, dd, time);
793 0 : if(dd->max_dofs(EDGE))
794 0 : adjust_jacobian<Edge, TUserData>(vUserData, si, J, u, dd, time);
795 0 : if(dd->max_dofs(FACE))
796 0 : adjust_jacobian<Face, TUserData>(vUserData, si, J, u, dd, time);
797 0 : if(dd->max_dofs(VOLUME))
798 0 : adjust_jacobian<Volume, TUserData>(vUserData, si, J, u, dd, time);
799 : }
800 0 : UG_CATCH_THROW("DirichletBoundary::adjust_jacobian:"
801 : " While calling 'adapt_jacobian' for TUserData, aborting.");
802 : }
803 0 : }
804 :
805 : template <typename TDomain, typename TAlgebra>
806 : template <typename TBaseElem, typename TUserData>
807 0 : void DirichletBoundary<TDomain, TAlgebra>::
808 : adjust_jacobian(const std::vector<TUserData*>& vUserData, int si,
809 : matrix_type& J, const vector_type& u,
810 : ConstSmartPtr<DoFDistribution> dd, number time)
811 : {
812 : // create Multiindex
813 : std::vector<DoFIndex> multInd;
814 :
815 : // dummy for readin
816 : typename TUserData::value_type val;
817 :
818 : // position of dofs
819 : std::vector<position_type> vPos;
820 :
821 : // save all dirichlet degree of freedom indices.
822 : std::set<size_t> dirichletDoFIndices;
823 :
824 :
825 : // iterators
826 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
827 0 : iter = dd->begin<TBaseElem>(si);
828 0 : iterEnd = dd->end<TBaseElem>(si);
829 :
830 : // loop elements
831 0 : for( ; iter != iterEnd; iter++)
832 : {
833 : // get vertex
834 : TBaseElem* elem = *iter;
835 :
836 : // loop dirichlet functions on this segment
837 0 : for(size_t i = 0; i < vUserData.size(); ++i)
838 : {
839 0 : for(size_t f = 0; f < TUserData::numFct; ++f)
840 : {
841 : // get function index
842 0 : const size_t fct = vUserData[i]->fct[f];
843 :
844 : // get local finite element id
845 : const LFEID& lfeID = dd->local_finite_element_id(fct);
846 :
847 : // get multi indices
848 0 : dd->inner_dof_indices(elem, fct, multInd);
849 :
850 : // get dof position
851 : if(TUserData::isConditional){
852 0 : InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
853 : UG_ASSERT(multInd.size() == vPos.size(), "Size mismatch");
854 : }
855 :
856 : // loop dofs on element
857 0 : for(size_t j = 0; j < multInd.size(); ++j)
858 : {
859 : // check if function is dirichlet
860 : if(TUserData::isConditional){
861 0 : if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
862 : }
863 :
864 0 : this->m_spAssTuner->set_dirichlet_row(J, multInd[j]);
865 0 : if(m_bDirichletColumns)
866 : dirichletDoFIndices.insert(multInd[j][0]);
867 : }
868 : }
869 : }
870 : }
871 :
872 :
873 0 : if(m_bDirichletColumns){
874 : // UG_LOG("adjust jacobian\n")
875 :
876 : // number of rows
877 : size_t nr = J.num_rows();
878 :
879 : // run over all rows of the local matrix J and save the colums
880 : // entries for the Dirichlet indices in the map
881 :
882 : typename std::set<size_t>::iterator currentDIndex;
883 :
884 0 : for(size_t i = 0; i<nr; i++)
885 : {
886 0 : for(typename matrix_type::row_iterator it = J.begin_row(i); it!=J.end_row(i); ++it){
887 :
888 : // look if the current index is a dirichlet index
889 : // if it.index is a dirichlet index
890 : // the iterator currentDIndex is delivered otherwise set::end()
891 : currentDIndex = dirichletDoFIndices.find(it.index());
892 :
893 : // fill dirichletMap & set corresponding entry to zero
894 0 : if(currentDIndex != dirichletDoFIndices.end()){
895 : // the dirichlet-dof-index it.index is assigned
896 : // the row i and the matrix entry it.value().
897 : // if necessary for defect remove comment
898 :
899 : // m_dirichletMap[it.index()][i] = it.value();
900 :
901 : // the corresponding entry at column it.index() is set zero
902 : // this corresponds to a dirichlet column.
903 : // diagonal stays unchanged
904 0 : if(i!=it.index())
905 0 : it.value() = 0.0;
906 : }
907 :
908 : }
909 : }
910 : }
911 :
912 0 : }
913 :
914 : ////////////////////////////////////////////////////////////////////////////////
915 : // adjust DEFECT
916 : ////////////////////////////////////////////////////////////////////////////////
917 :
918 : template <typename TDomain, typename TAlgebra>
919 0 : void DirichletBoundary<TDomain, TAlgebra>::
920 : adjust_defect(vector_type& d, const vector_type& u,
921 : ConstSmartPtr<DoFDistribution> dd, int type, number time,
922 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
923 : const std::vector<number>* vScaleMass,
924 : const std::vector<number>* vScaleStiff)
925 : {
926 0 : extract_data();
927 :
928 0 : adjust_defect<CondNumberData>(m_mBNDNumberBndSegment, d, u, dd, time);
929 0 : adjust_defect<NumberData>(m_mNumberBndSegment, d, u, dd, time);
930 0 : adjust_defect<ConstNumberData>(m_mConstNumberBndSegment, d, u, dd, time);
931 :
932 0 : adjust_defect<VectorData>(m_mVectorBndSegment, d, u, dd, time);
933 :
934 0 : adjust_defect<OldNumberData>(m_mOldNumberBndSegment, d, u, dd, time);
935 0 : }
936 :
937 :
938 : template <typename TDomain, typename TAlgebra>
939 : template <typename TUserData>
940 0 : void DirichletBoundary<TDomain, TAlgebra>::
941 : adjust_defect(const std::map<int, std::vector<TUserData*> >& mvUserData,
942 : vector_type& d, const vector_type& u,
943 : ConstSmartPtr<DoFDistribution> dd, number time)
944 : {
945 : // loop boundary subsets
946 : typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
947 0 : for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
948 : {
949 : // get subset index
950 0 : const int si = (*iter).first;
951 :
952 : // get vector of scheduled dirichlet data on this subset
953 0 : const std::vector<TUserData*>& vUserData = (*iter).second;
954 :
955 : // adapt jacobian for dofs in each base element type
956 : try
957 : {
958 0 : if(dd->max_dofs(VERTEX))
959 0 : adjust_defect<RegularVertex, TUserData>(vUserData, si, d, u, dd, time);
960 0 : if(dd->max_dofs(EDGE))
961 0 : adjust_defect<Edge, TUserData>(vUserData, si, d, u, dd, time);
962 0 : if(dd->max_dofs(FACE))
963 0 : adjust_defect<Face, TUserData>(vUserData, si, d, u, dd, time);
964 0 : if(dd->max_dofs(VOLUME))
965 0 : adjust_defect<Volume, TUserData>(vUserData, si, d, u, dd, time);
966 : }
967 0 : UG_CATCH_THROW("DirichletBoundary::adjust_defect:"
968 : " While calling 'adjust_defect' for TUserData, aborting.");
969 : }
970 0 : }
971 :
972 : template <typename TDomain, typename TAlgebra>
973 : template <typename TBaseElem, typename TUserData>
974 0 : void DirichletBoundary<TDomain, TAlgebra>::
975 : adjust_defect(const std::vector<TUserData*>& vUserData, int si,
976 : vector_type& d, const vector_type& u,
977 : ConstSmartPtr<DoFDistribution> dd, number time)
978 : {
979 : // create Multiindex
980 : std::vector<DoFIndex> multInd;
981 :
982 : // dummy for readin
983 : typename TUserData::value_type val;
984 :
985 : // position of dofs
986 : std::vector<position_type> vPos;
987 :
988 : // iterators
989 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
990 0 : iter = dd->begin<TBaseElem>(si);
991 0 : iterEnd = dd->end<TBaseElem>(si);
992 :
993 : // loop elements
994 0 : for( ; iter != iterEnd; iter++)
995 : {
996 : // get vertex
997 : TBaseElem* elem = *iter;
998 :
999 : // loop dirichlet functions on this segment
1000 0 : for(size_t i = 0; i < vUserData.size(); ++i)
1001 : {
1002 0 : for(size_t f = 0; f < TUserData::numFct; ++f)
1003 : {
1004 : // get function index
1005 0 : const size_t fct = vUserData[i]->fct[f];
1006 :
1007 : // get local finite element id
1008 : const LFEID& lfeID = dd->local_finite_element_id(fct);
1009 :
1010 : // get multi indices
1011 0 : dd->inner_dof_indices(elem, fct, multInd);
1012 :
1013 : // get dof position
1014 : if(TUserData::isConditional){
1015 0 : InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1016 : UG_ASSERT(multInd.size() == vPos.size(), "Size mismatch. (multInd.size()="<<
1017 : multInd.size()<<", vPos.size()="<<vPos.size()<<")");
1018 : }
1019 :
1020 : // loop dofs on element
1021 0 : for(size_t j = 0; j < multInd.size(); ++j)
1022 : {
1023 : // check if function is dirichlet
1024 : if(TUserData::isConditional){
1025 0 : if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
1026 : }
1027 :
1028 : // set zero for dirichlet values
1029 0 : this->m_spAssTuner->set_dirichlet_val(d, multInd[j], 0.0);
1030 : }
1031 : }
1032 : }
1033 : }
1034 0 : }
1035 :
1036 : ////////////////////////////////////////////////////////////////////////////////
1037 : // adjust SOLUTION
1038 : ////////////////////////////////////////////////////////////////////////////////
1039 :
1040 : template <typename TDomain, typename TAlgebra>
1041 0 : void DirichletBoundary<TDomain, TAlgebra>::
1042 : adjust_solution(vector_type& u, ConstSmartPtr<DoFDistribution> dd, int type, number time)
1043 : {
1044 0 : extract_data();
1045 :
1046 0 : adjust_solution<CondNumberData>(m_mBNDNumberBndSegment, u, dd, time);
1047 0 : adjust_solution<NumberData>(m_mNumberBndSegment, u, dd, time);
1048 0 : adjust_solution<ConstNumberData>(m_mConstNumberBndSegment, u, dd, time);
1049 :
1050 0 : adjust_solution<VectorData>(m_mVectorBndSegment, u, dd, time);
1051 :
1052 0 : adjust_solution<OldNumberData>(m_mOldNumberBndSegment, u, dd, time);
1053 0 : }
1054 :
1055 :
1056 : template <typename TDomain, typename TAlgebra>
1057 : template <typename TUserData>
1058 0 : void DirichletBoundary<TDomain, TAlgebra>::
1059 : adjust_solution(const std::map<int, std::vector<TUserData*> >& mvUserData,
1060 : vector_type& u, ConstSmartPtr<DoFDistribution> dd, number time)
1061 : {
1062 : // loop boundary subsets
1063 : typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1064 0 : for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1065 : {
1066 : // get subset index
1067 0 : const int si = (*iter).first;
1068 :
1069 : // get vector of scheduled dirichlet data on this subset
1070 0 : const std::vector<TUserData*>& vUserData = (*iter).second;
1071 :
1072 : // adapt jacobian for dofs in each base element type
1073 : try
1074 : {
1075 0 : if(dd->max_dofs(VERTEX))
1076 0 : adjust_solution<RegularVertex, TUserData>(vUserData, si, u, dd, time);
1077 0 : if(dd->max_dofs(EDGE))
1078 0 : adjust_solution<Edge, TUserData>(vUserData, si, u, dd, time);
1079 0 : if(dd->max_dofs(FACE))
1080 0 : adjust_solution<Face, TUserData>(vUserData, si, u, dd, time);
1081 0 : if(dd->max_dofs(VOLUME))
1082 0 : adjust_solution<Volume, TUserData>(vUserData, si, u, dd, time);
1083 : }
1084 0 : UG_CATCH_THROW("DirichletBoundary::adjust_solution:"
1085 : " While calling 'adjust_solution' for TUserData, aborting.");
1086 : }
1087 0 : }
1088 :
1089 : template <typename TDomain, typename TAlgebra>
1090 : template <typename TBaseElem, typename TUserData>
1091 0 : void DirichletBoundary<TDomain, TAlgebra>::
1092 : adjust_solution(const std::vector<TUserData*>& vUserData, int si,
1093 : vector_type& u, ConstSmartPtr<DoFDistribution> dd, number time)
1094 : {
1095 : // check if the solution is to be adjusted
1096 : if (! TUserData::setSolValue)
1097 : return;
1098 :
1099 : // create Multiindex
1100 : std::vector<DoFIndex> multInd;
1101 :
1102 : // value readin
1103 : typename TUserData::value_type val;
1104 :
1105 : // position of dofs
1106 : std::vector<position_type> vPos;
1107 :
1108 : // iterators
1109 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
1110 0 : iter = dd->begin<TBaseElem>(si);
1111 0 : iterEnd = dd->end<TBaseElem>(si);
1112 :
1113 : // loop elements
1114 0 : for( ; iter != iterEnd; iter++)
1115 : {
1116 : // get vertex
1117 : TBaseElem* elem = *iter;
1118 :
1119 : // loop dirichlet functions on this segment
1120 0 : for(size_t i = 0; i < vUserData.size(); ++i)
1121 : {
1122 0 : for(size_t f = 0; f < TUserData::numFct; ++f)
1123 : {
1124 : // get function index
1125 0 : const size_t fct = vUserData[i]->fct[f];
1126 :
1127 : // get local finite element id
1128 : const LFEID& lfeID = dd->local_finite_element_id(fct);
1129 :
1130 : // get dof position
1131 0 : InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1132 :
1133 : // get multi indices
1134 0 : dd->inner_dof_indices(elem, fct, multInd);
1135 :
1136 : UG_ASSERT(multInd.size() == vPos.size(), "Size mismatch");
1137 :
1138 : // loop dofs on element
1139 0 : for(size_t j = 0; j < multInd.size(); ++j)
1140 : {
1141 : // get dirichlet value
1142 0 : if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
1143 :
1144 0 : this->m_spAssTuner->set_dirichlet_val(u, multInd[j], val[f]);
1145 : }
1146 : }
1147 : }
1148 : }
1149 0 : }
1150 :
1151 :
1152 :
1153 : ////////////////////////////////////////////////////////////////////////////////
1154 : // adjust CORRECTION
1155 : ////////////////////////////////////////////////////////////////////////////////
1156 :
1157 : template <typename TDomain, typename TAlgebra>
1158 0 : void DirichletBoundary<TDomain, TAlgebra>::adjust_correction
1159 : (
1160 : vector_type& c,
1161 : ConstSmartPtr<DoFDistribution> dd,
1162 : int type,
1163 : number time
1164 : )
1165 : {
1166 0 : extract_data();
1167 :
1168 0 : adjust_correction<CondNumberData>(m_mBNDNumberBndSegment, c, dd, time);
1169 0 : adjust_correction<NumberData>(m_mNumberBndSegment, c, dd, time);
1170 0 : adjust_correction<ConstNumberData>(m_mConstNumberBndSegment, c, dd, time);
1171 :
1172 0 : adjust_correction<VectorData>(m_mVectorBndSegment, c, dd, time);
1173 :
1174 0 : adjust_correction<OldNumberData>(m_mOldNumberBndSegment, c, dd, time);
1175 0 : }
1176 :
1177 : template <typename TDomain, typename TAlgebra>
1178 : template <typename TUserData>
1179 0 : void DirichletBoundary<TDomain, TAlgebra>::
1180 : adjust_correction(const std::map<int, std::vector<TUserData*> >& mvUserData,
1181 : vector_type& c, ConstSmartPtr<DoFDistribution> dd, number time)
1182 : {
1183 : // loop boundary subsets
1184 : typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1185 0 : for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1186 : {
1187 : // get subset index
1188 0 : const int si = (*iter).first;
1189 :
1190 : // get vector of scheduled dirichlet data on this subset
1191 0 : const std::vector<TUserData*>& vUserData = (*iter).second;
1192 :
1193 : // adapt correction for dofs in each base element type
1194 : try
1195 : {
1196 0 : if(dd->max_dofs(VERTEX))
1197 0 : adjust_correction<RegularVertex, TUserData>(vUserData, si, c, dd, time);
1198 0 : if(dd->max_dofs(EDGE))
1199 0 : adjust_correction<Edge, TUserData>(vUserData, si,c, dd, time);
1200 0 : if(dd->max_dofs(FACE))
1201 0 : adjust_correction<Face, TUserData>(vUserData, si, c, dd, time);
1202 0 : if(dd->max_dofs(VOLUME))
1203 0 : adjust_correction<Volume, TUserData>(vUserData, si, c, dd, time);
1204 : }
1205 0 : UG_CATCH_THROW("DirichletBoundary::adjust_correction:"
1206 : " While calling 'adjust_correction' for TUserData, aborting.");
1207 : }
1208 0 : }
1209 :
1210 : template <typename TDomain, typename TAlgebra>
1211 : template <typename TBaseElem, typename TUserData>
1212 0 : void DirichletBoundary<TDomain, TAlgebra>::
1213 : adjust_correction(const std::vector<TUserData*>& vUserData, int si,
1214 : vector_type& c, ConstSmartPtr<DoFDistribution> dd, number time)
1215 : {
1216 : // create Multiindex
1217 : std::vector<DoFIndex> multInd;
1218 :
1219 : // value readin
1220 : typename TUserData::value_type val;
1221 :
1222 : // position of dofs
1223 : std::vector<position_type> vPos;
1224 :
1225 : // iterators
1226 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
1227 0 : iter = dd->begin<TBaseElem>(si);
1228 0 : iterEnd = dd->end<TBaseElem>(si);
1229 :
1230 : // loop elements
1231 0 : for( ; iter != iterEnd; iter++)
1232 : {
1233 : // get vertex
1234 : TBaseElem* elem = *iter;
1235 :
1236 : // loop dirichlet functions on this segment
1237 0 : for(size_t i = 0; i < vUserData.size(); ++i)
1238 : {
1239 0 : for(size_t f = 0; f < TUserData::numFct; ++f)
1240 : {
1241 : // get function index
1242 0 : const size_t fct = vUserData[i]->fct[f];
1243 :
1244 : // get local finite element id
1245 : const LFEID& lfeID = dd->local_finite_element_id(fct);
1246 :
1247 : // get dof position
1248 0 : InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1249 :
1250 : // get multi indices
1251 0 : dd->inner_dof_indices(elem, fct, multInd);
1252 :
1253 : UG_ASSERT(multInd.size() == vPos.size(), "Size mismatch");
1254 :
1255 : // loop dofs on element
1256 0 : for(size_t j = 0; j < multInd.size(); ++j)
1257 : {
1258 : // find out whether to use dirichlet value; concrete value is of no consequence
1259 0 : if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
1260 :
1261 0 : this->m_spAssTuner->set_dirichlet_val(c, multInd[j], 0.0);
1262 : }
1263 : }
1264 : }
1265 : }
1266 0 : }
1267 :
1268 :
1269 : ////////////////////////////////////////////////////////////////////////////////
1270 : // adjust LINEAR
1271 : ////////////////////////////////////////////////////////////////////////////////
1272 :
1273 : template <typename TDomain, typename TAlgebra>
1274 0 : void DirichletBoundary<TDomain, TAlgebra>::
1275 : adjust_linear(matrix_type& A, vector_type& b,
1276 : ConstSmartPtr<DoFDistribution> dd, int type, number time)
1277 : {
1278 0 : extract_data();
1279 :
1280 0 : adjust_linear<CondNumberData>(m_mBNDNumberBndSegment, A, b, dd, time);
1281 0 : adjust_linear<NumberData>(m_mNumberBndSegment, A, b, dd, time);
1282 0 : adjust_linear<ConstNumberData>(m_mConstNumberBndSegment, A, b, dd, time);
1283 :
1284 0 : adjust_linear<VectorData>(m_mVectorBndSegment, A, b, dd, time);
1285 :
1286 0 : adjust_linear<OldNumberData>(m_mOldNumberBndSegment, A, b, dd, time);
1287 0 : }
1288 :
1289 :
1290 : template <typename TDomain, typename TAlgebra>
1291 : template <typename TUserData>
1292 0 : void DirichletBoundary<TDomain, TAlgebra>::
1293 : adjust_linear(const std::map<int, std::vector<TUserData*> >& mvUserData,
1294 : matrix_type& A, vector_type& b,
1295 : ConstSmartPtr<DoFDistribution> dd, number time)
1296 : {
1297 : // loop boundary subsets
1298 : typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1299 0 : for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1300 : {
1301 : // get subset index
1302 0 : const int si = (*iter).first;
1303 :
1304 : // get vector of scheduled dirichlet data on this subset
1305 0 : const std::vector<TUserData*>& vUserData = (*iter).second;
1306 :
1307 : // adapt jacobian for dofs in each base element type
1308 : try
1309 : {
1310 0 : if(dd->max_dofs(VERTEX))
1311 0 : adjust_linear<RegularVertex, TUserData>(vUserData, si, A, b, dd, time);
1312 0 : if(dd->max_dofs(EDGE))
1313 0 : adjust_linear<Edge, TUserData>(vUserData, si, A, b, dd, time);
1314 0 : if(dd->max_dofs(FACE))
1315 0 : adjust_linear<Face, TUserData>(vUserData, si, A, b, dd, time);
1316 0 : if(dd->max_dofs(VOLUME))
1317 0 : adjust_linear<Volume, TUserData>(vUserData, si, A, b, dd, time);
1318 : }
1319 0 : UG_CATCH_THROW("DirichletBoundary::adjust_linear:"
1320 : " While calling 'adjust_linear' for TUserData, aborting.");
1321 : }
1322 0 : }
1323 :
1324 : template <typename TDomain, typename TAlgebra>
1325 : template <typename TBaseElem, typename TUserData>
1326 0 : void DirichletBoundary<TDomain, TAlgebra>::
1327 : adjust_linear(const std::vector<TUserData*>& vUserData, int si,
1328 : matrix_type& A, vector_type& b,
1329 : ConstSmartPtr<DoFDistribution> dd, number time)
1330 : {
1331 : // create Multiindex
1332 : std::vector<DoFIndex> multInd;
1333 :
1334 : // readin value
1335 : typename TUserData::value_type val;
1336 :
1337 : // save all dirichlet degree of freedom indices.
1338 : std::set<size_t> dirichletDoFIndices;
1339 :
1340 : // position of dofs
1341 : std::vector<position_type> vPos;
1342 :
1343 : // iterators
1344 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
1345 0 : iter = dd->begin<TBaseElem>(si);
1346 0 : iterEnd = dd->end<TBaseElem>(si);
1347 :
1348 : // loop elements
1349 0 : for( ; iter != iterEnd; iter++)
1350 : {
1351 : // get vertex
1352 : TBaseElem* elem = *iter;
1353 :
1354 : // loop dirichlet functions on this segment
1355 0 : for(size_t i = 0; i < vUserData.size(); ++i)
1356 : {
1357 0 : for(size_t f = 0; f < TUserData::numFct; ++f)
1358 : {
1359 : // get function index
1360 0 : const size_t fct = vUserData[i]->fct[f];
1361 :
1362 : // get local finite element id
1363 : const LFEID& lfeID = dd->local_finite_element_id(fct);
1364 :
1365 : // get dof position
1366 0 : InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1367 :
1368 : // get multi indices
1369 0 : dd->inner_dof_indices(elem, fct, multInd);
1370 :
1371 : UG_ASSERT(multInd.size() == vPos.size(),
1372 : "Mismatch: numInd="<<multInd.size()<<", numPos="
1373 : <<vPos.size()<<" on "<<elem->reference_object_id());
1374 :
1375 : // loop dofs on element
1376 0 : for(size_t j = 0; j < multInd.size(); ++j)
1377 : {
1378 : // check if function is dirichlet and read value
1379 0 : if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
1380 :
1381 0 : this->m_spAssTuner->set_dirichlet_row(A, multInd[j]);
1382 :
1383 0 : if(m_bDirichletColumns)
1384 : {
1385 : // FIXME: Beware, this is dangerous!
1386 : // It will not work for blocked algebras.
1387 0 : UG_COND_THROW(multInd[j][1] != 0,
1388 : "adjust_linear() is not implemented for block matrices and the symmetric case!");
1389 : dirichletDoFIndices.insert(multInd[j][0]);
1390 : }
1391 :
1392 : if (TUserData::setSolValue)
1393 0 : this->m_spAssTuner->set_dirichlet_val(b, multInd[j], val[f]);
1394 : }
1395 : }
1396 : }
1397 : }
1398 :
1399 :
1400 :
1401 0 : if(m_bDirichletColumns){
1402 : // UG_LOG("adjust linear\n")
1403 0 : m_A = &A;
1404 : // number of rows
1405 : size_t nr = A.num_rows();
1406 :
1407 : typename std::set<size_t>::iterator currentDIndex;
1408 :
1409 : // run over all rows of the local matrix J and save the column
1410 : // entries for the Dirichlet indices in the map
1411 :
1412 0 : for(size_t i = 0; i<nr; i++)
1413 : {
1414 : // do not save any entries in Dirichlet rows!
1415 0 : if (dirichletDoFIndices.find(i) != dirichletDoFIndices.end())
1416 0 : continue;
1417 :
1418 0 : for(typename matrix_type::row_iterator it = A.begin_row(i); it!=A.end_row(i); ++it)
1419 : {
1420 : currentDIndex = dirichletDoFIndices.find(it.index());
1421 :
1422 : // fill dirichletMap & set corresponding entry to zero
1423 0 : if(currentDIndex != dirichletDoFIndices.end()){
1424 :
1425 : // the dirichlet-dof-index it.index is assigned
1426 : // the row i and the matrix entry it.value().
1427 0 : m_dirichletMap[si][i][it.index()] = it.value();
1428 0 : it.value() = 0.0;
1429 : }
1430 : }
1431 : }
1432 :
1433 : // TODO: for better efficiency use vectors instead of maps (rows and columns are ordered!)
1434 : typename std::map<int, std::map<int, value_type> >::iterator itdirichletMap;
1435 : typename std::map<int, value_type>::iterator itdirichletRowMap;
1436 0 : for(size_t i = 0; i<nr; i++)
1437 : {
1438 : // step over if this row did not contain any connections to Dirichlet nodes
1439 0 : if ((itdirichletMap = m_dirichletMap[si].find(i)) == m_dirichletMap[si].end())
1440 0 : continue;
1441 :
1442 0 : for(typename matrix_type::row_iterator it = m_A->begin_row(i); it!=m_A->end_row(i); ++it){
1443 :
1444 : // current column index is a dirichlet index
1445 0 : if ((itdirichletRowMap = itdirichletMap->second.find(it.index())) != itdirichletMap->second.end())
1446 0 : b[i] -= itdirichletRowMap->second*b[it.index()];
1447 : }
1448 : }
1449 : }
1450 0 : }
1451 :
1452 : ////////////////////////////////////////////////////////////////////////////////
1453 : // adjust RHS
1454 : ////////////////////////////////////////////////////////////////////////////////
1455 :
1456 : template <typename TDomain, typename TAlgebra>
1457 0 : void DirichletBoundary<TDomain, TAlgebra>::
1458 : adjust_rhs(vector_type& b, const vector_type& u,
1459 : ConstSmartPtr<DoFDistribution> dd, int type, number time)
1460 : {
1461 0 : extract_data();
1462 :
1463 0 : adjust_rhs<CondNumberData>(m_mBNDNumberBndSegment, b, u, dd, time);
1464 0 : adjust_rhs<NumberData>(m_mNumberBndSegment, b, u, dd, time);
1465 0 : adjust_rhs<ConstNumberData>(m_mConstNumberBndSegment, b, u, dd, time);
1466 :
1467 0 : adjust_rhs<VectorData>(m_mVectorBndSegment, b, u, dd, time);
1468 :
1469 0 : adjust_rhs<OldNumberData>(m_mOldNumberBndSegment, b, u, dd, time);
1470 0 : }
1471 :
1472 :
1473 : template <typename TDomain, typename TAlgebra>
1474 : template <typename TUserData>
1475 0 : void DirichletBoundary<TDomain, TAlgebra>::
1476 : adjust_rhs(const std::map<int, std::vector<TUserData*> >& mvUserData,
1477 : vector_type& b, const vector_type& u,
1478 : ConstSmartPtr<DoFDistribution> dd, number time)
1479 : {
1480 : // loop boundary subsets
1481 : typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1482 0 : for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1483 : {
1484 : // get subset index
1485 0 : const int si = (*iter).first;
1486 :
1487 : // get vector of scheduled dirichlet data on this subset
1488 0 : const std::vector<TUserData*>& vUserData = (*iter).second;
1489 :
1490 : // adapt jacobian for dofs in each base element type
1491 : try
1492 : {
1493 0 : if(dd->max_dofs(VERTEX))
1494 0 : adjust_rhs<RegularVertex, TUserData>(vUserData, si, b, u, dd, time);
1495 0 : if(dd->max_dofs(EDGE))
1496 0 : adjust_rhs<Edge, TUserData>(vUserData, si, b, u, dd, time);
1497 0 : if(dd->max_dofs(FACE))
1498 0 : adjust_rhs<Face, TUserData>(vUserData, si, b, u, dd, time);
1499 0 : if(dd->max_dofs(VOLUME))
1500 0 : adjust_rhs<Volume, TUserData>(vUserData, si, b, u, dd, time);
1501 : }
1502 0 : UG_CATCH_THROW("DirichletBoundary::adjust_rhs:"
1503 : " While calling 'adjust_rhs' for TUserData, aborting.");
1504 : }
1505 0 : }
1506 :
1507 : template <typename TDomain, typename TAlgebra>
1508 : template <typename TBaseElem, typename TUserData>
1509 0 : void DirichletBoundary<TDomain, TAlgebra>::
1510 : adjust_rhs(const std::vector<TUserData*>& vUserData, int si,
1511 : vector_type& b, const vector_type& u,
1512 : ConstSmartPtr<DoFDistribution> dd, number time)
1513 : {
1514 : // create Multiindex
1515 : std::vector<DoFIndex> multInd;
1516 :
1517 : // readin value
1518 : typename TUserData::value_type val;
1519 :
1520 : // position of dofs
1521 : std::vector<position_type> vPos;
1522 :
1523 : // iterators
1524 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
1525 0 : iter = dd->begin<TBaseElem>(si);
1526 0 : iterEnd = dd->end<TBaseElem>(si);
1527 :
1528 : // loop elements
1529 0 : for( ; iter != iterEnd; iter++)
1530 : {
1531 : // get vertex
1532 : TBaseElem* elem = *iter;
1533 :
1534 : // loop dirichlet functions on this segment
1535 0 : for(size_t i = 0; i < vUserData.size(); ++i)
1536 : {
1537 0 : for(size_t f = 0; f < TUserData::numFct; ++f)
1538 : {
1539 : // get function index
1540 0 : const size_t fct = vUserData[i]->fct[f];
1541 :
1542 : // get local finite element id
1543 : const LFEID& lfeID = dd->local_finite_element_id(fct);
1544 :
1545 : // get dof position
1546 0 : InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1547 :
1548 : // get multi indices
1549 0 : dd->inner_dof_indices(elem, fct, multInd);
1550 :
1551 : UG_ASSERT(multInd.size() == vPos.size(), "Size mismatch");
1552 :
1553 : // loop dofs on element
1554 0 : for(size_t j = 0; j < multInd.size(); ++j)
1555 : {
1556 : // check if function is dirichlet and read value
1557 0 : if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
1558 :
1559 : if (TUserData::setSolValue)
1560 0 : this->m_spAssTuner->set_dirichlet_val(b, multInd[j], val[f]);
1561 : else
1562 0 : this->m_spAssTuner->set_dirichlet_val(b, multInd[j], DoFRef(u, multInd[j]));
1563 : }
1564 : }
1565 : }
1566 :
1567 : }
1568 :
1569 :
1570 : // adjust the right hand side
1571 0 : if(m_bDirichletColumns){
1572 : typename std::map<int, std::map<int, value_type> >::iterator itdirichletMap;
1573 : typename std::map<int, value_type>::iterator itdirichletRowMap;
1574 0 : size_t nr = m_A->num_rows();
1575 0 : for(size_t i = 0; i<nr; i++)
1576 : {
1577 : // step over if this row did not contain any connections to Dirichlet nodes
1578 0 : if ((itdirichletMap = m_dirichletMap[si].find(i)) == m_dirichletMap[si].end())
1579 0 : continue;
1580 :
1581 0 : for(typename matrix_type::row_iterator it = m_A->begin_row(i); it!=m_A->end_row(i); ++it){
1582 :
1583 : // current column index is a dirichlet index
1584 0 : if ((itdirichletRowMap = itdirichletMap->second.find(it.index())) != itdirichletMap->second.end())
1585 0 : b[i] -= itdirichletRowMap->second*b[it.index()];
1586 : }
1587 : }
1588 : }
1589 0 : }
1590 :
1591 :
1592 : // //////////////////////////////////////////////////////////////////////////////
1593 : // adjust error
1594 : // //////////////////////////////////////////////////////////////////////////////
1595 :
1596 : template <typename TDomain, typename TAlgebra>
1597 0 : void DirichletBoundary<TDomain, TAlgebra>::
1598 : adjust_error
1599 : (
1600 : const vector_type& u,
1601 : ConstSmartPtr<DoFDistribution> dd,
1602 : int type,
1603 : number time,
1604 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
1605 : const std::vector<number>* vScaleMass,
1606 : const std::vector<number>* vScaleStiff
1607 : )
1608 : {
1609 : // get the error estimator data object and check that it is of the right type
1610 0 : if (this->m_spErrEstData.get() == NULL)
1611 : {
1612 0 : UG_THROW("No ErrEstData object has been given to this constraint!");
1613 : }
1614 :
1615 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
1616 :
1617 0 : if (!err_est_data)
1618 : {
1619 0 : UG_THROW("Dynamic cast to MultipleSideAndElemErrEstData failed."
1620 : << std::endl << "Make sure you handed the correct type of ErrEstData to this discretization.");
1621 : }
1622 :
1623 :
1624 0 : extract_data();
1625 :
1626 0 : adjust_error<CondNumberData>(m_mBNDNumberBndSegment, u, dd, time);
1627 0 : adjust_error<NumberData>(m_mNumberBndSegment, u, dd, time);
1628 0 : adjust_error<ConstNumberData>(m_mConstNumberBndSegment, u, dd, time);
1629 :
1630 0 : adjust_error<VectorData>(m_mVectorBndSegment, u, dd, time);
1631 :
1632 0 : adjust_error<OldNumberData>(m_mOldNumberBndSegment, u, dd, time);
1633 0 : }
1634 :
1635 : template <typename TDomain, typename TAlgebra>
1636 : template <typename TUserData>
1637 0 : void DirichletBoundary<TDomain, TAlgebra>::
1638 : adjust_error
1639 : (
1640 : const std::map<int, std::vector<TUserData*> >& mvUserData,
1641 : const vector_type& u,
1642 : ConstSmartPtr<DoFDistribution> dd,
1643 : number time
1644 : )
1645 : {
1646 : // cast error estimator data object to the right type
1647 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
1648 :
1649 : typedef typename err_est_type::side_type side_type;
1650 :
1651 : // loop boundary subsets
1652 : typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1653 0 : for (iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1654 : {
1655 : // get subset index
1656 0 : const int si = (*iter).first;
1657 :
1658 : // get vector of scheduled dirichlet data on this subset
1659 : const std::vector<TUserData*>& vUserData = (*iter).second;
1660 :
1661 : try
1662 : {
1663 : // create multi-index
1664 : std::vector<DoFIndex> multInd;
1665 :
1666 : // dummy for readin
1667 : typename TUserData::value_type val;
1668 :
1669 : // position of dofs
1670 : std::vector<position_type> vPos;
1671 :
1672 : // iterators
1673 : typename DoFDistribution::traits<side_type>::const_iterator iter, iterEnd;
1674 0 : iter = dd->begin<side_type>(si);
1675 0 : iterEnd = dd->end<side_type>(si);
1676 :
1677 : // loop elements of side_type (only!)
1678 : // elements of measure 0 in the boundary are ignored.
1679 0 : for( ; iter != iterEnd; iter++)
1680 : {
1681 : // get vertex
1682 : side_type* elem = *iter;
1683 :
1684 : // get reference object id
1685 0 : ReferenceObjectID roid = elem->reference_object_id();
1686 :
1687 : // get corner coords (for later use in calculating global IPs)
1688 : std::vector<typename TDomain::position_type> vCoCo;
1689 : CollectCornerCoordinates(vCoCo, elem, *m_spDomain, false);
1690 :
1691 : // loop dirichlet functions on this segment
1692 0 : for(size_t i = 0; i < vUserData.size(); ++i)
1693 : {
1694 0 : for(size_t f = 0; f < TUserData::numFct; ++f)
1695 : {
1696 : // get function index
1697 0 : const size_t fct = vUserData[i]->fct[f];
1698 :
1699 : // get lfeID for function
1700 0 : LFEID lfeID = dd->local_finite_element_id(fct);
1701 :
1702 : // get local and global IPs
1703 : size_t numSideIPs;
1704 : const MathVector<side_type::dim>* sideLocIPs;
1705 : const MathVector<dim>* sideGlobIPs;
1706 :
1707 : try
1708 : {
1709 0 : numSideIPs = err_est_data->get(fct)->num_side_ips(roid);
1710 0 : sideLocIPs = err_est_data->get(fct)->template side_local_ips<side_type::dim>(roid);
1711 0 : sideGlobIPs = err_est_data->get(fct)->side_global_ips(elem, &vCoCo[0]);
1712 : }
1713 0 : UG_CATCH_THROW("Global integration points for error estimator cannot be determined.");
1714 :
1715 : // loop IPs
1716 0 : for (size_t ip = 0; ip < numSideIPs; ++ip)
1717 : {
1718 : // get Dirichlet value (and do nothing, if conditional D value is false)
1719 0 : if (!(*vUserData[i])(val, sideGlobIPs[ip], time, si)) continue;
1720 :
1721 : // check if we take the values directly from the solution
1722 : if (! TUserData::setSolValue)
1723 : {
1724 0 : (*err_est_data->get(fct))(elem,ip) = 0;
1725 : continue;
1726 : }
1727 :
1728 : // evaluate shapes at ip
1729 0 : LFEID new_lfeID(lfeID.type(), lfeID.dim()-1, lfeID.order());
1730 : const LocalShapeFunctionSet<side_type::dim>& rTrialSpace =
1731 : LocalFiniteElementProvider::get<side_type::dim>(roid, new_lfeID);
1732 : std::vector<number> vShape;
1733 0 : rTrialSpace.shapes(vShape, sideLocIPs[ip]);
1734 :
1735 : // get multiindices of element
1736 : std::vector<DoFIndex> ind;
1737 0 : dd->dof_indices(elem, fct, ind);
1738 :
1739 : // compute solution at integration point
1740 : number uAtIP = 0.0;
1741 0 : for (size_t sh = 0; sh < vShape.size(); ++sh)
1742 0 : uAtIP += DoFRef(u, ind[sh]) * vShape[sh];
1743 :
1744 : // set error integrand value at IP
1745 0 : (*err_est_data->get(fct))(elem,ip) = val[f] - uAtIP;
1746 : }
1747 : }
1748 : }
1749 : }
1750 0 : }
1751 0 : UG_CATCH_THROW("DirichletBoundary::adjust_error:"
1752 : " While calling 'adjust_error' for TUserData, aborting.");
1753 : }
1754 0 : }
1755 :
1756 :
1757 : } // end namespace ug
1758 :
1759 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__CONSTRAINTS__LAGRANGE_DIRICHLET_BOUNDARY_IMPL__ */
|
